diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/CoordinateStuff.c | 165 | ||||
-rw-r--r-- | src/WaveBinary.F77 | 289 | ||||
-rw-r--r-- | src/make.code.defn | 2 |
3 files changed, 455 insertions, 1 deletions
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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#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_value<cmin)||(coord_value>cmax)) { + 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<lgs) + index_low = -1; + + 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_value<cmin)||(coord_value>cmax)) { + 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<lgs) + index_up = -1; + + 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 = |