new setup.cc

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

1 files changed, 451 insertions, 0 deletions

diff --git a/src/setup.cc b/src/setup.cc
new file mode 100644
index 0000000..f7aa0be
--- /dev/null
+++ b/src/setup.cc
@@ -0,0 +1,451 @@
+/* $Header$ */
+
+#include <assert.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
+
+
+#define round(x) (x<0?ceil((x)-0.5):floor((x)+0.5))
+
+#define INITIAL_VAL 10000
+
+
+//TODO:  overlap with finer levels still seems to be incorrect!
+
+
+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;
+  
+  int n;
+  int ierr;
+  
+  
+  
+  if (nsurfaces == 0) return;
+  
+  
+  
+  group = CCTK_GroupIndex ("SphericalSurface::sf_radius");
+  assert (group>=0);
+  ierr = CCTK_GroupData (group, &groupinfo);
+  assert (!ierr);
+  ierr = CCTK_GroupDynamicData (cctkGH, group, &groupdata);
+  assert (!ierr);
+  
+  
+  
+  for (n=0; n<nsurfaces; ++n) {
+    
+    if (!auto_res[n])
+    {
+      // set resolution according to given parameters
+    
+      /* internal consistency checks */
+      assert (groupdata.dim == 2);
+      assert (groupdata.gsh[0] >= ntheta[n]);
+      assert (groupdata.gsh[1] >= nphi[n]);
+      
+      assert (ntheta[n] >= 3*nghoststheta[n] && ntheta[n] <= maxntheta);
+      assert (nphi[n] >= 3*nghostsphi[n] && nphi[n] <= maxnphi);
+      
+      
+      
+      /* copy parameters into grid functions */
+      sf_ntheta[n] = ntheta[n];
+      sf_nphi[n] = nphi[n];
+      sf_nghoststheta[n] = nghoststheta[n];
+      sf_nghostsphi[n] = nghostsphi[n];
+    }
+    else
+    {
+      // 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;
+      
+      if (symmetric_x[n])
+         phi_range /= 2;
+         
+      if (symmetric_y[n])
+         phi_range /= 2;
+      
+      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];
+      
+      // ...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)
+      sf_nphi[n]   += 4 - (sf_nphi[n]-2*nghostsphi[n]) % 4 + 2*nghostsphi[n];
+      
+      // consistency check
+      assert (sf_ntheta[n] >= 3*nghoststheta[n] && sf_ntheta[n] <= maxntheta);
+      assert (sf_nphi[n] >= 3*nghostsphi[n] && sf_nphi[n] <= maxnphi);
+      
+      // check, if there are enough maxntheta,maxnphi gridpoints to remove symmetries
+      // (some thorns require that)
+      if (symmetric_z[n])
+      {
+         if (2*(sf_ntheta[n]-2*nghoststheta[n])+2*nghoststheta[n] > maxntheta)
+	    CCTK_WARN(0, "maxntheta is not big enough to remove z-symmetry (some thorns may require that)!");
+      }
+      
+      if (symmetric_x[n])
+      {
+         if (symmetric_y[n])
+         {
+            if (4*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
+            {
+               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
+            }
+         }
+         else
+         {
+            if (2*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
+            {
+               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
+            }
+         }
+      }
+      else
+      {
+         if (symmetric_y[n])
+         {
+            if (2*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
+            {
+               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
+            }
+         }
+      }
+
+      
+      sf_nghoststheta[n] = nghoststheta[n];
+      sf_nghostsphi[n] = nghostsphi[n];
+      
+      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]);
+      }
+    }
+    
+    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]);
+    }
+    
+    /* coordinates in the theta direction */
+    /* avoid_sf_origin_theta = 1 */
+    if (symmetric_z[n]) {
+      
+      /* 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) */
+      sf_delta_theta[n] = pi / (sf_ntheta[n] - 2*nghoststheta[n]);
+      sf_origin_theta[n] = - (nghoststheta[n] - 0.5) * sf_delta_theta[n];
+      
+    }
+    
+    
+    
+    /* 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] */
+        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] */
+        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];
+        
+      }
+    } else {
+      if (symmetric_y[n]) {
+        
+        /* 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) */
+        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];
+        
+      }
+    }
+    
+#ifdef DEBUG
+      printf("SphericalSurface[%d]: sf_delta_theta = %.6f, sf_delta_phi = %.6f\n", n, sf_delta_theta[n], sf_delta_phi[n]);
+#endif
+    
+    /* mark surface as uninitialised */
+    sf_active[n] = 0;
+    sf_valid[n] = 0;
+    
+  } /* for n */
+}
+
+
+
+extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  static int first_call = 1;
+  int n;
+  
+  //printf("test\n");
+  
+  for (n=0; n<nsurfaces; ++n) 
+  {
+    
+    // try to find out with which Carpet::reflevel completely contains
+    // the current surface and will intersect and and store the
+    // result so that other thorns such as the kick-thorn
+    // can decide which reflevel to use
+     
+    // set resolution according to Radius and Cartesian resolution
+      
+    CCTK_REAL my_radius;
+      
+    if (!set_spherical[n] && !set_elliptic[n])
+    {
+       // we have a problem...without a radius we cannot estimate a proper resolution
+       // Here, we simply assume r=1 if radius[n] == 0.
+       // If you want a better estimate of the resolution, set the radius parameter to an initial guess!
+       
+       if (radius[n] == 0)
+          my_radius = 1;
+       else
+	  my_radius = radius[n];
+    }
+    else
+    {
+       my_radius = radius[n];
+    }
+      
+#ifdef DEBUG
+    printf("SphericalSurface: myradius = %f\n", my_radius);
+#endif
+    // get theta/phi delta-spacing by looking at Cartesian grid
+      
+    // first test, if we have an overlap with a certain reflevel
+    // we do this by calculating a sphere based on the corners of the Cartesian local patch
+    // on this processor.
+    unsigned long l[2];
+      
+    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
+      
+    // 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_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
+
+    // set origin of Cartesian grid
+    cart_origin[0] += x[l[0]];
+    cart_origin[1] += y[l[0]];
+    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])
+                      + (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
+        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;
+         
+    // get resolution
+    if (is_overlapping)
+    {
+       if (CCTK_Equals(auto_res_grid[n], "overlap"))
+       {
+          sf_delta_theta_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
+                                                                          + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
+                                                                          + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
+         
+          sf_delta_phi_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
+                                                                        + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
+                                                                        + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
+       }
+       
+       // store reflevel
+       sf_minreflevel[n] = Carpet::reflevel;
+    }
+    else if (first_call)
+    {
+       // set to very high value initially (otherwise would be zero
+       // and reported as minimum value in following MPI communication)
+       sf_delta_theta_estimate[n] = 10000;
+       sf_delta_phi_estimate[n] = 10000; 
+    }
+    
+    if (CCTK_Equals(auto_res_grid[n], "overlap"))
+    {
+      // find minimum value of "sf_delta_phi/theta" over all processors
+      double dummy = 0;
+      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      sf_delta_theta_estimate[n] = dummy;
+         
+      dummy = 0;
+      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+      sf_delta_phi_estimate[n] = dummy;
+    }
+    
+    // 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);
+    sf_minreflevel[n] = idummy;
+
+
+
+    // find reflevel that completely contains this surface
+    // FIXME: This algorithm does not take into account multiple
+    //        disconnected refinement regions.
+    
+    
+    //safety margin of 6 points around sphere
+    CCTK_INT np = 6;
+   
+    CCTK_REAL max_radius = my_radius + ( 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) ) ); 
+
+#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) 
+                                                             + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1)
+                                                             + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ));*/
+#endif
+
+   
+    unsigned long k[8];
+      
+    k[0] = CCTK_GFINDEX3D(cctkGH, 0,             0,             0);
+    k[1] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, 0,             0);
+    k[2] = CCTK_GFINDEX3D(cctkGH, 0,             cctk_lsh[1]-1, 0);
+    k[3] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, 0);
+    k[4] = CCTK_GFINDEX3D(cctkGH, 0,             0,             cctk_lsh[2]-1);
+    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);
+   
+    int i;
+    int is_contained = 0;
+    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;
+    }
+   
+    int max_rl = 0;
+    if (is_contained)
+    {
+       if (CCTK_Equals(auto_res_grid[n], "fully contained"))
+       {
+          sf_delta_theta_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
+                                                                          + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
+                                                                          + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
+         
+          sf_delta_phi_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
+                                                                        + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
+                                                                        + 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; 
+    }
+    else if (first_call)
+    {
+       // set to zero initially
+       sf_delta_theta_estimate[n] = 0;
+       sf_delta_phi_estimate[n] = 0; 
+    }
+    
+    
+    // find minimum value of "sf_maxreflevel" over all processors
+    MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+    if (idummy != 0 || first_call)
+	sf_maxreflevel[n] = idummy;
+    
+    // this condition can become true because of this approximate algorithm
+    // which does not take into account disconnected ref. regions
+    if (sf_maxreflevel[n] > sf_minreflevel[n])
+        sf_maxreflevel[n] = sf_minreflevel[n];   // can still be incorrect!
+	
+    
+    if (CCTK_Equals(auto_res_grid[n], "fully contained"))
+    {
+      // find maximum value of "sf_delta_phi/theta" over all processors
+      double dummy = 0;
+      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      sf_delta_theta_estimate[n] = dummy;
+         
+      dummy = 0;
+      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      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
+
+  }
+  
+  first_call = 0;
+  
+}
+