AMR boundaries

Add support for using boundaries with AMR. Seems to work with symmetry
boundaries but gets a bit confused if you try refining outer
boundaries. Requires that you setup the domain using CoordBase.

darcs-hash:20050211135704-58c7f-6ed6ce1449de17599dd8b23abc909f5ba313d048.gz

4 files changed, 257 insertions, 12 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par
index 359c35d88..588d5de82 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par
+++ b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par
@@ -13,15 +13,25 @@ mol::initial_data_is_crap = "yes"
 ##carpet::veryverbose = "yes"
 
 grid::domain = "full"
-grid::type = "byspacing"
+grid::type = "coordbase"
 grid::avoid_origin = "no"
-driver::global_nx = 51
-driver::global_ny = 51
-driver::global_nz = 51
-grid::dxyz = 0.02
+
+coordbase::xmin = -0.5
+coordbase::ymin = -0.5
+coordbase::zmin = -0.5
+coordbase::xmax =  0.5
+coordbase::ymax =  0.5
+coordbase::zmax =  0.5
+
+coordbase::dx = 0.02
+coordbase::dy = 0.02
+coordbase::dz = 0.02
+
 driver::ghost_size = 1
 time::dtfac = 0.5
 
+carpet::domain_from_coordbase = yes
+
 carpet::max_refinement_levels = 2
 carpet::init_3_timelevels = yes
 carpetadaptiveregrid::max_error = 1.e-6
@@ -33,6 +43,8 @@ carpetadaptiveregrid::regrid_every = 1
 
 #carpetadaptiveregrid::verbose = yes
 #carpetadaptiveregrid::veryverbose = yes
+#carpet::verbose = yes
+#carpet::veryverbose = yes
 
 cactus::terminate = "time"
 cactus::cctk_final_time = 0.5
diff --git a/CarpetDev/CarpetAdaptiveRegrid/par/AMR2.par b/CarpetDev/CarpetAdaptiveRegrid/par/AMR2.par
new file mode 100644
index 000000000..e53527333
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/par/AMR2.par
@@ -0,0 +1,76 @@
+
+ActiveThorns = "coordbase SymBase NaNChecker carpetReduce CartGrid3D carpet carpetlib carpetadaptiveregrid Boundary IOBasic IOUtil carpetIOASCII IDWaveMoL carpetSlab WaveMoL Time MoL reflectionsymmetry"
+
+
+IDWaveMoL::initial_data = "gaussian"
+
+wavemol::bound = "radiation"
+
+mol::initial_data_is_crap = "yes"
+
+##carpet::adaptive_stepsize = "yes"
+
+##carpet::veryverbose = "yes"
+
+grid::domain = "full"
+grid::type = "coordbase"
+grid::avoid_origin = "no"
+
+coordbase::xmin = -0.5
+coordbase::ymin = -0.5
+coordbase::zmin =  0.0
+coordbase::xmax =  0.5
+coordbase::ymax =  0.5
+coordbase::zmax =  0.5
+
+coordbase::dx = 0.02
+coordbase::dy = 0.02
+coordbase::dz = 0.02
+
+reflectionsymmetry::reflection_z = yes
+reflectionsymmetry::avoid_origin_z = no
+coordbase::boundary_size_z_lower = 2
+coordbase::boundary_shiftout_z_lower = 1
+
+carpet::prolongation_order_space = 3
+carpet::prolongation_order_time = 2
+driver::ghost_size = 2
+time::dtfac = 0.5
+
+carpet::domain_from_coordbase = yes
+
+carpet::max_refinement_levels = 2
+carpet::init_3_timelevels = yes
+carpetadaptiveregrid::max_error = 1.e-6
+carpetadaptiveregrid::pad = 2
+carpetadaptiveregrid::min_width = 6
+carpetadaptiveregrid::error_var = "mol::errorestimate[0]"
+#carpetadaptiveregrid::error_var = "wavemol::phi"
+carpetadaptiveregrid::regrid_every = 1
+
+#carpetadaptiveregrid::verbose = yes
+#carpetadaptiveregrid::veryverbose = yes
+#carpet::verbose = yes
+#carpet::veryverbose = yes
+
+cactus::terminate = "time"
+cactus::cctk_final_time = 0.5
+
+iobasic::outScalar_every = 1
+iobasic::outScalar_vars = "wavemol::phi mol::errorestimate"
+
+iobasic::outInfo_every = 1
+iobasic::outInfo_vars = "wavemol::phi mol::errorestimate[0]"
+
+ioascii::out1D_every = 1
+ioascii::out1D_vars = "wavemol::scalarevolvemol_scalar mol::errorestimate"
+
+IO::out_dir = $.
+
+#mol::ode_method = "Generic"
+#mol::mol_intermediate_steps = 4
+#mol::mol_num_scratch_levels = 3
+mol::ode_method = "RK45"
+mol::mol_intermediate_steps = 6
+mol::mol_num_scratch_levels = 6
+mol::adaptive_stepsize = "no"
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index 5f198ba26..ce59cb533 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -44,6 +44,11 @@ namespace CarpetAdaptiveRegrid {
        const CCTK_INT& min_width, const CCTK_REAL& min_fraction,
        CCTK_INT& didit);
   }
