diff options
Diffstat (limited to 'src/WaveBinary.c')
-rw-r--r-- | src/WaveBinary.c | 210 |
1 files changed, 51 insertions, 159 deletions
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 */ } |