Remove run-time dependency on Carpet.

Use CCTK_VInfo instead of printf.
Improve schedule during initialisation.
Some cleanup.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/SphericalSurface/trunk@31 40f6ab95-0e4f-0410-8daa-ee8d7420be1d

3 files changed, 98 insertions, 96 deletions

diff --git a/interface.ccl b/interface.ccl
index c9dd373..bb66f2c 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -4,7 +4,13 @@
 IMPLEMENTS: SphericalSurface
 inherits: grid
 
-USES INCLUDE HEADER: carpet.hh
+
+
+# Get current refinement level
+CCTK_INT FUNCTION GetRefinementLevel (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION GetRefinementLevel
+
+
 
 PUBLIC:
 
diff --git a/schedule.ccl b/schedule.ccl
index 034bdd9..53827c4 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -12,24 +12,28 @@ STORAGE: sf_minreflevel
 STORAGE: sf_maxreflevel
 
 
+
+SCHEDULE SphericalSurface_SetupRes AT basegrid BEFORE SphericalSurface_Setup AFTER SpatialCoordinates AFTER CorrectCoordinates
+{
+   LANG: C
+} "Set surface resolution automatically"
+
 SCHEDULE SphericalSurface_Setup AT basegrid
 {
   LANG: C
-  #OPTIONS: global
-  OPTIONS: local
 } "Calculate surface coordinate descriptors"
 
 
 
-SCHEDULE GROUP SphericalSurface_HasBeenSet AT poststep
+SCHEDULE SphericalSurface_Set AT basegrid BEFORE SphericalSurface_HasBeenSet
 {
-} "Set the spherical surfaces before this group, and use it afterwards"
+  LANG: C
+  OPTIONS: global
+} "Set surface radii to be used for initial setup in other thorns"
 
-  SCHEDULE SphericalSurface_CheckState IN SphericalSurface_HasBeenSet
-  {
-    LANG: C
-    OPTIONS: global
-  } "Test the state of the spherical surfaces"
+SCHEDULE GROUP SphericalSurface_HasBeenSet AT basegrid
+{
+} "Set the spherical surfaces before this group, and use it afterwards"
 
 
 
@@ -39,16 +43,14 @@ SCHEDULE SphericalSurface_Set AT poststep BEFORE SphericalSurface_HasBeenSet
   OPTIONS: global
 } "Set surface radii"
 
-
-SCHEDULE SphericalSurface_Set AT basegrid AFTER SphericalSurface_Setup
+SCHEDULE GROUP SphericalSurface_HasBeenSet AT poststep
 {
-  LANG: C
-  OPTIONS: local
-} "Set surface radii to be used for initial setup in other thorns"
+} "Set the spherical surfaces before this group, and use it afterwards"
 
 
-SCHEDULE SphericalSurface_SetupRes AT basegrid BEFORE SphericalSurface_Setup AFTER CartGrid3D_SetCoordinates
+
+SCHEDULE SphericalSurface_CheckState IN SphericalSurface_HasBeenSet
 {
-   LANG: C
-   OPTIONS: local
-} "Set surface resolution automatically"
+  LANG: C
+  OPTIONS: global
+} "Test the state of the spherical surfaces"
diff --git a/src/setup.cc b/src/setup.cc
index f7aa0be..dc24e82 100644
--- a/src/setup.cc
+++ b/src/setup.cc
@@ -1,24 +1,25 @@
-/* $Header$ */
+// $Header$
 
 #include <assert.h>
+#include <math.h>
 #include <stdio.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#include "mpi.h"
-#include <math.h>
-
-#include "carpet.hh"
-
-
-//#define DEBUG
+#include <mpi.h>
 
 
-#define round(x) (x<0?ceil((x)-0.5):floor((x)+0.5))
 
-#define INITIAL_VAL 10000
+static int get_reflevel (cGH const * const cctkGH)
+{
+  if (CCTK_IsFunctionAliased ("GetRefinementLevel")) {
+    return GetRefinementLevel (cctkGH);
+  } else {
+    return 0;
+  }
+}
 
 
 //TODO:  overlap with finer levels still seems to be incorrect!
@@ -29,13 +30,8 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  CCTK_REAL const pi = 3.1415926535897932384626433832795028841971693993751;
-  
-  int group;
-  cGroup groupinfo;
-  cGroupDynamicData groupdata;
+  CCTK_REAL const pi = 4 * atan(CCTK_REAL(1.0));
   
-  int n;
   int ierr;
   
   
@@ -44,22 +40,24 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
   
   
   
