c Using Cactus infrastructure #include "cctk.h" c Using Cactus parameters #include "cctk_parameters.h" c Using Cactus arguments lists #include "cctk_arguments.h" /*@@ @routine WaveBinarySlow @date Fri Jan 7 09:10:28 2000 @author Gerd Lanfermann @desc Provides the sources for rotating binary charges. It does this in a straight forward, slow approach, by looping over all grid points. @enddesc @calls CCTK_CoordRange @calledby @@*/ subroutine WaveBinarySlow(CCTK_FARGUMENTS) implicit none c Declare variables in argument list DECLARE_CCTK_FARGUMENTS c Declare parameters DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS INTEGER i,j,k,ierr INTEGER istart, jstart, kstart, iend, jend, kend CCTK_REAL dx,dy,dz,dt CCTK_REAL xs,ys,zs,rad CCTK_REAL zmin,zmax,charge_factor charge_factor = 3.0d0*binary_charge/ $ (4.0d0*3.1415*binary_size*binary_size*binary_size) call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z"); xs = binary_radius * cos(binary_omega*(cctk_time-dt)) ys = binary_radius * sin(binary_omega*(cctk_time-dt)) zs = (zmax-zmin)/2 +zmin c Set up shorthands c ----------------- dx = CCTK_DELTA_SPACE(1) dy = CCTK_DELTA_SPACE(2) dz = CCTK_DELTA_SPACE(3) dt = CCTK_DELTA_TIME istart = 2 jstart = 2 kstart = 2 iend = cctk_lsh(1)-1 jend = cctk_lsh(2)-1 kend = cctk_lsh(3)-1 do k = kstart, kend do j = jstart, jend do i = istart, iend rad =((x(i,j,k)-xs)**2)+ $ ((y(i,j,k)-ys)**2)+ $ ((z(i,j,k)-zs)**2) if (rad.le.(binary_size**2)) then phi(i,j,k) =phi(i,j,k)+ charge_factor end if rad =((x(i,j,k)+xs)**2)+ $ ((y(i,j,k)+ys)**2)+ $ ((z(i,j,k)-zs)**2) if (rad.le.(binary_size**2)) then phi(i,j,k) =phi(i,j,k)+ charge_factor end if end do end do end do return end /*@@ @routine WaveBinary @date Fri Jan 7 09:12:02 2000 @author Gerd Lanfermann @desc Provides the sources for rotating binary charges. It does this in a more intelligent approach by looping over the only gridpoints that are covered by charge. @enddesc @calls CCTK_CoordRange IndexFloor IndexCeil IndexFloor IndexCeil IndexFloor IndexCeil @calledby @history @@*/ subroutine WaveBinary(CCTK_FARGUMENTS) implicit none c Declare variables in argument list DECLARE_CCTK_FARGUMENTS c Declare parameters DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS INTEGER i,j,k,ierr,d,f c the start/end points in local index space of the binary charge INTEGER lowerloc(3) INTEGER upperloc(3) c lower/upper grid points, that we own (no ghostzones) INTEGER mingp(3),maxgp(3) c some flags INTEGER nothere,sign CCTK_REAL dx,dy,dz,dt CCTK_REAL xs,ys,zs,rad CCTK_REAL zmin,zmax,charge_factor c localgridstart/end: global physical coordinates where the c patch of grid starts/ends CCTK_REAL locgridstart(3),locgridend(3) c Things to do on first iteration INTEGER firstcall DATA firstcall /1/ SAVE firstcall,charge_factor if (firstcall.eq.1) then charge_factor = 3.0d0*binary_charge/ $ (4.0d0*3.1415*binary_size*binary_size*binary_size) end if c Initialize the range arrays do d=1,3 if (cctk_bbox(2*d-1).eq.1) then mingp(d) = 1 else mingp(d) = cctk_nghostzones(d) end if if (cctk_bbox(2*d).eq.1) then maxgp(d) = cctk_lsh(d) else maxgp(d) = cctk_lsh(d)-cctk_nghostzones(d) end if end do c we have two charges opposite, origin is the center c sign mulitplies the charge xy-positions by +1 and -1 do sign=1,-1,-2 c default flag: the charge is not here nothere = 1 c if we need to update phiold on first call c calculate the position of the binary sources call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z"); xs = sign*binary_radius * cos(binary_omega*cctk_time) ys = sign*binary_radius * sin(binary_omega*cctk_time) zs = (zmax-zmin)/2 +zmin c get the local indices for the extension of the binary source call IndexFloor(cctkGH, xs-binary_size, lowerloc(1), 1) call IndexCeil (cctkGH, xs+binary_size, upperloc(1), 1) call IndexFloor(cctkGH, ys-binary_size, lowerloc(2), 2) call IndexCeil (cctkGH, ys+binary_size, upperloc(2), 2) call IndexFloor(cctkGH, zs-binary_size, lowerloc(3), 3) call IndexCeil (cctkGH, zs+binary_size, upperloc(3), 3) do d=1,3 c Find out if the charge border can be found on our grid patch c (return value in lowerloc >= 0) if ((lowerloc(d).ge.0).or.(upperloc(d).ge.0)) then nothere = 0 if ((lowerloc(d).ge.0).and.(upperloc(d).lt.0)) then upperloc(d) = cctk_lsh(d) end if if ((upperloc(d).ge.0).and.(lowerloc(d).lt.0)) then lowerloc(d) = 1 end if else c else check if the region of charge is spread across c more than one processor: we have to make sure that the processor in c between knows that it has to loop as well! While lower/upperloc c have the values (-1), they should rather be 1 and cctk_lsh. c We need to find out if the size of the binary is covering the c physical size of the grid. To do that we look up the xyz value with c the min/max index we actually own (no ghostzones). if (d.eq.1) then if ((xs-binary_size.lt.x(mingp(1),mingp(2),mingp(3))).and. $ (xs+binary_size.gt.x(maxgp(1),maxgp(2),maxgp(3)))) then lowerloc(d)=1 upperloc(d)=cctk_lsh(d) nothere = 0 end if else if (d.eq.2) then if ((ys-binary_size.lt.y(mingp(1),mingp(2),mingp(3))).and. $ (ys+binary_size.gt.y(maxgp(1),maxgp(2),maxgp(3)))) then lowerloc(d)=1 upperloc(d)=cctk_lsh(d) nothere = 0 end if else if (d.eq.3) then if ((zs-binary_size.lt.z(mingp(1),mingp(2),mingp(3))).and. $ (zs+binary_size.gt.z(maxgp(1),maxgp(2),maxgp(3)))) then lowerloc(d)=1 upperloc(d)=cctk_lsh(d) nothere = 0 end if end if end if end do if (((firstcall.eq.1).and.(CCTK_EQUALS(binary_verbose,"yes"))) $ .or.(CCTK_EQUALS(binary_verbose,"debug"))) then write (*,*) write (*,*) "Charge: ",charge_factor write (*,*) "Charge center: xs,ys,zs ", xs,ys,zs write (*,*) "Charge extension: " write (*,*) " x-extension: ",xs-binary_size, xs+binary_size write (*,*) " y-extension: ",ys-binary_size, ys+binary_size write (*,*) " z-extension: ",zs-binary_size, zs+binary_size write (*,*) "Charge local index range" write (*,*) " x-index", lowerloc(1),upperloc(1) write (*,*) " y-index", lowerloc(2),upperloc(2) write (*,*) " z-index", lowerloc(3),upperloc(3) write (*,*) end if c Now loop over the grid points, that are covered by the charge if (nothere.eq.0) then do i=lowerloc(1),upperloc(1) do j=lowerloc(2),upperloc(2) do k=lowerloc(3),upperloc(3) rad =((x(i,j,k)-xs)**2)+ $ ((y(i,j,k)-ys)**2)+ $ ((z(i,j,k)-zs)**2) if (rad.le.(binary_size**2)) then phi(i,j,k) = phi(i,j,k)+charge_factor end if end do end do end do end if c end of the sign-loop end do c we reset firstcall to zero firstcall=0 c Note, that we do not need to sync anything, since each grid c patch has filled out its ghostzones. return end