From 147be29173d4c68db8771ee8813ffe0a9b018c1a Mon Sep 17 00:00:00 2001 From: schnetter Date: Sat, 5 Dec 2009 05:10:48 +0000 Subject: Correct and clean up WaveBinarySource: Multiply the RHS with the time step size. Use the correct time when evaluating the RHS. Use the correct time step size on refined levels. Place the sources at z=0 instead of at the center of the domain. This is relevant for bitant mode, where the sources should still be at z=0, not at the center of the simulated domain. Remove complex logic to determine which processors should apply the source term; instead, apply it on all processors. Small other cleanups. git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/WaveBinarySource/trunk@57 b9286e40-80fe-41ab-903a-d6b447012e1e --- interface.ccl | 2 +- schedule.ccl | 2 +- src/CoordinateStuff.c | 165 --------------------------------------- src/WaveBinary.c | 210 ++++++++++++-------------------------------------- src/make.code.defn | 2 +- 5 files changed, 54 insertions(+), 327 deletions(-) diff --git a/interface.ccl b/interface.ccl index 624520a..6971639 100644 --- a/interface.ccl +++ b/interface.ccl @@ -2,4 +2,4 @@ # $Header$ implements: binarysource -inherits: wavetoy grid idscalarwave +inherits: wavetoy grid diff --git a/schedule.ccl b/schedule.ccl index 9341491..3340956 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -3,6 +3,6 @@ schedule WaveBinaryC at EVOL after WaveToy_Evolution { - STORAGE: wavetoy::scalarevolve[3] + STORAGE: wavetoy::scalarevolve[1] LANG: C } "Provide binary source during evolution (C)" diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c index ae43de0..e69de29 100644 --- a/src/CoordinateStuff.c +++ b/src/CoordinateStuff.c @@ -1,165 +0,0 @@ - /*@@ - @file CoordinateStuff.c - @date - @author Cactus Maintainers - @desc - Extra functions for coordinate information - @enddesc - @version $Header$ - @@*/ - -#include -#include -#include -#include - -#include "cctk.h" - -static const char *rcsid = "$Header$"; - -CCTK_FILEVERSION(CactusWave_WaveBinarySource_CoordinateStuff_c) - -/******************************************************************** - ********************* Local Data Types *********************** - ********************************************************************/ - - -/******************************************************************** - ********************* Local Routine Prototypes ********************* - ********************************************************************/ - - -/******************************************************************** - ********************* Local Data ***************************** - ********************************************************************/ - - -/******************************************************************** - ***************** External Routines Prototype ******************** - ********************************************************************/ - -int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d); -int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d); - - -/******************************************************************** - ********************** External Routines ************************ - ********************************************************************/ - - - /*@@ - @routine IndexFloorC - @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, - this includes ghostzones; - - Other return values: - -1 : gridindex not on this gridpatch - -2 : illegal dimension - -3 : coordinate value outside grid - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - -int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d) -{ - int index_low; - CCTK_REAL cmin,cmax; - - if (d>=0 || d<3) - { - CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d"); - - if ((coord_valuecmax)) - { - CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, - "IndexFloor: coordinate value outside grid [%f,%f]: %f \n", - cmin, cmax, coord_value); - index_low = -3; - } - else - { - index_low = floor((coord_value-cmin)/ GH->cctk_delta_space[d]) - - GH->cctk_lbnd[d]; - - if (index_low<0 || index_low>=GH->cctk_lsh[d]) - { - index_low = -1; - } - } - } - else - { - CCTK_WARN(1,"IndexFloorC: dimension is not valid"); - index_low = -2; - } - - return(index_low); -} - - - /*@@ - @routine IndexCeilC - @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 - including ghostzones. - Other return values: - -1 : grid index not on this gridpatch - -2 : illegal dimension - -3 : coordinate value outside grid - @enddesc - @calls - @calledby - @history - - @endhistory - -@@*/ - - -int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d) -{ - int index_up; - CCTK_REAL cmin,cmax; - - if (d>=0 || d<3) - { - CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d"); - - if ((coord_valuecmax)) - { - CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, - "IndexCeil: coordinate value outside grid [%f,%f]: %f \n", - cmin, cmax, coord_value); - index_up = -3; - } - else - { - index_up = ceil((coord_value-cmin)/ GH->cctk_delta_space[d]) - -GH->cctk_lbnd[d]; - - if (index_up<0 || index_up>=GH->cctk_lsh[d]) - { - index_up = -1; - } - } - } - else - { - CCTK_WARN(1,"IndexCeilC: dimension is not valid"); - index_up = -2; - } - - return(index_up); -} 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