-  group = CCTK_GroupIndex ("SphericalSurface::sf_radius");
+  int const group = CCTK_GroupIndex ("SphericalSurface::sf_radius");
   assert (group>=0);
+  cGroup groupinfo;
   ierr = CCTK_GroupData (group, &groupinfo);
   assert (!ierr);
+  cGroupDynamicData groupdata;
   ierr = CCTK_GroupDynamicData (cctkGH, group, &groupdata);
   assert (!ierr);
   
   
   
-  for (n=0; n<nsurfaces; ++n) {
+  for (int n=0; n<nsurfaces; ++n) {
     
     if (!auto_res[n])
     {
       // set resolution according to given parameters
     
-      /* internal consistency checks */
+      // internal consistency checks
       assert (groupdata.dim == 2);
       assert (groupdata.gsh[0] >= ntheta[n]);
       assert (groupdata.gsh[1] >= nphi[n]);
@@ -69,7 +67,7 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       
       
       
-      /* copy parameters into grid functions */
+      // copy parameters into grid functions
       sf_ntheta[n] = ntheta[n];
       sf_nphi[n] = nphi[n];
       sf_nghoststheta[n] = nghoststheta[n];
@@ -80,8 +78,8 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       // we already have an approximate delta-spacing calculated in "SetupRes"
       // based on that, we calculate the number of necessary gridpoints...
       
-      double theta_range = pi;
-      double phi_range = 2*pi;
+      CCTK_REAL theta_range = pi;
+      CCTK_REAL phi_range = 2*pi;
       
       if (symmetric_x[n])
          phi_range /= 2;
@@ -92,14 +90,14 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       if (symmetric_z[n])
          theta_range /= 2;
       
-      sf_ntheta[n] = round(theta_range/sf_delta_theta_estimate[n]) + 2*nghoststheta[n];
-      sf_nphi[n]   = round(phi_range/sf_delta_phi_estimate[n]) + 2*nghostsphi[n];
+      sf_ntheta[n] = floor(theta_range/sf_delta_theta_estimate[n]+0.5) + 2*nghoststheta[n];
+      sf_nphi[n]   = floor(phi_range  /sf_delta_phi_estimate[n]  +0.5) + 2*nghostsphi[n];
       
       // ...make number of theta-gridpoints odd...
       if (sf_ntheta[n] % 2 != 1)
         sf_ntheta[n] += 1;
       
-      // ... and make inner phi-gridpoints devidable by 4 since some thorns require that (eg IsolatedHorizon)
+      // ... and make inner phi-gridpoints divisible by 4 since some thorns require that (e.g. IsolatedHorizon)
       sf_nphi[n]   += 4 - (sf_nphi[n]-2*nghostsphi[n]) % 4 + 2*nghostsphi[n];
       
       // consistency check
@@ -148,28 +146,28 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       
       if (verbose)
       {
-         printf("SphericalSurface[%d]: setting sf_ntheta = %d\n", n, sf_ntheta[n]);
-         printf("SphericalSurface[%d]: setting sf_nphi   = %d\n", n, sf_nphi[n]);
+        CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_ntheta = %d\n", n, sf_ntheta[n]);
+         CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_nphi   = %d\n", n, sf_nphi[n]);
       }
     }
     
     if (verbose)
     {
-       printf("SphericalSurface[%d]: finest refinement-level that fully contains this surface:  sf_maxreflevel = %d\n", n, sf_maxreflevel[n]);
-       printf("SphericalSurface[%d]: finest refinement-level that overlaps with this surface :  sf_minreflevel = %d\n", n, sf_minreflevel[n]);
+       CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that fully contains this surface:  sf_maxreflevel = %d\n", n, sf_maxreflevel[n]);
+       CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that overlaps with this surface :  sf_minreflevel = %d\n", n, sf_minreflevel[n]);
     }
     
-    /* coordinates in the theta direction */
-    /* avoid_sf_origin_theta = 1 */
+    // coordinates in the theta direction
+    // avoid_sf_origin_theta = 1
     if (symmetric_z[n]) {
       
-      /* upper hemisphere: z>=0, theta in (0, pi/2) */
+      // upper hemisphere: z>=0, theta in (0, pi/2)
       sf_delta_theta[n] = pi/2 / (sf_ntheta[n] - 2*nghoststheta[n] - 0.5);
       sf_origin_theta[n] = - (nghoststheta[n] - 0.5) * sf_delta_theta[n];
       
     } else {
       
-      /* both hemispheres: theta in (0, pi) */
+      // both hemispheres: theta in (0, pi)
       sf_delta_theta[n] = pi / (sf_ntheta[n] - 2*nghoststheta[n]);
       sf_origin_theta[n] = - (nghoststheta[n] - 0.5) * sf_delta_theta[n];
       
@@ -177,19 +175,19 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
     
     
     
-    /* coordinates in the phi direction */
-    /* avoid_sf_origin_phi = 0 */
+    // coordinates in the phi direction
+    // avoid_sf_origin_phi = 0
     if (symmetric_x[n]) {
       if (symmetric_y[n]) {
         
-        /* one quadrant: x>=0, y>=0, phi in [0, pi/2] */
+        // one quadrant: x>=0, y>=0, phi in [0, pi/2]
         assert (sf_nphi[n] - 2*nghostsphi[n] >= 1);
         sf_delta_phi[n] = pi/2 / (sf_nphi[n] - 2*nghostsphi[n] - 1);
         sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
         
       } else {
         
-        /* two quadrants: x>=0, phi in [-pi/2, pi/2] */
+        // two quadrants: x>=0, phi in [-pi/2, pi/2]
         assert (sf_nphi[n] - 2*nghostsphi[n] >= 2);
         sf_delta_phi[n] = pi / (sf_nphi[n] - 2*nghostsphi[n] - 1);
         sf_origin_phi[n] = - pi/2 - nghostsphi[n] * sf_delta_phi[n];
@@ -198,14 +196,14 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
     } else {
       if (symmetric_y[n]) {
         
-        /* two quadrants: y>=0, phi in [0, pi] */
+        // two quadrants: y>=0, phi in [0, pi]
         assert (sf_nphi[n] - 2*nghostsphi[n] >= 2);
         sf_delta_phi[n] = pi / (sf_nphi[n] - 2*nghostsphi[n] - 1);
         sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
         
       } else {
         
-        /* all quadrants: phi in [0, 2pi) */
+        // all quadrants: phi in [0, 2pi)
         assert (sf_nphi[n] - 2*nghostsphi[n] >= 4);
         sf_delta_phi[n] = 2*pi / (sf_nphi[n] - 2*nghostsphi[n]);
         sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
@@ -213,15 +211,15 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       }
     }
     
-#ifdef DEBUG
-      printf("SphericalSurface[%d]: sf_delta_theta = %.6f, sf_delta_phi = %.6f\n", n, sf_delta_theta[n], sf_delta_phi[n]);
-#endif
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: sf_delta_theta = %.6f, sf_delta_phi = %.6f\n", n, sf_delta_theta[n], sf_delta_phi[n]);
+    }
     
-    /* mark surface as uninitialised */
+    // mark surface as uninitialised
     sf_active[n] = 0;
     sf_valid[n] = 0;
     
-  } /* for n */
+  } // for n
 }
 
 
