- fixed "surface not valid" bug

- replaced MPI calls by ReduceLocScalar calls


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

2 files changed, 65 insertions, 29 deletions

diff --git a/schedule.ccl b/schedule.ccl
index 53827c4..89aec32 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -16,11 +16,13 @@ STORAGE: sf_maxreflevel
 SCHEDULE SphericalSurface_SetupRes AT basegrid BEFORE SphericalSurface_Setup AFTER SpatialCoordinates AFTER CorrectCoordinates
 {
    LANG: C
+     OPTIONS: GLOBAL LOOP-LOCAL
 } "Set surface resolution automatically"
 
 SCHEDULE SphericalSurface_Setup AT basegrid
 {
   LANG: C
+    OPTIONS: GLOBAL
 } "Calculate surface coordinate descriptors"
 
 
diff --git a/src/setup.cc b/src/setup.cc
index 897efb3..fa22c5d 100644
--- a/src/setup.cc
+++ b/src/setup.cc
@@ -164,15 +164,15 @@ extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
       
       if (verbose)
       {
-        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]);
+        CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_ntheta = %d", n, sf_ntheta[n]);
+        CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: setting sf_nphi   = %d", n, sf_nphi[n]);
       }
     }
     
     if (verbose)
     {
-       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]);
+       CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that fully contains this surface:  sf_maxreflevel = %d", n, sf_maxreflevel[n]);
+       CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface[%d]: finest refinement-level that overlaps with this surface :  sf_minreflevel = %d", n, sf_minreflevel[n]);
     }
     
     // coordinates in the theta direction
@@ -249,6 +249,14 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
   
   static bool first_call = true;
   
+  int min_reduction_handle = CCTK_ReductionArrayHandle("minimum");
+  int max_reduction_handle = CCTK_ReductionArrayHandle("maximum");
+
+  if (min_reduction_handle < 0)
+     CCTK_WARN(0, "Cannot get reduction handle for minimum operation.");
+  if (max_reduction_handle < 0)
+     CCTK_WARN(0, "Cannot get reduction handle for maximum operation.");
+  
   for (int n=0; n<nsurfaces; ++n) 
   {
     
@@ -278,7 +286,7 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     }
       
     if (verbose) {
-      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: myradius = %f\n", double(my_radius));
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: myradius = %f", double(my_radius));
     }
     // get theta/phi delta-spacing by looking at Cartesian grid
       
@@ -293,27 +301,29 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     bool is_overlapping = false;
     CCTK_REAL cart_radius = 0;    // radius of Cart. grid
     CCTK_REAL cart_origin[3];     // origin of Cart. grid
-      
+    
+    //CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0]) * CCTK_DELTA_SPACE(0) - (CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0] + cctk_lsh[0]) * CCTK_DELTA_SPACE(0))
+    
     // calculate radius of Cartesian grid
-    cart_origin[0] = (x[l[1]] - x[l[0]]) / 2.0;
-    cart_origin[1] = (y[l[1]] - y[l[0]]) / 2.0;
-    cart_origin[2] = (z[l[1]] - z[l[0]]) / 2.0;
+    cart_origin[0] = 0.5 * cctk_lsh[0] * CCTK_DELTA_SPACE(0); //(x[l[1]] - x[l[0]]) / 2.0;
+    cart_origin[1] = 0.5 * cctk_lsh[1] * CCTK_DELTA_SPACE(1); //(y[l[1]] - y[l[0]]) / 2.0;
+    cart_origin[2] = 0.5 * cctk_lsh[2] * CCTK_DELTA_SPACE(2); //(z[l[1]] - z[l[0]]) / 2.0;
       
     cart_radius = sqrt(cart_origin[0]*cart_origin[0] + cart_origin[1]*cart_origin[1] + cart_origin[2]*cart_origin[2]);
       
     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]));
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f", 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]];
-    cart_origin[1] += y[l[0]];
-    cart_origin[2] += z[l[0]];
+    cart_origin[0] += CCTK_ORIGIN_SPACE(0) + cctk_lbnd[0] * CCTK_DELTA_SPACE(0); //x[l[0]];
+    cart_origin[1] += CCTK_ORIGIN_SPACE(1) + cctk_lbnd[1] * CCTK_DELTA_SPACE(1); //y[l[0]];
+    cart_origin[2] += CCTK_ORIGIN_SPACE(2) + cctk_lbnd[2] * CCTK_DELTA_SPACE(2); //z[l[0]];
       
     // norm of displacement vector between origin of Cart. grid and origin of SphericalSurface
     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]));
