BinarySources in C: IDBinarySourceC

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

13 files changed, 492 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..6f74b98
--- /dev/null
+++ b/README
@@ -0,0 +1,7 @@
+Cactus Code Thorn IDBinarySourceC
+Authors    : Gerd Lanfermann
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn: provides a source term to the scalar field evolution
+for two rotating binary charges.
diff --git a/doc/README b/doc/README
new file mode 100644
index 0000000..666abab
--- /dev/null
+++ b/doc/README
@@ -0,0 +1,7 @@
+
+
+This thorn adds a source term to the evolution equation. It is
+scheduled to run after wavetoy evolution and adds data corresponding
+to two rotating binary charges. It simply overwrites the scalarfield
+grid variable where the binary sources are located. The binary sources 
+are not evolved. they are rotated based on the evolution time.
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..6971639
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn IDBinarySourceC
+# $Header$
+
+implements: binarysource
+inherits: wavetoy grid
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..e081016
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,39 @@
+# Parameter definitions for thorn IDBinarySourceC
+# $Header$
+
+shares: idscalarwave
+
+EXTENDS KEYWORD initial_data ""
+{
+ "binarysource" ::  "Rotating binary source (fast'n fancy); set binary_[radius/size/omega]"
+}
+
+private:
+
+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"
+
+REAL binary_size "Radial extension of the binary source" STEERABLE = ALWAYS
+{
+  0.0: :: "Some positive value"
+} 0.5
+
+REAL binary_omega "Frequency of the circular binary orbit" STEERABLE = ALWAYS
+{
+ 0.0: :: "Some positive value"
+} 2.0
+
+REAL binary_charge "Charge of source" STEERABLE = ALWAYS
+{
+ : :: "No restriction"
+} 0.1
+
+REAL binary_radius "Radius of the circular binary orbit" STEERABLE = ALWAYS
+{
+ 0.0: :: "Some positive value"
+} 2.0
+
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..f8efc49
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,12 @@
+# Schedule definitions for thorn IDBinarySourceC
+# $Header$
+
+if (CCTK_Equals(initial_data,"binarysource"))
+{ 
+  schedule WaveBinaryC at EVOL after WaveToy_Evolution
+  {  
+    STORAGE: wavetoy::scalarevolve
+    LANG:    C
+  } "Provide binary source during evolution (C)"
+}
+
diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c
new file mode 100644
index 0000000..377b141
--- /dev/null
+++ b/src/CoordinateStuff.c
@@ -0,0 +1,150 @@
+#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 IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d)
+{ 
+  int index_low,idummy,ugs,lgs;
+  char *message;
+  CCTK_REAL cmin,cmax;
+
+  if (d==0)
+    CCTK_CoordRange(GH,&cmin,&cmax,"x");
+  else if (d==1)
+    CCTK_CoordRange(GH,&cmin,&cmax,"y");
+  else if (d==2)
+    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 IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d)
+{  
+  int index_up,idummy,ugs,lgs;
+  char *message;
+  CCTK_REAL cmin,cmax;
+
+  if (d==0)
+    CCTK_CoordRange(GH,&cmin,&cmax,"x");
+  else if (d==1)
+    CCTK_CoordRange(GH,&cmin,&cmax,"y");
+  else if (d==2)
+    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);
+}
+
+ 
diff --git a/src/WaveBinary.c b/src/WaveBinary.c
new file mode 100644
index 0000000..6a8ca79
--- /dev/null
+++ b/src/WaveBinary.c
@@ -0,0 +1,204 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/* Using Cactus infrastructure */
+#include "cctk.h" 
+
+/* Using Cactus parameters */
+#include "cctk_Parameters.h"
+
+/* Using Cactus arguments lists */
+#include "cctk_Arguments.h"
+
+int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);
+int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);
+
+void WaveBinaryC(CCTK_ARGUMENTS)
+{
+  
+  /* Declare variables in argument list */
+  DECLARE_CCTK_ARGUMENTS
+    
+  /* Declare parameters */
+  DECLARE_CCTK_PARAMETERS
+
+  int i,j,k,d,index;
+  
+  /* the start/end points in local index 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];
+
+  /* some flags */
+  int onthisproc, sign;
+  static int firstcall;
+  
+  CCTK_REAL xs,ys,zs,rad;
+  CCTK_REAL zmin,zmax,charge_factor;
+
+  /* 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*3.1415*binary_size*binary_size*binary_size);
+
+  /* Initialize the range arrays */
+  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. */
+  index   = CCTK_GFINDEX3D(cctkGH,0,0,0);
+  minc[0] = x[index];
+  minc[1] = y[index];
+  minc[2] = z[index];
+  index   = CCTK_GFINDEX3D(cctkGH,cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1);
+  maxc[0] = x[index];
+  maxc[1] = y[index];
+  maxc[2] = z[index];
+        
+
+  /*  we have two charges opposite, origin is the center
+      sign mulitplies the charge xy-positions by +1 and -1 */
+  for (sign=1; sign>=-1; sign-=2)
+  {
+
+    /* default flag: the charge is on this proc. */
+    onthisproc = 0;
+
+    /* calculate the center of a single binary source */
+    CCTK_CoordRange(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;
+
+    /* 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 binaries are:
+	 - the 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 completetly 
+	 
+	 Get the floor-index of the lower bound and 
+	 ceiling-index of the upper bound.
+
+	 Increase onthisproc : if equal to three, we have the source
+      */
+      if ((lbin[d]>=minc[d]) && (ubin[d]<=minc[d]))
+      {   
+	lowerloc[d]=IndexFloorC(cctkGH, lbin[d],d);
+	upperloc[d]=IndexCeilC (cctkGH, ubin[d],d);
+	onthisproc += 1;
+      }
+      else if ((lbin[d]>=minc[d]) && (lbin[d]<=maxc[d]))
+      {
+	lowerloc[d]=IndexFloorC(cctkGH, lbin[d],d);
+	upperloc[d]=maxgp[d];
+	onthisproc += 1;
+      }
+      else if ((ubin[d]>=minc[d]) && (ubin[d]<=maxc[d]))
+      {
+	lowerloc[d]=mingp[d];
+	upperloc[d]=IndexCeilC (cctkGH, ubin[d],d);
+	onthisproc += 1;
+      }
+      else if ((ubin[d]>=maxc[d]) && (lbin[d]<=minc[d]))
+      {
+	lowerloc[d]=mingp[d];
+	upperloc[d]=maxgp[d];
+	onthisproc += 1;
+      }
+      else {
+	lowerloc[d] = -1;
+	upperloc[d] = -1;
+      }
+    }
+
+    /* 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 index range: \n");
+      printf("     x-index   %d %d \n", lowerloc[0],upperloc[0]);
+      printf("     y-index   %d %d \n", lowerloc[1],upperloc[1]);
+      printf("     z-index   %d %d \n", lowerloc[2],upperloc[2]);
+      index = CCTK_GFINDEX3D(cctkGH, lowerloc[0],lowerloc[1],lowerloc[2]);
+      printf("Charge corner  LOW: [%f %f %f] \n",x[index],y[index],z[index]);
+      index = CCTK_GFINDEX3D(cctkGH, upperloc[0],upperloc[1],upperloc[2]);
+      printf("Charge corner  UP : [%f %f %f] \n",x[index],y[index],z[index]);
+      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, no loop over the index boudary */
+    if (onthisproc==3)
+    {
+      for (k=lowerloc[2];k<=upperloc[2];k++)
+	for (j=lowerloc[1];j<=upperloc[1];j++)
+	  for (i=lowerloc[0];i<=upperloc[0];i++)
+	  {
+	    index =  CCTK_GFINDEX3D(cctkGH,i,j,k);
+	    rad = 
+	      ((x[index]-xs)*(x[index]-xs)) +
+	      ((y[index]-ys)*(y[index]-ys)) +
+	      ((z[index]-zs)*(z[index]-zs));
+
+	    if (rad<binary_size*binary_size) 
+	      phi[index] = phi[index] + 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. */
+
+  return;
+}
+    
+
+      
+	
+	
+
+    
+
+  
+  
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..781200f
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn IDBinarySourceC
+# $Header$
+
+# Source files in this directory
+SRCS = CoordinateStuff.c WaveBinary.c
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/test/test_binary_1.par b/test/test_binary_1.par
new file mode 100644
index 0000000..f4c2913
--- /dev/null
+++ b/test/test_binary_1.par
@@ -0,0 +1,31 @@
+!DESC "Wavetoy scalarfield evolution with rotating binary sources (C)"
+ActiveThorns = "time wavetoyc pugh cartgrid3d ioutil ioascii iobasic idscalarwaveC idbinarysourceC"
+
+time::dtfac = 0.25
+
+idbinarysourceC::binary_omega  = 26 
+idbinarysourceC::binary_charge = 0.00001
+idbinarysourceC::binary_radius = 0.07
+idbinarysourceC::binary_size   = 0.04
+
+wavetoy::bound         = "radiation"
+
+idscalarwave::initial_data= "binarysource"
+
+grid::type   = "BySpacing"
+grid::domain = "full"
+grid::dxyz   = 0.01
+
+driver::global_nx = 20
+driver::global_ny = 20
+driver::global_nz = 20
+
+cactus::cctk_itlast = 5
+
+IOBasic::outScalar_every = 1
+IOBasic::outScalar_vars  = "wavetoy::phi"
+
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars  = "wavetoy::phi"
+
+IO::outdir = "test_binary_1"
diff --git a/test/test_binary_1/phi_max.tl b/test/test_binary_1/phi_max.tl
new file mode 100644
index 0000000..6b944db
--- /dev/null
+++ b/test/test_binary_1/phi_max.tl
@@ -0,0 +1,7 @@
+"phi max v time
+0.000000           0.0000000000000
+0.002500           0.0581927140823
+0.005000           0.1124229518062
+0.007500           0.1485896005588
+0.010000           0.1603854727938
+0.012500           0.1466981949482
diff --git a/test/test_binary_1/phi_min.tl b/test/test_binary_1/phi_min.tl
new file mode 100644
index 0000000..cd0653c
--- /dev/null
+++ b/test/test_binary_1/phi_min.tl
@@ -0,0 +1,7 @@
+"phi min v time
+0.000000           0.0000000000000
+0.002500          -0.0000000000000
+0.005000          -0.0373030399491
+0.007500          -0.1119091198472
+0.010000          -0.2238182396944
+0.012500          -0.3730204043529
diff --git a/test/test_binary_1/phi_nm1.tl b/test/test_binary_1/phi_nm1.tl
new file mode 100644
index 0000000..068e992
--- /dev/null
+++ b/test/test_binary_1/phi_nm1.tl
@@ -0,0 +1,7 @@
+"phi norm1 v time
+0.000000           0.0000000000000
+0.002500           0.0005706410186
+0.005000           0.0032716809104
+0.007500           0.0078138538369
+0.010000           0.0141102471339
+0.012500           0.0221190898545
diff --git a/test/test_binary_1/phi_nm2.tl b/test/test_binary_1/phi_nm2.tl
new file mode 100644
index 0000000..4644ac4
--- /dev/null
+++ b/test/test_binary_1/phi_nm2.tl
@@ -0,0 +1,7 @@
+"phi norm2 v time
+0.000000           0.0000000000000
+0.002500           0.0055341917418
+0.005000           0.0137955970905
+0.007500           0.0288227405794
+0.010000           0.0507403802788
+0.012500           0.0791238991074