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

5 files changed, 1 insertions, 514 deletions

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 <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 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_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 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_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 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 =