diff options
Diffstat (limited to 'Carpet/Carpet/src/Recompose.cc')
| -rw-r--r-- | Carpet/Carpet/src/Recompose.cc | 1037 |
1 files changed, 657 insertions, 380 deletions
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index 958a5d058..f1f0b8834 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -15,21 +15,25 @@ #include <sys/stat.h> #include <sys/types.h> -#include "cctk.h" -#include "cctk_Parameters.h" +#include <mpi.h> -#include "bbox.hh" -#include "bboxset.hh" -#include "defs.hh" -#include "dh.hh" -#include "gh.hh" -#include "region.hh" -#include "vect.hh" +#include <cctk.h> +#include <cctk_Parameters.h> -#include "carpet.hh" -#include "modes.hh" +#include <loopcontrol.h> -#define DEBUG false // false or true +#include <bbox.hh> +#include <bboxset.hh> +#include <defs.hh> +#include <dh.hh> +#include <gh.hh> +#include <region.hh> +#include <vect.hh> + +#include <carpet.hh> +#include <modes.hh> +#include <variables.hh> +#include <Timers.hh> @@ -61,6 +65,11 @@ namespace Carpet { static void + ClassifyPoints (cGH const * cctkGH, int rl); + + + + static void SplitRegionsMaps_Automatic_Recursively (bvect const & dims, int const firstproc, int const nprocs, @@ -76,57 +85,68 @@ namespace Carpet { void CheckRegions (gh::mregs const & regsss) { + char const * const where = "Carpet::CheckRegions"; + static Timer timer (where); + timer.start(); + // At least one multigrid level if (regsss.size() == 0) { CCTK_WARN (0, "I cannot set up a grid hierarchy with zero multigrid levels."); } assert (regsss.size() > 0); // At least one refinement level - if (regsss.at(0).size() == 0) { + if (regsss.AT(0).size() == 0) { CCTK_WARN (0, "I cannot set up a grid hierarchy with zero refinement levels."); } - assert (regsss.at(0).size() > 0); + assert (regsss.AT(0).size() > 0); // At most maxreflevels levels - if ((int)regsss.at(0).size() > maxreflevels) { + if ((int)regsss.AT(0).size() > maxreflevels) { CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, "I cannot set up a grid hierarchy with more than Carpet::max_refinement_levels refinement levels. I found Carpet::max_refinement_levels=%d, while %d levels were requested.", - (int)maxreflevels, (int)regsss.at(0).size()); + (int)maxreflevels, (int)regsss.AT(0).size()); } - assert ((int)regsss.at(0).size() <= maxreflevels); + assert ((int)regsss.AT(0).size() <= maxreflevels); for (int ml=0; ml<(int)regsss.size(); ++ml) { int num_regions = 0; - for (int rl=0; rl<(int)regsss.at(0).size(); ++rl) { + for (int rl=0; rl<(int)regsss.AT(0).size(); ++rl) { // No empty levels // (but allow some empty maps) - // assert (regsss.at(ml).at(rl).size() > 0); - num_regions += regsss.at(ml).at(rl).size(); - for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) { + // assert (regsss.AT(ml).AT(rl).size() > 0); + num_regions += regsss.AT(ml).AT(rl).size(); + for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) { // Check sizes // (but allow processors with zero grid points) - // assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <= - // regsss.at(rl).at(c).at(ml).extent.upper())); + // assert (all(regsss.AT(rl).AT(c).AT(ml).extent.lower() <= + // regsss.AT(rl).AT(c).AT(ml).extent.upper())); // Check strides ivect const str = - (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml)); - assert (all(regsss.at(ml).at(rl).at(c).extent.stride() % str == 0)); + (maxspacereflevelfact / spacereffacts.AT(rl) * ipow(mgfact, ml)); + assert (all(regsss.AT(ml).AT(rl).AT(c).extent.stride() % str == 0)); // Check alignments - assert (all(regsss.at(ml).at(rl).at(c).extent.lower() % str == 0)); - assert (all(regsss.at(ml).at(rl).at(c).extent.upper() % str == 0)); + assert (all(regsss.AT(ml).AT(rl).AT(c).extent.lower() % str == 0)); + assert (all(regsss.AT(ml).AT(rl).AT(c).extent.upper() % str == 0)); } } // No empty levels assert (num_regions > 0); } + + timer.stop(); } bool Regrid (cGH const * const cctkGH, - bool const force_recompose) + bool const force_recompose, + bool const do_init) { DECLARE_CCTK_PARAMETERS; + char const * const where = "Carpet::Regrid"; + static Timer timer (where); + timer.start(); + Checkpoint ("Regridding level %d...", reflevel); assert (is_level_mode()); @@ -141,6 +161,7 @@ namespace Carpet { CCTK_WARN (2, "The regridding routine Carpet_Regrid has not been provided. There will be no regridding. Maybe you forgot to activate a regridding thorn?"); didtell = true; } + timer.stop(); return false; } @@ -150,7 +171,7 @@ namespace Carpet { CCTK_WARN (2, "The regridding routine Carpet_RegridMaps has not been provided. Regridding will be performed in singlemap mode instead of level mode."); didtell = true; } - } + } @@ -160,7 +181,7 @@ namespace Carpet { BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - gh::rregs superregss = vhh.at(map)->superregions; + gh::rregs superregss = vhh.AT(map)->superregions; gh::mregs regsss; // Check whether to recompose @@ -170,7 +191,7 @@ namespace Carpet { did_change = did_change or do_recompose; if (do_recompose) { - RegridMap (cctkGH, map, superregss, regsss); + RegridMap (cctkGH, map, superregss, regsss, do_init); } } END_MAP_LOOP; @@ -180,20 +201,28 @@ namespace Carpet { vector<gh::rregs> superregsss (maps); vector<gh::mregs> regssss (maps); BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - superregsss.at(map) = vhh.at(map)->superregions; + superregsss.AT(map) = vhh.AT(map)->superregions; } END_MAP_LOOP; // Check whether to recompose CCTK_INT const do_recompose = Carpet_RegridMaps (cctkGH, & superregsss, & regssss, force_recompose); assert (do_recompose >= 0); +#warning "TODO" +#if 1 // #ifdef CARPET_DEBUG + { + int ival = do_recompose; + MPI_Bcast (& ival, 1, MPI_INT, 0, dist::comm()); + assert (ival == do_recompose); + } +#endif did_change = did_change or do_recompose; if (do_recompose) { BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - gh::rregs const & superregss = superregsss.at(map); - gh::mregs const & regsss = regssss.at(map); - RegridMap (cctkGH, map, superregss, regsss); + gh::rregs const & superregss = superregsss.AT(map); + gh::mregs const & regsss = regssss.AT(map); + RegridMap (cctkGH, map, superregss, regsss, do_init); } END_MAP_LOOP; } @@ -208,6 +237,14 @@ namespace Carpet { } // if did change + + timer.stop(); + + if (enable_no_storage) { + CCTK_WARN (CCTK_WARN_ALERT, + "Carpet completed its internal setup, and would now normally go on to allocate memory. Since the parameter Carpet::enable_no_storage has been set, Carpet will exit instead."); + CCTK_Exit (cctkGH, 0); + } return did_change; } @@ -218,12 +255,16 @@ namespace Carpet { RegridMap (cGH const * const cctkGH, int const m, gh::rregs const & superregss, - gh::mregs const & regsss) + gh::mregs const & regsss, + bool const do_init) { DECLARE_CCTK_PARAMETERS; + char const * const where = "Carpet::RegridMap"; + static Timer timer (where); + timer.start(); + Waypoint ("Regridding map %d...", m); - #warning "TODO: keep levels fixed here" #if 0 @@ -239,15 +280,15 @@ namespace Carpet { int const do_every = use_tapered_grids ? - maxtimereflevelfact / timereffacts.at(max(0,rl-1)): - maxtimereflevelfact / timereffacts.at( rl ); + maxtimereflevelfact / timereffacts.AT(max(0,rl-1)): + maxtimereflevelfact / timereffacts.AT( rl ); bool const regrid_this_level = (cctkGH->cctk_iteration - 1) % do_every == 0; if (not regrid_this_level) { // Set regions from current grid structure - regions.at(rl) = ...; + regions.AT(rl) = ...; } } #endif @@ -258,18 +299,18 @@ namespace Carpet { // not change // Regrid - vhh.at(m)->regrid (superregss, regsss); + vhh.AT(m)->regrid (superregss, regsss, do_init); + + // Output grid structure to screen + OutputSuperregions (cctkGH, m, * vhh.AT(m), superregss); + OutputGrids (cctkGH, m, * vhh.AT(m), * vdd.AT(m)); // Write grid structure to file OutputGridStructure (cctkGH, m, regsss); OutputGridCoordinates (cctkGH, m, regsss); -#warning "TODO: output superregss" - - if (verbose or veryverbose) { - OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m)); - } Waypoint ("Done regridding map %d.", m); + timer.stop(); } @@ -279,13 +320,27 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; + char const * const where = "Carpet::PostRegrid"; + static Timer timer (where); + timer.start(); + // Calculate new number of levels int const oldreflevels = reflevels; - reflevels = vhh.at(0)->reflevels(); + reflevels = vhh.AT(0)->reflevels(); for (int m=0; m<maps; ++m) { - assert (vhh.at(m)->reflevels() == reflevels); + assert (vhh.AT(m)->reflevels() == reflevels); } +#warning "TODO" +#if 1 // #ifdef CARPET_DEBUG + { + // All processes must use the same number of levels + int ival = reflevels; + MPI_Bcast (& ival, 1, MPI_INT, 0, dist::comm()); + assert (ival == reflevels); + } +#endif + // One cannot switch off the current level assert (reflevels > reflevel); @@ -294,19 +349,23 @@ namespace Carpet { int const grouptype = CCTK_GroupTypeI (n); if (grouptype == CCTK_GF) { for (int ml=0; ml<mglevels; ++ml) { - groupdata.at(n).activetimelevels.at(ml).resize - (reflevels, groupdata.at(n).activetimelevels.at(ml).at(0)); + groupdata.AT(n).activetimelevels.AT(ml).resize + (reflevels, groupdata.AT(n).activetimelevels.AT(ml).AT(0)); } } } // Set new level times for (int ml=0; ml<mglevels; ++ml) { - leveltimes.at(ml).resize - (reflevels, leveltimes.at(ml).at(oldreflevels-1)); + leveltimes.AT(ml).resize + (reflevels, leveltimes.AT(ml).AT(oldreflevels-1)); } + ++ regridding_epoch; + OutputGridStatistics (cctkGH); + + timer.stop(); } @@ -316,200 +375,291 @@ namespace Carpet { int const rl, bool const do_init) { + char const * const where = "Carpet::Recompose"; + static Timer timer (where); + timer.start(); + bool did_recompose = false; for (int m=0; m<maps; ++m) { Waypoint ("Recomposing the grid hierarchy for map %d level %d...", m, rl); - assert (rl>=0 and rl<vhh.at(m)->reflevels()); - did_recompose |= vhh.at(m)->recompose (rl, do_init); + assert (rl>=0 and rl<vhh.AT(m)->reflevels()); + did_recompose |= vhh.AT(m)->recompose (rl, do_init); Waypoint ("Done recomposing the grid hierarchy for map %d level %d.", m, rl); } + + ClassifyPoints (cctkGH, rl); + + timer.stop(); return did_recompose; } void - OutputGrids (cGH const * const cctkGH, - int const m, - gh const & hh, - dh const & dd) + RegridFree (cGH const * const cctkGH, + bool const do_init) { - CCTK_INFO ("Grid structure (grid points):"); - for (int ml=0; ml<hh.mglevels(); ++ml) { - for (int rl=0; rl<hh.reflevels(); ++rl) { - for (int c=0; c<hh.components(rl); ++c) { - const int convfact = ipow(mgfact, ml); - const ivect levfact = spacereffacts.at(rl); - const ibbox ext = hh.extent(ml,rl,c); - const ivect & lower = ext.lower(); - const ivect & upper = ext.upper(); - const ivect & stride = ext.stride(); - assert (all(lower * levfact % maxspacereflevelfact == 0)); - assert (all(upper * levfact % maxspacereflevelfact == 0)); - assert (all(((upper - lower) * levfact / maxspacereflevelfact) - % convfact == 0)); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " exterior: " - << "proc " - << hh.processor(rl,c) - << " " - << lower / stride - << " : " - << upper / stride - << " (" - << (upper - lower) / stride + 1 - << ") " - << prod ((upper - lower) / stride + 1) - << endl; - } - } + char const * const where = "Carpet::RegridFree"; + static Timer timer (where); + timer.start(); + + Checkpoint ("Freeing after regridding level %d...", reflevel); + + assert (is_level_mode()); + + // Free old grid structure + for (int m=0; m<maps; ++m) { + vhh.AT(m)->regrid_free (do_init); } - CCTK_INFO ("Grid structure (boundaries):"); - for (int rl=0; rl<hh.reflevels(); ++rl) { - for (int c=0; c<hh.components(rl); ++c) { + timer.stop(); + } + + + + void + OutputSuperregions (cGH const * const cctkGH, + int const m, + gh const & hh, + gh::rregs const & superregss) + { + CCTK_INFO ("Grid structure (superregions, grid points):"); + for (int rl=0; rl<(int)superregss.size(); ++rl) { + assert (rl < hh.reflevels()); + for (int c=0; c<(int)superregss.AT(rl).size(); ++c) { + const ivect & levfact = spacereffacts.AT(rl); + const ibbox & ext = superregss.AT(rl).AT(c).extent; + const ivect & lower = ext.lower(); + const ivect & upper = ext.upper(); + const ivect & stride = ext.stride(); + assert (all(lower * levfact % maxspacereflevelfact == 0)); + assert (all(upper * levfact % maxspacereflevelfact == 0)); cout << " [" << rl << "][" << m << "][" << c << "]" - << " bbox: " - << hh.outer_boundaries(rl,c) - << endl; + << " exterior: " + << lower / stride + << " : " + << upper / stride + << " (" + << (upper - lower) / stride + 1 + << ") " + << prod ((upper - lower) / stride + 1) + << eol; } } - CCTK_INFO ("Grid structure (coordinates):"); - for (int ml=0; ml<hh.mglevels(); ++ml) { - for (int rl=0; rl<hh.reflevels(); ++rl) { - for (int c=0; c<hh.components(rl); ++c) { - const rvect origin = domainspecs.at(m).exterior_min; - const rvect delta = (domainspecs.at(m).exterior_max - domainspecs.at(m).exterior_min) / rvect (domainspecs.at(m).npoints - 1); - const ibbox ext = hh.extent(ml,rl,c); - const ivect & lower = ext.lower(); - const ivect & upper = ext.upper(); - const int convfact = ipow(mgfact, ml); - const ivect levfact = spacereffacts.at(rl); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " exterior: " - << origin + delta * rvect(lower) / rvect(maxspacereflevelfact) - << " : " - << origin + delta * rvect(upper) / rvect(maxspacereflevelfact) - << " : " - << delta * rvect(convfact) / rvect(levfact) << endl; - } + CCTK_INFO ("Grid structure (superregions, coordinates):"); + for (int rl=0; rl<(int)superregss.size(); ++rl) { + assert (rl < hh.reflevels()); + for (int c=0; c<(int)superregss.AT(rl).size(); ++c) { + const rvect origin = domainspecs.AT(m).exterior_min; + const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1); + const ibbox & ext = superregss.AT(rl).AT(c).extent; + const ivect & lower = ext.lower(); + const ivect & upper = ext.upper(); + const ivect & levfact = spacereffacts.AT(rl); + cout << " [" << rl << "][" << m << "][" << c << "]" + << " exterior: " + << origin + delta * rvect(lower) / rvect(maxspacereflevelfact) + << " : " + << origin + delta * rvect(upper) / rvect(maxspacereflevelfact) + << " : " + << delta / rvect(levfact) << eol; } } - CCTK_INFO ("Grid structure (coordinates, including ghosts):"); - for (int ml=0; ml<hh.mglevels(); ++ml) { - for (int rl=0; rl<hh.reflevels(); ++rl) { - for (int c=0; c<hh.components(rl); ++c) { - const rvect origin = domainspecs.at(m).exterior_min; - const rvect delta = (domainspecs.at(m).exterior_max - domainspecs.at(m).exterior_min) / rvect (domainspecs.at(m).npoints - 1); - const ivect lower = dd.boxes.at(ml).at(rl).at(c).exterior.lower(); - const ivect upper = dd.boxes.at(ml).at(rl).at(c).exterior.upper(); - const int convfact = ipow(mgfact, ml); - const ivect levfact = spacereffacts.at(rl); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " exterior: " - << origin + delta * rvect(lower) / rvect(maxspacereflevelfact) - << " : " - << origin + delta * rvect(upper) / rvect(maxspacereflevelfact) - << " : " - << delta * rvect(convfact) / rvect(levfact) << endl; + fflush (stdout); + } + + + + void + OutputGrids (cGH const * const cctkGH, + int const m, + gh const & hh, + dh const & dd) + { + DECLARE_CCTK_PARAMETERS; + + if (verbose or veryverbose) { + + CCTK_INFO ("Grid structure (grid points):"); + for (int ml=0; ml<hh.mglevels(); ++ml) { + for (int rl=0; rl<hh.reflevels(); ++rl) { + for (int c=0; c<hh.components(rl); ++c) { + const int convfact = ipow(mgfact, ml); + const ivect levfact = spacereffacts.AT(rl); + const ibbox ext = hh.extent(ml,rl,c); + const ivect & lower = ext.lower(); + const ivect & upper = ext.upper(); + const ivect & stride = ext.stride(); + assert (all(lower * levfact % maxspacereflevelfact == 0)); + assert (all(upper * levfact % maxspacereflevelfact == 0)); + assert (all(((upper - lower) * levfact / maxspacereflevelfact) + % convfact == 0)); + cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" + << " exterior: " + << "proc " + << hh.processor(rl,c) + << " " + << lower / stride + << " : " + << upper / stride + << " (" + << (upper - lower) / stride + 1 + << ") " + << prod ((upper - lower) / stride + 1) + << eol; + } } } - } - - CCTK_INFO ("Grid statistics:"); - const int oldprecision = cout.precision(); - const ios_base::fmtflags oldflags = cout.flags(); - cout.setf (ios::fixed); - for (int ml=0; ml<hh.mglevels(); ++ml) { - CCTK_REAL coarsevolume = 0; + + CCTK_INFO ("Grid structure (boundaries):"); for (int rl=0; rl<hh.reflevels(); ++rl) { - - const CCTK_REAL basevolume = hh.baseextents.AT(0).AT(0).size(); - CCTK_REAL countvolume = 0; - CCTK_REAL totalvolume = 0; - CCTK_REAL totalvolume2 = 0; - for (int c=0; c<hh.components(rl); ++c) { - const CCTK_REAL volume = hh.extent(ml,rl,c).size(); - ++ countvolume; - totalvolume += volume; - totalvolume2 += ipow(volume, 2); + cout << " [" << rl << "][" << m << "][" << c << "]" + << " bbox: " + << hh.outer_boundaries(rl,c) + << eol; } - - const CCTK_REAL avgvolume = totalvolume / countvolume; - const CCTK_REAL stddevvolume - = sqrt (max (CCTK_REAL(0), - totalvolume2 / countvolume - ipow (avgvolume, 2))); - - for (int c=0; c<hh.components(rl); ++c) { - const CCTK_REAL volume = hh.extent(ml,rl,c).size(); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " volume: " << setprecision(0) << volume - << " of parent: " << setprecision(1) << 100 * volume / totalvolume << "%" - << " of domain: " << setprecision(1) << 100 * volume / basevolume << "%" - << endl; + } + + CCTK_INFO ("Grid structure (coordinates):"); + for (int ml=0; ml<hh.mglevels(); ++ml) { + for (int rl=0; rl<hh.reflevels(); ++rl) { + for (int c=0; c<hh.components(rl); ++c) { + const rvect origin = domainspecs.AT(m).exterior_min; + const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1); + const ibbox ext = hh.extent(ml,rl,c); + const ivect & lower = ext.lower(); + const ivect & upper = ext.upper(); + const int convfact = ipow(mgfact, ml); + const ivect levfact = spacereffacts.AT(rl); + cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" + << " exterior: " + << origin + delta * rvect(lower) / rvect(maxspacereflevelfact) + << " : " + << origin + delta * rvect(upper) / rvect(maxspacereflevelfact) + << " : " + << delta * rvect(convfact) / rvect(levfact) << eol; + } } - - cout << " [" << ml << "][" << rl << "][" << m << "]" - << " average volume: " << setprecision(0) << avgvolume - << " of parent: " << setprecision(1) << 100 * avgvolume / totalvolume << "%" - << " of domain: " << setprecision(1) << 100 * avgvolume / basevolume << "%" - << endl; - cout << " [" << ml << "][" << rl << "][" << m << "]" - << " standard deviation: " << setprecision(0) << stddevvolume - << " of parent: " << setprecision(1) << 100 * stddevvolume / totalvolume << "%" - << " of domain: " << setprecision(1) << 100 * stddevvolume / basevolume << "%" - << endl; - - // TODO: ghost points vs. interior points (and boundary - // points) - - CCTK_REAL countquadrupole = 0; - CCTK_REAL minquadrupole = 1; - CCTK_REAL totalquadrupole = 0; - CCTK_REAL totalquadrupole2 = 0; - - for (int c=0; c<hh.components(rl); ++c) { - const ibbox ext = hh.extent(ml,rl,c); - const ivect shape = ext.shape(); - const ivect stride = ext.stride(); - const CCTK_REAL minlength = minval (rvect (shape) / rvect (stride)); - const CCTK_REAL maxlength = maxval (rvect (shape) / rvect (stride)); - const CCTK_REAL quadrupole = minlength / maxlength; - ++ countquadrupole; - minquadrupole = min (minquadrupole, quadrupole); - totalquadrupole += quadrupole; - totalquadrupole2 += ipow (quadrupole, 2); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " length ratio: " << setprecision(3) << quadrupole - << endl; + } + + CCTK_INFO ("Grid structure (coordinates, including ghosts):"); + for (int ml=0; ml<hh.mglevels(); ++ml) { + for (int rl=0; rl<hh.reflevels(); ++rl) { + for (int c=0; c<hh.components(rl); ++c) { + const rvect origin = domainspecs.AT(m).exterior_min; + const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1); + const ivect lower = dd.boxes.AT(ml).AT(rl).AT(c).exterior.lower(); + const ivect upper = dd.boxes.AT(ml).AT(rl).AT(c).exterior.upper(); + const int convfact = ipow(mgfact, ml); + const ivect levfact = spacereffacts.AT(rl); + cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" + << " exterior: " + << origin + delta * rvect(lower) / rvect(maxspacereflevelfact) + << " : " + << origin + delta * rvect(upper) / rvect(maxspacereflevelfact) + << " : " + << delta * rvect(convfact) / rvect(levfact) << eol; + } + } + } + + CCTK_INFO ("Grid statistics:"); + const int oldprecision = cout.precision(); + const ios_base::fmtflags oldflags = cout.flags(); + cout.setf (ios::fixed); + for (int ml=0; ml<hh.mglevels(); ++ml) { + CCTK_REAL coarsevolume = 0; + for (int rl=0; rl<hh.reflevels(); ++rl) { + + const CCTK_REAL basevolume = hh.baseextents.AT(0).AT(0).size(); + CCTK_REAL countvolume = 0; + CCTK_REAL totalvolume = 0; + CCTK_REAL totalvolume2 = 0; + + for (int c=0; c<hh.components(rl); ++c) { + const CCTK_REAL volume = hh.extent(ml,rl,c).size(); + ++ countvolume; + totalvolume += volume; + totalvolume2 += ipow(volume, 2); + } + + const CCTK_REAL avgvolume = totalvolume / countvolume; + const CCTK_REAL stddevvolume + = sqrt (max (CCTK_REAL(0), + totalvolume2 / countvolume - ipow (avgvolume, 2))); + + for (int c=0; c<hh.components(rl); ++c) { + const CCTK_REAL volume = hh.extent(ml,rl,c).size(); + cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" + << " volume: " << setprecision(0) << volume + << " of parent: " << setprecision(1) << 100 * volume / totalvolume << "%" + << " of domain: " << setprecision(1) << 100 * volume / basevolume << "%" + << eol; + } + + cout << " [" << ml << "][" << rl << "][" << m << "]" + << " average volume: " << setprecision(0) << avgvolume + << " of parent: " << setprecision(1) << 100 * avgvolume / totalvolume << "%" + << " of domain: " << setprecision(1) << 100 * avgvolume / basevolume << "%" + << eol; + cout << " [" << ml << "][" << rl << "][" << m << "]" + << " standard deviation: " << setprecision(0) << stddevvolume + << " of parent: " << setprecision(1) << 100 * stddevvolume / totalvolume << "%" + << " of domain: " << setprecision(1) << 100 * stddevvolume / basevolume << "%" + << eol; + + // TODO: ghost points vs. interior points (and boundary + // points) + + CCTK_REAL countquadrupole = 0; + CCTK_REAL minquadrupole = 1; + CCTK_REAL totalquadrupole = 0; + CCTK_REAL totalquadrupole2 = 0; + + for (int c=0; c<hh.components(rl); ++c) { + const ibbox ext = hh.extent(ml,rl,c); + const ivect shape = ext.shape(); + const ivect stride = ext.stride(); + const CCTK_REAL minlength = minval (rvect (shape) / rvect (stride)); + const CCTK_REAL maxlength = maxval (rvect (shape) / rvect (stride)); + const CCTK_REAL quadrupole = minlength / maxlength; + ++ countquadrupole; + minquadrupole = min (minquadrupole, quadrupole); + totalquadrupole += quadrupole; + totalquadrupole2 += ipow (quadrupole, 2); + cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" + << " length ratio: " << setprecision(3) << quadrupole + << eol; + } + + const CCTK_REAL avgquadrupole = totalquadrupole / countquadrupole; + const CCTK_REAL stddevquadrupole + = sqrt (max (CCTK_REAL(0), + (totalquadrupole2 / countquadrupole + - ipow (avgquadrupole, 2)))); + + cout << " [" << ml << "][" << rl << "][" << m << "]" + << " average length ratio: " << setprecision(3) << avgquadrupole + << " standard deviation: " << setprecision(3) << stddevquadrupole + << eol; + + // TODO: processor distribution, average load, std deviation + + coarsevolume = totalvolume * prod (rvect (spacereflevelfact)); } - - const CCTK_REAL avgquadrupole = totalquadrupole / countquadrupole; - const CCTK_REAL stddevquadrupole - = sqrt (max (CCTK_REAL(0), - (totalquadrupole2 / countquadrupole - - ipow (avgquadrupole, 2)))); - - cout << " [" << ml << "][" << rl << "][" << m << "]" - << " average length ratio: " << setprecision(3) << avgquadrupole - << " standard deviation: " << setprecision(3) << stddevquadrupole - << endl; - - // TODO: processor distribution, average load, std deviation - - coarsevolume = totalvolume * prod (rvect (spacereflevelfact)); } + cout.precision (oldprecision); + cout.setf (oldflags); + + fflush (stdout); } - cout.precision (oldprecision); - cout.setf (oldflags); - } @@ -561,14 +711,14 @@ namespace Carpet { file << "maps " << maps << eol; file << m << " mglevels " << regsss.size() << eol; for (int ml=0; ml<(int)regsss.size(); ++ml) { - file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol; - for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) { - file << m << " " << ml << " " << rl << " components " << regsss.at(ml).at(rl).size() << eol; - for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) { + file << m << " " << ml << " reflevels " << regsss.AT(ml).size() << eol; + for (int rl=0; rl<(int)regsss.AT(ml).size(); ++rl) { + file << m << " " << ml << " " << rl << " components " << regsss.AT(ml).AT(rl).size() << eol; + for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) { file << m << " " << ml << " " << rl << " " << c << " " - << regsss.at(ml).at(rl).at(c).processor << " " - << regsss.at(ml).at(rl).at(c).extent << " " - << regsss.at(ml).at(rl).at(c).outer_boundaries << eol; + << regsss.AT(ml).AT(rl).AT(c).processor << " " + << regsss.AT(ml).AT(rl).AT(c).extent << " " + << regsss.AT(ml).AT(rl).AT(c).outer_boundaries << eol; } } } @@ -668,7 +818,7 @@ namespace Carpet { // Affine transformation between index space and coordinate space rvect const origin = exterior_lower; rvect const scale = - rvect (vhh.at(m)->baseextents.at(0).at(0).stride()) / spacing; + rvect (vhh.AT(m)->baseextents.AT(0).AT(0).stride()) / spacing; @@ -676,11 +826,11 @@ namespace Carpet { file << "maps " << maps << eol; file << m << " mglevels " << regsss.size() << eol; for (int ml=0; ml<(int)regsss.size(); ++ml) { - file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol; - for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) { + file << m << " " << ml << " reflevels " << regsss.AT(ml).size() << eol; + for (int rl=0; rl<(int)regsss.AT(ml).size(); ++rl) { ibset extents; - for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) { - extents += regsss.at(ml).at(rl).at(c).extent; + for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) { + extents += regsss.AT(ml).AT(rl).AT(c).extent; } extents.normalize(); file << m << " " << ml << " " << rl << " regions " << extents.setsize() << eol; @@ -690,12 +840,12 @@ namespace Carpet { { #if 0 ibbox const & ext = * bi; - ibbox const & baseext = vhh.at(m)->baseextents.at(ml).at(rl); - ibbox const & coarseext = vhh.at(m)->baseextents.at(ml).at(0 ); + ibbox const & baseext = vhh.AT(m)->baseextents.AT(ml).AT(rl); + ibbox const & coarseext = vhh.AT(m)->baseextents.AT(ml).AT(0 ); // This is nice, but is wrong since CartGrid3D has not yet // initialised the coordinates - ivect const cctk_levfac = spacereffacts.at(rl); + ivect const cctk_levfac = spacereffacts.AT(rl); ivect const cctk_levoff = baseext.lower() - coarseext.lower(); ivect const cctk_levoffdenom = baseext.stride(); @@ -705,9 +855,9 @@ namespace Carpet { (ext.upper() - baseext.lower()) / ext.stride(); rvect const cctk_origin_space = - origin_space.at(m).at(ml); + origin_space.AT(m).AT(ml); rvect const cctk_delta_space = - delta_space.at(m) * rvect (ipow (mgfact, ml)); + delta_space.AT(m) * rvect (ipow (mgfact, ml)); rvect const CCTK_ORIGIN_SPACE = cctk_origin_space + @@ -773,28 +923,35 @@ namespace Carpet { OutputGridStatistics (cGH const * const cctkGH) { // Grid array statistics - int num_gfs = 0; int num_arrays = 0; + int num_gfs = 0; + int size_gfs = 0; CCTK_REAL num_active_array_points = 0; CCTK_REAL num_total_array_points = 0; + CCTK_REAL size_total_array_points = 0; for (int g=0; g<CCTK_NumGroups(); ++g) { + cGroup gdata; + check (not CCTK_GroupData (g, &gdata)); int const num_tl = CCTK_ActiveTimeLevelsGI (cctkGH, g); - switch (CCTK_GroupTypeI(g)) { + int const num_vars = gdata.numvars; + int const size_vars = gdata.numvars * CCTK_VarTypeSize (gdata.vartype); + switch (gdata.grouptype) { case CCTK_SCALAR: case CCTK_ARRAY: { - int const num_vars = CCTK_NumVarsInGroupI (g); num_arrays += num_tl * num_vars; gh const * const hh = arrdata.AT(g).AT(0).hh; dh const * const dd = arrdata.AT(g).AT(0).dd; for (int c=0; c<hh->components(0); ++c) { dh::dboxes const & b = dd->boxes.AT(0).AT(0).AT(c); - num_active_array_points += num_tl * num_vars * b.active.size(); - num_total_array_points += num_tl * num_vars * b.exterior.size(); + num_active_array_points += num_tl * num_vars * b.active_size; + num_total_array_points += num_tl * num_vars * b.exterior_size; + size_total_array_points += num_tl * size_vars * b.exterior_size; } break; } case CCTK_GF: - num_gfs += num_tl; + num_gfs += num_tl * num_vars; + size_gfs += num_tl * size_vars; break; default: assert (0); @@ -807,6 +964,7 @@ namespace Carpet { CCTK_REAL num_active_mem_points = 0; CCTK_REAL num_owned_mem_points = 0; CCTK_REAL num_total_mem_points = 0; + CCTK_REAL size_total_mem_points = 0; CCTK_REAL num_active_cpu_points = 0; CCTK_REAL num_owned_cpu_points = 0; CCTK_REAL num_total_cpu_points = 0; @@ -822,12 +980,13 @@ namespace Carpet { for (int c=0; c<hh->components(rl); ++c) { ++ num_comps; dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c); - num_active_mem_points += num_gfs * b.active.size(); - num_owned_mem_points += num_gfs * b.owned.size(); - num_total_mem_points += num_gfs * b.exterior.size(); - num_active_cpu_points += trf * b.active.size(); - num_owned_cpu_points += trf * b.owned.size(); - num_total_cpu_points += trf * b.exterior.size(); + num_active_mem_points += num_gfs * b.active_size; + num_owned_mem_points += num_gfs * b.owned_size; + num_total_mem_points += num_gfs * b.exterior_size; + size_total_mem_points += size_gfs * b.exterior_size; + num_active_cpu_points += trf * b.active_size; + num_owned_cpu_points += trf * b.owned_size; + num_total_cpu_points += trf * b.exterior_size; } } } @@ -841,27 +1000,30 @@ namespace Carpet { "Grid structure statistics:"); CCTK_VInfo (CCTK_THORNSTRING, "GF: rhs: %.0fk active, %.0fk owned (+%.0f%%), %.0fk total (+%.0f%%), %.3g steps/time", - double (num_active_cpu_points / 1000), - double (num_owned_cpu_points / 1000), + double (num_active_cpu_points / 1e+3), + double (num_owned_cpu_points / 1e+3), double (num_owned_cpu_points / num_active_cpu_points * 100 - 100), - double (num_total_cpu_points / 1000), + double (num_total_cpu_points / 1e+3), double (num_total_cpu_points / num_owned_cpu_points * 100 - 100), double (num_steps / delta_time)); CCTK_VInfo (CCTK_THORNSTRING, "GF: vars: %d, pts: %.0fM active, %.0fM owned (+%.0f%%), %.0fM total (+%.0f%%), %.1f comp/proc", num_gfs, - double (num_active_mem_points / 1000000), - double (num_owned_mem_points / 1000000), + double (num_active_mem_points / 1e+6), + double (num_owned_mem_points / 1e+6), double (num_owned_mem_points / num_active_mem_points * 100 - 100), - double (num_total_mem_points / 1000000), + double (num_total_mem_points / 1e+6), double (num_total_mem_points / num_owned_mem_points * 100 - 100), double (num_comps / (reflevels * CCTK_nProcs (cctkGH)))); CCTK_VInfo (CCTK_THORNSTRING, "GA: vars: %d, pts: %.0fM active, %.0fM total (+%.0f%%)", num_arrays, - double (num_active_array_points / 1000000), - double (num_total_array_points / 1000000), + double (num_active_array_points / 1e+6), + double (num_total_array_points / 1e+6), double (num_total_array_points / num_active_array_points * 100 - 100)); + CCTK_VInfo (CCTK_THORNSTRING, + "Total required memory: %.3f GByte (for GAs and currently active GFs)", + double ((size_total_array_points + size_total_mem_points) / 1e+9)); // After this, we will begin to allocate memory for the grid // structure. If we run out of memory, ensure that this output @@ -902,16 +1064,16 @@ namespace Carpet { regs = superregs; if (nprocs == 1) { - regs.at(0).processor = 0; - pseudoregion_t const preg (regs.at(0).extent, regs.at(0).processor); - assert (superregs.at(0).processors == NULL); - superregs.at(0).processors = new ipfulltree (preg); + regs.AT(0).processor = 0; + pseudoregion_t const preg (regs.AT(0).extent, regs.AT(0).processor); + assert (superregs.AT(0).processors == NULL); + superregs.AT(0).processors = new ipfulltree (preg); return; } assert (dir>=0 and dir<dim); - region_t const & reg0 = regs.at(0); + region_t const & reg0 = regs.AT(0); const ivect rstr0 = reg0.extent.stride(); const ivect rlb0 = reg0.extent.lower(); const ivect rub0 = reg0.extent.upper() + rstr0; @@ -933,7 +1095,7 @@ namespace Carpet { if (cub[dir] > rub0[dir]) cub[dir] = rub0[dir]; assert (clb[dir] <= cub[dir]); assert (cub[dir] <= rub0[dir]); - region_t & reg = regs.at(c); + region_t & reg = regs.AT(c); ibbox & ext = reg.extent; b2vect & obnd = reg.outer_boundaries; int & proc = reg.processor; @@ -943,16 +1105,16 @@ namespace Carpet { obnd[1][dir] &= cub[dir] == rub0[dir]; proc = c; pseudoregion_t const preg (reg.extent, c); - bounds.at(c) = reg.extent.lower()[dir]; - subtrees.at(c) = new ipfulltree (preg); + bounds.AT(c) = reg.extent.lower()[dir]; + subtrees.AT(c) = new ipfulltree (preg); } - bounds.at(nprocs) = rub0[dir] + rstr0[dir]; + bounds.AT(nprocs) = rub0[dir] + rstr0[dir]; - assert (superregs.at(0).processors == NULL); - superregs.at(0).processors = new ipfulltree (dir, bounds, subtrees); + assert (superregs.AT(0).processors == NULL); + superregs.AT(0).processors = new ipfulltree (dir, bounds, subtrees); for (int c=0; c<(int)regs.size(); ++c) { - assert (regs.at(c).processor == c); + assert (regs.AT(c).processor == c); } } @@ -968,9 +1130,9 @@ namespace Carpet { vector<vector<region_t> > regss (1); SplitRegionsMaps_Automatic (cctkGH, superregss, regss); assert (superregss.size() == 1); - superregs = regss.at(0); + superregs = superregss.AT(0); assert (regss.size() == 1); - regs = regss.at(0); + regs = regss.AT(0); } @@ -993,7 +1155,7 @@ namespace Carpet { assert (superregs.size() == 1); - region_t const & reg0 = superregs.at(0); + region_t const & reg0 = superregs.AT(0); const ivect rstr0 = reg0.extent.stride(); const ivect rlb0 = reg0.extent.lower(); const ivect rub0 = reg0.extent.upper() + rstr0; @@ -1052,7 +1214,7 @@ namespace Carpet { assert (all (cub <= rub0)); assert (all (not (ipos==0) or clb==rlb0)); assert (all (not (ipos==nprocs_dir-1) or cub==rub0)); - region_t & reg = regs.at(c); + region_t & reg = regs.AT(c); ibbox & ext = reg.extent; b2vect & obnd = reg.outer_boundaries; int & proc = reg.processor; @@ -1063,21 +1225,21 @@ namespace Carpet { proc = c; pseudoregion_t preg (reg.extent, c); - subtreesx.at(i) = new ipfulltree (preg); + subtreesx.AT(i) = new ipfulltree (preg); } - boundsx.at(nprocs_dir[0]) = rub0[0] + rstr0[0]; - subtreesy.at(j) = new ipfulltree (0, boundsx, subtreesx); + boundsx.AT(nprocs_dir[0]) = rub0[0] + rstr0[0]; + subtreesy.AT(j) = new ipfulltree (0, boundsx, subtreesx); } - boundsy.at(nprocs_dir[1]) = rub0[1] + rstr0[1]; - subtreesz.at(k) = new ipfulltree (1, boundsy, subtreesy); + boundsy.AT(nprocs_dir[1]) = rub0[1] + rstr0[1]; + subtreesz.AT(k) = new ipfulltree (1, boundsy, subtreesy); } - boundsz.at(nprocs_dir[2]) = rub0[2] + rstr0[2]; + boundsz.AT(nprocs_dir[2]) = rub0[2] + rstr0[2]; - assert (superregs.at(0).processors == NULL); - superregs.at(0).processors = new ipfulltree (2, boundsz, subtreesz); + assert (superregs.AT(0).processors == NULL); + superregs.AT(0).processors = new ipfulltree (2, boundsz, subtreesz); for (int c=0; c<(int)regs.size(); ++c) { - assert (regs.at(c).processor == c); + assert (regs.AT(c).processor == c); } } @@ -1121,7 +1283,9 @@ namespace Carpet { region_t & superreg, vector<region_t> & newregs) { - if (DEBUG) cout << "SRMAR enter" << endl; + DECLARE_CCTK_PARAMETERS; + + if (recompose_verbose) cout << "SRMAR enter" << endl; // Check preconditions assert (nprocs >= 1); @@ -1130,7 +1294,7 @@ namespace Carpet { // Are we done? if (all(dims)) { - if (DEBUG) cout << "SRMAR bottom" << endl; + if (recompose_verbose) cout << "SRMAR bottom" << endl; // Check precondition assert (nprocs == 1); @@ -1145,30 +1309,57 @@ namespace Carpet { // Check postcondition assert (newregs.size() == oldsize + nprocs); - if (DEBUG) cout << "SRMAR exit" << endl; + if (recompose_verbose) cout << "SRMAR exit" << endl; return; } // Special case if (superreg.extent.empty()) { - if (DEBUG) cout << "SRMAR empty" << endl; + if (recompose_verbose) cout << "SRMAR empty" << endl; + + // Choose a direction + int mydim = -1; + for (int d=0; d<dim; ++d) { + if (not dims[d]) { + mydim = d; + break; + } + } + assert (mydim>=0 and mydim<dim); // Create a new region region_t newreg (superreg); newreg.outer_boundaries = b2vect(false); - if (DEBUG) cout << "SRMAR newreg " << newreg << endl; + if (recompose_verbose) cout << "SRMAR newreg " << newreg << endl; +#if 0 // Store for (int p=0; p<nprocs; ++p) { newreg.processor = firstproc + p; newregs.push_back (newreg); } + assert (superreg.processors == NULL); superreg.processors = new ipfulltree (); +#endif + + // Store + vector<int> bounds (nprocs+1); + vector<ipfulltree *> subtrees (nprocs); + for (int p=0; p<nprocs; ++p) { + newreg.processor = firstproc + p; + newregs.push_back (newreg); + bounds.AT(p) = newreg.extent.lower()[mydim]; + pseudoregion_t const preg (newreg.extent, newreg.processor); + subtrees.AT(p) = new ipfulltree (preg); + } + bounds.AT(nprocs) = newreg.extent.lower()[mydim]; + assert (superreg.processors == NULL); + superreg.processors = new ipfulltree (mydim, bounds, subtrees); // Check postcondition assert (newregs.size() == oldsize + nprocs); - if (DEBUG) cout << "SRMAR exit" << endl; + if (recompose_verbose) cout << "SRMAR exit" << endl; return; } @@ -1177,7 +1368,7 @@ namespace Carpet { CCTK_REAL const rfact = pow (nprocs / prod(rcost), CCTK_REAL(1)/dim); rcost *= rfact; assert (abs (prod (rcost) - nprocs) <= 1.0e-6); - if (DEBUG) cout << "SRMA shapes " << rcost << endl; + if (recompose_verbose) cout << "SRMA shapes " << rcost << endl; // Choose a direction int mydim = -1; @@ -1198,22 +1389,20 @@ namespace Carpet { } assert (mydim>=0 and mydim<dim); assert (mycost>=0); - if (DEBUG) cout << "SRMAR mydim " << mydim << endl; - if (DEBUG) cout << "SRMAR mycost " << mycost << endl; + if (recompose_verbose) cout << "SRMAR mydim " << mydim << endl; + if (recompose_verbose) cout << "SRMAR mycost " << mycost << endl; // Mark this direction as done assert (not dims[mydim]); bvect const newdims = dims.replace(mydim, true); // Choose a number of slices for this direction - int const nslices - = min (nprocs, - int (floor (mycost * pow(nprocs / totalcost, - CCTK_REAL(1) / alldims) + - CCTK_REAL(0.5)))); + CCTK_REAL const mycost1 = + mycost * pow(nprocs / totalcost, CCTK_REAL(1) / alldims); + int const nslices = min (nprocs, int (floor (mycost1 + CCTK_REAL(0.5)))); assert (nslices <= nprocs); - if (DEBUG) cout << "SRMAR " << mydim << " nprocs " << nprocs << endl; - if (DEBUG) cout << "SRMAR " << mydim << " nslices " << nslices << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " nprocs " << nprocs << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " nslices " << nslices << endl; // Split the remaining processors vector<int> mynprocs(nslices); @@ -1221,11 +1410,11 @@ namespace Carpet { int const mynprocs_left = nprocs - nslices * mynprocs_base; int sum_mynprocs = 0; for (int n=0; n<nslices; ++n) { - mynprocs.at(n) = mynprocs_base + int (n < mynprocs_left); - sum_mynprocs += mynprocs.at(n); + mynprocs.AT(n) = mynprocs_base + int (n < mynprocs_left); + sum_mynprocs += mynprocs.AT(n); } assert (sum_mynprocs == nprocs); - if (DEBUG) cout << "SRMAR " << mydim << " mynprocs " << mynprocs << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " mynprocs " << mynprocs << endl; // Split the region vector<int> mynpoints(nslices); @@ -1237,20 +1426,22 @@ namespace Carpet { int npoints_left = npoints; int nprocs_left = nprocs; for (int n=0; n<nslices; ++n) { - mynpoints.at(n) = int (floor (CCTK_REAL(1) * npoints_left * mynprocs.at(n) - / nprocs_left + CCTK_REAL(0.5))); - assert (mynprocs .at(n) > 0); - npoints_left -= mynpoints.at(n); - nprocs_left -= mynprocs.at(n); + assert (nprocs_left > 0); + CCTK_REAL const npoints1 = + CCTK_REAL(1) * npoints_left * mynprocs.AT(n) / nprocs_left; + mynpoints.AT(n) = int (floor (npoints1 + CCTK_REAL(0.5))); + assert (mynpoints.AT(n) >= 0 and mynpoints.AT(n) <= npoints_left); + npoints_left -= mynpoints.AT(n); + nprocs_left -= mynprocs.AT(n); assert (npoints_left >= 0); assert (nprocs_left >= 0); } assert (npoints_left == 0); assert (nprocs_left == 0); - if (DEBUG) cout << "SRMAR " << mydim << " mynpoints " << mynpoints << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " mynpoints " << mynpoints << endl; // Create the regions and recurse - if (DEBUG) cout << "SRMAR " << mydim << ": create bboxes" << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << ": create bboxes" << endl; ivect lo = superreg.extent.lower(); ivect up = superreg.extent.upper(); ivect const str = superreg.extent.stride(); @@ -1258,41 +1449,41 @@ namespace Carpet { vector<int> bounds (nslices+1); vector<ipfulltree *> subtrees (nslices); for (int n=0; n<nslices; ++n) { - if (DEBUG) cout << "SRMAR " << mydim << " n " << n << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " n " << n << endl; // Create a new region - up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim]; + up[mydim] = lo[mydim] + (mynpoints.AT(n) - 1) * str[mydim]; b2vect newob (superreg.outer_boundaries); if (n > 0) { newob[0][mydim] = false; } if (n < nslices-1) { - up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim]; + up[mydim] = lo[mydim] + (mynpoints.AT(n) - 1) * str[mydim]; newob[1][mydim] = false; } region_t newreg (superreg); newreg.extent = ibbox(lo, up, str); newreg.outer_boundaries = newob; - if (DEBUG) cout << "SRMAR " << mydim << " newreg " << newreg << endl; + if (recompose_verbose) cout << "SRMAR " << mydim << " newreg " << newreg << endl; // Recurse - bounds.at(n) = lo[mydim]; + bounds.AT(n) = lo[mydim]; SplitRegionsMaps_Automatic_Recursively - (newdims, p, mynprocs.at(n), newreg, newregs); - if (DEBUG) cout << "SRMAR " << mydim << " newregs " << newregs << endl; + (newdims, p, mynprocs.AT(n), newreg, newregs); + if (recompose_verbose) cout << "SRMAR " << mydim << " newregs " << newregs << endl; assert (newreg.processors != NULL); - subtrees.at(n) = newreg.processors; + subtrees.AT(n) = newreg.processors; newreg.processors = NULL; newreg.processor = p; // Next slice lo[mydim] = up[mydim] + str[mydim]; - p += mynprocs.at(n); + p += mynprocs.AT(n); } assert (up[mydim] == superreg.extent.upper()[mydim]); assert (p == firstproc + nprocs); - bounds.at(nslices) = up[mydim] + str[mydim]; + bounds.AT(nslices) = up[mydim] + str[mydim]; // Create tree assert (superreg.processors == NULL); @@ -1301,7 +1492,7 @@ namespace Carpet { // Check postcondition assert (newregs.size() == oldsize + nprocs); - if (DEBUG) cout << "SRMAR exit" << endl; + if (recompose_verbose) cout << "SRMAR exit" << endl; } @@ -1313,20 +1504,20 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; - if (DEBUG) cout << "SRMA enter" << endl; + if (recompose_verbose) cout << "SRMA enter" << endl; int const nmaps = superregss.size(); int minmap = 1000000000; for (int m=0; m<nmaps; ++m) { - for (int r=0; r<int(superregss.at(m).size()); ++r) { - minmap = min (minmap, superregss.at(m).at(r).map); + for (int r=0; r<int(superregss.AT(m).size()); ++r) { + minmap = min (minmap, superregss.AT(m).AT(r).map); } } int nregs = 0; for (int m=0; m<nmaps; ++m) { - nregs += superregss.at(m).size(); + nregs += superregss.AT(m).size(); } - if (DEBUG) cout << "SRMA nregs " << nregs << endl; + if (recompose_verbose) cout << "SRMA nregs " << nregs << endl; // Something to do? if (nregs == 0) { @@ -1337,7 +1528,7 @@ namespace Carpet { vector<region_t> superregs; { for (int m=0; m<nmaps; ++m) { - combine_regions (superregss.at(m), superregs); + combine_regions (superregss.AT(m), superregs); } nregs = superregs.size(); @@ -1351,14 +1542,13 @@ namespace Carpet { reg.extent = ibbox (ivect (0), ivect (-1), ivect (1)); reg.outer_boundaries = b2vect (bvect (true), bvect (true)); reg.map = 0; - reg.processor = -1; superregs.push_back (reg); nregs = superregs.size(); } } int const real_nprocs = CCTK_nProcs (cctkGH); - if (DEBUG) cout << "SRMA real_nprocs " << real_nprocs << endl; + if (recompose_verbose) cout << "SRMA real_nprocs " << real_nprocs << endl; // Deactivate some processors if there are too many int nprocs; @@ -1367,109 +1557,155 @@ namespace Carpet { } else { CCTK_REAL mycost = 0; for (int r=0; r<nregs; ++r) { - mycost += prod (cost (superregs.at(r))); + mycost += prod (cost (superregs.AT(r))); } int const goodnprocs = int (floor (mycost / min_points_per_proc)); nprocs = max (1, min (real_nprocs, goodnprocs)); } - if (DEBUG) cout << "SRMA nprocs " << nprocs << endl; + if (recompose_verbose) cout << "SRMA nprocs " << nprocs << endl; // ncomps: number of components per processor int const ncomps = (nregs + nprocs - 1) / nprocs; - if (DEBUG) cout << "SRMA ncomps " << ncomps << endl; + if (recompose_verbose) cout << "SRMA ncomps " << ncomps << endl; assert (ncomps > 0); int const newnregs = nprocs * ncomps; - if (DEBUG) cout << "SRMA newnregs " << newnregs << endl; + if (recompose_verbose) cout << "SRMA newnregs " << newnregs << endl; - if (DEBUG) cout << "SRMA: distributing processors to regions" << endl; + if (recompose_verbose) cout << "SRMA: distributing processors to regions" << endl; vector<CCTK_REAL> mycosts(nregs); for (int r=0; r<nregs; ++r) { - mycosts.at(r) = prod (cost (superregs.at(r))); + mycosts.AT(r) = prod (cost (superregs.AT(r))); } int nregs_left = newnregs; vector<int> mynprocs(nregs); for (int r=0; r<nregs; ++r) { - mynprocs.at(r) = 1; + mynprocs.AT(r) = 1; -- nregs_left; } #warning "TODO: split regions if necessary" while (nregs_left > 0) { - if (DEBUG) cout << "SRMA nregs_left " << nregs_left << endl; + if (recompose_verbose) cout << "SRMA nregs_left " << nregs_left << endl; int maxr = -1; CCTK_REAL maxratio = -1; for (int r=0; r<nregs; ++r) { - CCTK_REAL const ratio = mycosts.at(r) / mynprocs.at(r); + CCTK_REAL const ratio = mycosts.AT(r) / mynprocs.AT(r); if (ratio > maxratio) { maxr=r; maxratio=ratio; } } assert (maxr>=0 and maxr<nregs); - ++ mynprocs.at(maxr); - if (DEBUG) cout << "SRMA maxr " << maxr << endl; - if (DEBUG) cout << "SRMA mynprocs[maxr] " << mynprocs.at(maxr) << endl; + ++ mynprocs.AT(maxr); + if (recompose_verbose) cout << "SRMA maxr " << maxr << endl; + if (recompose_verbose) cout << "SRMA mynprocs[maxr] " << mynprocs.AT(maxr) << endl; -- nregs_left; } assert (nregs_left == 0); int sum_nprocs = 0; for (int r=0; r<nregs; ++r) { - sum_nprocs += mynprocs.at(r); + sum_nprocs += mynprocs.AT(r); } assert (sum_nprocs == newnregs); - if (DEBUG) cout << "SRMA mynprocs " << mynprocs << endl; + if (recompose_verbose) cout << "SRMA mynprocs " << mynprocs << endl; - if (DEBUG) cout << "SRMA: splitting work units" << endl; + if (recompose_verbose) cout << "SRMA: splitting work units" << endl; #warning "TODO: rename newregs to regs" #warning "TODO: rename nregs to nsuperregs" #warning "TODO: rename newnregs to nregs" vector<region_t> newregs; newregs.reserve (newnregs); - for (int r=0, p=0; r<nregs; p+=mynprocs.at(r), ++r) { - if (DEBUG) cout << "SRMA superreg[" << r << "] " << superregs.at(r) << endl; + for (int r=0, p=0; r<nregs; p+=mynprocs.AT(r), ++r) { + if (recompose_verbose) cout << "SRMA superreg[" << r << "] " << superregs.AT(r) << endl; bvect const dims = false; SplitRegionsMaps_Automatic_Recursively - (dims, p, mynprocs.at(r), superregs.at(r), newregs); + (dims, p, mynprocs.AT(r), superregs.AT(r), newregs); } // for r - if (DEBUG) cout << "SRMA newregs " << newregs << endl; + if (recompose_verbose) cout << "SRMA newregs " << newregs << endl; assert (int(newregs.size()) == newnregs); // Count components per map vector<int> myncomps(nmaps, 0); +#if 0 vector<int> empty_comps(nmaps, 0); +#endif for (int r=0; r<newnregs; ++r) { - int const m = newregs.at(r).map - minmap; + int const m = newregs.AT(r).map - minmap; assert (m>=0 and m<nmaps); - if (not newregs.at(r).extent.empty()) { +#if 0 + if (not newregs.AT(r).extent.empty()) { // Ignore empty regions, which may come from empty grid arrays - ++ myncomps.at(m); + ++ myncomps.AT(m); } else { - ++ empty_comps.at(m); + ++ empty_comps.AT(m); } +#endif + ++ myncomps.AT(m); } vector<int> mynregs(nmaps, 0); +#if 0 + vector<int> empty_regs(nmaps, 0); +#endif for (int r=0; r<nregs; ++r) { - int const m = superregs.at(r).map - minmap; + int const m = superregs.AT(r).map - minmap; assert (m>=0 and m<nmaps); - ++ mynregs.at(m); + ++ mynregs.AT(m); +#if 0 + if (not superregs.AT(r).extent.empty()) { + // Ignore empty regions, which may come from empty grid arrays + ++ mynregs.AT(m); + } else { + ++ empty_regs.AT(m); + } +#endif } // Convert virtual to real processors for (int r=0; r<newnregs; ++r) { - newregs.at(r).processor /= ncomps; - assert (newregs.at(r).processor >= 0 and - newregs.at(r).processor < nprocs); + newregs.AT(r).processor /= ncomps; + assert (newregs.AT(r).processor >= 0 and + newregs.AT(r).processor < nprocs); } { vector<int> tmpncomps(nmaps, 0); for (int r=0; r<nregs; ++r) { - int const m = superregs.at(r).map - minmap; - ipfulltree * const regf = superregs.at(r).processors; + int const m = superregs.AT(r).map - minmap; + ipfulltree * const regf = superregs.AT(r).processors; assert (regf != NULL); for (ipfulltree::iterator fti (* regf); not fti.done(); ++ fti) { pseudoregion_t & preg = (* fti).payload(); - preg.component = tmpncomps.at(m)++; + preg.component = tmpncomps.AT(m)++; } } for (int m=0; m<nmaps; ++m) { - assert (tmpncomps.at(m) == myncomps.at(m)); +#warning "TODO" + if (not (tmpncomps.AT(m) == myncomps.AT(m))) { + cout << "Recompose.cc" << endl + << "superregss=" << superregss << endl + << "regss=" << regss << endl + << "nregs=" << nregs << endl + << "newregs=" << newregs << endl + << "newnregs=" << newnregs << endl + << "m=" << m << endl + << "tmpncomps.AT(m)=" << tmpncomps.AT(m) << endl + << "myncomps.AT(m)=" << myncomps.AT(m) << endl; + cout << "newregs:" << endl; + for (int r=0; r<newnregs; ++r) { + int const m = newregs.AT(r).map - minmap; + assert (m>=0 and m<nmaps); + cout << "r=" << r << endl + << "newregs=" << newregs.at(r) << endl + << "newregs.extent.size()=" << newregs.at(r).extent.size() << endl + << "newregs.extent.empty()=" << newregs.at(r).extent.empty() << endl; + } + cout << "*regf:" << endl; + for (int r=0; r<nregs; ++r) { + int const m = superregs.AT(r).map - minmap; + ipfulltree * const regf = superregs.AT(r).processors; + cout << "r=" << r << endl + << "m=" << m << endl; + cout << "*regf=" << *regf << endl; + } + cout.flush(); + } + assert (tmpncomps.AT(m) == myncomps.AT(m)); } } @@ -1477,34 +1713,37 @@ namespace Carpet { // Allocate regions assert ((int)regss.size() == nmaps); for (int m=0; m<nmaps; ++m) { - assert (regss.at(m).empty()); - regss.at(m).reserve (myncomps.at(m)); - superregss.at(m).clear(); - superregss.at(m).reserve (mynregs.at(m)); + assert (regss.AT(m).empty()); + regss.AT(m).reserve (myncomps.AT(m)); + superregss.AT(m).clear(); + superregss.AT(m).reserve (mynregs.AT(m)); } // Assign regions for (int r=0; r<newnregs; ++r) { - int const m = newregs.at(r).map - minmap; + int const m = newregs.AT(r).map - minmap; assert (m>=0 and m<nmaps); - regss.at(m).push_back (newregs.at(r)); + regss.AT(m).push_back (newregs.AT(r)); } for (int r=0; r<nregs; ++r) { - int const m = superregs.at(r).map - minmap; + int const m = superregs.AT(r).map - minmap; assert (m>=0 and m<nmaps); - superregss.at(m).push_back (superregs.at(r)); + superregss.AT(m).push_back (superregs.AT(r)); } // Output regions - if (DEBUG) { + if (recompose_verbose) { cout << "SRMA superregss " << superregss << endl; cout << "SRMA regss " << regss << endl; } // Check sizes for (int m=0; m<nmaps; ++m) { - assert (int(regss.at(m).size()) == myncomps.at(m) + empty_comps.at(m)); - assert (int(superregss.at(m).size()) == mynregs.at(m)); +#if 0 + assert (int(regss.AT(m).size()) == myncomps.AT(m) + empty_comps.AT(m)); +#endif + assert (int(regss.AT(m).size()) == myncomps.AT(m)); + assert (int(superregss.AT(m).size()) == mynregs.AT(m)); } - if (DEBUG) cout << "SRMA exit" << endl; + if (recompose_verbose) cout << "SRMA exit" << endl; } @@ -1547,30 +1786,30 @@ namespace Carpet { } } vector<ibbox> bases(mglevels); - bases.at(0) = base; + bases.AT(0) = base; for (int ml=1; ml<mglevels; ++ml) { // next finer base - ivect const fbaselo = bases.at(ml-1).lower(); - ivect const fbasehi = bases.at(ml-1).upper(); - ivect const fbasestr = bases.at(ml-1).stride(); + ivect const fbaselo = bases.AT(ml-1).lower(); + ivect const fbasehi = bases.AT(ml-1).upper(); + ivect const fbasestr = bases.AT(ml-1).stride(); // this base ivect const basestr = fbasestr * mgfact; ivect const baselo = fbaselo + (xpose(offset)[0] - ivect(mgfact) * xpose(offset)[0]) * fbasestr; ivect const basehi = fbasehi + (xpose(offset)[1] - ivect(mgfact) * xpose(offset)[1]) * fbasestr; ivect const baselo1 = baselo; ivect const basehi1 = baselo1 + (basehi - baselo1) / basestr * basestr; - bases.at(ml) = ibbox(baselo1, basehi1, basestr); + bases.AT(ml) = ibbox(baselo1, basehi1, basestr); // next finer grid - ivect const flo = regs.at(ml-1).extent.lower(); - ivect const fhi = regs.at(ml-1).extent.upper(); - ivect const fstr = regs.at(ml-1).extent.stride(); + ivect const flo = regs.AT(ml-1).extent.lower(); + ivect const fhi = regs.AT(ml-1).extent.upper(); + ivect const fstr = regs.AT(ml-1).extent.stride(); // this grid ivect const str = fstr * mgfact; ivect const lo = flo + either (reg.outer_boundaries[0], (xpose(offset)[0] - ivect(mgfact) * xpose(offset)[0]) * fstr, ivect(0)); ivect const hi = fhi + either (reg.outer_boundaries[1], - (xpose(offset)[1] - ivect(mgfact) * xpose(offset)[1]) * fstr, ivect(0)); ivect const lo1 = baselo1 + (lo - baselo1 + str - 1) / str * str; ivect const hi1 = lo1 + (hi - lo1) / str * str; - regs.at(ml).extent = ibbox(lo1, hi1, str); + regs.AT(ml).extent = ibbox(lo1, hi1, str); } } } @@ -1583,25 +1822,25 @@ namespace Carpet { { regsss.resize (mglevels); for (int ml=0; ml<mglevels; ++ml) { - regsss.at(ml).resize (regss.size()); + regsss.AT(ml).resize (regss.size()); for (int rl=0; rl<(int)regss.size(); ++rl) { - regsss.at(ml).at(rl).resize (regss.at(rl).size()); + regsss.AT(ml).AT(rl).resize (regss.AT(rl).size()); } } for (int rl=0; rl<(int)regss.size(); ++rl) { ibbox base; - for (int c=0; c<(int)regss.at(rl).size(); ++c) { - base = base.expanded_containing(regss.at(rl).at(c).extent); + for (int c=0; c<(int)regss.AT(rl).size(); ++c) { + base = base.expanded_containing(regss.AT(rl).AT(c).extent); } - for (int c=0; c<(int)regss.at(rl).size(); ++c) { + for (int c=0; c<(int)regss.AT(rl).size(); ++c) { vector<region_t> mg_regs; - MakeMultigridBoxes (cctkGH, m, base, regss.at(rl).at(c), mg_regs); + MakeMultigridBoxes (cctkGH, m, base, regss.AT(rl).AT(c), mg_regs); assert ((int)mg_regs.size() == mglevels); for (int ml=0; ml<mglevels; ++ml) { - regsss.at(ml).at(rl).at(c) = mg_regs.at(ml); + regsss.AT(ml).AT(rl).AT(c) = mg_regs.AT(ml); } } // for c @@ -1615,8 +1854,46 @@ namespace Carpet { vector<vector<vector<vector<region_t> > > > & regssss) { for (int m = 0; m < maps; ++m) { - MakeMultigridBoxes (cctkGH, m, regsss.at(m), regssss.at(m)); + MakeMultigridBoxes (cctkGH, m, regsss.AT(m), regssss.AT(m)); } // for m } + + + static + void + ClassifyPoints (cGH const * const cctkGH, int const rl) + { + // negative: needs to be set explicitly (e.g. boundary) + // zero: unused (e.g. ghost) + // positive: needs to be evolved + // -1: boundary point (needs to be set explicitly) + // 0: unused (e.g. ghost point, or restriction target) + // n=1..N: evolved, used for integrator substeps i<=n + // (i=N..1, counting backwards; see MoL documentation) + // i.e.: n=1: used between time steps (i.e., should be visualised) + // n>1: used only while time stepping (e.g. buffer zones) + + BEGIN_META_MODE(cctkGH) { + BEGIN_MGLEVEL_LOOP (cctkGH) { + ENTER_LEVEL_MODE (cctkGH, rl) { + BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; +#pragma omp parallel + LC_LOOP3 (CarpetClassifyPoints, + i,j,k, + 0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + point_class[ind] = 1; + } LC_ENDLOOP3(CarpetClassifyPoints); + } END_LOCAL_COMPONENT_LOOP; + } END_LOCAL_MAP_LOOP; + } LEAVE_LEVEL_MODE; + } END_MGLEVEL_LOOP; + } END_META_MODE; + } + } // namespace Carpet |