+
+  ivect pos2int (const cGH* const cctkGH, const gh& hh,
+                 const rvect & rpos, const int rl);
+  rvect int2pos (const cGH* const cctkGH, const gh& hh,
+                 const ivect & ipos, const int rl);
   
   CCTK_INT CarpetAdaptiveRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
                                 CCTK_POINTER const bbsss_,
@@ -163,6 +168,16 @@ namespace CarpetAdaptiveRegrid {
     ivect imin = (bb.lower() - baseext.lower())/bb.stride(), 
       imax = (bb.upper() - baseext.lower())/bb.stride();
 
+    rvect physical_min, physical_max;
+    rvect interior_min, interior_max;
+    rvect exterior_min, exterior_max;
+    rvect base_spacing;
+    int 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);
+
     BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) 
     {
       const CCTK_REAL *error_var_ptr = 
@@ -430,14 +445,104 @@ namespace CarpetAdaptiveRegrid {
       vector<bbvect> obs;
       while (! final.empty()) {
         ibbox bb = final.top(); final.pop();
-        ivect lo = bb.lower();
-        ivect hi = bb.upper();
-        ivect str = bb.stride() / reffact;
-        newbbs.push_back (ibbox(lo,hi,str));
+        ivect ilo = bb.lower();
+        ivect ihi = bb.upper();
+        rvect lo = int2pos(cctkGH, hh, ilo, reflevel);
+        rvect hi = int2pos(cctkGH, hh, ihi, reflevel);
+        rvect str = base_spacing * 
+            ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel);
+        rbbox newbbcoord(lo, hi, str);
         
+        if (veryverbose) {
+          ostringstream buf;
+          buf << "Dealing with boundaries. Coord box is:"
+              << endl << newbbcoord;
+          CCTK_INFO(buf.str().c_str());
+        }
+
         // FIXME: Set the correct ob here.
 
         bbvect ob(false);
+        for (int d=0; d<dim; ++d) {
+          assert (mglevel==0);
+
+          // Find the size of the physical domain
+
+          rvect const spacing = base_spacing * 
+            ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel+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);
+
+          // If need be clip the domain
+
+          rvect lo = newbbcoord.lower();
+          if (newbbcoord.lower()[d] < physical_min[d]) {
+            lo[d] = exterior_min[d];
+          }
+          rvect up = newbbcoord.upper();
+          if (newbbcoord.upper()[d] > physical_max[d]) {
+            up[d] = exterior_max[d];
+          }
+          rvect str = newbbcoord.stride();
+
+          // Set the ob if outside the physical domain
+
+          ob[d][0] = 
+            abs(lo[d] - exterior_min[d]) < 1.0e-6 * spacing[d];
+          ob[d][1] = 
+            abs(up[d] - exterior_max[d]) < 1.0e-6 * spacing[d];
+
+          newbbcoord = rbbox(lo, up, str);
+        }
+        if (veryverbose) {
+          ostringstream buf;
+          buf << "Done dealing with boundaries. Coord box is:"
+              << endl << newbbcoord << endl
+              << "obox is:" << endl << ob;
+          CCTK_INFO(buf.str().c_str());
+        }
+
+        // Convert back to integer coordinates
+
+        ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel);
+        ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel);
+        ivect istr = bb.stride() / reffact;
+
+        // Check that the width is sufficient
+        // This can only be too small if the domain was clipped
+        for (int d=0; d < dim; ++d) {
+          if (ihi[d] - ilo[d] < min_width * istr[d]) {
+            if (ob[d][0]) {
+              if (ob[d][1]) {
+                CCTK_WARN(0, "The domain is too small?!");
+              }
+              ihi[d] = ilo[d] + min_width * istr[d];
+            }
+            else if (ob[d][1]) {
+              if (ob[d][0]) {
+                CCTK_WARN(0, "The domain is too small?!");
+              }
+              ilo[d] = ihi[d] - min_width * istr[d];
+            }
+            else {
+              CCTK_WARN(0, "The grid is unclipped and too small?!");
+            }
+          }
+        }
+
+        ibbox newbb(ilo, ihi, istr);          
+
+        if (veryverbose) {
+          ostringstream buf;
+          buf << "After dealing with boundaries. Final box is:"
+              << endl << newbb;
+          CCTK_INFO(buf.str().c_str());
+        }
+
+        newbbs.push_back (newbb);
         obs.push_back(ob);
       }
       // FIXME: check if the newbbs is really different from the current bbs
