From 88a0a1799137bd7ff3506301b3024fc9d3fcc366 Mon Sep 17 00:00:00 2001 From: lanfer Date: Fri, 7 Jan 2000 09:52:35 +0000 Subject: rotating binary charges git-svn-id: http://svn.cactuscode.org/arrangements/CactusWave/IDScalarWave/trunk@19 f5a6acaf-da7d-456b-b0a8-35edbc60b392 --- param.ccl | 44 ++++++++ schedule.ccl | 21 ++++ src/CoordinateStuff.c | 165 ++++++++++++++++++++++++++++ src/WaveBinary.F77 | 289 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 2 +- 5 files changed, 520 insertions(+), 1 deletion(-) create mode 100644 src/CoordinateStuff.c create mode 100644 src/WaveBinary.F77 diff --git a/param.ccl b/param.ccl index df5a216..8b325c4 100644 --- a/param.ccl +++ b/param.ccl @@ -9,6 +9,48 @@ USES KEYWORD type "" restricted: +## 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" +} "yes" + +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" + +CCTK_REAL binary_size "Radial extension of the binary source" +{ + 0.0: :: "Some positive value" +} 0.5 + +CCTK_REAL binary_omega "Frequency of the circular binary orbit" +{ + 0.0: :: "Some positive value" +} 2.0 + +CCTK_REAL binary_charge "Charge of source" +{ + : :: "No restriction" +} 0.1 + +CCTK_REAL binary_radius "Radius of the circular binary orbit" +{ + 0.0: :: "Some positive value" +} 2.0 + + + + +## Parameter for initial wavepulses + KEYWORD initial_data "Type of initial data" { "plane" :: "Plane wave" @@ -33,10 +75,12 @@ REAL kx "The wave number in the x-direction" { *:* :: "No restriction" } 4.0 + REAL ky "The wave number in the y-direction" { *:* :: "No restriction" } 0.0 + REAL kz "The wave number in the z-direction" { *:* :: "No restriction" diff --git a/schedule.ccl b/schedule.ccl index d5d569f..9725db8 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,6 +1,26 @@ # Schedule definitions for thorn IDScalarWave # $Header$ + +if (CCTK_Equals(binary_source,"yes")||CCTK_Equals(binary_source,"fast")) +{ + schedule WaveBinary at EVOL after WaveToyF77_Evolution + { + STORAGE: wavetoy::scalarevolve + LANG: Fortran + } "Provide binary source during evolution" +} + +if (CCTK_Equals(binary_source,"slow")) +{ + schedule WaveBinarySlow at EVOL after WaveToyF77_Evolution + { + STORAGE: wavetoy::scalarevolve + LANG: Fortran + } "Provide binary source during evolution" +} + + schedule IDScalarWave_CheckParameters at CCTK_PARAMCHECK { LANG: Fortran @@ -12,3 +32,4 @@ schedule IDScalarWave_InitialData at CCTK_INITIAL COMMUNICATION: wavetoy::scalarevolve LANG: Fortran } "Initial data for 3D wave equation" + diff --git a/src/CoordinateStuff.c b/src/CoordinateStuff.c new file mode 100644 index 0000000..59e99a1 --- /dev/null +++ b/src/CoordinateStuff.c @@ -0,0 +1,165 @@ +#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 ierr,index_low,idummy,ugs,lgs; + char *message; + CCTK_REAL cmin,cmax; + + if (d==0) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x"); + else if (d==1) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y"); + else if (d==2) + ierr= 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 ierr,index_up,idummy,ugs,lgs; + char *message; + CCTK_REAL cmin,cmax; + + if (d==0) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"x"); + else if (d==1) + ierr= CCTK_CoordRange(GH,&cmin,&cmax,"y"); + else if (d==2) + ierr= 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 FMODIFIER FORTRAN_NAME(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 FMODIFIER FORTRAN_NAME(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 new file mode 100644 index 0000000..980e5f7 --- /dev/null +++ b/src/WaveBinary.F77 @@ -0,0 +1,289 @@ +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_FARGUMENTS) + + implicit none + +c Declare variables in argument list + DECLARE_CCTK_FARGUMENTS + +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 + + + 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 + +c Set up shorthands +c ----------------- + dx = CCTK_DELTA_SPACE(1) + dy = CCTK_DELTA_SPACE(2) + dz = CCTK_DELTA_SPACE(3) + dt = CCTK_DELTA_TIME + + 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_FARGUMENTS) + + implicit none + +c Declare variables in argument list + DECLARE_CCTK_FARGUMENTS + +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 + + + if (firstcall.eq.1) then + charge_factor = 3.0d0*binary_charge/ + $ (4.0d0*3.1415*binary_size*binary_size*binary_size) + write(*,*) "Charge: ",charge_factor + 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 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 2fabbf2..3816f66 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 +SRCS = InitialData.F77 CheckParameters.F77 CoordinateStuff.c WaveBinary.F77 # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3