CarpetRegrid2: Enforce regridding when the number of levels has changed

When using the movement_threshold parameters to regrid less often,
still regrid when the number of levels has changed, or when a region
has become active or inactive, or when the radius has changed
significantly.

darcs-hash:20070927164001-dae7b-d70720cb48753f0554f2301c4da2407d52ccb2f4.gz

5 files changed, 167 insertions, 20 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index dc688d24d..7b55ed523 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -78,16 +78,16 @@ CCTK_INT last_map "Last map at which there was a regridding at last_iteration"
 
 PUBLIC:
 
-CCTK_INT num_levels[10] TYPE=scalar
-{
-  num_levels
-} "Number of refinement levels"
-
 CCTK_INT active[10] TYPE=scalar
 {
   active
 } "Whether this centre is active"
 
+CCTK_INT num_levels[10] TYPE=scalar
+{
+  num_levels
+} "Number of refinement levels"
+
 CCTK_REAL positions[10] TYPE=scalar
 {
   position_x position_y position_z
@@ -100,7 +100,22 @@ CCTK_REAL radii[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant
 
 PRIVATE:
 
+CCTK_INT old_active[10] TYPE=scalar
+{
+  old_active
+} "Old whether this centre is active"
+
+CCTK_INT old_num_levels[10] TYPE=scalar
+{
+  old_num_levels
+} "Old number of refinement levels"
+
 CCTK_REAL old_positions[10] TYPE=scalar
 {
   old_position_x old_position_y old_position_z
 } "Old positions of refined regions"
+
+CCTK_REAL old_radii[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant
+{
+  old_radius
+} "Old radii of refined regions for each level"
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index 087fd9fb7..7212a40b9 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -96,6 +96,11 @@ CCTK_REAL movement_threshold_1 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_1 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -140,6 +145,11 @@ CCTK_REAL movement_threshold_2 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_2 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -184,6 +194,11 @@ CCTK_REAL movement_threshold_3 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_3 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -228,6 +243,11 @@ CCTK_REAL movement_threshold_4 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_4 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -272,6 +292,11 @@ CCTK_REAL movement_threshold_5 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_5 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -316,6 +341,11 @@ CCTK_REAL movement_threshold_6 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_6 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -360,6 +390,11 @@ CCTK_REAL movement_threshold_7 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_7 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -404,6 +439,11 @@ CCTK_REAL movement_threshold_8 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_8 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -448,6 +488,11 @@ CCTK_REAL movement_threshold_9 "Minimum movement to trigger a regridding" STEERA
   0:* :: ""
 } 0.0
 
+CCTK_REAL radius_change_threshold_9 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
+
 
 
 ################################################################################
@@ -491,3 +536,8 @@ CCTK_REAL movement_threshold_10 "Minimum movement to trigger a regridding" STEER
 {
   0:* :: ""
 } 0.0
+
+CCTK_REAL radius_change_threshold_10 "Minimum change in radius to trigger a regridding" STEERABLE=always
+{
+  0:* :: ""
+} 0.0
diff --git a/Carpet/CarpetRegrid2/schedule.ccl b/Carpet/CarpetRegrid2/schedule.ccl
index 27aee5200..3f53b6d38 100644
--- a/Carpet/CarpetRegrid2/schedule.ccl
+++ b/Carpet/CarpetRegrid2/schedule.ccl
@@ -1,7 +1,8 @@
 # Schedule definitions for thorn CarpetRegrid2
 
 STORAGE: last_iteration last_map
-STORAGE: num_levels active positions radii old_positions
+STORAGE: active num_levels positions radii
+STORAGE: old_active old_positions old_num_levels old_radii
 
 
 
diff --git a/Carpet/CarpetRegrid2/src/paramcheck.cc b/Carpet/CarpetRegrid2/src/paramcheck.cc
index fcc845016..5b0b08c31 100644
--- a/Carpet/CarpetRegrid2/src/paramcheck.cc
+++ b/Carpet/CarpetRegrid2/src/paramcheck.cc
@@ -25,16 +25,16 @@ namespace CarpetRegrid2 {
     DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
-    if (num_centres >= 1 and num_levels_1 > maxreflevels or
-        num_centres >= 2 and num_levels_2 > maxreflevels or
-        num_centres >= 3 and num_levels_3 > maxreflevels or
-        num_centres >= 4 and num_levels_4 > maxreflevels or
-        num_centres >= 5 and num_levels_5 > maxreflevels or
-        num_centres >= 6 and num_levels_6 > maxreflevels or
-        num_centres >= 7 and num_levels_7 > maxreflevels or
-        num_centres >= 8 and num_levels_8 > maxreflevels or
-        num_centres >= 9 and num_levels_9 > maxreflevels or
-        num_centres >= 10 and num_levels_10 > maxreflevels)
+    if ((num_centres >= 1 and num_levels_1 > maxreflevels) or
+        (num_centres >= 2 and num_levels_2 > maxreflevels) or
+        (num_centres >= 3 and num_levels_3 > maxreflevels) or
+        (num_centres >= 4 and num_levels_4 > maxreflevels) or
+        (num_centres >= 5 and num_levels_5 > maxreflevels) or
+        (num_centres >= 6 and num_levels_6 > maxreflevels) or
+        (num_centres >= 7 and num_levels_7 > maxreflevels) or
+        (num_centres >= 8 and num_levels_8 > maxreflevels) or
+        (num_centres >= 9 and num_levels_9 > maxreflevels) or
+        (num_centres >= 10 and num_levels_10 > maxreflevels))
     {
       CCTK_PARAMWARN ("The number of requested refinement levels is larger than the maximum number of levels specified by Carpet::max_refinement_levels");
     }
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index eea342919..4e6efe469 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -1,3 +1,4 @@
+#include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
@@ -720,9 +721,23 @@ namespace CarpetRegrid2 {
     }
     
     if (do_recompose and * last_iteration != -1) {
-      // Regrid only if the positions have changed sufficiently
+      
+      // Regrid only if the regions have changed sufficiently
       do_recompose = false;
       for (int n = 0; n < num_centres; ++ n) {
+        
+        // Regrid if a region became active or inactive
+        do_recompose = active[n] != old_active[n];
+        if (do_recompose) break;
+        
+        // Check only active regions
+        if (! active[n]) continue;
+        
+        // Regrid if the number of levels changed
+        do_recompose = num_levels[n] != old_num_levels[n];
+        if (do_recompose) break;
+        
+        // Regrid if the positions have changed sufficiently
         CCTK_REAL const dist2 =
           pow (position_x[n] - old_position_x[n], 2) +
           pow (position_y[n] - old_position_y[n], 2) +
@@ -743,11 +758,31 @@ namespace CarpetRegrid2 {
         }
         do_recompose = dist2 >= pow (mindist, 2);
         if (do_recompose) break;
+        
+        // Regrid if the radii have changed sufficiently
+        CCTK_REAL const dr = abs (radius[n] - old_radius[n]);
+        CCTK_REAL mindr;
+        switch (n) {
+        case 0: mindr = radius_change_threshold_1; break;
+        case 1: mindr = radius_change_threshold_2; break;
+        case 2: mindr = radius_change_threshold_3; break;
+        case 3: mindr = radius_change_threshold_4; break;
+        case 4: mindr = radius_change_threshold_5; break;
+        case 5: mindr = radius_change_threshold_6; break;
+        case 6: mindr = radius_change_threshold_7; break;
+        case 7: mindr = radius_change_threshold_8; break;
+        case 8: mindr = radius_change_threshold_9; break;
+        case 9: mindr = radius_change_threshold_10; break;
+        default: assert (0);
+        }
+        do_recompose = dr >= mindr;
+        if (do_recompose) break;
+        
       } // for n
       if (verbose) {
         if (not do_recompose) {
           CCTK_INFO
-            ("Centres have not moved sufficiently; skipping regridding");
+            ("Refined regions have not changed sufficiently; skipping regridding");
         }
       }
     }
@@ -776,9 +811,15 @@ namespace CarpetRegrid2 {
       
       // Remember current positions
       for (int n = 0; n < num_centres; ++ n) {
+        old_active[n] = active[n];
+        
+        old_num_levels[n] = num_levels[n];
+        
         old_position_x[n] = position_x[n];
         old_position_y[n] = position_y[n];
         old_position_z[n] = position_z[n];
+        
+        old_radius[n] = radius[n];
       }
       
     } // if do_recompose
@@ -839,9 +880,23 @@ namespace CarpetRegrid2 {
     }
     
     if (do_recompose and * last_iteration != -1) {
-      // Regrid only if the positions have changed sufficiently
+      
+      // Regrid only if the regions have changed sufficiently
       do_recompose = false;
       for (int n = 0; n < num_centres; ++ n) {
+        
+        // Regrid if a region became active or inactive
+        do_recompose = active[n] != old_active[n];
+        if (do_recompose) break;
+        
+        // Check only active regions
+        if (! active[n]) continue;
+        
+        // Regrid if the number of levels changed
+        do_recompose = num_levels[n] != old_num_levels[n];
+        if (do_recompose) break;
+        
+        // Regrid if the positions have changed sufficiently
         CCTK_REAL const dist2 =
           pow (position_x[n] - old_position_x[n], 2) +
           pow (position_y[n] - old_position_y[n], 2) +
@@ -862,11 +917,31 @@ namespace CarpetRegrid2 {
         }
         do_recompose = dist2 >= pow (mindist, 2);
         if (do_recompose) break;
+        
+        // Regrid if the radii have changed sufficiently
+        CCTK_REAL const dr = abs (radius[n] - old_radius[n]);
+        CCTK_REAL mindr;
+        switch (n) {
+        case 0: mindr = radius_change_threshold_1; break;
+        case 1: mindr = radius_change_threshold_2; break;
+        case 2: mindr = radius_change_threshold_3; break;
+        case 3: mindr = radius_change_threshold_4; break;
+        case 4: mindr = radius_change_threshold_5; break;
+        case 5: mindr = radius_change_threshold_6; break;
+        case 6: mindr = radius_change_threshold_7; break;
+        case 7: mindr = radius_change_threshold_8; break;
+        case 8: mindr = radius_change_threshold_9; break;
+        case 9: mindr = radius_change_threshold_10; break;
+        default: assert (0);
+        }
+        do_recompose = dr >= mindr;
+        if (do_recompose) break;
+        
       } // for n
       if (verbose) {
         if (not do_recompose) {
           CCTK_INFO
-            ("Centres have not moved sufficiently; skipping regridding");
+            ("Refined regions have not changed sufficiently; skipping regridding");
         }
       }
     }
@@ -921,9 +996,15 @@ namespace CarpetRegrid2 {
       
       // Remember current positions
       for (int n = 0; n < num_centres; ++ n) {
+        old_active[n] = active[n];
+        
+        old_num_levels[n] = num_levels[n];
+        
         old_position_x[n] = position_x[n];
         old_position_y[n] = position_y[n];
         old_position_z[n] = position_z[n];
+        
+        old_radius[n] = radius[n];
       }
       
     } // if do_recompose