@@ -231,12 +229,9 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
-  static int first_call = 1;
-  int n;
-  
-  //printf("test\n");
+  static bool first_call = true;
   
-  for (n=0; n<nsurfaces; ++n) 
+  for (int n=0; n<nsurfaces; ++n) 
   {
     
     // try to find out with which Carpet::reflevel completely contains
@@ -264,9 +259,9 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
        my_radius = radius[n];
     }
       
-#ifdef DEBUG
-    printf("SphericalSurface: myradius = %f\n", my_radius);
-#endif
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: myradius = %f\n", double(my_radius));
+    }
     // get theta/phi delta-spacing by looking at Cartesian grid
       
     // first test, if we have an overlap with a certain reflevel
@@ -277,9 +272,9 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     l[0] = CCTK_GFINDEX3D(cctkGH, 0,             0,             0);
     l[1] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1);
       
-    int is_overlapping = 0;
-    double cart_radius = 0;    // radius of Cart. grid
-    double cart_origin[3];     // origin of Cart. grid
+    bool is_overlapping = false;
+    CCTK_REAL cart_radius = 0;    // radius of Cart. grid
+    CCTK_REAL cart_origin[3];     // origin of Cart. grid
       
     // calculate radius of Cartesian grid
     cart_origin[0] = (x[l[1]] - x[l[0]]) / 2.0;
@@ -288,9 +283,9 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
       
     cart_radius = sqrt(cart_origin[0]*cart_origin[0] + cart_origin[1]*cart_origin[1] + cart_origin[2]*cart_origin[2]);
       
