diff options
author | Erik Schnetter <schnetter@cct.lsu.edu> | 2007-08-21 19:10:00 +0000 |
---|---|---|
committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2007-08-21 19:10:00 +0000 |
commit | f3dfa3429fef576e91f851cd8965fcf7aad6e208 (patch) | |
tree | 4d36ce643a0c37efe95fd9d9b409bce7c65db14c | |
parent | 9f33cb74c49834fc8de304019793162b8bf00e5b (diff) |
CarpetRegrid2: Add new parameter "active"
Add new parameters and grid scalars "active" for refined regions,
which specify whether a region is active. This allows switching
refined regions on and off, e.g. when a common horizon is found.
darcs-hash:20070821191037-dae7b-8f5c0ec10326b5eb876fa8ccb7e3147b6fa2ff2f.gz
-rw-r--r-- | Carpet/CarpetRegrid2/interface.ccl | 5 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/param.ccl | 40 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/schedule.ccl | 2 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/src/initialise.cc | 130 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/src/regrid.cc | 102 |
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); + } } } } |