aboutsummaryrefslogtreecommitdiff
path: root/src/WaveBinary.F77
diff options
context:
space:
mode:
Diffstat (limited to 'src/WaveBinary.F77')
-rw-r--r--src/WaveBinary.F77291
1 files changed, 0 insertions, 291 deletions
diff --git a/src/WaveBinary.F77 b/src/WaveBinary.F77
deleted file mode 100644
index b457d1a..0000000
--- a/src/WaveBinary.F77
+++ /dev/null
@@ -1,291 +0,0 @@
-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_ARGUMENTS)
-
- implicit none
-
-c Declare variables in argument list
- DECLARE_CCTK_ARGUMENTS
-
-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
-
-c Set up shorthands
-c -----------------
- dx = CCTK_DELTA_SPACE(1)
- dy = CCTK_DELTA_SPACE(2)
- dz = CCTK_DELTA_SPACE(3)
- dt = CCTK_DELTA_TIME
-
- 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
-
-
- 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_ARGUMENTS)
-
- implicit none
-
-c Declare variables in argument list
- DECLARE_CCTK_ARGUMENTS
-
-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
-
-
-c Because binary_charge and binary_size are a steerable parameters now,
-c charge_factor needs to be recomputed at every iteration.
-c if (firstcall.eq.1) then
- charge_factor = 3.0d0*binary_charge/
- $ (4.0d0*3.1415*binary_size*binary_size*binary_size)
-c 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
-
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 06:51:46 +0000