aboutsummaryrefslogtreecommitdiff
path: root/src/WaveBinary.F77
diff options
context:
space:
mode:
Diffstat (limited to 'src/WaveBinary.F77')
-rw-r--r--src/WaveBinary.F77289
1 files changed, 289 insertions, 0 deletions
diff --git a/src/WaveBinary.F77 b/src/WaveBinary.F77
new file mode 100644
index 0000000..980e5f7
--- /dev/null
+++ b/src/WaveBinary.F77
@@ -0,0 +1,289 @@
+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)
+ write(*,*) "Charge: ",charge_factor
+ 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 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
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 03:58:38 +0000