CarpetAdaptiveRegrid: Allow non-trivial initial grids

Use the CarpetRegrid::manual-coordinate-list notation to allow the
setting up of initially refined grids. Literal cut 'n paste in many places.

darcs-hash:20050223235501-0ff1f-4be0e5c479aadabc7e1c5682476f95a23c00ec71.gz

5 files changed, 296 insertions, 8 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/param.ccl b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
index 9660d454c..0ce964b32 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/param.ccl
+++ b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
@@ -42,3 +42,14 @@ CCTK_BOOLEAN verbose "Should we be verbose"
 CCTK_BOOLEAN veryverbose "Should we be very verbose"
 {
 } "no"
+
+CCTK_INT refinement_levels "Number of initial refinement levels (including the base level)" STEERABLE=always
+{
+  1:* :: "must be positive, and must not be larger than Carpet::max_refinement_levels"
+} 1
+
+CCTK_STRING coordinates "List of bounding box coordinates for the initial data" STEERABLE=always
+{
+  "^$" :: "leave empty for no refinement"
+  ".*" :: "[ [ ([<xmin>,<ymin>,<zmin>]:[<xmax>,<ymax>,<zmax>]:[<xstride>,<ystride>,<zstride>]), ... ], ... ]"
+} ""
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index 80ccb1e9e..331ca504b 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -24,13 +24,14 @@ extern "C" {
 
 
 namespace CarpetAdaptiveRegrid {
+
+
+  static gh::mexts local_bbsss;
+  static CCTK_INT last_iteration = -1;  
   
   using namespace std;
   using namespace Carpet;
 
-  static gh::mexts local_bbsss;
-  static CCTK_INT last_iteration = -1;
-
   extern "C" {
     void CCTK_FCALL CCTK_FNAME(copy_mask)
       (const CCTK_INT& snx, const CCTK_INT& sny, const CCTK_INT& snz,
@@ -75,6 +76,18 @@ namespace CarpetAdaptiveRegrid {
       // Is this really the right thing to do on 
       // multiprocessors?
       local_bbsss = bbsss;
+      last_iteration = cctkGH->cctk_iteration;
+      CCTK_INT do_recompose = 
+        ManualCoordinateList (cctkGH, hh, bbsss, obss, pss, local_bbsss);
+
+      if (verbose) {
+        ostringstream buf;
+        buf << "Done with manual coordinate list. Total list is:"
+            << endl << local_bbsss;
+        CCTK_INFO(buf.str().c_str());
+      }      
+
+      return do_recompose;
     }
 
     // FIXME: We should check that the local reflevel "agrees"
@@ -103,9 +116,11 @@ namespace CarpetAdaptiveRegrid {
       }
 
       // Return if it's initial data as we can't handle that yet.
-      if (cctkGH->cctk_iteration ==0) {
-        return 0;
-      }
+      // Actually don't as initial data should now be handled
+      // by the manualcoordinatelist above.
+      //       if (cctkGH->cctk_iteration == 0) {
+      //         return 0;
+      //       }
       
     }
 
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
index d5371068c..ca91b9511 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
@@ -23,7 +23,6 @@ namespace CarpetAdaptiveRegrid {
   using namespace Carpet;
   
   
-  
   extern "C" {
     
     /* Scheduled functions */
@@ -41,6 +40,39 @@ namespace CarpetAdaptiveRegrid {
 				  CCTK_INT force);
   }
 
+  int ManualCoordinateList (cGH const * const cctkGH,
+                            gh const & hh,
+                            gh::mexts  & bbsss,
+                            gh::rbnds  & obss,
+                            gh::rprocs & pss,
+                            gh::mexts  & local_bbsss);
+    
+  void ManualCoordinates_OneLevel (const cGH * const cctkGH,
+                                   const gh & hh,
+                                   const int rl,
+                                   const int numrl,
+                                   const rvect lower,
+                                   const rvect upper,
+                                   const bbvect obound,
+                                   vector<ibbox> & bbs,
+                                   vector<bbvect> & obs);
+  
+  void ManualGridpoints_OneLevel (const cGH * const cctkGH,
+                                  const gh & hh,
+                                  const int rl,
+                                  const int numrl,
+                                  const ivect ilower,
+                                  const ivect iupper,
+                                  const bbvect obound,
+                                  vector<ibbox> & bbs,
+                                  vector<bbvect> & obs);
+
+  rvect int2pos (const cGH* const cctkGH, const gh& hh,
+                 const ivect & ipos, const int rl);
+
+  ivect pos2int (const cGH* const cctkGH, const gh& hh,
+                 const rvect & rpos, const int rl);
+
 } // namespace CarpetAdaptiveregrid
 
 #endif // !defined(CARPETADAPTIVEREGRID_HH)
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
index 8f36e4674..a818be026 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
@@ -4,7 +4,8 @@
 # Source files in this directory
 SRCS = 	CAR_Paramcheck.cc \
 	CAR.cc \
-	CAR_utils.F90
+	CAR_utils.F90 \
+	manualcoordinatelist.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
new file mode 100644
index 000000000..20e5dfc0e
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
@@ -0,0 +1,229 @@
+#include <cassert>
+#include <cmath>
+#include <cstring>
+#include <sstream>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "gh.hh"
+
+#include "carpet.hh"
+#include "CAR.hh"
+
+
+
+namespace CarpetAdaptiveRegrid {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  int ManualCoordinateList (cGH const * const cctkGH,
+                            gh const & hh,
+                            gh::mexts  & bbsss,
+                            gh::rbnds  & obss,
+                            gh::rprocs & pss,
+                            gh::mexts  & local_bbsss)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    int ierr;
+    
+    assert (refinement_levels >= 1);
+    
+    // do nothing if the levels already exist
+    if (reflevel == refinement_levels) return 0;
+    
+    assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
+    
+    jjvect nboundaryzones, is_internal, is_staggered, shiftout;
+    ierr = GetBoundarySpecification
+      (2*dim, &nboundaryzones[0][0], &is_internal[0][0],
+       &is_staggered[0][0], &shiftout[0][0]);
+    assert (!ierr);
+    rvect physical_min, physical_max;
+    rvect interior_min, interior_max;
+    rvect exterior_min, exterior_max;
+    rvect base_spacing;
+    ierr = GetDomainSpecification
+      (dim, &physical_min[0], &physical_max[0],
+       &interior_min[0], &interior_max[0],
+       &exterior_min[0], &exterior_max[0], &base_spacing[0]);
+    assert (!ierr);
+    
+    bbss.resize (refinement_levels);
+    obss.resize (refinement_levels);
+    pss.resize (refinement_levels);
+    
+    vector<vector<rbbox> > newbbss;
+    if (strcmp(coordinates, "") != 0) {
+      istringstream gp_str(coordinates);
+      try {
+        gp_str >> newbbss;
+      } catch (input_error) {
+        CCTK_WARN (0, "Could not parse parameter \"coordinates\"");
+      }
+    }
+
+    vector<vector<bbvect> > newobss;
+
+    newobss.resize(newbbss.size());
+    for (size_t rl=0; rl<newobss.size(); ++rl) {
+      newobss.at(rl).resize(newbbss.at(rl).size());
+      for (size_t c=0; c<newobss.at(rl).size(); ++c) {
+        for (int d=0; d<dim; ++d) {
+          assert (mglevel==0);
+          rvect const spacing = base_spacing * ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, rl+1);
+          ierr = ConvertFromPhysicalBoundary
+            (dim, &physical_min[0], &physical_max[0],
+             &interior_min[0], &interior_max[0],
+             &exterior_min[0], &exterior_max[0], &spacing[0]);
+          assert (!ierr);
+          newobss.at(rl).at(c)[d][0] = abs(newbbss.at(rl).at(c).lower()[d] - physical_min[d]) < 1.0e-6 * spacing[d];
+          if (newobss.at(rl).at(c)[d][0]) {
+            rvect lo = newbbss.at(rl).at(c).lower();
+            rvect up = newbbss.at(rl).at(c).upper();
+            rvect str = newbbss.at(rl).at(c).stride();
+            lo[d] = exterior_min[d];
+            newbbss.at(rl).at(c) = rbbox(lo, up, str);
+          }
+          newobss.at(rl).at(c)[d][1] = abs(newbbss.at(rl).at(c).upper()[d] - physical_max[d]) < 1.0e-6 * base_spacing[d] / ipow(reffact, rl);
+          if (newobss.at(rl).at(c)[d][1]) {
+            rvect lo = newbbss.at(rl).at(c).lower();
+            rvect up = newbbss.at(rl).at(c).upper();
+            rvect str = newbbss.at(rl).at(c).stride();
+            up[d] = exterior_max[d];
+            newbbss.at(rl).at(c) = rbbox(lo, up, str);
+          }
+        }
+      }
+    }
+        
+    if (newbbss.size() < refinement_levels-1) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The parameter \"coordinates\" must contain at least \"refinement_levels-1\" (here: %d) levels", int(refinement_levels-1));
+    }
+    
+    for (size_t rl=1; rl<refinement_levels; ++rl) {
+      
+      vector<ibbox> bbs;
+      gh::cbnds obs;
+      
+      bbs.reserve (newbbss.at(rl-1).size());
+      obs.reserve (newbbss.at(rl-1).size());
+      
+      for (size_t c=0; c<newbbss.at(rl-1).size(); ++c) {
+        rbbox const & ext = newbbss.at(rl-1).at(c);
+        bbvect const & ob = newobss.at(rl-1).at(c);
+        // TODO: why can basemglevel not be used here?
+        // rvect const spacing = base_spacing * ipow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl);
+        rvect const spacing = base_spacing / ipow(reffact, rl);
+        if (! all(abs(ext.stride() - spacing) < spacing * 1.0e-10)) {
+          assert (dim==3);
+          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "The grid spacing on refinement level %d is incorrect.  I expected [%g,%g,%g], but I found [%g,%g,%g].",
+                      int(rl),
+                      double(spacing[0]), double(spacing[1]), double(spacing[2]),
+                      double(ext.stride()[0]), double(ext.stride()[1]), double(ext.stride()[2]));
+        }
+        assert (all(abs(ext.stride() - spacing) < spacing * 1.0e-10));
+        
+        ManualCoordinates_OneLevel
+          (cctkGH, hh, rl, refinement_levels,
+           ext.lower(), ext.upper(), ob, bbs, obs);
+      }
+
+      bbss.at(rl) = bbs;
+      obss.at(rl) = obs;
+      
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+
+      if (verbose) {
+        ostringstream buf;
+        buf << "Doing manual coordinate list, level " << rl 
+            << ". Total list is:"
+            << endl << local_bbsss;
+        CCTK_INFO(buf.str().c_str());
+      }      
+      
+      // make multiprocessor aware
+      gh::cprocs ps;
+      SplitRegions (cctkGH, bbs, obs, ps);
+
+      bbss.at(rl) = bbs;
+      obss.at(rl) = obs;
+      pss.at(rl) = ps;
+              
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+
+    } // for rl
+    
+    return 1;
+  }
+  
+  void ManualCoordinates_OneLevel (const cGH * const cctkGH,
+                                   const gh & hh,
+                                   const int rl,
+                                   const int numrl,
+                                   const rvect lower,
+                                   const rvect upper,
+                                   const bbvect obound,
+                                   vector<ibbox> & bbs,
+                                   vector<bbvect> & obs)
+  {
+    if (rl >= numrl) return;
+    
+    jvect const ilower = pos2int (cctkGH, hh, lower, rl);
+    jvect const iupper = pos2int (cctkGH, hh, upper, rl);
+    
+    ManualGridpoints_OneLevel
+      (cctkGH, hh, rl, numrl, ilower, iupper, obound, bbs, obs);
+  }
+
+  void ManualGridpoints_OneLevel (const cGH * const cctkGH,
+                                  const gh & hh,
+                                  const int rl,
+                                  const int numrl,
+                                  const ivect ilower,
+                                  const ivect iupper,
+                                  const bbvect obound,
+                                  vector<ibbox> & bbs,
+                                  vector<bbvect> & obs)
+  {
+    const ivect rstr = hh.baseextent.stride();
+    const ivect rlb  = hh.baseextent.lower();
+    const ivect rub  = hh.baseextent.upper();
+    
+    const int levfac = ipow(hh.reffact, rl);
+    assert (all (rstr % levfac == 0));
+    const ivect str (rstr / levfac);
+    const ivect lb  (ilower);
+    const ivect ub  (iupper);
+    if (! all(lb>=rlb && ub<=rub)) {
+      ostringstream buf;
+      buf << "The refinement region boundaries for refinement level #" << rl << " are not within the main grid.  Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
+      CCTK_WARN (0, buf.str().c_str());
+    }
+    if (! all(lb<=ub)) {
+      ostringstream buf;
+      buf << "The refinement region boundaries for refinement level #" << rl << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends;
+      CCTK_WARN (0, buf.str().c_str());
+    }
+    if (! all(lb%str==0 && ub%str==0)) {
+      CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
+    }
+    assert (all(lb>=rlb && ub<=rub));
+    assert (all(lb<=ub));
+    assert (all(lb%str==0 && ub%str==0));
+    
+    bbs.push_back (ibbox(lb, ub, str));
+    obs.push_back (obound);
+  }
+  
+} // namespace CarpetRegrid