From ad2cfd1860cd602045d171e17b30b4d12509c989 Mon Sep 17 00:00:00 2001 From: lanfer Date: Fri, 14 Jul 2000 13:13:50 +0000 Subject: BinarySources in C: IDBinarySourceC git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/WaveBinarySource/trunk@2 b9286e40-80fe-41ab-903a-d6b447012e1e --- src/CoordinateStuff.c | 150 +++++++++++++++++++++++++++++++++++++ src/WaveBinary.c | 204 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 9 +++ 3 files changed, 363 insertions(+) create mode 100644 src/CoordinateStuff.c create mode 100644 src/WaveBinary.c create mode 100644 src/make.code.defn (limited to 'src') 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 +#include +#include +#include + +#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_valuecmax)) { + 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=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_valuecmax)) { + 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=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 +#include +#include + +/* 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