Merged changes

5 files changed, 93 insertions, 82 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index a49f122b5..7ecba6785 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -997,7 +997,8 @@ regrid (bool const do_init)
       static Carpet::Timer timer_mask ("CarpetLib::dh::regrid::mask");
       timer_mask.start();
       
-      // Declare this here to save it for 'unused-mask'
+      // Declare this here to save it for later use. Contains all the boxes
+      // which are active minus the boundary
       ibset all_refined;
 
       if (rl > 0) {
@@ -1092,8 +1093,6 @@ regrid (bool const do_init)
           full_dboxes const& box = full_level.AT(c);
           local_dboxes & local_box = local_level.AT(lc);
           
-          // Subtract the boundaries from the refined region
-          all_refined = allactive;
           // Set prolongation information for current level
           for (int d=0; d<dim; ++d) {
             local_box.prolongation_boundaries[d] =
@@ -1125,7 +1124,7 @@ regrid (bool const do_init)
         if (all (reffact == 2)) {
           // use the already computed 'all_refined' to get region from where
           // no information will be used later (overwritten)
-          // First: get the region which will get restricted
+          // First: get the region which will get restricted, on the coarse level
           ibset restricted_region = all_refined.contracted_for(h.baseextent(ml,orl));
           // This is too big - during MoL-substeps information within this
           // region will be used to update points outside -> need to
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index b281031a4..23e823257 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -405,7 +405,7 @@ rpos2ipos1 (rvect const & rpos,
 // Find the refinement level and component to which a grid point
 // belongs.  This uses a tree search over the superregions in the grid
 // struction, which should scale reasonably (O(n log n)) better with
-// the number of componets components.
+// the number of components.
 void
 gh::
 locate_position (rvect const & rpos,
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index 46bd9f6c2..5f87e9a9a 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -132,9 +132,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_1 "Minimum RELATIVE change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -200,9 +200,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_2 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -268,9 +268,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_3 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -336,9 +336,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_4 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -404,9 +404,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_5 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -472,9 +472,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_6 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -540,9 +540,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_7 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -608,9 +608,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_8 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -676,9 +676,9 @@ 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
+CCTK_REAL radius_rel_change_threshold_9 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
 
 
@@ -744,7 +744,7 @@ 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
+CCTK_REAL radius_rel_change_threshold_10 "Minimum change in radius to trigger a regridding" STEERABLE=always
 {
-  0:* :: ""
+  0.0:* :: ""
 } 0.0
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 35b048a4b..95561d8c1 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -33,10 +33,10 @@ namespace CarpetRegrid2 {
   
   
   struct centre_description {
-    int            num_levels;
-    int            active;
-    rvect          position;
-    vector <rvect> radius;
+    int            _num_levels;
+    int            _active;
+    rvect          _position;
+    vector <rvect> _radius;
     
     centre_description (cGH const * cctkGH, int n);
   };
@@ -56,21 +56,21 @@ namespace CarpetRegrid2 {
     int lsh[2];
     getvectorindex2 (cctkGH, "CarpetRegrid2::radii", lsh);
     
-    this->num_levels = num_levels[n];
-    this->active = active[n];
-    this->position = rvect (position_x[n], position_y[n], position_z[n]);
-    if (any (not isfinite(this->position))) {
+    this->_num_levels = num_levels[n];
+    this->_active = active[n];
+    this->_position = rvect (position_x[n], position_y[n], position_z[n]);
+    if (any (not isfinite(this->_position))) {
       CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "The position of region %d is [%g,%g,%g], which is not finite",
                   n + 1,
-                  double(this->position[0]),
-                  double(this->position[1]),
-                  double(this->position[2]));
+                  double(this->_position[0]),
+                  double(this->_position[1]),
+                  double(this->_position[2]));
       found_error = true;
     }
-    this->radius.resize (this->num_levels);
-    this->radius.at(0) = rvect(-1.0, -1.0, -1.0); // unused
-    for (int rl = 1; rl < this->num_levels; ++ rl) {
+    this->_radius.resize (this->_num_levels);
+    this->_radius.at(0) = rvect(-1.0, -1.0, -1.0); // unused
+    for (int rl = 1; rl < this->_num_levels; ++ rl) {
       int const ind = index2 (lsh, rl, n);
       CCTK_REAL const rx = radius_x[ind] < 0.0 ? radius[ind] : radius_x[ind];
       CCTK_REAL const ry = radius_y[ind] < 0.0 ? radius[ind] : radius_y[ind];
@@ -80,27 +80,27 @@ namespace CarpetRegrid2 {
         CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "The radius of refinement level %d of region %d is [%g,%g,%g], which is not finite",
                     rl, n + 1,
-                    double(this->radius.at(rl)[0]),
-                    double(this->radius.at(rl)[1]),
-                    double(this->radius.at(rl)[2]));
+                    double(this->_radius.at(rl)[0]),
+                    double(this->_radius.at(rl)[1]),
+                    double(this->_radius.at(rl)[2]));
         found_error = true;
       }
       if (any (rad < 0.0)) {
         CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                     "The radius of refinement level %d of region %d is [%g,%g,%g], which is non-negative",
                     rl, n + 1,
-                    double(this->radius.at(rl)[0]),
-                    double(this->radius.at(rl)[1]),
-                    double(this->radius.at(rl)[2]));
+                    double(this->_radius.at(rl)[0]),
+                    double(this->_radius.at(rl)[1]),
+                    double(this->_radius.at(rl)[2]));
         found_error = true;
       }
-      this->radius.at(rl) = rad;
+      this->_radius.at(rl) = rad;
     }
     
-    if (this->num_levels > maxreflevels) {
+    if (this->_num_levels > maxreflevels) {
       CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "Region %d has %d levels active, which is larger than the maximum number of refinement levels %d",
-                  n + 1, this->num_levels, maxreflevels);
+                  n + 1, this->_num_levels, maxreflevels);
       found_error = true;
     }
     
@@ -439,18 +439,18 @@ namespace CarpetRegrid2 {
       // Loop over all centres
       for (int n = 0; n < num_centres; ++ n) {
         centre_description centre (cctkGH, n);
-        if (centre.active) {
+        if (centre._active) {
           
           // Loop over all levels for this centre
-          for (int rl = min_rl; rl < centre.num_levels; ++ rl) {
+          for (int rl = min_rl; rl < centre._num_levels; ++ rl) {
             
             if (veryverbose) {
               cout << "Refinement level " << rl << ": importing refined region settings..." << endl;
             }
             
             // Calculate a bbox for this region
-            rvect const rmin = centre.position - centre.radius.at(rl);
-            rvect const rmax = centre.position + centre.radius.at(rl);
+            rvect const rmin = centre._position - centre._radius.at(rl);
+            rvect const rmax = centre._position + centre._radius.at(rl);
             
             // Convert to an integer bbox
             ivect const istride = hh.baseextents.at(0).at(rl).stride();
@@ -1172,22 +1172,23 @@ namespace CarpetRegrid2 {
           rvect const rad (rx, ry, rz);
           rvect const oldrad
             (old_radius_x[ind], old_radius_y[ind], old_radius_z[ind]);
-          CCTK_REAL const dr = sqrt (sum (ipow (rad - oldrad, 2)));
-          CCTK_REAL mindr;
+          CCTK_REAL const drfac = 
+	    (sqrt (sum (ipow (rad - oldrad, 2))))/(sqrt (sum (ipow (oldrad, 2))));
+          CCTK_REAL mindrfac;
           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;
+          case 0: mindrfac = radius_rel_change_threshold_1; break;
+          case 1: mindrfac = radius_rel_change_threshold_2; break;
+          case 2: mindrfac = radius_rel_change_threshold_3; break;
+          case 3: mindrfac = radius_rel_change_threshold_4; break;
+          case 4: mindrfac = radius_rel_change_threshold_5; break;
+          case 5: mindrfac = radius_rel_change_threshold_6; break;
+          case 6: mindrfac = radius_rel_change_threshold_7; break;
+          case 7: mindrfac = radius_rel_change_threshold_8; break;
+          case 8: mindrfac = radius_rel_change_threshold_9; break;
+          case 9: mindrfac = radius_rel_change_threshold_10; break;
           default: assert (0);
           }
-          do_recompose = dr > mindr;
+          do_recompose = drfac > mindrfac;
           if (do_recompose) break;
         } // for rl
         if (do_recompose) break;
@@ -1352,22 +1353,23 @@ namespace CarpetRegrid2 {
           rvect const rad (rx, ry, rz);
           rvect const oldrad
             (old_radius_x[ind], old_radius_y[ind], old_radius_z[ind]);
-          CCTK_REAL const dr = sqrt (sum (ipow (rad - oldrad, 2)));
-          CCTK_REAL mindr;
+          CCTK_REAL const drfac = 
+	    (sqrt (sum (ipow (rad - oldrad, 2))))/(sqrt (sum (ipow (oldrad, 2))));
+          CCTK_REAL mindrfac;
           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;
+          case 0: mindrfac = radius_rel_change_threshold_1; break;
+          case 1: mindrfac = radius_rel_change_threshold_2; break;
+          case 2: mindrfac = radius_rel_change_threshold_3; break;
+          case 3: mindrfac = radius_rel_change_threshold_4; break;
+          case 4: mindrfac = radius_rel_change_threshold_5; break;
+          case 5: mindrfac = radius_rel_change_threshold_6; break;
+          case 6: mindrfac = radius_rel_change_threshold_7; break;
+          case 7: mindrfac = radius_rel_change_threshold_8; break;
+          case 8: mindrfac = radius_rel_change_threshold_9; break;
+          case 9: mindrfac = radius_rel_change_threshold_10; break;
           default: assert (0);
           }
-          do_recompose = dr > mindr;
+          do_recompose = drfac > mindrfac;
           if (do_recompose) break;
         } // for rl
         if (do_recompose) break;
diff --git a/Carpet/LoopControl/src/loopcontrol.c b/Carpet/LoopControl/src/loopcontrol.c
index 875c6fc60..a0bbe564b 100644
--- a/Carpet/LoopControl/src/loopcontrol.c
+++ b/Carpet/LoopControl/src/loopcontrol.c
@@ -406,9 +406,10 @@ lc_statset_init (lc_statset_t * restrict const ls,
   /*** Threads ****************************************************************/
   
   static int saved_maxthreads = -1;
-  static lc_topology_t * restrict * saved_topologies = NULL;
+  static lc_topology_t * * saved_topologies = NULL;
   static int * restrict saved_ntopologies = NULL;
-  
+
+  // Allocate memory the first time this function is called
   if (saved_maxthreads < 0) {
     saved_maxthreads = omp_get_max_threads();
     saved_topologies  = malloc (saved_maxthreads * sizeof * saved_topologies );
@@ -420,9 +421,18 @@ lc_statset_init (lc_statset_t * restrict const ls,
       saved_ntopologies[n] = -1;
     }
   }
-  
+  // Reallocate memory in case we need more
   if (num_threads > saved_maxthreads) {
-    CCTK_WARN (CCTK_WARN_ABORT, "Thread count inconsistency");
+    int old_saved_maxthreads = saved_maxthreads;
+    saved_maxthreads = num_threads;
+    saved_topologies  = realloc (saved_topologies,  saved_maxthreads * sizeof * saved_topologies );
+    saved_ntopologies = realloc (saved_ntopologies, saved_maxthreads * sizeof * saved_ntopologies);
+    assert (saved_topologies );
+    assert (saved_ntopologies);
+    for (int n=old_saved_maxthreads; n<saved_maxthreads; ++n) {
+      saved_topologies [n] = NULL;
+      saved_ntopologies[n] = -1;
+    }
   }
   
   if (! saved_topologies[num_threads-1]) {