From 0c8e8432a8fdc9f2cdf3925b076e29f28abe3acc Mon Sep 17 00:00:00 2001 From: allen Date: Wed, 7 Jun 2000 20:05:52 +0000 Subject: Removing the binary source initial data to its own thorn git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWave/trunk@38 f5a6acaf-da7d-456b-b0a8-35edbc60b392 --- param.ccl | 37 ------- schedule.ccl | 20 ---- src/CoordinateStuff.c | 165 ---------------------------- src/WaveBinary.F77 | 291 -------------------------------------------------- src/make.code.defn | 2 +- 5 files changed, 1 insertion(+), 514 deletions(-) delete mode 100644 src/CoordinateStuff.c delete mode 100644 src/WaveBinary.F77 diff --git a/param.ccl b/param.ccl index 6df8a22..05d5335 100644 --- a/param.ccl +++ b/param.ccl @@ -19,43 +19,6 @@ KEYWORD initial_data "Type of initial data" private: -## Parameter for rotating binary charges - -KEYWORD binary_source "Rotating binary source" -{ - "yes" :: "Rotating binary source (fast'n fancy); set binary_[radius/size/omega]" - "fast":: "Rotating binary source (fast'n fancy); set binary_[radius/size/omega]" - "slow":: "Rotating binary source (slown'n easy); set binary_[radius/size/omega]" - "no" :: "No rotating binary source" -} "no" - -KEYWORD binary_verbose "Rotating binary source verbose" -{ - "yes" :: "Info on charge location/extension on first iteration" - "debug":: "Info on charge location/extension on all iterations" - "no" :: "no output" -} "no" - -REAL binary_size "Radial extension of the binary source" STEERABLE = ALWAYS -{ - 0.0: :: "Some positive value" -} 0.5 - -REAL binary_omega "Frequency of the circular binary orbit" STEERABLE = ALWAYS -{ - 0.0: :: "Some positive value" -} 2.0 - -REAL binary_charge "Charge of source" STEERABLE = ALWAYS -{ - : :: "No restriction" -} 0.1 - -REAL binary_radius "Radius of the circular binary orbit" STEERABLE = ALWAYS -{ - 0.0: :: "Some positive value" -} 2.0 - ## Parameter for initial wavepulses REAL radius "The radius of the gaussian wave" diff --git a/schedule.ccl b/schedule.ccl index ddc4857..7c128ca 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,26 +1,6 @@ # Schedule definitions for thorn IDScalarWave # $Header$ - -if (CCTK_Equals(binary_source,"yes")||CCTK_Equals(binary_source,"fast")) -{ - schedule WaveBinary at EVOL after WaveToy_Evolution - { - STORAGE: wavetoy::scalarevolve - LANG: Fortran - } "Provide binary source during evolution" -} - -if (CCTK_Equals(binary_source,"slow")) -{ - schedule WaveBinarySlow at EVOL after WaveToy_Evolution - { - STORAGE: wavetoy::scalarevolve - LANG: Fortran - } "Provide binary source during evolution" -} - - schedule IDScalarWave_CheckParameters at CCTK_PARAMCHECK { LANG: Fortran diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c deleted file mode 100644 index c630e4e..0000000 --- a/src/CoordinateStuff.c +++ /dev/null @@ -1,165 +0,0 @@ -#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 IndexFloor(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 IndexCeil(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); -} - -void CCTK_FCALL CCTK_FNAME(IndexFloor) - (cGH *GH, CCTK_REAL *coord_value, int *index_low, int *fdim) { - - /* note the increase/decrease of indeces for fortran compatability */ - *index_low = IndexFloor(GH, *coord_value, (*fdim)-1); - if ((*index_low)>=0) (*index_low)++; -} - -void CCTK_FCALL CCTK_FNAME(IndexCeil) - (cGH *GH, CCTK_REAL *coord_value, int *index_up, int *fdim) { - /* note the increase/decrease of indeces for fortran compatability */ - *index_up = IndexCeil(GH, *coord_value, (*fdim)-1); - if ((*index_up)>=0) (*index_up)++; -} - - diff --git a/src/WaveBinary.F77 b/src/WaveBinary.F77 deleted file mode 100644 index b457d1a..0000000 --- a/src/WaveBinary.F77 +++ /dev/null @@ -1,291 +0,0 @@ -c Using Cactus infrastructure -#include "cctk.h" - -c Using Cactus parameters -#include "cctk_Parameters.h" - -c Using Cactus arguments lists -#include "cctk_Arguments.h" - - - - /*@@ - @routine WaveBinarySlow - @date Fri Jan 7 09:10:28 2000 - @author Gerd Lanfermann - @desc - Provides the sources for rotating binary charges. - It does this in a straight forward, slow approach, - by looping over all grid points. - @enddesc - @calls CCTK_CoordRange - @calledby - @@*/ - - subroutine WaveBinarySlow(CCTK_ARGUMENTS) - - implicit none - -c Declare variables in argument list - DECLARE_CCTK_ARGUMENTS - -c Declare parameters - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - INTEGER i,j,k,ierr - INTEGER istart, jstart, kstart, iend, jend, kend - CCTK_REAL dx,dy,dz,dt - CCTK_REAL xs,ys,zs,rad - CCTK_REAL zmin,zmax,charge_factor - -c Set up shorthands -c ----------------- - dx = CCTK_DELTA_SPACE(1) - dy = CCTK_DELTA_SPACE(2) - dz = CCTK_DELTA_SPACE(3) - dt = CCTK_DELTA_TIME - - charge_factor = 3.0d0*binary_charge/ - $ (4.0d0*3.1415*binary_size*binary_size*binary_size) - - call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z"); - xs = binary_radius * cos(binary_omega*(cctk_time-dt)) - ys = binary_radius * sin(binary_omega*(cctk_time-dt)) - zs = (zmax-zmin)/2 +zmin - - - istart = 2 - jstart = 2 - kstart = 2 - - iend = cctk_lsh(1)-1 - jend = cctk_lsh(2)-1 - kend = cctk_lsh(3)-1 - - do k = kstart, kend - do j = jstart, jend - do i = istart, iend - - - - rad =((x(i,j,k)-xs)**2)+ - $ ((y(i,j,k)-ys)**2)+ - $ ((z(i,j,k)-zs)**2) - - if (rad.le.(binary_size**2)) then - phi(i,j,k) =phi(i,j,k)+ charge_factor - end if - - rad =((x(i,j,k)+xs)**2)+ - $ ((y(i,j,k)+ys)**2)+ - $ ((z(i,j,k)-zs)**2) - - if (rad.le.(binary_size**2)) then - phi(i,j,k) =phi(i,j,k)+ charge_factor - end if - - - end do - end do - end do - - return - end - - - /*@@ - @routine WaveBinary - @date Fri Jan 7 09:12:02 2000 - @author Gerd Lanfermann - @desc - Provides the sources for rotating binary charges. - It does this in a more intelligent approach - by looping over the only gridpoints that are covered by - charge. - @enddesc - @calls CCTK_CoordRange IndexFloor IndexCeil IndexFloor IndexCeil IndexFloor IndexCeil - @calledby - @history - @@*/ - - - - subroutine WaveBinary(CCTK_ARGUMENTS) - - implicit none - -c Declare variables in argument list - DECLARE_CCTK_ARGUMENTS - -c Declare parameters - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - INTEGER i,j,k,ierr,d,f - -c the start/end points in local index space of the binary charge - INTEGER lowerloc(3) - INTEGER upperloc(3) - -c lower/upper grid points, that we own (no ghostzones) - INTEGER mingp(3),maxgp(3) - -c some flags - INTEGER nothere,sign - - CCTK_REAL dx,dy,dz,dt - CCTK_REAL xs,ys,zs,rad - CCTK_REAL zmin,zmax,charge_factor - -c localgridstart/end: global physical coordinates where the -c patch of grid starts/ends - CCTK_REAL locgridstart(3),locgridend(3) - -c Things to do on first iteration - INTEGER firstcall - DATA firstcall /1/ - SAVE firstcall,charge_factor - - -c Because binary_charge and binary_size are a steerable parameters now, -c charge_factor needs to be recomputed at every iteration. -c if (firstcall.eq.1) then - charge_factor = 3.0d0*binary_charge/ - $ (4.0d0*3.1415*binary_size*binary_size*binary_size) -c end if - - - -c Initialize the range arrays - do d=1,3 - if (cctk_bbox(2*d-1).eq.1) then - mingp(d) = 1 - else - mingp(d) = cctk_nghostzones(d) - end if - if (cctk_bbox(2*d).eq.1) then - maxgp(d) = cctk_lsh(d) - else - maxgp(d) = cctk_lsh(d)-cctk_nghostzones(d) - end if - end do - - - -c we have two charges opposite, origin is the center -c sign mulitplies the charge xy-positions by +1 and -1 - do sign=1,-1,-2 - -c default flag: the charge is not here - nothere = 1 - -c if we need to update phiold on first call - -c calculate the position of the binary sources - call CCTK_CoordRange(ierr,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 - -c get the local indices for the extension of the binary source - call IndexFloor(cctkGH, xs-binary_size, lowerloc(1), 1) - call IndexCeil (cctkGH, xs+binary_size, upperloc(1), 1) - - call IndexFloor(cctkGH, ys-binary_size, lowerloc(2), 2) - call IndexCeil (cctkGH, ys+binary_size, upperloc(2), 2) - - call IndexFloor(cctkGH, zs-binary_size, lowerloc(3), 3) - call IndexCeil (cctkGH, zs+binary_size, upperloc(3), 3) - - - do d=1,3 - -c Find out if the charge border can be found on our grid patch -c (return value in lowerloc >= 0) - if ((lowerloc(d).ge.0).or.(upperloc(d).ge.0)) then - nothere = 0 - if ((lowerloc(d).ge.0).and.(upperloc(d).lt.0)) then - upperloc(d) = cctk_lsh(d) - end if - if ((upperloc(d).ge.0).and.(lowerloc(d).lt.0)) then - lowerloc(d) = 1 - end if - else -c else check if the region of charge is spread across -c more than one processor: we have to make sure that the processor in -c between knows that it has to loop as well! While lower/upperloc -c have the values (-1), they should rather be 1 and cctk_lsh. -c We need to find out if the size of the binary is covering the -c physical size of the grid. To do that we look up the xyz value with -c the min/max index we actually own (no ghostzones). - - if (d.eq.1) then - if ((xs-binary_size.lt.x(mingp(1),mingp(2),mingp(3))).and. - $ (xs+binary_size.gt.x(maxgp(1),maxgp(2),maxgp(3)))) then - lowerloc(d)=1 - upperloc(d)=cctk_lsh(d) - nothere = 0 - end if - else if (d.eq.2) then - if ((ys-binary_size.lt.y(mingp(1),mingp(2),mingp(3))).and. - $ (ys+binary_size.gt.y(maxgp(1),maxgp(2),maxgp(3)))) then - lowerloc(d)=1 - upperloc(d)=cctk_lsh(d) - nothere = 0 - end if - else if (d.eq.3) then - if ((zs-binary_size.lt.z(mingp(1),mingp(2),mingp(3))).and. - $ (zs+binary_size.gt.z(maxgp(1),maxgp(2),maxgp(3)))) then - lowerloc(d)=1 - upperloc(d)=cctk_lsh(d) - nothere = 0 - end if - end if - end if - end do - - if (((firstcall.eq.1).and.(CCTK_EQUALS(binary_verbose,"yes"))) - $ .or.(CCTK_EQUALS(binary_verbose,"debug"))) then - write (*,*) - write (*,*) "Charge: ",charge_factor - write (*,*) "Charge center: xs,ys,zs ", xs,ys,zs - write (*,*) "Charge extension: " - write (*,*) " x-extension: ",xs-binary_size, xs+binary_size - write (*,*) " y-extension: ",ys-binary_size, ys+binary_size - write (*,*) " z-extension: ",zs-binary_size, zs+binary_size - write (*,*) "Charge local index range" - write (*,*) " x-index", lowerloc(1),upperloc(1) - write (*,*) " y-index", lowerloc(2),upperloc(2) - write (*,*) " z-index", lowerloc(3),upperloc(3) - write (*,*) - end if - -c Now loop over the grid points, that are covered by the charge - if (nothere.eq.0) then - do i=lowerloc(1),upperloc(1) - do j=lowerloc(2),upperloc(2) - do k=lowerloc(3),upperloc(3) - rad =((x(i,j,k)-xs)**2)+ - $ ((y(i,j,k)-ys)**2)+ - $ ((z(i,j,k)-zs)**2) - if (rad.le.(binary_size**2)) then - phi(i,j,k) = phi(i,j,k)+charge_factor - end if - end do - end do - end do - end if - -c end of the sign-loop - end do - - -c we reset firstcall to zero - firstcall=0 - -c Note, that we do not need to sync anything, since each grid -c patch has filled out its ghostzones. - - return - end - diff --git a/src/make.code.defn b/src/make.code.defn index 3816f66..ef4c3d5 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = InitialData.F77 CheckParameters.F77 CoordinateStuff.c WaveBinary.F77 +SRCS = InitialData.F77 CheckParameters.F77 # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3