-#ifdef DEBUG
-    printf("SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f\n", cart_radius, cart_origin[0], cart_origin[1], cart_origin[2]);
-#endif
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f\n", double(cart_radius), double(cart_origin[0]), double(cart_origin[1]), double(cart_origin[2]));
+    }
 
     // set origin of Cartesian grid
     cart_origin[0] += x[l[0]];
@@ -298,7 +293,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     cart_origin[2] += z[l[0]];
       
     // norm of displacement vector between origin of Cart. grid and origin of SphericalSurface
-    double dist = sqrt( (cart_origin[0]-origin_x[n]) * (cart_origin[0]-origin_x[n])
+    CCTK_REAL dist = sqrt( (cart_origin[0]-origin_x[n]) * (cart_origin[0]-origin_x[n])
                       + (cart_origin[1]-origin_y[n]) * (cart_origin[1]-origin_y[n])
                       + (cart_origin[2]-origin_z[n]) * (cart_origin[2]-origin_z[n]));
       
@@ -306,7 +301,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     if (my_radius+cart_radius >= dist &&  // Cart. grid and SphericalSurface overlap
         cart_radius+dist > my_radius  &&  // Cart. grid patch is not fully contained in this SphericalSurface
         -cart_radius+dist < my_radius)    // Cart. grid patch is not outside of the spherical shell
-       is_overlapping = 1;
+       is_overlapping = true;
          
     // get resolution
     if (is_overlapping)
@@ -323,7 +318,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
        }
        
        // store reflevel
-       sf_minreflevel[n] = Carpet::reflevel;
+       sf_minreflevel[n] = get_reflevel(cctkGH);
     }
     else if (first_call)
     {
@@ -336,7 +331,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     if (CCTK_Equals(auto_res_grid[n], "overlap"))
     {
       // find minimum value of "sf_delta_phi/theta" over all processors
-      double dummy = 0;
+      CCTK_REAL dummy = 0;
       MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
       sf_delta_theta_estimate[n] = dummy;
          
@@ -364,12 +359,12 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
                                               + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1)
                                               + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ) ); 
 
-#ifdef DEBUG						 
-    /*printf("SphericalSurface: calculated max_radius = %.12f\n", max_radius);
-    printf("SphericalSurface: buffer_radius = %.12f\n", sqrt(  np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
+    if (verbose) {
+    //CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: calculated max_radius = %.12f\n", max_radius);
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: buffer_radius = %.12f\n", double (sqrt(  np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                              + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1)
-                                                             + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ));*/
-#endif
+                                                                                                + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) )));
+    }
 
    
     unsigned long k[8];
@@ -384,11 +379,11 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     k[7] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1);
    
     int i;
-    int is_contained = 0;
+    bool is_contained = false;
     for (i=0; i < 8; i++)
     {
        if (max_radius < sqrt(x[k[i]]*x[k[i]] + y[k[i]]*y[k[i]] + z[k[i]]*z[k[i]]))
-          is_contained = 1;
+          is_contained = true;
     }
    
     int max_rl = 0;
@@ -405,7 +400,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
                                                                         + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
        }
     
-       max_rl = Carpet::reflevel; //(int) floor(log((double) (cctk_levfac[0])/log(2.0))+0.5)+1; 
+       max_rl = get_reflevel(cctkGH);
     }
     else if (first_call)
     {
@@ -429,7 +424,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     if (CCTK_Equals(auto_res_grid[n], "fully contained"))
     {
       // find maximum value of "sf_delta_phi/theta" over all processors
-      double dummy = 0;
+      CCTK_REAL dummy = 0;
       MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
       sf_delta_theta_estimate[n] = dummy;
          
@@ -438,14 +433,13 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
       sf_delta_phi_estimate[n] = dummy;
     }
     
-#ifdef DEBUG
-    printf("SphericalSurface: sf_maxreflevel[%d] = %d,  sf_minreflevel[%d] = %d \n", n, sf_maxreflevel[n], n, sf_minreflevel[n]);
-    printf("SphericalSurface: sf_delta_theta[%d] = %.6f, sf_delta_phi[%d] = %.6f,  \n", n, sf_delta_theta_estimate[n], n, sf_delta_phi_estimate[n]);
-#endif
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_maxreflevel[%d] = %d,  sf_minreflevel[%d] = %d \n", n, sf_maxreflevel[n], n, sf_minreflevel[n]);
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_delta_theta[%d] = %.6f, sf_delta_phi[%d] = %.6f,  \n", n, sf_delta_theta_estimate[n], n, sf_delta_phi_estimate[n]);
+    }
 
   }
   
-  first_call = 0;
+  first_call = false;
   
 }
-