rotating binary charges

git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWave/trunk@19 f5a6acaf-da7d-456b-b0a8-35edbc60b392

5 files changed, 520 insertions, 1 deletions

diff --git a/param.ccl b/param.ccl
index df5a216..8b325c4 100644
--- a/param.ccl
+++ b/param.ccl
@@ -9,6 +9,48 @@ USES KEYWORD type ""
 
 restricted:
 
+## Parameter for rotating binary charges
+
+KEYWORD binary_source "Rotating binary source"
+{
+  "yes" :: "Rotating binary source (fast'n fancy); set binary_[radius/size/omega]"
+  "fast":: "Rotating binary source (fast'n fancy); set binary_[radius/size/omega]"
+  "slow":: "Rotating binary source (slown'n easy); set binary_[radius/size/omega]"
+  "no"  :: "No rotating binary source"
+} "yes"
+
+KEYWORD binary_verbose "Rotating binary source verbose"
+{
+  "yes"  :: "Info on charge location/extension on first iteration"
+  "debug":: "Info on charge location/extension on all iterations"
+  "no"   :: "no output"
+} "no"
+
+CCTK_REAL binary_size "Radial extension of the binary source"
+{
+  0.0: :: "Some positive value"
+} 0.5
+
+CCTK_REAL binary_omega "Frequency of the circular binary orbit"
+{
+ 0.0: :: "Some positive value"
+} 2.0
+
+CCTK_REAL binary_charge "Charge of source"
+{
+ : :: "No restriction"
+} 0.1
+
+CCTK_REAL binary_radius "Radius of the circular binary orbit"
+{
+ 0.0: :: "Some positive value"
+} 2.0
+
+
+
+
+## Parameter for initial wavepulses
+
 KEYWORD initial_data "Type of initial data"
 {
   "plane"      :: "Plane wave"
@@ -33,10 +75,12 @@ REAL kx "The wave number in the x-direction"
 {
  *:* :: "No restriction"
 } 4.0
+
 REAL ky "The wave number in the y-direction"
 {
  *:* :: "No restriction"
 } 0.0
+
 REAL kz "The wave number in the z-direction"
 {
  *:* :: "No restriction"
diff --git a/schedule.ccl b/schedule.ccl
index d5d569f..9725db8 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,6 +1,26 @@
 # Schedule definitions for thorn IDScalarWave
 # $Header$
 
+
+if (CCTK_Equals(binary_source,"yes")||CCTK_Equals(binary_source,"fast"))
+{ 
+  schedule WaveBinary at EVOL after WaveToyF77_Evolution
+  {  
+    STORAGE: wavetoy::scalarevolve
+    LANG:    Fortran
+  } "Provide binary source during evolution"
+}
+
+if (CCTK_Equals(binary_source,"slow"))
+{ 
+  schedule WaveBinarySlow at EVOL after WaveToyF77_Evolution
+  {  
+    STORAGE: wavetoy::scalarevolve
+    LANG:    Fortran
+  } "Provide binary source during evolution"
+}
+
+
 schedule IDScalarWave_CheckParameters at CCTK_PARAMCHECK
 {
   LANG: Fortran
@@ -12,3 +32,4 @@ schedule IDScalarWave_InitialData at CCTK_INITIAL
   COMMUNICATION: wavetoy::scalarevolve
   LANG:          Fortran
 } "Initial data for 3D wave equation"
+
diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c
new file mode 100644
index 0000000..59e99a1
--- /dev/null
+++ b/src/CoordinateStuff.c
@@ -0,0 +1,165 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "cctk.h"
+
+ /*@@
+   @routine    IndexFloor
+   @date       Fri Jan  7 10:34:29 2000
+   @author     Gerd Lanfermann
+   @desc 
+     For a given physical, global coordinate, IndexFloor returns
+     the closest lower grid index *locally* in the direction d, 
+     which is owned by the grid path. (Ghostzone are not owned!)
+     Other return values: 
+     -1 : gridindex not on this gridpatch
+     -3 : coordinate is outside global grid
+
+     This routine and the Fortran wrapper are Fortran/C index aware.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+int IndexFloor(cGH *GH, CCTK_REAL coord_value, int d)
+{ 
+  int ierr,index_low,idummy,ugs,lgs;
+  char *message;
+  CCTK_REAL cmin,cmax;
+
+  if (d==0)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x");
+  else if (d==1)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y");
+  else if (d==2)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"z");
+ else 
+   CCTK_WARN(0,"IndexFloor: dimension is not valid. Fortran:1..3  C:0..2");
+  
+  if ((coord_value<cmin)||(coord_value>cmax)) {
+    message = (char*)malloc(128*sizeof(char));
+    sprintf(message,
+	    "IndexFloor: coordinate value outside grid [%f,%f]:  %f \n", 
+	    cmin, cmax, coord_value); 
+    CCTK_WARN(2,message);
+    free(message);
+    return(-3);
+  }
+  
+  idummy    = floor((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]);
+  index_low = idummy-GH->cctk_lbnd[d];
+
+  if (GH->cctk_bbox[2*d]==1) 
+    lgs = 0;
+  else
+    lgs = GH->cctk_nghostzones[d];
+
+  if (GH->cctk_bbox[2*d+1]==1)
+    ugs = 0;
+  else
+    ugs = GH->cctk_nghostzones[d];
+
+  if (index_low<lgs)
+    index_low = -1;
+
+  if (index_low>=GH->cctk_lsh[d]-ugs)
+    index_low = -1;
+   
+  return(index_low);
+}
+
+
+ /*@@
+   @routine    IndexCeil
+   @date       Fri Jan  7 10:34:29 2000
+   @author     Gerd Lanfermann
+   @desc 
+     For a given physical, global coordinate, IndexCeil returns
+     the closest upper grid index *locally* in the direction d, 
+     which is owned by the grid path. (Ghostzone are not owned!)
+     Other return values: 
+     -1 : grid index not on this gridpatch
+     -3 : coordinate is outside global grid
+
+     This routine and the Fortran wrapper are Fortran/C index aware.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+
+int IndexCeil(cGH *GH, CCTK_REAL coord_value, int d)
+{  
+  int ierr,index_up,idummy,ugs,lgs;
+  char *message;
+  CCTK_REAL cmin,cmax;
+
+  if (d==0)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x");
+  else if (d==1)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y");
+  else if (d==2)
+    ierr= CCTK_CoordRange(GH,&cmin,&cmax,"z");
+  else 
+    CCTK_WARN(0,"IndexCeil: dimension is not valid. Fortran:1..3  C:0..2");
+ 
+
+  if ((coord_value<cmin)||(coord_value>cmax)) { 
+    message = (char*) malloc(256*sizeof(char));
+    sprintf(message,
+	    "IndexCeil: coordinate value outside grid [%f,%f]:  %f \n", 
+	    cmin,cmax, coord_value);
+    CCTK_WARN(2,message);
+    free(message);
+    return(-3);
+  }
+  
+  idummy   = ceil((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]);
+  index_up = idummy-GH->cctk_lbnd[d];
+
+  if (GH->cctk_bbox[2*d]==1) 
+    lgs = 0;
+  else
+    lgs = GH->cctk_nghostzones[d];
+
+  if (GH->cctk_bbox[2*d+1]==1)
+    ugs = 0;
+  else
+    ugs = GH->cctk_nghostzones[d];
+
+  if (index_up<lgs)
+    index_up = -1;
+
+  if (index_up>=GH->cctk_lsh[d]-ugs)
+    index_up = -1;
+   
+  return(index_up);
+}
+
+void FMODIFIER FORTRAN_NAME(IndexFloor)
+  (cGH *GH, CCTK_REAL *coord_value, int *index_low, int *fdim) {
+  
+  /* note the increase/decrease of indeces for fortran compatability */
+  *index_low = IndexFloor(GH, *coord_value,  (*fdim)-1);
+  if ((*index_low)>=0) (*index_low)++; 
+}
+
+void FMODIFIER FORTRAN_NAME(IndexCeil)
+  (cGH *GH, CCTK_REAL *coord_value, int *index_up, int *fdim) {
+  /* note the increase/decrease of indeces for fortran compatability */
+  *index_up = IndexCeil(GH, *coord_value, (*fdim)-1);
+  if ((*index_up)>=0) (*index_up)++; 
+}
+
+  
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
+      
diff --git a/src/make.code.defn b/src/make.code.defn
index 2fabbf2..3816f66 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = InitialData.F77 CheckParameters.F77
+SRCS = InitialData.F77 CheckParameters.F77 CoordinateStuff.c WaveBinary.F77
 
 # Subdirectories containing source files
 SUBDIRS =