CarpetRegrid2: Implement true AMR

5 files changed, 256 insertions, 1 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 8202f4072..ee2bb307d 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -139,6 +139,8 @@ CCTK_REAL radiixyz[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant
   radius_x radius_y radius_z
 } "Radii of refined regions for each level"
 
+CCTK_REAL level_mask TYPE=gf TAGS='checkpoint="no" prolongation="none"' "Requested refinement level for this grid point"
+
 PRIVATE:
 
 CCTK_INT old_active[10] TYPE=scalar "Old whether this centre is active"
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index c45241de9..70b1be5f1 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -64,6 +64,33 @@ BOOLEAN symmetry_rotating180 "Ensure a 180 degree rotating symmetry about the z
 
 
 
+BOOLEAN adaptive_refinement "Use level_mask for adaptive refinement" STEERABLE=always
+{
+} "no"
+
+CCTK_INT adaptive_block_size "Block size for adaptive refinement" STEERABLE=always
+{
+  1:* :: ""
+} 8
+
+CCTK_INT adaptive_block_size_x "Block size in x direction for adaptive refinement" STEERABLE=always
+{
+  -1 :: "use adaptive_block_size"
+  1:* :: ""
+} -1
+CCTK_INT adaptive_block_size_y "Block size in y direction for adaptive refinement" STEERABLE=always
+{
+  -1 :: "use adaptive_block_size"
+  1:* :: ""
+} -1
+CCTK_INT adaptive_block_size_z "Block size in z direction for adaptive refinement" STEERABLE=always
+{
+  -1 :: "use adaptive_block_size"
+  1:* :: ""
+} -1
+
+
+
 CCTK_INT num_centres "Number of refinement centres"
 {
   0:10 :: ""
diff --git a/Carpet/CarpetRegrid2/src/amr.cc b/Carpet/CarpetRegrid2/src/amr.cc
new file mode 100644
index 000000000..572654f6d
--- /dev/null
+++ b/Carpet/CarpetRegrid2/src/amr.cc
@@ -0,0 +1,210 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <cmath>
+
+#include <carpet.hh>
+#include <mpi_string.hh>
+
+#include "boundary.hh"
+
+
+
+namespace CarpetRegrid2 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  void evaluate_level_mask (cGH const* restrict const cctkGH,
+                            vector<ibset>& regions, int const rl)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (is_singlemap_mode());
+    int const m = Carpet::map;
+    assert (rl > 0);
+    BEGIN_GLOBAL_MODE (cctkGH) {
+      ENTER_LEVEL_MODE (cctkGH, rl-1) {
+        ENTER_SINGLEMAP_MODE (cctkGH, m, CCTK_GF) {
+          if (verbose or veryverbose) {
+            cout << "Refinement level " << rl << ":\n";
+          }
+          
+          gh const& hh = *vhh.AT(Carpet::map);
+          
+          ivect const gsh         = ivect::ref(cctkGH->cctk_gsh);
+          ivect const nghostzones = ivect::ref(cctkGH->cctk_nghostzones);
+          if (veryverbose) {
+            cout << "   gsh: " << gsh << "\n"
+                 << "   nghostzones: " << nghostzones << "\n";
+          }
+          
+          // Determine the block size
+          ivect block_size = adaptive_block_size;
+          if (adaptive_block_size_x > 0) block_size[0] = adaptive_block_size_x;
+          if (adaptive_block_size_y > 0) block_size[1] = adaptive_block_size_y;
+          if (adaptive_block_size_z > 0) block_size[2] = adaptive_block_size_z;
+          assert (all (block_size > 0));
+          
+          // Determine the block offset, i.e. the grid point index of
+          // the coordinate origin. We don't start blocks at index
+          // zero, so that one can change the location of the outer
+          // boundary without influencing the refinement mechanism.
+          // Instead, we align the block such that one of the blocks
+          // begins at the coordinate origin.
+          ivect origin;
+          rvect rorigin;
+          {
+            DECLARE_CCTK_ARGUMENTS;
+            for (int d=0; d<dim; ++d) {
+              // We round up, because the origin is located between
+              // grid points for cell-centred grids, and we want the
+              // next upper grid point in this case. We subtract 0.25
+              // just for proper rounding.
+              origin[d] =
+                ceil(- CCTK_ORIGIN_SPACE(d) / CCTK_DELTA_SPACE(d) - 0.25);
+              rorigin[d] =
+                CCTK_ORIGIN_SPACE(d) + origin[d] * CCTK_DELTA_SPACE(d);
+            }
+          }
+          // block_index := ([0...gsh] + block_offset) / block_size
+          // [0...gsh] = block_index * block_size - block_offset
+          ivect const block_offset =
+            (block_size - origin % block_size) % block_size;
+          ivect const num_blocks =
+            (gsh + block_offset + block_size - 1) / block_size;
+          if (veryverbose) {
+            cout << "   origin: " << origin << "\n"
+                 << "   coordinates of origin: " << rorigin << "\n"
+                 << "   offset: " << block_offset << "\n"
+                 << "   num blocks: " << num_blocks << "\n";
+          }
+          assert (all (num_blocks > 0));
+          
+          // Mask containing all refined block indices (in arbitrary
+          // order, and potentially even duplicated). This is really a
+          // set, but we need to communicate this data structure, and
+          // therefore a vector is a better choice.
+          vector<ivect> mask;
+            
+          // For vertex centred grids, the blocks need to overlap by
+          // one grid point, so that e.g. an 8^3 block is turned into
+          // a 9^3 block. This will ensure symmetry of refined regions
+          // in the simulation domain.
+          int const overlap = hh.refcent == vertex_centered ? 1 : 0;
+          
+          BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+            DECLARE_CCTK_ARGUMENTS;
+            
+            if (veryverbose) {
+              cout << "   component " << component << ":\n";
+            }
+            
+            ivect const lbnd = ivect::ref(cctk_lbnd);
+            ivect const lssh (CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2));
+            ivect const bboxlo (cctk_bbox[0], cctk_bbox[2], cctk_bbox[4]);
+            ivect const bboxhi (cctk_bbox[1], cctk_bbox[3], cctk_bbox[5]);
+            
+            ivect const imin = 0    + either(bboxlo, 0, nghostzones);
+            ivect const imax = lssh - either(bboxhi, 0, nghostzones);
+            
+            ivect const bmin =
+              max (0,
+                   ((lbnd + imin + block_offset + block_size - overlap) /
+                    block_size) -
+                   1);
+            ivect const bmax = (lbnd + imax + block_offset) / block_size;
+            
+            // Loop over all blocks
+            int nblocks = 0;
+            int nrefined = 0;
+            for (int bk=bmin[2]; bk<bmax[2]; ++bk) {
+              for (int bj=bmin[1]; bj<bmax[1]; ++bj) {
+                for (int bi=bmin[0]; bi<bmax[0]; ++bi) {
+                  ivect const bind (bi,bj,bk);
+                  
+                  ivect const bimin =
+                    max(imin,
+                        (bind  ) * block_size - block_offset - lbnd          );
+                  ivect const bimax =
+                    min(imax,
+                        (bind+1) * block_size - block_offset - lbnd + overlap);
+                  
+                  bool refine = false;
+                  
+                  // Loop over all points in this block
+                  for (int k=bimin[2]; k<bimax[2]; ++k) {
+                    for (int j=bimin[1]; j<bimax[1]; ++j) {
+                      for (int i=bimin[0]; i<bimax[0]; ++i) {
+                        int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
+                        refine = refine or level_mask[ind3d] >= rl;
+                      }
+                    }
+                  }
+                  
+                  // Refine this block if any point in this block requires
+                  // refinement
+                  if (refine) {
+                    if (veryverbose) {
+                      cout << "      refining block " << bind << "\n";
+                    }
+                    mask.push_back(bind);
+                    ++ nrefined;
+                  }
+                  ++ nblocks;
+                  
+                }
+              }
+            }
+            
+            if (veryverbose) {
+              cout << "   refining " << nrefined << " "
+                   << "of " << nblocks << " blocks on this component\n";
+            }
+            
+          } END_LOCAL_COMPONENT_LOOP;
+          
+          // Combine this mask from all processes, and to all processes
+          vector<ivect> const fullmask = allgatherv1 (dist::comm(), mask);
+          
+          // Convert vector into ibset
+          const ibbox& cbaseext = hh.baseextent(mglevel, rl-1);
+          const ibbox& baseext = hh.baseextent(mglevel, rl);
+          
+          ibset& region = regions.at(rl);
+          for (vector<ivect>::const_iterator
+                 it = fullmask.begin(); it != fullmask.end(); ++it)
+          {
+            ivect const& bind = *it;
+            ivect const cstr = cbaseext.stride();
+            ivect const clo  =
+              cbaseext.lower() + (bind * block_size - block_offset) * cstr;
+            ivect const cup  = clo + (block_size + overlap - 1) * cstr;
+            ibbox const cblock (clo, cup, cstr);
+            ibbox const block = cblock.expanded_for(baseext);
+            region |= block;
+          }
+          
+          if (verbose or veryverbose) {
+            if (veryverbose) {
+              cout << "   refined region is " << region << "\n"
+                   << "   domain is " << baseext << "\n"
+                   << "   next coarser level is " << regions.at(rl-1) << "\n";
+            }
+            ibbox::size_type const nrefined = region.size();
+            ibbox::size_type const npoints = baseext.size();
+            ibbox::size_type const ncoarse =
+              regions.at(rl-1).expanded_for(baseext).size();
+            cout << "   refining "
+                 << (100.0*nrefined/ncoarse) << "% of the next coarser level\n"
+                 << "   refining "
+                 << (100.0*nrefined/npoints) << "% of the domain\n";
+          }
+          
+        } LEAVE_SINGLEMAP_MODE;
+      } LEAVE_LEVEL_MODE;
+    } END_GLOBAL_MODE;
+  }
+  
+} // namespace CarpetRegrid2
diff --git a/Carpet/CarpetRegrid2/src/amr.hh b/Carpet/CarpetRegrid2/src/amr.hh
new file mode 100644
index 000000000..27abb6349
--- /dev/null
+++ b/Carpet/CarpetRegrid2/src/amr.hh
@@ -0,0 +1,15 @@
+#ifndef AMR_HH
+#define AMR_HH
+
+#include <carpet.hh>
+
+
+
+namespace CarpetRegrid2 {
+  
+  void evaluate_level_mask (cGH const* restrict cctkGH,
+                            vector<ibset>& regions, int rl);
+  
+} // namespace CarpetRegrid2
+
+#endif  // #ifndef AMR_HH
diff --git a/Carpet/CarpetRegrid2/src/make.code.defn b/Carpet/CarpetRegrid2/src/make.code.defn
index e8bd87f0b..db068520d 100644
--- a/Carpet/CarpetRegrid2/src/make.code.defn
+++ b/Carpet/CarpetRegrid2/src/make.code.defn
@@ -1,7 +1,8 @@
 # Main make.code.defn file for thorn CarpetRegrid2
 
 # Source files in this directory
-SRCS =	boundary.cc				\
+SRCS =	amr.cc					\
+	boundary.cc				\
 	indexing.cc				\
 	initialise.cc				\
 	paramcheck.cc				\