Correct and clean up WaveBinarySource:

Multiply the RHS with the time step size.
Use the correct time when evaluating the RHS.
Use the correct time step size on refined levels.

Place the sources at z=0 instead of at the center of the domain.  This
is relevant for bitant mode, where the sources should still be at z=0,
not at the center of the simulated domain.

Remove complex logic to determine which processors should apply the
source term; instead, apply it on all processors.

Small other cleanups.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/WaveBinarySource/trunk@57 b9286e40-80fe-41ab-903a-d6b447012e1e

5 files changed, 54 insertions, 327 deletions

diff --git a/interface.ccl b/interface.ccl
index 624520a..6971639 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,4 +2,4 @@
 # $Header$
 
 implements: binarysource
-inherits: wavetoy grid idscalarwave
+inherits: wavetoy grid
diff --git a/schedule.ccl b/schedule.ccl
index 9341491..3340956 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -3,6 +3,6 @@
 
 schedule WaveBinaryC at EVOL after WaveToy_Evolution
 {  
-  STORAGE: wavetoy::scalarevolve[3]
+  STORAGE: wavetoy::scalarevolve[1]
   LANG:    C
 } "Provide binary source during evolution (C)"
diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c
index ae43de0..e69de29 100644
--- a/src/CoordinateStuff.c
+++ b/src/CoordinateStuff.c
@@ -1,165 +0,0 @@
- /*@@
-   @file      CoordinateStuff.c
-   @date      
-   @author    Cactus Maintainers
-   @desc 
-   Extra functions for coordinate information
-   @enddesc 
-   @version $Header$
- @@*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-
-#include "cctk.h"
-
-static const char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusWave_WaveBinarySource_CoordinateStuff_c)
-
-/********************************************************************
- *********************     Local Data Types   ***********************
- ********************************************************************/
-
-
-/********************************************************************
- ********************* Local Routine Prototypes *********************
- ********************************************************************/
-
-
-/********************************************************************
- *********************     Local Data   *****************************
- ********************************************************************/
-
-
-/********************************************************************
- *****************   External Routines Prototype ********************
- ********************************************************************/
-
-int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);
-int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);
-
-
-/********************************************************************
- **********************   External Routines  ************************
- ********************************************************************/
-
-
- /*@@
-   @routine    IndexFloorC
-   @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,
-     this includes ghostzones;
-
-     Other return values: 
-     -1 : gridindex not on this gridpatch
-     -2 : illegal dimension
-     -3 : coordinate value outside grid
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d)
-{ 
-  int index_low;
-  CCTK_REAL cmin,cmax;
-
-  if (d>=0 || d<3)
-  {
-    CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d");
-
-    if ((coord_value<cmin)||(coord_value>cmax)) 
-    {
-      CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, 
-                 "IndexFloor: coordinate value outside grid [%f,%f]:  %f \n", 
-                 cmin, cmax, coord_value); 
-      index_low = -3;
-    }
-    else
-    {
-      index_low    = floor((coord_value-cmin)/ GH->cctk_delta_space[d])
-          - GH->cctk_lbnd[d];
-
-      if (index_low<0 || index_low>=GH->cctk_lsh[d])
-      {
-        index_low = -1;
-      }
-    }
-  }
-  else 
-  {
-    CCTK_WARN(1,"IndexFloorC: dimension is not valid");
-    index_low = -2;
-  }
- 
-  return(index_low);
-}
-
-
- /*@@
-   @routine    IndexCeilC
-   @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
-     including ghostzones.
-     Other return values: 
-     -1 : grid index not on this gridpatch
-     -2 : illegal dimension
-     -3 : coordinate value outside grid
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
-@@*/
-
-
-int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d)
-{  
-  int index_up;
-  CCTK_REAL cmin,cmax;
-
-  if (d>=0 || d<3)
-  {
-    CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d");
-
-    if ((coord_value<cmin)||(coord_value>cmax)) 
-    {
-      CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, 
-                 "IndexCeil: coordinate value outside grid [%f,%f]:  %f \n", 
-                 cmin, cmax, coord_value); 
-      index_up = -3;
-    }
-    else
-    {
-      index_up   = ceil((coord_value-cmin)/ GH->cctk_delta_space[d])
-          -GH->cctk_lbnd[d];
-  
-      if (index_up<0 || index_up>=GH->cctk_lsh[d])
-      {
-        index_up = -1;
-      }
-    }
-  }
-  else 
-  {
-    CCTK_WARN(1,"IndexCeilC: dimension is not valid");
-    index_up = -2;
-  }
-
-  return(index_up);
-}
diff --git a/src/WaveBinary.c b/src/WaveBinary.c
index c02befd..f89f934 100644
--- a/src/WaveBinary.c
+++ b/src/WaveBinary.c
@@ -29,9 +29,6 @@ CCTK_FILEVERSION(CactusWave_WaveBinarySource_WaveBinary_c)
  ********************* Local Routine Prototypes *********************
  ********************************************************************/
 
