CarpetLib: Do not interpolate from buffer zones, if requested

Introduce parameter CarpetLib::interpolate_from_buffer_zones that allows disabling interpolation from buffer zones.

2 files changed, 36 insertions, 7 deletions

diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index 25e87dded..deca16868 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -33,6 +33,12 @@ BOOLEAN use_dgfe "Use DGFE operators instead of Lagrange operators"
 
 
 
+BOOLEAN interpolate_from_buffer_zones "Use buffer points for interpolation" STEERABLE=always
+{
+} "yes"
+
+
+
 restricted:
 
 BOOLEAN use_higher_order_restriction "Use third order cell centered restriction operators instead of first order"
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 7e0dba448..6bb0adb1d 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -443,6 +443,8 @@ locate_position (rvect const & rpos,
                  int const minrl, int const maxrl,
                  int & rl, int & c, ivect & aligned_ipos) const
 {
+  DECLARE_CCTK_PARAMETERS;
+  
   assert (ml>=0 and ml<mglevels());
   assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
   
@@ -452,6 +454,10 @@ locate_position (rvect const & rpos,
     return;
   }
   
+  // Find associated dh
+  assert (dhs.size() == 1);
+  dh const& dd = *dhs.front();
+  
   // Try finer levels first
   for (rl = maxrl-1; rl >= minrl; --rl) {
     
@@ -460,14 +466,17 @@ locate_position (rvect const & rpos,
     ivect const str = baseextent(ml,rl).stride();
     aligned_ipos = ivect(floor(rpos / rvect(str) + rvect(0.5))) * str;
     
-    gh::cregs const & regs = superregions.AT(rl);
+    // Ignore this level if this point is not in the active region
+    // (i.e. if it is a buffer point or similar)
+    if (not interpolate_from_buffer_zones and 
+        not dd.level_boxes.AT(ml).AT(rl).active.contains(aligned_ipos))
+    {
+      continue;
+    }
     
-    // TODO: Ignore this level if this point is not in the active
-    // region (i.e. if it is a buffer point or similar)
-    // TODO: check the level_bbox
-    // TODO: (assume this gh has only one dh)
+    gh::cregs const & regs = superregions.AT(rl);
     
-    // Search all superregions linearly.  Each superregion corresponds
+    // Search all superregions linearly. Each superregion corresponds
     // to a "refined region", and the number of superregions is thus
     // presumably independent of the number of processors.
     for (size_t r = 0; r < regs.size(); ++r) {
@@ -501,9 +510,15 @@ locate_position (ivect const & ipos,
                  int const minrl, int const maxrl,
                  int & rl, int & c, ivect & aligned_ipos) const
 {
+  DECLARE_CCTK_PARAMETERS;
+  
   assert (ml>=0 and ml<mglevels());
   assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
   
+  // Find associated dh
+  assert (dhs.size() == 1);
+  dh const& dd = *dhs.front();
+  
   // Try finer levels first
   for (rl = maxrl-1; rl >= minrl; --rl) {
     
@@ -512,9 +527,17 @@ locate_position (ivect const & ipos,
     ivect const str = baseextent(ml, rl).stride();
     aligned_ipos = ivect(floor(rvect(ipos) / rvect(str) + rvect(0.5))) * str;
     
+    // Ignore this level if this point is not in the active region
+    // (i.e. if it is a buffer point or similar)
+    if (not interpolate_from_buffer_zones and 
+        not dd.level_boxes.AT(ml).AT(rl).active.contains(aligned_ipos))
+    {
+      continue;
+    }
+    
     gh::cregs const & regs = superregions.AT(rl);
     
-    // Search all superregions linearly.  Each superregion corresponds
+    // Search all superregions linearly. Each superregion corresponds
     // to a "refined region", and the number of superregions is thus
     // presumably independent of the number of processors.
     for (size_t r = 0; r < regs.size(); ++r) {