@@ -487,4 +592,56 @@ namespace CarpetAdaptiveRegrid {
     return do_recompose;
   }
   
+  ivect pos2int (const cGH* const cctkGH, const gh& hh,
+                 const rvect & rpos, const int rl)
+  {
+    rvect global_lower, global_upper;
+    for (int d=0; d<dim; ++d) {
+      const int ierr = CCTK_CoordRange
+	(cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
+      if (ierr<0) {
+	global_lower[d] = 0;
+	global_upper[d] = 1;
+      }
+    }
+    const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
+    
+    const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
+    const int levfac = ipow(hh.reffact, rl);
+    assert (all (hh.baseextent.stride() % levfac == 0));
+    const ivect istride = hh.baseextent.stride() / levfac;
+    
+    const ivect ipos
+      = (ivect(floor((rpos - global_lower) * scale / rvect(istride) + 0.5))
+         * istride);
+    
+    return ipos;
+  }
+  
+  rvect int2pos (const cGH* const cctkGH, const gh& hh,
+                 const ivect & ipos, const int rl)
+  {
+    rvect global_lower, global_upper;
+    for (int d=0; d<dim; ++d) {
+      const int ierr = CCTK_CoordRange
+	(cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
+      
+      if (ierr<0) {
+	global_lower[d] = 0;
+	global_upper[d] = 1;
+      }
+    }
+    const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
+    
+    const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
+    const int levfac = ipow(hh.reffact, rl);
+    assert (all (hh.baseextent.stride() % levfac == 0));
+    const ivect istride = hh.baseextent.stride() / levfac;    
+    
+    const rvect rpos
+      = rvect(ipos) / scale + global_lower;
+    
+    return rpos;
+  }
+  
 } // namespace CarpetAdaptiveRegrid
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
index 28fe821ed..5a6970565 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
@@ -26,9 +26,9 @@ namespace CarpetAdaptiveRegrid {
       = (const CCTK_INT *) CCTK_ParameterGet ("domain_from_coordbase", "Carpet", &type);
     assert (domain_from_coordbase);
     assert (type == PARAMETER_BOOLEAN);
-//     if (! *domain_from_coordbase) {
-//       CCTK_PARAMWARN ("CarpetAdaptiveRegrid requires that Carpet::domain_from_coordbase=yes");
-//     }
+    if (! *domain_from_coordbase) {
+      CCTK_PARAMWARN ("CarpetAdaptiveRegrid requires that Carpet::domain_from_coordbase=yes");
+    }
     
     return 0;
   }