-int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);
-int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);
-
 
 /********************************************************************
  *********************     Local Data   *****************************
@@ -52,176 +49,71 @@ void WaveBinaryC(CCTK_ARGUMENTS);
 void WaveBinaryC(CCTK_ARGUMENTS)
 {
   
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-
-  int i,j,k,d,vindex;
-  
-  /* the start/end points in local vindex space of the binary charge */
-  int lowerloc[3];
-  int upperloc[3];
-  
-  /* lower/upper physical coord. of one binary charge */
-  CCTK_REAL lbin[3],ubin[3];
-  CCTK_REAL minc[3],maxc[3];
-  
-  /*  lower/upper grid points, that we own (no ghostzones) */ 
-  int mingp[3], maxgp[3];
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
 
-  /* some flags */
-  int onthisproc, xproc, sign;
-  static int firstcall;
+  static int firstcall=1;
 
-  CCTK_REAL xs,ys,zs,rad;
-  CCTK_REAL zmin,zmax,charge_factor;
   const CCTK_REAL binary_pi=3.141592653589793;
+  CCTK_REAL charge_factor;
+  CCTK_REAL xsm,ysm,zsm,rad2m;
+  CCTK_REAL xsp,ysp,zsp,rad2p;
 
-  /* Because binary_charge and binary_size are a steerable parameters now,
-     charge_factor needs to be recomputed at every iteration. */
-  charge_factor = 3.0*binary_charge/
-    (4.0*binary_pi*binary_size*binary_size*binary_size);
+  int i,j,k,vindex;
 
-  /* Initialize the range arrays */
+  charge_factor = 3.0*binary_charge/(4.0*binary_pi*pow(binary_size,3));
 
-  for (d=0;d<3;d++)
-  {
-    mingp[d] = 0;
-    maxgp[d] = cctk_lsh[d]-1;
-  } 
-
-
-  /* The physical minimum/maximum coordinates of our grid chunk;
-     will use those to check if the charges are on our proc. */
-  vindex   = CCTK_GFINDEX3D(cctkGH,0,0,0);
-  minc[0] = x[vindex];
-  minc[1] = y[vindex];
-  minc[2] = z[vindex];
-  vindex   = CCTK_GFINDEX3D(cctkGH,cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1);
-  maxc[0] = x[vindex];
-  maxc[1] = y[vindex];
-  maxc[2] = z[vindex];
-        
-
-  /*  we have two charges opposite, origin is the center
-      sign multiplies the charge xy-positions by +1 and -1 */
-  for (sign=1; sign>=-1; sign-=2)
-  {
+  /* Calculate the centers of the binary sources  */
+  xsm = - binary_radius * cos(binary_omega*cctk_time);
+  ysm = - binary_radius * sin(binary_omega*cctk_time);
+  zsm = 0;
+  xsp = + binary_radius * cos(binary_omega*cctk_time);
+  ysp = + binary_radius * sin(binary_omega*cctk_time);
+  zsp = 0;
 
-    /* default flag: the charge is on this proc. */
-    onthisproc = 0;
-
-    /* calculate the center of a single binary source */
-    CCTK_CoordRange(cctkGH, &zmin, &zmax, -1, "z", "cart3d");
-    xs = sign * binary_radius * cos(binary_omega*(cctk_time-cctk_delta_time));
-    ys = sign * binary_radius * sin(binary_omega*(cctk_time-cctk_delta_time));
-    zs = (zmax-zmin)/2 + zmin;
-
-    /* shortcuts for the physical boundaries of 
-       a single charge */
-    lbin[0] = xs - binary_size;
-    lbin[1] = ys - binary_size;
-    lbin[2] = zs - binary_size;
-
-    ubin[0] = xs + binary_size;
-    ubin[1] = ys + binary_size;
-    ubin[2] = zs + binary_size;
-
-   
-    /* loop over all three dimensions */
-    for (d=0;d<3;d++)
-    {   
-     
-      /* check where our binary sources are:
-         - they can be on our grid chunk completly
-         - we have a lower part of the source on our grid chunk
-         - we have a upper part
-         - the source covers a grid chunk completely 
-         
-         Get the floor-vindex of the lower bound and 
-         ceiling-vindex of the upper bound.
-
-         Increase onthisproc : if equal to three, we have the source
-      */
-      
-      xproc = 0;
-      lowerloc[d]=mingp[d];
-      upperloc[d]=maxgp[d];
-
-      if ((lbin[d]<=minc[d]) && (ubin[d]>=maxc[d]))
-      {    
-        xproc = 1;
-      } 
-      else 
-      {
-        if ((minc[d]<=lbin[d]) && (lbin[d]<=maxc[d]))
-        {   
-          lowerloc[d]=IndexFloorC(cctkGH, lbin[d],d);
-          xproc = 1;
-        }
-        if ((minc[d]<=ubin[d]) && (ubin[d]<=maxc[d]))
-        {
-          upperloc[d]=IndexCeilC (cctkGH, ubin[d],d);
-          xproc = 1;
-        }        
-      }
-      onthisproc+=xproc;
-    }
+  if ((CCTK_EQUALS(binary_verbose, "yes") && firstcall) ||
+      CCTK_EQUALS(binary_verbose, "debug"))
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "Charge factor: %g\n",(double)charge_factor);
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "Charge center (-): %g %g %g\n",
+               (double)xsm,(double)ysm,(double)zsm);
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "Charge center (+): %g %g %g\n",
+               (double)xsp,(double)ysp,(double)zsp);
+    CCTK_VInfo(CCTK_THORNSTRING,
+               "Charge extent: %g\n",(double)binary_size);
+  }
+  firstcall=0;
 
