aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Carpet/CarpetRegrid2/interface.ccl5
-rw-r--r--Carpet/CarpetRegrid2/param.ccl40
-rw-r--r--Carpet/CarpetRegrid2/schedule.ccl2
-rw-r--r--Carpet/CarpetRegrid2/src/initialise.cc130
-rw-r--r--Carpet/CarpetRegrid2/src/regrid.cc102
5 files changed, 139 insertions, 140 deletions
diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 7c3a4e90f..dc688d24d 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -83,6 +83,11 @@ 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_REAL positions[10] TYPE=scalar
{
position_x position_y position_z
diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl
index c97e71810..087fd9fb7 100644
--- a/Carpet/CarpetRegrid2/param.ccl
+++ b/Carpet/CarpetRegrid2/param.ccl
@@ -61,6 +61,10 @@ CCTK_INT num_levels_1 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_1 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_1 "Position of this centre" STEERABLE=always
@@ -101,6 +105,10 @@ CCTK_INT num_levels_2 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_2 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_2 "Position of this centre" STEERABLE=always
@@ -141,6 +149,10 @@ CCTK_INT num_levels_3 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_3 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_3 "Position of this centre" STEERABLE=always
@@ -181,6 +193,10 @@ CCTK_INT num_levels_4 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_4 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_4 "Position of this centre" STEERABLE=always
@@ -221,6 +237,10 @@ CCTK_INT num_levels_5 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_5 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_5 "Position of this centre" STEERABLE=always
@@ -261,6 +281,10 @@ CCTK_INT num_levels_6 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_6 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_6 "Position of this centre" STEERABLE=always
@@ -301,6 +325,10 @@ CCTK_INT num_levels_7 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_7 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_7 "Position of this centre" STEERABLE=always
@@ -341,6 +369,10 @@ CCTK_INT num_levels_8 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_8 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_8 "Position of this centre" STEERABLE=always
@@ -381,6 +413,10 @@ CCTK_INT num_levels_9 "Number of refinement levels for this centre" STEERABLE=al
0:30 :: ""
} 0
+BOOLEAN active_9 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_9 "Position of this centre" STEERABLE=always
@@ -421,6 +457,10 @@ CCTK_INT num_levels_10 "Number of refinement levels for this centre" STEERABLE=a
0:30 :: ""
} 0
+BOOLEAN active_10 "Is this region active?" STEERABLE=never
+{
+} "yes"
+
CCTK_REAL position_x_10 "Position of this centre" STEERABLE=always
diff --git a/Carpet/CarpetRegrid2/schedule.ccl b/Carpet/CarpetRegrid2/schedule.ccl
index 431e7d7d9..27aee5200 100644
--- a/Carpet/CarpetRegrid2/schedule.ccl
+++ b/Carpet/CarpetRegrid2/schedule.ccl
@@ -1,7 +1,7 @@
# Schedule definitions for thorn CarpetRegrid2
STORAGE: last_iteration last_map
-STORAGE: num_levels positions radii old_positions
+STORAGE: num_levels active positions radii old_positions
diff --git a/Carpet/CarpetRegrid2/src/initialise.cc b/Carpet/CarpetRegrid2/src/initialise.cc
index 8d1a32ce4..9b3c77694 100644
--- a/Carpet/CarpetRegrid2/src/initialise.cc
+++ b/Carpet/CarpetRegrid2/src/initialise.cc
@@ -35,111 +35,35 @@ namespace CarpetRegrid2 {
int lsh[2];
getvectorindex2 (cctkGH, "CarpetRegrid2::radii", lsh);
- if (num_centres >= 1) {
- num_levels[0] = num_levels_1;
- position_x[0] = position_x_1;
- position_y[0] = position_y_1;
- position_z[0] = position_z_1;
- for (int rl = 0; rl < num_levels[0]; ++ rl) {
- radius[index2 (lsh, rl, 0)] = radius_1[rl];
- }
- }
-
- if (num_centres >= 2) {
- num_levels[1] = num_levels_2;
- position_x[1] = position_x_2;
- position_y[1] = position_y_2;
- position_z[1] = position_z_2;
- for (int rl = 0; rl < num_levels[1]; ++ rl) {
- radius[index2 (lsh, rl, 1)] = radius_2[rl];
- }
- }
-
- if (num_centres >= 3) {
- num_levels[2] = num_levels_3;
- position_x[2] = position_x_3;
- position_y[2] = position_y_3;
- position_z[2] = position_z_3;
- for (int rl = 0; rl < num_levels[2]; ++ rl) {
- radius[index2 (lsh, rl, 2)] = radius_3[rl];
- }
- }
-
- if (num_centres >= 4) {
- num_levels[3] = num_levels_4;
- position_x[3] = position_x_4;
- position_y[3] = position_y_4;
- position_z[3] = position_z_4;
- for (int rl = 0; rl < num_levels[3]; ++ rl) {
- radius[index2 (lsh, rl, 3)] = radius_4[rl];
- }
- }
-
- if (num_centres >= 5) {
- num_levels[4] = num_levels_5;
- position_x[4] = position_x_5;
- position_y[4] = position_y_5;
- position_z[4] = position_z_5;
- for (int rl = 0; rl < num_levels[4]; ++ rl) {
- radius[index2 (lsh, rl, 4)] = radius_5[rl];
- }
- }
-
- if (num_centres >= 6) {
- num_levels[5] = num_levels_6;
- position_x[5] = position_x_6;
- position_y[5] = position_y_6;
- position_z[5] = position_z_6;
- for (int rl = 0; rl < num_levels[5]; ++ rl) {
- radius[index2 (lsh, rl, 5)] = radius_6[rl];
- }
- }
-
- if (num_centres >= 7) {
- num_levels[6] = num_levels_7;
- position_x[6] = position_x_7;
- position_y[6] = position_y_7;
- position_z[6] = position_z_7;
- for (int rl = 0; rl < num_levels[6]; ++ rl) {
- radius[index2 (lsh, rl, 6)] = radius_7[rl];
- }
- }
+#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 < num_levels[N-1]; ++ 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]; \
+ } \
+ } while (0)
- if (num_centres >= 8) {
- num_levels[7] = num_levels_8;
- position_x[7] = position_x_8;
- position_y[7] = position_y_8;
- position_z[7] = position_z_8;
- for (int rl = 0; rl < num_levels[7]; ++ rl) {
- radius[index2 (lsh, rl, 7)] = radius_8[rl];
- }
- }
+ INIT_CENTRE( 1);
+ INIT_CENTRE( 2);
+ INIT_CENTRE( 3);
+ INIT_CENTRE( 4);
+ INIT_CENTRE( 5);
+ INIT_CENTRE( 6);
+ INIT_CENTRE( 7);
+ INIT_CENTRE( 8);
+ INIT_CENTRE( 9);
+ INIT_CENTRE(10);
- if (num_centres >= 9) {
- num_levels[8] = num_levels_9;
- position_x[8] = position_x_9;
- position_y[8] = position_y_9;
- position_z[8] = position_z_9;
- for (int rl = 0; rl < num_levels[8]; ++ rl) {
- radius[index2 (lsh, rl, 8)] = radius_9[rl];
- }
- }
-
- if (num_centres >= 10) {
- num_levels[9] = num_levels_10;
- position_x[9] = position_x_10;
- position_y[9] = position_y_10;
- position_z[9] = position_z_10;
- for (int rl = 0; rl < num_levels[9]; ++ rl) {
- radius[index2 (lsh, rl, 9)] = radius_10[rl];
- }
- }
-
- for (int n = 0; n < 10; ++ n) {
- old_position_x[n] = position_x[n];
- old_position_y[n] = position_y[n];
- old_position_z[n] = position_z[n];
- }
+#undef INIT_CENTRE
if (verbose) {
for (int n = 0; n < num_centres; ++ n) {
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 6d885f5cc..da4023f6e 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -38,6 +38,7 @@ namespace CarpetRegrid2 {
struct centre_description {
int num_levels;
+ int active;
rvect position;
vector <CCTK_REAL> radius;
@@ -58,6 +59,7 @@ namespace CarpetRegrid2 {
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]);
this->radius.resize (this->num_levels);
for (int rl = 0; rl < this->num_levels; ++ rl) {
@@ -274,32 +276,34 @@ namespace CarpetRegrid2 {
// Loop over all centres
for (int n = 0; n < num_centres; ++ n) {
centre_description centre (cctkGH, n);
-
- // Loop over all levels for this centre
- for (int rl = 1; rl < centre.num_levels; ++ rl) {
-
- // Calculate a bbox for this region
- 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 imin =
- rpos2ipos (rmin, origin, scale, hh, rl);
- ivect const imax =
- rpos2ipos1 (rmax, origin, scale, hh, rl);
-
- ivect const istride = hh.baseextents.at(0).at(rl).stride();
-
- ibbox const region (imin, imax, istride);
+ if (centre.active) {
- // Add this region to the list of regions
- if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
- regions.at(rl) |= region;
- regions.at(rl).normalize();
+ // Loop over all levels for this centre
+ for (int rl = 1; rl < centre.num_levels; ++ rl) {
+
+ // Calculate a bbox for this region
+ 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 imin =
+ rpos2ipos (rmin, origin, scale, hh, rl);
+ ivect const imax =
+ rpos2ipos1 (rmax, origin, scale, hh, rl);
+
+ ivect const istride = hh.baseextents.at(0).at(rl).stride();
+
+ ibbox const region (imin, imax, istride);
+
+ // Add this region to the list of regions
+ if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
+ regions.at(rl) |= region;
+ regions.at(rl).normalize();
+
+ } // for rl
- } // for rl
-
- } // for n
+ } // if centre is active
+ } // for n
@@ -510,6 +514,22 @@ namespace CarpetRegrid2 {
min_bnd_dist_incl[0] * bb.stride()),
bb.upper())),
bb.stride());
+ if (not clipped_bb.is_contained_in (hh.baseextents.at(0).at(rl))) {
+ ostringstream msg;
+ msg << "Level " << rl << " of the refinement hierarchy is not contained in the simulation domain."
+ << " (There may be too many ghost of buffer zones.)"
+ << " One bbox is " << clipped_bb << "."
+ << " lower_is_outside_lower=" << lower_is_outside_lower
+ << " upper_is_outside_upper=" << upper_is_outside_upper
+ << " lower_is_almost_outside_upper=" << lower_is_almost_outside_upper
+ << " upper_is_almost_outside_lower=" << upper_is_almost_outside_lower
+ << " level_exterior_ilower=" << level_exterior_ilower
+ << " level_exterior_iupper=" << level_exterior_iupper
+ << " physical_ilower=" << physical_ilower
+ << " physical_iupper=" << physical_iupper
+ << " baseextent=" << hh.baseextents.at(0).at(rl);
+ CCTK_WARN (CCTK_WARN_ABORT, msg.str().c_str());
+ }
assert (clipped_bb.is_contained_in (hh.baseextents.at(0).at(rl)));
clipped |= clipped_bb;
@@ -665,12 +685,17 @@ namespace CarpetRegrid2 {
if (verbose) {
if (do_recompose) {
for (int n = 0; n < num_centres; ++ n) {
- CCTK_VInfo (CCTK_THORNSTRING,
- "Centre %d is at position [%g,%g,%g]",
- n,
- static_cast <double> (position_x[n]),
- static_cast <double> (position_y[n]),
- static_cast <double> (position_z[n]));
+ if (active[n]) {
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Centre %d is at position [%g,%g,%g]",
+ n,
+ static_cast <double> (position_x[n]),
+ static_cast <double> (position_y[n]),
+ static_cast <double> (position_z[n]));
+ } else {
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Centre %d is not active", n);
+ }
}
}
}
@@ -779,12 +804,17 @@ namespace CarpetRegrid2 {
if (verbose) {
if (do_recompose) {
for (int n = 0; n < num_centres; ++ n) {
- CCTK_VInfo (CCTK_THORNSTRING,
- "Centre %d is at position [%g,%g,%g]",
- n,
- static_cast <double> (position_x[n]),
- static_cast <double> (position_y[n]),
- static_cast <double> (position_z[n]));
+ if (active[n]) {
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Centre %d is at position [%g,%g,%g]",
+ n,
+ static_cast <double> (position_x[n]),
+ static_cast <double> (position_y[n]),
+ static_cast <double> (position_z[n]));
+ } else {
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Centre %d is not active", n);
+ }
}
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-03 02:51:02 +0000