BinarySources in C: IDBinarySourceC

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

3 files changed, 363 insertions, 0 deletions

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 = 
+