-    /* Debugging debugging */
-    if (((firstcall==1)&&(CCTK_EQUALS(binary_verbose,"yes"))) ||
-        (CCTK_EQUALS(binary_verbose,"debug")))
-    {
-      printf("Charge: %f \n",charge_factor);
-      printf("On this proc: %d \n",onthisproc);
-      printf("cctk_lsh: %d %d %d\n",cctk_lsh[0],cctk_lsh[0],cctk_lsh[0]);
-      printf("Charge center: %f %f %f \n", xs,ys,zs);
-      printf("Charge extension: \n");
-      printf("     x-extension: %f -> %f \n",lbin[0],ubin[0]);
-      printf("     y-extension: %f -> %f \n",lbin[1],ubin[1]);
-      printf("     z-extension: %f -> %f \n",lbin[2],ubin[2]);
-      printf("Charge local vindex range: \n");
-      printf("     x-vindex   %d %d \n", lowerloc[0],upperloc[0]);
-      printf("     y-vindex   %d %d \n", lowerloc[1],upperloc[1]);
-      printf("     z-vindex   %d %d \n", lowerloc[2],upperloc[2]);
-      vindex = CCTK_GFINDEX3D(cctkGH, lowerloc[0],lowerloc[1],lowerloc[2]);
-      printf("Charge corner  LOW: [%f %f %f] \n",x[vindex],y[vindex],z[vindex]);
-      vindex = CCTK_GFINDEX3D(cctkGH, upperloc[0],upperloc[1],upperloc[2]);
-      printf("Charge corner  UP : [%f %f %f] \n",x[vindex],y[vindex],z[vindex]);
-      printf(" MINC/MAXC  [%f %f %f] [%f %f %f] \n\n\n",
-             minc[0],minc[1],minc[2],
-             maxc[0],maxc[1],maxc[2]);
-    }
-      
-    /* we have the source, now loop over the vindex boundary */
-    if (onthisproc==3)
+  for (k=0;k<cctk_lsh[2];k++)
+  {
+    for (j=0;j<=cctk_lsh[1];j++)
     {
-      for (k=lowerloc[2];k<=upperloc[2];k++)
+      for (i=0;i<=cctk_lsh[0];i++)
       {
-        for (j=lowerloc[1];j<=upperloc[1];j++)
+        vindex = CCTK_GFINDEX3D(cctkGH,i,j,k);
+
+        rad2m =
+          pow(x[vindex]-xsm,2) + pow(y[vindex]-ysm,2) + pow(z[vindex]-zsm,2);
+        rad2p =
+          pow(x[vindex]-xsp,2) + pow(y[vindex]-ysp,2) + pow(z[vindex]-zsp,2);
+
+        /* Note that both sources are positive, leading to a net
+           monopole moment  */
+        if (rad2m<pow(binary_size,2))
+        {
+          phi[vindex] += 0.5*pow(CCTK_DELTA_TIME,2)*charge_factor;
+        }
+        if (rad2p<pow(binary_size,2))
         {
-          for (i=lowerloc[0];i<=upperloc[0];i++)
-          {
-            vindex =  CCTK_GFINDEX3D(cctkGH,i,j,k);
-            rad = 
-              ((x[vindex]-xs)*(x[vindex]-xs)) +
-              ((y[vindex]-ys)*(y[vindex]-ys)) +
-              ((z[vindex]-zs)*(z[vindex]-zs));
-
-            if (rad<binary_size*binary_size) 
-            {
-              phi[vindex] = phi[vindex] + charge_factor;
-            }
-          }
+          phi[vindex] += 0.5*pow(CCTK_DELTA_TIME,2)*charge_factor;
         }
       }
     }
-
-    /* end of the sign loop */
   }
 
-  /* reset firstcall to zero */
-  firstcall = 0;
-  
-  /*  Note, that we do not need to sync anything, since each grid 
-      patch has filled out its ghostzones. */
+  /* Note that we do not need to sync anything, since each grid patch
+     has filled out its ghostzones  */
 }
diff --git a/src/make.code.defn b/src/make.code.defn
index 781200f..4c9617c 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = CoordinateStuff.c WaveBinary.c
+SRCS = WaveBinary.c
 
 # Subdirectories containing source files
 SUBDIRS =