CarpetRegrid2: Allow rectangular (cuboid, non-cubic) refined regions

5 files changed, 236 insertions, 28 deletions

diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 30fa070d9..41e235a46 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -133,6 +133,11 @@ CCTK_REAL radii[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant
   radius
 } "Radii of refined regions for each level"
 
+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"
+
 PRIVATE:
 
 CCTK_INT old_active[10] TYPE=scalar "Old whether this centre is active"
@@ -144,7 +149,7 @@ 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
+CCTK_REAL old_radiixyz[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant
 {
-  old_radius
+  old_radius_x old_radius_y old_radius_z
 } "Old radii of refined regions for each level"
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index cacdd6d2e..8a7e53baa 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -102,6 +102,24 @@ CCTK_REAL radius_1[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_1[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_1[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_1[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_1 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -151,6 +169,24 @@ CCTK_REAL radius_2[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_2[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_2[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_2[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_2 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -200,6 +236,24 @@ CCTK_REAL radius_3[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_3[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_3[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_3[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_3 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -249,6 +303,24 @@ CCTK_REAL radius_4[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_4[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_4[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_4[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_4 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -298,6 +370,24 @@ CCTK_REAL radius_5[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_5[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_5[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_5[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_5 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -347,6 +437,24 @@ CCTK_REAL radius_6[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_6[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_6[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_6[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_6 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -396,6 +504,24 @@ CCTK_REAL radius_7[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_7[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_7[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_7[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_7 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -445,6 +571,24 @@ CCTK_REAL radius_8[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_8[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_8[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_8[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_8 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -494,6 +638,24 @@ CCTK_REAL radius_9[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_9[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_9[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_9[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_9 "Minimum movement to trigger a regridding" STEERABLE=always
@@ -543,6 +705,24 @@ CCTK_REAL radius_10[30] "Radius of refined region for this level of this centre"
   0:* :: ""
 } 1.0
 
+CCTK_REAL radius_x_10[30] "x-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_y_10[30] "y-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
+CCTK_REAL radius_z_10[30] "z-radius of refined region for this level of this centre"
+{
+  -1  :: "ignore this radius"
+  0:* :: ""
+} -1.0
+
 
 
 CCTK_REAL movement_threshold_10 "Minimum movement to trigger a regridding" STEERABLE=always
diff --git a/Carpet/CarpetRegrid2/schedule.ccl b/Carpet/CarpetRegrid2/schedule.ccl
index 3f53b6d38..2d93afae5 100644
--- a/Carpet/CarpetRegrid2/schedule.ccl
+++ b/Carpet/CarpetRegrid2/schedule.ccl
@@ -1,8 +1,8 @@
 # Schedule definitions for thorn CarpetRegrid2
 
 STORAGE: last_iteration last_map
-STORAGE: active num_levels positions radii
-STORAGE: old_active old_positions old_num_levels old_radii
+STORAGE: active num_levels positions radii radiixyz
+STORAGE: old_active old_positions old_num_levels old_radiixyz
 
 
 
diff --git a/Carpet/CarpetRegrid2/src/initialise.cc b/Carpet/CarpetRegrid2/src/initialise.cc
index 09bf84fee..789f168e5 100644
--- a/Carpet/CarpetRegrid2/src/initialise.cc
+++ b/Carpet/CarpetRegrid2/src/initialise.cc
@@ -35,21 +35,25 @@ namespace CarpetRegrid2 {
     int lsh[2];
     getvectorindex2 (cctkGH, "CarpetRegrid2::radii", lsh);
     
-#define INIT_CENTRE(N)                                          \
-    do {                                                        \
-      if (num_centres >= N) {                                   \
-        num_levels[N-1] = num_levels_##N;                       \
-        active    [N-1] = active_##N;                           \
-        position_x[N-1] = position_x_##N;                       \
-        position_y[N-1] = position_y_##N;                       \
-        position_z[N-1] = position_z_##N;                       \
-        for (int rl = 0; rl < 30; ++ rl) {                      \
-          radius[index2 (lsh, rl, N-1)] = radius_##N[rl];       \
-        }                                                       \
-        old_position_x[N-1] = position_x[N-1];                  \
-        old_position_y[N-1] = position_y[N-1];                  \
-        old_position_z[N-1] = position_z[N-1];                  \
-      }                                                         \
+#define INIT_CENTRE(N)                                  \
+    do {                                                \
+      if (num_centres >= N) {                           \
+        num_levels[N-1] = num_levels_##N;               \
+        active    [N-1] = active_##N;                   \
+        position_x[N-1] = position_x_##N;               \
+        position_y[N-1] = position_y_##N;               \
+        position_z[N-1] = position_z_##N;               \
+        for (int rl = 0; rl < 30; ++ rl) {              \
+          int const ind = index2 (lsh, rl, N-1);        \
+          radius[ind] = radius_##N[rl];                 \
+          radius_x[ind] = radius_x_##N[rl];             \
+          radius_y[ind] = radius_y_##N[rl];             \
+          radius_z[ind] = radius_z_##N[rl];             \
+        }                                               \
+        old_position_x[N-1] = position_x[N-1];          \
+        old_position_y[N-1] = position_y[N-1];          \
+        old_position_z[N-1] = position_z[N-1];          \
+      }                                                 \
     } while (0)
     
     INIT_CENTRE( 1);
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 79af7c84b..2712b2f32 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -40,10 +40,10 @@ namespace CarpetRegrid2 {
   
   
   struct centre_description {
-    int                num_levels;
-    int                active;
-    rvect              position;
-    vector <CCTK_REAL> radius;
+    int            num_levels;
+    int            active;
+    rvect          position;
+    vector <rvect> radius;
     
     centre_description (cGH const * cctkGH, int n);
   };
@@ -66,7 +66,12 @@ namespace CarpetRegrid2 {
     this->position = rvect (position_x[n], position_y[n], position_z[n]);
     this->radius.resize (this->num_levels);
     for (int rl = 0; rl < this->num_levels; ++ rl) {
-      this->radius.at(rl) = radius[index2 (lsh, rl, n)];
+      int const ind = index2 (lsh, rl, n);
+      CCTK_REAL const rx = radius_x[ind] < 0 ? radius[ind] : radius_x[ind];
+      CCTK_REAL const ry = radius_y[ind] < 0 ? radius[ind] : radius_y[ind];
+      CCTK_REAL const rz = radius_z[ind] < 0 ? radius[ind] : radius_z[ind];
+      rvect const rad (rx, ry, rz);
+      this->radius.at(rl) = rad;
     }
     
     assert (this->num_levels <= maxreflevels);
@@ -904,7 +909,12 @@ namespace CarpetRegrid2 {
         if (do_recompose) break;
         
         // Regrid if the radii have changed sufficiently
-        CCTK_REAL const dr = abs (radius[n] - old_radius[n]);
+        CCTK_REAL const rx = radius_x[n] < 0 ? radius[n] : radius_x[n];
+        CCTK_REAL const ry = radius_y[n] < 0 ? radius[n] : radius_y[n];
+        CCTK_REAL const rz = radius_z[n] < 0 ? radius[n] : radius_z[n];
+        rvect const rad (rx, ry, rz);
+        rvect const oldrad (old_radius_x[n], old_radius_y[n], old_radius_z[n]);
+        CCTK_REAL const dr = sqrt (sum (ipow (rad - oldrad, 2)));
         CCTK_REAL mindr;
         switch (n) {
         case 0: mindr = radius_change_threshold_1; break;
@@ -962,7 +972,9 @@ namespace CarpetRegrid2 {
         old_position_y[n] = position_y[n];
         old_position_z[n] = position_z[n];
         
-        old_radius[n] = radius[n];
+        old_radius_x[n] = radius_x[n] < 0 ? radius[n] : radius_x[n];
+        old_radius_y[n] = radius_y[n] < 0 ? radius[n] : radius_y[n];
+        old_radius_z[n] = radius_z[n] < 0 ? radius[n] : radius_z[n];
       }
       
     } // if do_recompose
@@ -1069,7 +1081,12 @@ namespace CarpetRegrid2 {
         if (do_recompose) break;
         
         // Regrid if the radii have changed sufficiently
-        CCTK_REAL const dr = abs (radius[n] - old_radius[n]);
+        CCTK_REAL const rx = radius_x[n] < 0 ? radius[n] : radius_x[n];
+        CCTK_REAL const ry = radius_y[n] < 0 ? radius[n] : radius_y[n];
+        CCTK_REAL const rz = radius_z[n] < 0 ? radius[n] : radius_z[n];
+        rvect const rad (rx, ry, rz);
+        rvect const oldrad (old_radius_x[n], old_radius_y[n], old_radius_z[n]);
+        CCTK_REAL const dr = sqrt (sum (ipow (rad - oldrad, 2)));
         CCTK_REAL mindr;
         switch (n) {
         case 0: mindr = radius_change_threshold_1; break;
@@ -1156,7 +1173,9 @@ namespace CarpetRegrid2 {
         old_position_y[n] = position_y[n];
         old_position_z[n] = position_z[n];
         
-        old_radius[n] = radius[n];
+        old_radius_x[n] = radius_x[n] < 0 ? radius[n] : radius_x[n];
+        old_radius_y[n] = radius_y[n] < 0 ? radius[n] : radius_y[n];
+        old_radius_z[n] = radius_z[n] < 0 ? radius[n] : radius_z[n];
       }
       
     } // if do_recompose