diff options
Diffstat (limited to 'src/WaveBinary.c')
| -rw-r--r-- | src/WaveBinary.c | 204 |
1 files changed, 204 insertions, 0 deletions
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; +} + + + + + + + + + + |