+                         + (cart_origin[1]-origin_y[n]) * (cart_origin[1]-origin_y[n])
+                         + (cart_origin[2]-origin_z[n]) * (cart_origin[2]-origin_z[n]));
       
     // does our spherical shell overlap with the processor's local grid patch? 
     if (my_radius+cart_radius >= dist &&  // Cart. grid and SphericalSurface overlap
@@ -350,11 +360,13 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     {
       // find minimum value of "sf_delta_phi/theta" over all processors
       CCTK_REAL dummy = 0;
-      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      //MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &sf_delta_theta_estimate[n], &dummy, CCTK_VARIABLE_REAL);
       sf_delta_theta_estimate[n] = dummy;
          
       dummy = 0;
-      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      //MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &sf_delta_phi_estimate[n], &dummy, CCTK_VARIABLE_REAL);
       sf_delta_phi_estimate[n] = dummy;
       
       /*
@@ -379,7 +391,8 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     
     // find maximum value of "sf_minreflevel" over all processors
     int idummy;
-    MPI_Allreduce(&sf_minreflevel[n], &idummy, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+    //MPI_Allreduce(&sf_minreflevel[n], &idummy, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+    CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_minreflevel[n], &idummy, CCTK_VARIABLE_INT);
     sf_minreflevel[n] = idummy;
 
 
@@ -398,12 +411,12 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
 
     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) )));
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: buffer_radius = %.12f", 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) )));
     }
 
-   
+   /*
     unsigned long k[8];
       
     k[0] = CCTK_GFINDEX3D(cctkGH, 0,             0,             0);
@@ -414,12 +427,30 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     k[5] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, 0,             cctk_lsh[2]-1);
     k[6] = CCTK_GFINDEX3D(cctkGH, 0,             cctk_lsh[1]-1, cctk_lsh[2]-1);
     k[7] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1);
-   
+    */
+    
+    unsigned long kx[8], ky[8], kz[8];
+      
+    kx[0] = 0;             ky[0] = 0;             kz[0] = 0;
+    kx[1] = cctk_lsh[0]-1, ky[1] = 0;             kz[1] = 0;
+    kx[2] = 0;             ky[2] = cctk_lsh[1]-1; kz[2] = 0;
+    kx[3] = cctk_lsh[0]-1; ky[3] = cctk_lsh[1]-1; kz[3] = 0;
+    kx[4] = 0;             ky[4] = 0;             kz[4] = cctk_lsh[2]-1;
+    kx[5] = cctk_lsh[0]-1; ky[5] = 0;             kz[5] = cctk_lsh[2]-1;
+    kx[6] = 0;             ky[6] = cctk_lsh[1]-1; kz[6] = cctk_lsh[2]-1;
+    kx[7] = cctk_lsh[0]-1; ky[7] = cctk_lsh[1]-1; kz[7] = cctk_lsh[2]-1;
+    
+    
     int i;
     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]]))
+       //if (max_radius < sqrt(x[k[i]]*x[k[i]] + y[k[i]]*y[k[i]] + z[k[i]]*z[k[i]]))
+       //   is_contained = true;
+       CCTK_REAL xx = CCTK_ORIGIN_SPACE(0) + (cctk_lbnd[0] + kx[i]) * CCTK_DELTA_SPACE(0);
+       CCTK_REAL yy = CCTK_ORIGIN_SPACE(1) + (cctk_lbnd[1] + ky[i]) * CCTK_DELTA_SPACE(1);
+       CCTK_REAL zz = CCTK_ORIGIN_SPACE(2) + (cctk_lbnd[2] + kz[i]) * CCTK_DELTA_SPACE(2);
+       if (max_radius < sqrt(xx*xx + yy*yy + zz*zz))
           is_contained = true;
     }
    
@@ -448,7 +479,8 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     
     
     // find minimum value of "sf_maxreflevel" over all processors
-    MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+    //MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+    CCTK_ReduceLocScalar(cctkGH, -1, min_reduction_handle, &max_rl, &idummy, CCTK_VARIABLE_INT);
     if (idummy != 0 || first_call)
 	sf_maxreflevel[n] = idummy;
     
@@ -462,17 +494,19 @@ extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
     {
       // find maximum value of "sf_delta_phi/theta" over all processors
       CCTK_REAL dummy = 0;
-      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      //MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_delta_theta_estimate[n], &dummy, CCTK_VARIABLE_REAL);
       sf_delta_theta_estimate[n] = dummy;
          
       dummy = 0;
-      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      //MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      CCTK_ReduceLocScalar(cctkGH, -1, max_reduction_handle, &sf_delta_phi_estimate[n], &dummy, CCTK_VARIABLE_REAL);
       sf_delta_phi_estimate[n] = dummy;
     }
     
     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]);
+      CCTK_VInfo (CCTK_THORNSTRING, "SphericalSurface: sf_maxreflevel[%d] = %d,  sf_minreflevel[%d] = %d", 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]);
     }
 
   }