diff options
| author | Erik Schnetter <schnetter@cct.lsu.edu> | 2008-06-20 16:29:06 -0500 |
|---|---|---|
| committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2008-06-20 16:29:06 -0500 |
| commit | 4b16584382e52728dc658deed7b38ba78f41e865 (patch) | |
| tree | 80ec1ec1a21dc870b50dfd0eb818160f6bc7ec81 | |
| parent | bdfdbc1b611cd9339fb693b92d12ae3becc68d37 (diff) | |
Introduce a tree data structure to speed up domain decomposition
Introduce a tree data structure "fulltree", which decomposes a single,
rectangular region into a tree of non-overlapping, rectangular sub-regions.
Move the processor decomposition from the regridding thorns into Carpet.
Create such trees during processor decomposition.
Store these trees with the grid hierarchy.
30 files changed, 1256 insertions, 333 deletions
diff --git a/Carpet/Carpet/interface.ccl b/Carpet/Carpet/interface.ccl index 734020bae..cee9c1075 100644 --- a/Carpet/Carpet/interface.ccl +++ b/Carpet/Carpet/interface.ccl @@ -12,10 +12,12 @@ uses include header: dist.hh uses include header: bbox.hh uses include header: bboxset.hh +uses include header: fulltree.hh uses include header: region.hh uses include header: vect.hh uses include header: data.hh + uses include header: dh.hh uses include header: gf.hh uses include header: ggf.hh @@ -205,18 +207,22 @@ PROVIDES FUNCTION GetRefinementLevels \ # The true prototype of the routine below: # int Carpet_Regrid (const cGH * cctkGH, +# gh::rregs * superregss, # gh::mregs * regsss, # int force); CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_POINTER IN superregss, \ CCTK_POINTER IN regsss, \ CCTK_INT IN force) USES FUNCTION Carpet_Regrid # The true prototype of the routine below: # int Carpet_Regrid (const cGH * cctkGH, +# vector<gh::rregs> * superregsss, # vector<gh::mregs> * regssss, # int force); CCTK_INT FUNCTION Carpet_RegridMaps (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_POINTER IN superregsss, \ CCTK_POINTER IN regssss, \ CCTK_INT IN force) USES FUNCTION Carpet_RegridMaps diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index 9bd7d4066..204542f02 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -60,12 +60,14 @@ namespace Carpet { static void - SplitRegionsMaps_Automatic_Recursively (bvect const & dims, - int const nprocs, - region_t const & reg, + SplitRegionsMaps_Automatic_Recursively (bvect const & dims, + int const firstproc, + int const nprocs, + region_t & superreg, vector<region_t> & newregs); static void SplitRegions_AsSpecified (cGH const * cctkGH, + vector<region_t> & superregs, vector<region_t> & regs); @@ -124,6 +126,8 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; + Checkpoint ("Regridding level %d...", reflevel); + assert (is_level_mode()); bool const have_regrid = CCTK_IsFunctionAliased ("Carpet_Regrid"); @@ -155,37 +159,40 @@ namespace Carpet { BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - gh::mregs regsss = vhh.at(map)->regions; + gh::rregs superregss = vhh.at(map)->superregions; + gh::mregs regsss; // Check whether to recompose CCTK_INT const do_recompose = - Carpet_Regrid (cctkGH, & regsss, force_recompose); + Carpet_Regrid (cctkGH, & superregss, & regsss, force_recompose); assert (do_recompose >= 0); did_change = did_change or do_recompose; if (do_recompose) { - RegridMap (cctkGH, map, regsss); + RegridMap (cctkGH, map, superregss, regsss); } } END_MAP_LOOP; } else { + vector<gh::rregs> superregsss (maps); vector<gh::mregs> regssss (maps); BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - regssss.at(map) = vhh.at(map)->regions; + superregsss.at(map) = vhh.at(map)->superregions; } END_MAP_LOOP; // Check whether to recompose CCTK_INT const do_recompose = - Carpet_RegridMaps (cctkGH, & regssss, force_recompose); + Carpet_RegridMaps (cctkGH, & superregsss, & regssss, force_recompose); assert (do_recompose >= 0); 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, regsss); + RegridMap (cctkGH, map, superregss, regsss); } END_MAP_LOOP; } @@ -209,6 +216,7 @@ namespace Carpet { void RegridMap (cGH const * const cctkGH, int const m, + gh::rregs const & superregss, gh::mregs const & regsss) { DECLARE_CCTK_PARAMETERS; @@ -221,13 +229,15 @@ namespace Carpet { // not change // Regrid - vhh.at(m)->regrid (regsss); + vhh.at(m)->regrid (superregss, regsss); // Write grid structure to file OutputGridStructure (cctkGH, m, regsss); OutputGridCoordinates (cctkGH, m, regsss); - if (verbose) OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m)); + if (verbose or veryverbose) { + OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m)); + } Waypoint ("Done regridding map %d.", m); } @@ -717,18 +727,20 @@ namespace Carpet { // SplitRegions_AlongZ for grid arrays) void SplitRegions (cGH const * const cctkGH, + vector<region_t> & superregs, vector<region_t> & regs) { DECLARE_CCTK_PARAMETERS; + assert (regs.empty()); if (CCTK_EQUALS (processor_topology, "along-z")) { - SplitRegions_AlongZ (cctkGH, regs); + SplitRegions_AlongZ (cctkGH, superregs, regs); } else if (CCTK_EQUALS (processor_topology, "along-dir")) { - SplitRegions_AlongDir (cctkGH, regs, split_direction); + SplitRegions_AlongDir (cctkGH, superregs, regs, split_direction); } else if (CCTK_EQUALS (processor_topology, "automatic")) { - SplitRegions_Automatic (cctkGH, regs); + SplitRegions_Automatic (cctkGH, superregs, regs); } else if (CCTK_EQUALS (processor_topology, "manual")) { - SplitRegions_AsSpecified (cctkGH, regs); + SplitRegions_AsSpecified (cctkGH, superregs, regs); } else { assert (0); } @@ -835,29 +847,38 @@ namespace Carpet { void SplitRegions_AlongZ (cGH const * const cctkGH, + vector<region_t> & superregs, vector<region_t> & regs) { - SplitRegions_AlongDir (cctkGH, regs, dim - 1); + assert (regs.empty()); + SplitRegions_AlongDir (cctkGH, superregs, regs, dim - 1); } void SplitRegions_AlongDir (cGH const * const cctkGH, + vector<region_t> & superregs, vector<region_t> & regs, int const dir) { + assert (regs.empty()); + // Something to do? - if (regs.size() == 0) { + if (superregs.size() == 0) { return; } const int nprocs = CCTK_nProcs (cctkGH); - assert (regs.size() == 1); + assert (superregs.size() == 1); + 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); return; } @@ -870,6 +891,8 @@ namespace Carpet { const b2vect obnd0 = reg0.outer_boundaries; regs.resize (nprocs); + vector<int> bounds (nprocs+1); + vector<ipfulltree *> subtrees (nprocs); for (int c=0; c<nprocs; ++c) { ivect cstr = rstr0; ivect clb = rlb0; @@ -892,7 +915,14 @@ namespace Carpet { obnd[0][dir] &= clb[dir] == rlb0[dir]; 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(nprocs) = rub0[dir] + rstr0[dir]; + + 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); @@ -903,10 +933,15 @@ namespace Carpet { void SplitRegions_Automatic (cGH const * const cctkGH, + vector<region_t> & superregs, vector<region_t> & regs) { - vector<vector<region_t> > regss (1, regs); - SplitRegionsMaps_Automatic (cctkGH, regss); + assert (regs.empty()); + vector<vector<region_t> > superregss (1, superregs); + vector<vector<region_t> > regss (1); + SplitRegionsMaps_Automatic (cctkGH, superregss, regss); + assert (superregss.size() == 1); + superregs = regss.at(0); assert (regss.size() == 1); regs = regss.at(0); } @@ -915,20 +950,23 @@ namespace Carpet { static void SplitRegions_AsSpecified (cGH const * const cctkGH, + vector<region_t> & superregs, vector<region_t> & regs) { DECLARE_CCTK_PARAMETERS; + assert (regs.empty()); + // Something to do? - if (regs.size() == 0) { + if (superregs.size() == 0) { return; } const int nprocs = CCTK_nProcs (cctkGH); - assert (regs.size() == 1); + assert (superregs.size() == 1); - region_t const & reg0 = regs.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; @@ -956,9 +994,17 @@ namespace Carpet { const ivect rem = glonp % nprocs_dir; const ivect step = locnp * cstr; assert (dim==3); + + vector<int> boundsz(nprocs_dir[2]); + vector<ipfulltree *> subtreesz(nprocs_dir[2]+1); for (int k=0; k<nprocs_dir[2]; ++k) { + vector<int> boundsy(nprocs_dir[2]); + vector<ipfulltree *> subtreesy(nprocs_dir[2]+1); for (int j=0; j<nprocs_dir[1]; ++j) { + vector<int> boundsx(nprocs_dir[2]); + vector<ipfulltree *> subtreesx(nprocs_dir[2]+1); for (int i=0; i<nprocs_dir[0]; ++i) { + const int c = i + nprocs_dir[0] * (j + nprocs_dir[1] * k); const ivect ipos (i, j, k); ivect clb = rlb0 + step * ipos; @@ -988,9 +1034,20 @@ namespace Carpet { obnd[0] &= clb == rlb0; obnd[1] &= cub == rub0; proc = c; + + pseudoregion_t preg (reg.extent, c); + subtreesx.at(i) = new ipfulltree (preg); } + 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); } + boundsz.at(nprocs_dir[2]) = rub0[2] + rstr0[2]; + + 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); @@ -1003,21 +1060,26 @@ namespace Carpet { // SplitRegions_AlongZ for grid arrays) void SplitRegionsMaps (cGH const * const cctkGH, + vector<vector<region_t> > & superregss, vector<vector<region_t> > & regss) { DECLARE_CCTK_PARAMETERS; + assert (regss.size() == superregss.size()); + for (size_t m=0; m<regss.size(); ++m) { + assert (regss.AT(m).empty()); + } if (CCTK_EQUALS (processor_topology, "along-z")) { assert (0); -// SplitRegionsMaps_AlongZ (cctkGH, regss); +// SplitRegionsMaps_AlongZ (cctkGH, superregss, regss); } else if (CCTK_EQUALS (processor_topology, "along-dir")) { assert (0); -// SplitRegionsMaps_AlongDir (cctkGH, regss, split_direction); +// SplitRegionsMaps_AlongDir (cctkGH, superregss, regss, split_direction); } else if (CCTK_EQUALS (processor_topology, "automatic")) { - SplitRegionsMaps_Automatic (cctkGH, regss); + SplitRegionsMaps_Automatic (cctkGH, superregss, regss); } else if (CCTK_EQUALS (processor_topology, "manual")) { assert (0); -// SplitRegionsMaps_AsSpecified (cctkGH, regss); +// SplitRegionsMaps_AsSpecified (cctkGH, superregss, regss); } else { assert (0); } @@ -1026,9 +1088,10 @@ namespace Carpet { static void - SplitRegionsMaps_Automatic_Recursively (bvect const & dims, - int const nprocs, - region_t const & reg, + SplitRegionsMaps_Automatic_Recursively (bvect const & dims, + int const firstproc, + int const nprocs, + region_t & superreg, vector<region_t> & newregs) { if (DEBUG) cout << "SRMAR enter" << endl; @@ -1046,7 +1109,11 @@ namespace Carpet { assert (nprocs == 1); // Return argument - newregs.push_back (reg); + region_t newreg = superreg; + newreg.processor = firstproc; + newregs.push_back (newreg); + pseudoregion_t const preg (newreg.extent, newreg.processor); + superreg.processors = new ipfulltree (preg); // Check postcondition assert (newregs.size() == oldsize + nprocs); @@ -1056,19 +1123,21 @@ namespace Carpet { } // Special case - if (reg.extent.empty()) { + if (superreg.extent.empty()) { if (DEBUG) cout << "SRMAR empty" << endl; // Create a new region - region_t newreg (reg); - newreg.outer_boundaries = b2vect(false); + region_t newreg (superreg); + newreg.outer_boundaries = b2vect(false); if (DEBUG) cout << "SRMAR newreg " << newreg << endl; // Store for (int p=0; p<nprocs; ++p) { - newreg.processor = reg.processor + p; + newreg.processor = firstproc + p; newregs.push_back (newreg); } + // Do not set any superreg.processors information + // TODO: Maybe create an empty tree for this? // Check postcondition assert (newregs.size() == oldsize + nprocs); @@ -1078,10 +1147,10 @@ namespace Carpet { } // Calculate cost factors - rvect rcost = cost (reg); + rvect rcost = cost (superreg); CCTK_REAL const rfact = pow (nprocs / prod(rcost), CCTK_REAL(1)/dim); rcost *= rfact; - assert (abs (prod (rcost) - nprocs) < 1.0e-6); + assert (abs (prod (rcost) - nprocs) <= 1.0e-6); if (DEBUG) cout << "SRMA shapes " << rcost << endl; // Choose a direction @@ -1134,7 +1203,8 @@ namespace Carpet { // Split the region vector<int> mynpoints(nslices); - int const npoints = (reg.extent.shape() / reg.extent.stride())[mydim]; + int const npoints = + (superreg.extent.shape() / superreg.extent.stride())[mydim]; // Keep track of how many points and processors we have left to // distribute @@ -1155,16 +1225,18 @@ namespace Carpet { // Create the regions and recurse if (DEBUG) cout << "SRMAR " << mydim << ": create bboxes" << endl; - ivect lo = reg.extent.lower(); - ivect up = reg.extent.upper(); - ivect const str = reg.extent.stride(); - int p = reg.processor; + ivect lo = superreg.extent.lower(); + ivect up = superreg.extent.upper(); + ivect const str = superreg.extent.stride(); + int p = firstproc; + vector<int> bounds (nslices+1); + vector<ipfulltree *> subtrees (nslices); for (int n=0; n<nslices; ++n) { if (DEBUG) cout << "SRMAR " << mydim << " n " << n << endl; // Create a new region up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim]; - b2vect newob (reg.outer_boundaries); + b2vect newob (superreg.outer_boundaries); if (n > 0) { newob[0][mydim] = false; } @@ -1172,23 +1244,33 @@ namespace Carpet { up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim]; newob[1][mydim] = false; } - region_t newreg (reg); + region_t newreg (superreg); newreg.extent = ibbox(lo, up, str); - newreg.outer_boundaries = newob; - newreg.processor = p; + newreg.outer_boundaries = newob; if (DEBUG) cout << "SRMAR " << mydim << " newreg " << newreg << endl; // Recurse + bounds.at(n) = lo[mydim]; SplitRegionsMaps_Automatic_Recursively - (newdims, mynprocs.at(n), newreg, newregs); + (newdims, p, mynprocs.at(n), newreg, newregs); if (DEBUG) cout << "SRMAR " << mydim << " newregs " << newregs << endl; + assert (newreg.processors != NULL); + subtrees.at(n) = newreg.processors; + newreg.processors = NULL; + newreg.processor = p; + // Next slice lo[mydim] = up[mydim] + str[mydim]; p += mynprocs.at(n); } - assert (up[mydim] == reg.extent.upper()[mydim]); - assert (p == reg.processor + nprocs); + assert (up[mydim] == superreg.extent.upper()[mydim]); + assert (p == firstproc + nprocs); + bounds.at(nslices) = up[mydim] + str[mydim]; + + // Create tree + assert (superreg.processors == NULL); + superreg.processors = new ipfulltree (mydim, bounds, subtrees); // Check postcondition assert (newregs.size() == oldsize + nprocs); @@ -1200,14 +1282,15 @@ namespace Carpet { void SplitRegionsMaps_Automatic (cGH const * const cctkGH, + vector<vector<region_t> > & superregss, vector<vector<region_t> > & regss) { if (DEBUG) cout << "SRMA enter" << endl; - int const nmaps = regss.size(); + int const nmaps = superregss.size(); int nregs = 0; for (int m=0; m<nmaps; ++m) { - nregs += regss.at(m).size(); + nregs += superregss.at(m).size(); } if (DEBUG) cout << "SRMA nregs " << nregs << endl; @@ -1217,84 +1300,12 @@ namespace Carpet { } // Collect slices - vector<region_t> regs; + vector<region_t> superregs; { -#if 0 - regs.resize (nregs); - int r=0; for (int m=0; m<nmaps; ++m) { - for (int c=0; c<int(regss.at(m).size()); ++c, ++r) { - regs.at(r) = regss.at(m).at(c); - regs.at(r).map = m; - regs.at(r).processor = -1; - } + combine_regions (superregss.at(m), superregs); } - assert (r == nregs); -#endif - for (int m=0; m<nmaps; ++m) { - ibset comps; - ibset cobnds[2][dim]; - for (size_t c=0; c<regss.at(m).size(); ++c) { - region_t const & reg = regss.at(m).at(c); - comps += reg.extent; - for (int f = 0; f < 2; ++ f) { - for (int d = 0; d < dim; ++ d) { - if (reg.outer_boundaries[f][d]) { - ibbox bnd = reg.extent; - ivect lo = bnd.lower(); - ivect up = bnd.upper(); - if (f==0) { - up[d] = lo[d]; - } else { - lo[d] = up[d]; - } - bnd = ibbox (lo, up, bnd.stride()); - cobnds[f][d] += bnd; - } - } - } - } - comps.normalize(); - for (int f = 0; f < 2; ++ f) { - for (int d = 0; d < dim; ++ d) { - cobnds[f][d].normalize(); - } - } - size_t const needsize = regs.size() + comps.setsize(); - if (regs.capacity() < needsize) { - regs.reserve (1000 + 2 * needsize); - } - for (ibset::const_iterator ci = comps.begin(); ci != comps.end(); ++ ci) - { - ibbox const & c = * ci; - b2vect obnds; - for (int f = 0; f < 2; ++ f) { - for (int d = 0; d < dim; ++ d) { - obnds[f][d] = cobnds[f][d].intersects (c); - if (obnds[f][d]) { - ivect lo = c.lower(); - ivect up = c.upper(); - if (f) lo[d]=up[d]; else up[d]=lo[d]; - ibbox const cbnds (lo, up, c.stride()); - if (not ((cobnds[f][d] & ibset(c)) == ibset(cbnds))) { - cout << "cobnds[f][d] = " << cobnds[f][d] << endl - << "ibset(c) = " << ibset(c) << endl - << "(cobnds[f][d] & ibset(c)) = " << (cobnds[f][d] & ibset(c)) << endl - << "ibset(cbnds) = " << ibset(cbnds) << endl; - } - assert ((cobnds[f][d] & ibset(c)) == ibset(cbnds)); - } - } - } - region_t reg; - reg.extent = c; - reg.outer_boundaries = obnds; - reg.map = m; - reg.processor = -1; - regs.push_back (reg); - } - } // for m - nregs = regs.size(); + nregs = superregs.size(); // If the last region was removed, add a new empty region again. // A set of regions (corresponding to a refinement level or a @@ -1307,8 +1318,8 @@ namespace Carpet { reg.outer_boundaries = b2vect (bvect (true), bvect (true)); reg.map = 0; reg.processor = -1; - regs.push_back (reg); - nregs = regs.size(); + superregs.push_back (reg); + nregs = superregs.size(); } } @@ -1326,7 +1337,7 @@ namespace Carpet { if (DEBUG) cout << "SRMA: distributing processors to regions" << endl; vector<CCTK_REAL> mycosts(nregs); for (int r=0; r<nregs; ++r) { - mycosts.at(r) = prod (cost (regs.at(r))); + mycosts.at(r) = prod (cost (superregs.at(r))); } int nregs_left = newnregs; vector<int> mynprocs(nregs); @@ -1361,34 +1372,61 @@ namespace Carpet { vector<region_t> newregs; newregs.reserve (newnregs); for (int r=0, p=0; r<nregs; p+=mynprocs.at(r), ++r) { - regs.at(r).processor = p; - if (DEBUG) cout << "SRMA reg[" << r << "] " << regs.at(r) << endl; + if (DEBUG) cout << "SRMA superreg[" << r << "] " << superregs.at(r) << endl; bvect const dims = false; SplitRegionsMaps_Automatic_Recursively - (dims, mynprocs.at(r), regs.at(r), newregs); + (dims, p, mynprocs.at(r), superregs.at(r), newregs); } // for r if (DEBUG) cout << "SRMA newregs " << newregs << endl; assert (int(newregs.size()) == newnregs); + // Count components per map + vector<int> myncomps(nmaps, 0); + for (int r=0; r<newnregs; ++r) { + int const m = newregs.at(r).map; + assert (m>=0 and m<nmaps); + ++ myncomps.at(m); + } + vector<int> mynregs(nmaps, 0); + for (int r=0; r<nregs; ++r) { + int const m = superregs.at(r).map; + assert (m>=0 and m<nmaps); + ++ mynregs.at(m); + } + // 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); } + { + vector<int> tmpncomps(nmaps, 0); + for (int r=0; r<nregs; ++r) { + int const m = superregs.at(r).map; + ipfulltree * const regf = superregs.at(r).processors; + assert (regf != NULL); + for (ipfulltree::iterator + fti = regf->begin(); fti != regf->end(); ++ fti) + { + pseudoregion_t & preg = (* fti).payload(); + preg.component = tmpncomps.at(m)++; + // preg.processor /= ncomps; + } + } + for (int m=0; m<nmaps; ++m) { + assert (tmpncomps.at(m) == myncomps.at(m)); + } + } // Distribute regions - // Count components per map - vector<int> myncomps(nmaps); - for (int r=0; r<newnregs; ++r) { - int const m = newregs.at(r).map; - assert (m>=0 and m<nmaps); - ++ myncomps.at(m); - } // Allocate regions + assert ((int)regss.size() == nmaps); for (int m=0; m<nmaps; ++m) { - regss.at(m).resize (0); + 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) { @@ -1396,9 +1434,15 @@ namespace Carpet { assert (m>=0 and m<nmaps); regss.at(m).push_back (newregs.at(r)); } + for (int r=0; r<nregs; ++r) { + int const m = superregs.at(r).map; + assert (m>=0 and m<nmaps); + superregss.at(m).push_back (superregs.at(r)); + } // Check sizes for (int m=0; m<nmaps; ++m) { assert (int(regss.at(m).size()) == myncomps.at(m)); + assert (int(superregss.at(m).size()) == mynregs.at(m)); } if (DEBUG) cout << "SRMA exit" << endl; @@ -1477,7 +1521,7 @@ namespace Carpet { void MakeMultigridBoxes (cGH const * const cctkGH, int const m, - vector<vector<region_t> > const & regss, + vector<vector<region_t> > const & regss, vector<vector<vector<region_t> > > & regsss) { regsss.resize (mglevels); diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index 5e34fcd30..b2bf94b70 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -128,7 +128,7 @@ namespace Carpet { static void find_processor_decomposition (cGH const * cctkGH, - vector<vector<vector<region_t> > > & regsss, + vector<vector<vector<region_t> > > & superregsss, vector<vector<vector<vector<region_t> > > > & regssss); static void @@ -204,7 +204,7 @@ namespace Carpet { // Do not call Util_GetHostName. Certain InfiniBand libraries // do not allow calling fork or exec, and getting the host name // seems also not allowed. It leads to random crashes. - if (verbose) { + if (verbose or veryverbose) { // Collect host names char hostnamebuf[1000]; Util_GetHostName (hostnamebuf, sizeof hostnamebuf); @@ -572,13 +572,13 @@ namespace Carpet { { DECLARE_CCTK_PARAMETERS; - vector<vector<vector<region_t> > > regsss (maps); + vector<vector<vector<region_t> > > superregsss (maps); for (int m=0; m<maps; ++m) { - set_base_extent (m, regsss.at(m)); + set_base_extent (m, superregsss.at(m)); } vector<vector<vector<vector<region_t> > > > regssss; - find_processor_decomposition (cctkGH, regsss, regssss); + find_processor_decomposition (cctkGH, superregsss, regssss); for (int m=0; m<maps; ++m) { @@ -586,7 +586,7 @@ namespace Carpet { CheckRegions (regssss.at(m)); // Recompose grid hierarchy - vhh.at(m)->regrid (regssss.at(m)); + vhh.at(m)->regrid (superregsss.at(m), regssss.at(m)); int const rl = 0; vhh.at(m)->recompose (rl, false); @@ -827,21 +827,25 @@ namespace Carpet { ivect const & convoffsets) { // Set refinement structure for scalars and arrays - vector<region_t> regs(1); - int const m=0; - regs.at(0).extent = arrdata.at(group).at(m).hh->baseextents.at(0).at(0); - regs.at(0).outer_boundaries = b2vect(true); - regs.at(m).map = m; + vector<region_t> superregs(1); + int const m = 0; + int const rl = 0; + int const c = 0; + superregs.at(c).extent = + arrdata.at(group).at(m).hh->baseextents.at(rl).at(c); + superregs.at(c).outer_boundaries = b2vect(true); + superregs.at(c).map = m; + vector<region_t> regs; // Split it into components, one for each processor switch (gdata.disttype) { case CCTK_DISTRIB_DEFAULT: { - SplitRegions_Automatic (cctkGH, regs); + SplitRegions_Automatic (cctkGH, superregs, regs); break; } case CCTK_DISTRIB_CONSTANT: { int const d = gdata.dim==0 ? 0 : gdata.dim-1; - SplitRegions_AlongDir (cctkGH, regs, d); + SplitRegions_AlongDir (cctkGH, superregs, regs, d); break; } default: @@ -849,8 +853,10 @@ namespace Carpet { } // Only one refinement level + vector<vector<region_t> > superregss(1); + superregss.at(rl) = superregs; vector<vector<region_t> > regss(1); - regss.at(0) = regs; + regss.at(rl) = regs; // Create all multigrid levels vector<vector<vector<region_t> > > regsss (mglevels); @@ -895,7 +901,7 @@ namespace Carpet { char * const groupname = CCTK_GroupName (group); assert (groupname); Checkpoint ("Recomposing grid array group \"%s\"...", groupname); - arrdata.at(group).at(0).hh->regrid (regsss); + arrdata.at(group).at(0).hh->regrid (superregss, regsss); arrdata.at(group).at(0).hh->recompose (0, false); Checkpoint ("Done recomposing grid array group \"%s\".", groupname); free (groupname); @@ -1445,7 +1451,7 @@ namespace Carpet { void find_processor_decomposition (cGH const * const cctkGH, - vector<vector<vector<region_t> > > & regsss, + vector<vector<vector<region_t> > > & superregsss, vector<vector<vector<vector<region_t> > > > & regssss) { DECLARE_CCTK_PARAMETERS; @@ -1459,11 +1465,13 @@ namespace Carpet { for (int m=0; m<maps; ++m) { int const rl=0; + vector<vector<region_t> > regss(1); + // Distribute onto the processors - SplitRegions (cctkGH, regsss.at(m).at(rl)); + SplitRegions (cctkGH, superregsss.at(m).at(rl), regss.at(rl)); // Create all multigrid levels - MakeMultigridBoxes (cctkGH, m, regsss.at(m), regssss.at(m)); + MakeMultigridBoxes (cctkGH, m, regss, regssss.at(m)); } // for m } else { @@ -1471,14 +1479,18 @@ namespace Carpet { int const rl=0; - vector<vector<region_t> > regss(maps); + vector<vector<region_t> > superregss(maps); for (int m=0; m<maps; ++m) { - regss.at(m) = regsss.at(m).at(rl); + superregss.at(m) = superregsss.at(m).at(rl); } - SplitRegionsMaps (cctkGH, regss); + vector<vector<region_t> > regss(maps); + SplitRegionsMaps (cctkGH, superregss, regss); + vector<vector<vector<region_t> > > regsss(maps); for (int m=0; m<maps; ++m) { + superregsss.at(m).at(rl) = superregss.at(m); + regsss.at(m).resize(1); regsss.at(m).at(rl) = regss.at(m); } diff --git a/Carpet/Carpet/src/functions.hh b/Carpet/Carpet/src/functions.hh index bbe2ff724..ef86ca3c9 100644 --- a/Carpet/Carpet/src/functions.hh +++ b/Carpet/Carpet/src/functions.hh @@ -22,6 +22,7 @@ namespace Carpet { using namespace std; + using namespace CarpetLib; int SyncGroupsByDirI (const cGH* cctkGH, int num_groups, const int* groups, const int* directions); @@ -102,6 +103,7 @@ namespace Carpet { void RegridMap (cGH const * cctkGH, int m, + gh::rregs const & supeerregss, gh::mregs const & regsss); void PostRegrid (cGH const * cctkGH); @@ -137,28 +139,30 @@ namespace Carpet { // Functions for recomposing the grid hierarchy void SplitRegions (cGH const * cctkGH, + vector<region_t> & superregs, vector<region_t> & regs); void SplitRegions_AlongZ (cGH const * cctkGH, + vector<region_t> & superregs, vector<region_t> & regs); void SplitRegions_AlongDir (cGH const * cctkGH, + vector<region_t> & superregs, vector<region_t> & regs, - int const dir); + int dir); void SplitRegions_Automatic (cGH const * cctkGH, + vector<region_t> & superregs, vector<region_t> & regs); - void - SplitRegions (cGH const * cctkGH, - vector<region_t> & regss); - - void - SplitRegionsMaps_Automatic (cGH const * cctkGH, - vector<vector<region_t> > & regss); void SplitRegionsMaps (cGH const * cctkGH, + vector<vector<region_t> > & superregss, vector<vector<region_t> > & regss); + void + SplitRegionsMaps_Automatic (cGH const * cctkGH, + vector<vector<region_t> > & superregss, + vector<vector<region_t> > & regss); void MakeMultigridBoxes (cGH const * cctkGH, diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc index 5242daedc..7abdd1232 100644 --- a/Carpet/CarpetIOHDF5/src/Input.cc +++ b/Carpet/CarpetIOHDF5/src/Input.cc @@ -108,7 +108,7 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS) // checkpoint file, or if the number of maps is wrong assert (int(fileset.grid_structure.size()) == maps); - vector<vector<vector<region_t> > > regsss = fileset.grid_structure; + vector<vector<vector<region_t> > > superregsss = fileset.grid_structure; vector<vector<vector<vector<region_t> > > > regssss (maps); int type; @@ -122,26 +122,28 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS) // Distribute each map independently for (int m = 0; m < maps; ++ m) { - - vector <vector <region_t> > & regss = regsss.at(m); + vector <vector <region_t> > & superregss = superregsss.at(m); // Make multiprocessor aware - for (size_t rl = 0; rl < regss.size(); ++ rl) { - Carpet::SplitRegions (cctkGH, regss.at(rl)); + vector <vector <region_t> > regss(superregss.size()); + for (size_t rl = 0; rl < superregss.size(); ++ rl) { + Carpet::SplitRegions (cctkGH, superregss.at(rl), regss.at(rl)); } // for rl // Make multigrid aware - Carpet::MakeMultigridBoxes (cctkGH, Carpet::map, regss, regssss.at(m)); + Carpet::MakeMultigridBoxes (cctkGH, m, regss, regssss.at(m)); } // for m - } else { + } else { // if regrid_in_level_mode // Distribute all maps at the same time + vector<vector<vector<region_t> > > regsss(maps); + // Count levels vector <int> rls (maps); for (int m = 0; m < maps; ++ m) { - rls.at(m) = regsss.at(m).size(); + rls.at(m) = superregsss.at(m).size(); } int maxrl = 0; for (int m = 0; m < maps; ++ m) { @@ -149,30 +151,33 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS) } // All maps must have the same number of levels for (int m = 0; m < maps; ++ m) { + superregsss.at(m).resize (maxrl); regsss.at(m).resize (maxrl); } // Make multiprocessor aware for (int rl = 0; rl < maxrl; ++ rl) { - vector <vector <region_t> > regss (maps); + vector <vector <region_t> > superregss (maps); for (int m = 0; m < maps; ++ m) { - regss.at(m) = regsss.at(m).at(rl); + superregss.at(m) = superregsss.at(m).at(rl); } - Carpet::SplitRegionsMaps (cctkGH, regss); + vector <vector <region_t> > regss (maps); + Carpet::SplitRegionsMaps (cctkGH, superregss, regss); for (int m = 0; m < maps; ++ m) { + superregsss.at(m).at(rl) = superregss.at(m); regsss.at(m).at(rl) = regss.at(m); } } // for rl - // Make multigrid aware + // Make multigrid aware Carpet::MakeMultigridBoxesMaps (cctkGH, regsss, regssss); - } // if + } // if regrid_in_level_mode for (int m = 0; m < maps; ++ m) { // Regrid - RegridMap (cctkGH, m, regssss.at(m)); + RegridMap (cctkGH, m, superregsss.at(m), regssss.at(m)); // Set time hierarchy correctly after RegridMap created it for (int ml = 0; ml < mglevels; ++ ml) { diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl index 98a160bf9..44b7729d5 100644 --- a/Carpet/CarpetLib/interface.ccl +++ b/Carpet/CarpetLib/interface.ccl @@ -7,6 +7,7 @@ includes header: dist.hh in dist.hh includes header: bbox.hh in bbox.hh includes header: bboxset.hh in bboxset.hh +includes header: fulltree.hh in fulltree.hh includes header: region.hh in region.hh includes header: vect.hh in vect.hh diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index 13ef01301..357061ee2 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.cc @@ -224,6 +224,7 @@ ostream& output (ostream& os, const vector<T>& v) { #include "bbox.hh" #include "bboxset.hh" #include "dh.hh" +#include "fulltree.hh" #include "gdata.hh" #include "ggf.hh" #include "region.hh" @@ -245,6 +246,7 @@ template size_t memoryof (vector<int> const & v); template size_t memoryof (vector<CCTK_REAL> const & v); template size_t memoryof (vector<bbox<int,3> > const & v); template size_t memoryof (vector<vect<int,3> > const & v); +template size_t memoryof (vector<fulltree <int,3,pseudoregion_t> *> const & f); template size_t memoryof (vector<pseudoregion_t> const & v); template size_t memoryof (vector<region_t> const & v); template size_t memoryof (vector<sendrecv_pseudoregion_t> const & v); diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh index 59d22f4c8..1152c154a 100644 --- a/Carpet/CarpetLib/src/defs.hh +++ b/Carpet/CarpetLib/src/defs.hh @@ -57,14 +57,19 @@ const int dim = 3; // Some shortcuts for type names +template<typename T, int D> class vect; template<typename T, int D> class bbox; template<typename T, int D> class bboxset; -template<typename T, int D> class vect; +template<typename T, int D, typename P> class fulltree; + +struct pseudoregion_t; +struct region_t; typedef vect<bool,dim> bvect; typedef vect<int,dim> ivect; typedef bbox<int,dim> ibbox; typedef bboxset<int,dim> ibset; +typedef fulltree<int,dim,pseudoregion_t> ipfulltree; // (Try to replace these by b2vect and i2vect) typedef vect<vect<bool,2>,dim> bbvect; diff --git a/Carpet/CarpetLib/src/fulltree.cc b/Carpet/CarpetLib/src/fulltree.cc new file mode 100644 index 000000000..b114dee19 --- /dev/null +++ b/Carpet/CarpetLib/src/fulltree.cc @@ -0,0 +1,492 @@ +#include <algorithm> +#include <cmath> +#include <iostream> + +#include "defs.hh" +#include "region.hh" + +#include "fulltree.hh" + + + +#if 0 +// Create an empty tree +template <typename T, int D, typename P> +fulltree<T,D,P>::fulltree () + : type (type_empty) +{ + assert (invariant()); +} +#endif + + + +// Create a tree branch from a list of bounds and subtrees +template <typename T, int D, typename P> +fulltree<T,D,P>::fulltree (int dir_, vector <T> const & bounds_, + vector <fulltree *> const & subtrees_) + : type (type_branch), dir (dir_), bounds (bounds_), subtrees (subtrees_) +{ + assert (dir>=0 and dir<D); + assert (bounds.size() == subtrees.size() + 1); + assert (invariant()); +} + + + +// Create a tree leaf from a payload +template <typename T, int D, typename P> +fulltree<T,D,P>::fulltree (P const & p_) + : type (type_leaf), p (p_) +{ + assert (invariant()); +} + + + +// Create a tree as copy from another tree +template <typename T, int D, typename P> +fulltree<T,D,P>::fulltree (fulltree const & t) + : type (t.type) +{ + switch (type) { + case type_empty: + // do nothing + break; + case type_branch: + dir = t.dir; + bounds = t.bounds; + subtrees.resize (t.subtrees.size()); + for (size_t i=0; i<subtrees.size(); ++i) { + subtrees.AT(i) = new fulltree (*t.subtrees.AT(i)); + } + break; + case type_leaf: + p = t.p; + break; + default: + assert (0); + } + assert (invariant()); +} + + + +// Copy a tree +template <typename T, int D, typename P> +fulltree<T,D,P> & +fulltree<T,D,P>::operator= (fulltree const & t) +{ + assert (invariant()); + if (is_branch()) { + for (size_t i=0; i<subtrees.size(); ++i) { + delete subtrees.AT(i); + } + } + type = t.type; + switch (type) { + case type_empty: + // do nothing + break; + case type_branch: + dir = t.dir; + bounds = t.bounds; + subtrees.resize (t.subtrees.size()); + for (size_t i=0; i<subtrees.size(); ++i) { + subtrees.AT(i) = new fulltree (*t.subtrees.AT(i)); + } + break; + case type_leaf: + p = t.p; + break; + default: + assert (0); + } + assert (invariant()); + return *this; +} + + + +// Delete a tree (and its subtrees) +template <typename T, int D, typename P> +fulltree<T,D,P>::~fulltree () +{ + assert (invariant()); + if (is_branch()) { + for (size_t i=0; i<subtrees.size(); ++i) { + delete subtrees.AT(i); + } + } +} + + + +// Compare trees +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::operator== (fulltree const & t) const +{ + assert (invariant()); + assert (t.invariant()); + if (type != t.type) return false; + switch (type) { + case type_empty: + break; + case type_branch: + if (dir != t.dir) return false; + if (bounds != t.bounds) return false; + if (subtrees != t.subtrees) return false; + break; + case type_leaf: + return p == t.p; + default: + assert (0); + } + return true; +} + + + +// Invariant +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::invariant () const +{ + return empty() or is_branch() or is_leaf(); +} + + + +// Find the leaf payload corresponding to a position +template <typename T, int D, typename P> +P const * +fulltree<T,D,P>::search (tvect const & ipos) const +{ + assert (not empty()); + // cout << "fulltree::search ipos=" << ipos << endl; + if (is_leaf()) return & p; + int const i = asearch (ipos[dir], bounds); + // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl; + if (i<0 or i>=int(subtrees.size())) return NULL; // not found + return subtrees.AT(i)->search(ipos); +} + +template <typename T, int D, typename P> +P * +fulltree<T,D,P>::search (tvect const & ipos) +{ + assert (not empty()); + // cout << "fulltree::search ipos=" << ipos << endl; + if (is_leaf()) return & p; + int const i = asearch (ipos[dir], bounds); + // cout << "fulltree::search i=" << i << " size=" << subtrees.size() << endl; + if (i<0 or i>=int(subtrees.size())) return NULL; // not found + return subtrees.AT(i)->search(ipos); +} + + + +// Const iterator +template <typename T, int D, typename P> +fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_) + : f(f_), i(0), it(0) +{ + if (f.is_branch()) { + assert (f.subtrees.size() > 0); + it = new const_iterator (* f.subtrees.at(i)); + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P>::const_iterator::const_iterator (fulltree const & f_, int) + : f(f_), it(0) +{ + if (f.empty()) { + i = 0; + } else if (f.is_leaf()) { + i = 1; + } else { + i = f.subtrees.size(); + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P>::const_iterator::~const_iterator () +{ + if (it) { + delete it; + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P> const & +fulltree<T,D,P>::const_iterator::operator* () const +{ + assert (not f.empty()); + if (f.is_leaf()) { + return f; + } else { + return **it; + } +} + +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::const_iterator::operator== (const_iterator const & it2) const +{ + // assert (f == it2.f); + assert (&f == &it2.f); + if (i != it2.i) return false; + if (it == 0 and it2.it == 0) return true; + if (it == 0 or it2.it == 0) return false; + return *it == *it2.it; +} + +template <typename T, int D, typename P> +typename fulltree<T,D,P>::const_iterator & +fulltree<T,D,P>::const_iterator::operator++ () +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + ++ i; + } else { + ++ *it; + if ((*it).done()) { + delete it; + it = 0; + ++ i; + if (not done()) { + // to do: use a new function "reset iterator" instead + it = new const_iterator (* f.subtrees.at(i)); + } + } + } + return *this; +} + +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::const_iterator::done () const +{ + if (f.empty()) { + return true; + } else if (f.is_leaf()) { + return i > 0; + } else { + return i == f.subtrees.size(); + } +} + + + +// Non-const iterator +template <typename T, int D, typename P> +fulltree<T,D,P>::iterator::iterator (fulltree & f_) + : f(f_), i(0), it(0) +{ + if (f.is_branch()) { + assert (f.subtrees.size() > 0); + it = new iterator (* f.subtrees.at(i)); + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P>::iterator::iterator (fulltree & f_, int) + : f(f_), it(0) +{ + if (f.empty()) { + i = 0; + } else if (f.is_leaf()) { + i = 1; + } else { + i = f.subtrees.size(); + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P>::iterator::~iterator () +{ + if (it) { + delete it; + } +} + +template <typename T, int D, typename P> +fulltree<T,D,P> & +fulltree<T,D,P>::iterator::operator* () +{ + assert (not f.empty()); + if (f.is_leaf()) { + return f; + } else { + return **it; + } +} + +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::iterator::operator== (iterator const & it2) const +{ + // assert (f == it2.f); + assert (&f == &it2.f); + if (i != it2.i) return false; + if (it == 0 and it2.it == 0) return true; + if (it == 0 or it2.it == 0) return false; + return *it == *it2.it; +} + +template <typename T, int D, typename P> +typename fulltree<T,D,P>::iterator & +fulltree<T,D,P>::iterator::operator++ () +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + ++ i; + } else { + ++ *it; + if ((*it).done()) { + delete it; + it = 0; + ++ i; + if (not done()) { + // to do: use a new function "reset iterator" instead + it = new iterator (* f.subtrees.at(i)); + } + } + } + return *this; +} + +template <typename T, int D, typename P> +bool +fulltree<T,D,P>::iterator::done () const +{ + if (f.empty()) { + return true; + } else if (f.is_leaf()) { + return i > 0; + } else { + return i == f.subtrees.size(); + } +} + + + +// Memory usage +template <typename T, int D, typename P> +size_t +fulltree<T,D,P>::memory () const +{ + size_t size = sizeof *this; + if (is_branch()) { + size += memoryof (bounds) + memoryof (subtrees); + for (typename vector <fulltree *>::const_iterator + i = subtrees.begin(); i != subtrees.end(); ++ i) + { + size += memoryof(**i); + } + } + return size; +} + + + +// Output helper +template <typename T, int D, typename P> +void +fulltree<T,D,P>::output (ostream & os) const +{ + os << "fulltree{"; + if (empty()) { + os << "empty"; + } else if (is_branch()) { + os << "branch:" + << "dir=" << dir << "," + << "subtrees=["; + for (size_t i=0; i<subtrees.size(); ++i) { + os << bounds.at(i) << ":[" << i << "]=" << subtrees.at(i) << ":"; + } + os << bounds.at(subtrees.size()) << "]"; + } else { + os << "leaf:" + << "payload=" << p; + } + os << "}"; +} + + + +// Generic arithmetic search +template <typename T> +static +int asearch (T const t, vector <T> const & ts) +{ + // cout << "fulltree::asearch t=" << t << " ts=" << ts << endl; + int imin = 0; + int imax = int(ts.size()) - 1; + if (imax <= imin) { + // cout << "fulltree::asearch: no values" << endl; + return imin; + } + T tmin = ts.AT(imin); + // cout << "fulltree::asearch: imin=" << imin << " tmin=" << tmin << endl; + if (t < tmin) { + // cout << "fulltree::asearch: below minimum" << endl; + return -1; + } + T tmax = ts.AT(imax); + // cout << "fulltree::asearch: imax=" << imax << " tmax=" << tmax << endl; + if (t >= tmax) { + // cout << "fulltree::asearch: above maximum" << endl; + return imax; + } + int isize = imax - imin; + for (;;) { + // check loop invariants + assert (imin < imax); + assert (t >= tmin); + assert (t < tmax); + // cout << "fulltree::asearch t=" << t << " imin=" << imin << " imax=" << imax << endl; + if (imax == imin+1) { + // cout << "fulltree::asearch: found value" << endl; + return imin; + } + assert (tmax > tmin); // require that ts is strictly ordered + CCTK_REAL const rguess = + (imax - imin) * CCTK_REAL(t - tmin) / (tmax - tmin); + int const iguess = imin + max (1, int(floor(rguess))); + // handle round-off errors + if (iguess == imax) { + // cout << "fulltree::asearch: accidental hit after roundoff error" << endl; + return imax - 1; + } + assert (iguess > imin and iguess < imax); + T const tguess = ts.AT(iguess); + // cout << "fulltree::asearch: intersecting at iguess=" << iguess << " tguess=" << tguess << endl; + if (t < tguess) { + imax = iguess; + tmax = tguess; + // cout << "fulltree::asearch: new imax=" << imax << " tmax=" << tmax << endl; + } else { + imin = iguess; + tmin = tguess; + // cout << "fulltree::asearch: new imin=" << imin << " tmin=" << tmin << endl; + } + // check loop variant + int const newisize = imax - imin; + assert (newisize < isize); + isize = newisize; + } +} + + + +template class fulltree <int, dim, pseudoregion_t>; + +template size_t memoryof (fulltree <int, dim, pseudoregion_t> const & f); +template size_t memoryof (fulltree <int, dim, pseudoregion_t> * const & f); + +template ostream & operator<< (ostream & os, fulltree <int, dim, pseudoregion_t> const & f); diff --git a/Carpet/CarpetLib/src/fulltree.hh b/Carpet/CarpetLib/src/fulltree.hh new file mode 100644 index 000000000..6ffca209e --- /dev/null +++ b/Carpet/CarpetLib/src/fulltree.hh @@ -0,0 +1,181 @@ +#ifndef FULLTREE_HH +#define FULLTREE_HH + +#include <cassert> +#include <cmath> +#include <cstdlib> +#include <iostream> +#include <vector> + +#include <vect.hh> + +using namespace std; + + + +// This is a "full tree" data structure, i.e., a tree data structure +// which decomposes a cuboid domain into a set of non-overlapping +// cuboid subdomains. It is an n-ary tree, i.e., each tree node can +// have arbitrarily many subtrees. Each node splits a domain in +// exactly one direction. Subdomains cannot be empty. + +// All intervals are closed at the lower and open at the upper +// boundary. This makes it possible to combine adjacent such +// intervals, obtaining again an interval with this property. (In +// particular, if bboxes are stored, the upper bound of bboxes must be +// increased by the stride to arrive at such intervals.) + + +// Generic arithmetic search +template <typename T> +static +int asearch (T t, vector <T> const & ts); + + + +// T: index space (usually integer, or maybe real) +// D: number of dimensions (rank) +// P: payload (e.g. region_t) +template <typename T, int D, typename P> +class fulltree { + +public: + + // Short name for a small vector + typedef vect <T, D> tvect; + +private: + + enum tree_type_t { type_empty, type_branch, type_leaf } type; + + // Direction in which the node is split (0 <= dir < D) + int dir; + + // If this is a branch: + // n+1 bounds, splitting the domain + vector <T> bounds; // [n+1] + // n pointers to subtrees + vector <fulltree *> subtrees; // [n] + + // If this is a leaf: + // just the payload + P p; + +public: + +#if 0 + // Create an empty tree + fulltree (); +#endif + + // Create a tree branch from a list of bounds and subtrees + fulltree (int dir_, vector <T> const & bounds_, + vector <fulltree *> const & subtrees_); + + // Create a tree leaf from a payload + fulltree (P const & p_); + + // Create a tree as copy from another tree + fulltree (fulltree const & t); + + // Copy a tree + fulltree & operator= (fulltree const & t); + + // Delete a tree (and its subtrees) + ~fulltree (); + + // Inquire tree properties + bool empty() const { return type == type_empty; } + bool is_branch() const { return type == type_branch; } + bool is_leaf() const { return type == type_leaf; } + + // Compare trees + bool operator== (fulltree const & t) const; + bool operator!= (fulltree const & t) const + { return not (*this == t); } + + // Invariant + bool invariant () const; + + // Access the payload + P const & payload () const { assert (is_leaf()); return p; } + P & payload () { assert (is_leaf()); return p; } + + // Find the leaf payload corresponding to a position + P const * search (tvect const & ipos) const; + P * search (tvect const & ipos); + + class const_iterator { + fulltree const & f; + size_t i; + const_iterator * it; + public: + const_iterator (fulltree const & f_); + const_iterator (fulltree const & f_, int); + ~const_iterator (); + fulltree const & operator* () const; + bool operator== (const_iterator const & it2) const; + bool operator!= (const_iterator const & it2) const + { return not (*this == it2); } + const_iterator & operator++ (); + bool done () const; + }; + + class iterator { + fulltree & f; + size_t i; + iterator * it; + public: + iterator (fulltree & f_); + iterator (fulltree & f_, int); + ~iterator (); + fulltree & operator* (); + bool operator== (iterator const & it2) const; + bool operator!= (iterator const & it2) const + { return not (*this == it2); } + iterator & operator++ (); + bool done () const; + }; + + const_iterator begin() const + { return const_iterator (*this); } + + const_iterator end() const + { return const_iterator (*this, 0); } + + iterator begin() + { return iterator (*this); } + + iterator end() + { return iterator (*this, 0); } + + // Memory usage + size_t memory () const; + + // Output helper + void output (ostream & os) const; +}; + + + +// Memory usage +template <typename T, int D, typename P> +inline size_t memoryof (fulltree<T,D,P> const & f) { return f.memory(); } + +template <typename T, int D, typename P> +inline size_t memoryof (fulltree<T,D,P> * const & fp) { return sizeof fp; } + + + +// Output +template <typename T, int D, typename P> +ostream & +operator<< (ostream & os, fulltree<T,D,P> const & f) +{ + f.output (os); + return os; +} + + + +#endif // #ifndef FULLTREE_HH diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index b6eefa499..ae42668d8 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -512,8 +512,8 @@ transfer_from_all (comm_state & state, pseudoregion_t const & precv = (* ipsendrecv).recv; ibbox const & send = psend.extent; ibbox const & recv = precv.extent; - int const c2 = psend.processor; - int const c1 = precv.processor; + int const c2 = psend.component; + int const c1 = precv.component; // Source and destination data gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1); diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index e3c351b9d..b2c86b8db 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -62,7 +62,7 @@ protected: public: const int vectorlength; // vector length const int vectorindex; // index of *this - const ggf* vectorleader; // first vector element + ggf* const vectorleader; // first vector element private: mdata oldstorage; // temporary storage diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index fa9ff2125..9b4a5f443 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -15,6 +15,8 @@ using namespace std; +using namespace CarpetLib; + // Constructors @@ -76,10 +78,12 @@ gh:: // Modifiers void gh:: -regrid (mregs const & regs) +regrid (rregs const & superregs, mregs const & regs) { DECLARE_CCTK_PARAMETERS; + superregions = superregs; + // Save the grid hierarchy oldregions.clear (); swap (oldregions, regions); diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh index 0a283b4cf..ffac38804 100644 --- a/Carpet/CarpetLib/src/gh.hh +++ b/Carpet/CarpetLib/src/gh.hh @@ -47,6 +47,10 @@ public: // should be readonly vector<vector<ibbox> > baseextents; // [ml][rl] const i2vect boundary_width; + // Extents of the regions before distributing them over the + // processors + rregs superregions; + mregs regions; // extents and properties of all grids mregs oldregions; // a copy, used during regridding @@ -65,7 +69,7 @@ public: ~gh (); // Modifiers - void regrid (mregs const & regs); + void regrid (rregs const & superregs, mregs const & regs); bool recompose (int rl, bool do_prolongate); private: diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index aa9b399c0..88f6261ce 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -8,6 +8,7 @@ SRCS = bbox.cc \ defs.cc \ dh.cc \ dist.cc \ + fulltree.cc \ gdata.cc \ gf.cc \ ggf.cc \ diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc index caf535598..0230d373d 100644 --- a/Carpet/CarpetLib/src/region.cc +++ b/Carpet/CarpetLib/src/region.cc @@ -1,5 +1,8 @@ +#include <cassert> +#include <cstdlib> #include <iostream> +#include "bboxset.hh" #include "defs.hh" #include "region.hh" @@ -7,6 +10,69 @@ using namespace std; +region_t::region_t () + : processor (-1), processors (NULL) +{ + assert (invariant()); +} + +region_t::region_t (region_t const & a) +{ + assert (a.invariant()); + extent = a.extent; + outer_boundaries = a.outer_boundaries; + map = a.map; + processor = a.processor; + if (a.processors == NULL) { + processors = NULL; + } else { + processors = new ipfulltree (*a.processors); + } + assert (invariant()); +} + +region_t & +region_t::operator= (region_t const & a) +{ + assert (invariant()); + if (processors != NULL) { + delete processors; + } + assert (a.invariant()); + extent = a.extent; + outer_boundaries = a.outer_boundaries; + map = a.map; + processor = a.processor; + if (a.processors == NULL) { + processors = NULL; + } else { + processors = new ipfulltree (*a.processors); + } + assert (invariant()); + return *this; +} + + +region_t::~region_t () +{ + assert (invariant()); + if (processors != NULL) { + delete processors; + } +} + + + +bool +region_t::invariant () const +{ + if (processor >= 0 and processors != NULL) return false; + return true; +} + + + +// Compare two regions for equality. bool operator== (region_t const & a, region_t const & b) { @@ -14,7 +80,113 @@ operator== (region_t const & a, region_t const & b) a.extent == b.extent and all (all (a.outer_boundaries == b.outer_boundaries)) and a.map == b.map and - a.processor == b.processor; + a.processor == b.processor and + ((a.processors == NULL and b.processors == NULL) or + (a.processors != NULL and b.processors != NULL and + *a.processors == *b.processors)); +} + + + +// Combine a collection of regions. Regions can be combined if they +// abutt on boundaries which are not outer boundaries, ignoring the +// processor distribution. This should lead to a canonical +// representations of collections of regions. +// +// We use vectors to represent the collection, but we could also use +// other containers. oldregs is read, newregs is added-to. newregs +// is not cleared. +void +combine_regions (vector<region_t> const & oldregs, + vector<region_t> & newregs) +{ + // Find the union of all regions' bounding boxes, and the union of + // all regions' outer boundaries. Represent the boundaries as the + // outermost layer of grid points of the corresponding bounding + // boxes. + int const m = oldregs.empty() ? -1 : oldregs.AT(0).map; + ibset comps; + ibset cobnds[2][dim]; + for (vector<region_t>::const_iterator + ri = oldregs.begin(); ri != oldregs.end(); ++ ri) + { + region_t const & reg = * ri; + assert (reg.map == m); + comps += reg.extent; + for (int f = 0; f < 2; ++ f) { + for (int d = 0; d < dim; ++ d) { + if (reg.outer_boundaries[f][d]) { + ibbox bnd = reg.extent; + ivect lo = bnd.lower(); + ivect up = bnd.upper(); + if (f==0) { + up[d] = lo[d]; + } else { + lo[d] = up[d]; + } + bnd = ibbox (lo, up, bnd.stride()); + cobnds[f][d] += bnd; + } + } + } + } + comps.normalize(); + for (int f = 0; f < 2; ++ f) { + for (int d = 0; d < dim; ++ d) { + cobnds[f][d].normalize(); + } + } + + // Reserve (generous) memory for the result + size_t const needsize = newregs.size() + comps.setsize(); + if (newregs.capacity() < needsize) { + newregs.reserve (1000 + 2 * needsize); + } + + // Insert the regions + for (ibset::const_iterator ci = comps.begin(); ci != comps.end(); ++ ci) { + ibbox const & c = * ci; + b2vect obnds; + for (int f = 0; f < 2; ++ f) { + for (int d = 0; d < dim; ++ d) { + obnds[f][d] = cobnds[f][d].intersects (c); + if (obnds[f][d]) { + ivect lo = c.lower(); + ivect up = c.upper(); + if (f) lo[d]=up[d]; else up[d]=lo[d]; + ibbox const cbnds (lo, up, c.stride()); + if (not ((cobnds[f][d] & ibset(c)) == ibset(cbnds))) { + cout << "cobnds[f][d] = " << cobnds[f][d] << endl + << "ibset(c) = " << ibset(c) << endl + << "(cobnds[f][d] & ibset(c)) = " << (cobnds[f][d] & ibset(c)) << endl + << "ibset(cbnds) = " << ibset(cbnds) << endl; + } + assert ((cobnds[f][d] & ibset(c)) == ibset(cbnds)); + } + } + } + + region_t reg; + reg.extent = c; + reg.outer_boundaries = obnds; + reg.map = m; + reg.processor = -1; + reg.processors = NULL; + newregs.push_back (reg); + } +} + + + +size_t memoryof (region_t const & reg) +{ + return + memoryof (reg.extent) + + memoryof (reg.outer_boundaries) + + memoryof (reg.map) + + memoryof (reg.processor) + + memoryof (reg.processors) + + (reg.processors != NULL ? memoryof (*reg.processors) : 0); } @@ -59,6 +231,8 @@ operator>> (istream & is, region_t & reg) skipws (is); consume (is, ')'); + reg.processors = NULL; + return is; } @@ -77,9 +251,20 @@ operator<< (ostream & os, region_t const & reg) +// Compare two pseudoregions for equality. +bool +operator== (pseudoregion_t const & a, pseudoregion_t const & b) +{ + return + a.extent == b.extent and + a.component == b.component; +} + + + ostream & operator<< (ostream & os, pseudoregion_t const & p) { - return os << p.extent << "/p:" << p.processor; + return os << p.extent << "/c:" << p.component; } ostream & operator<< (ostream & os, sendrecv_pseudoregion_t const & srp) diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh index 033e73a05..66037bdc7 100644 --- a/Carpet/CarpetLib/src/region.hh +++ b/Carpet/CarpetLib/src/region.hh @@ -2,20 +2,29 @@ #define REGION_HH #include <iostream> +#include <vector> #include "defs.hh" #include "bbox.hh" +#include "fulltree.hh" #include "vect.hh" // Region description struct region_t { - ibbox extent; // extent - b2vect outer_boundaries; // outer boundaries - int map; // map to which this - // region belongs - int processor; // processor number + ibbox extent; // extent + b2vect outer_boundaries; // outer boundaries + int map; // map to which this region belongs + int processor; // processor number + ipfulltree * processors; // processor decomposition + + region_t (); + region_t (region_t const & a); + region_t & operator= (region_t const & a); + ~region_t (); + + bool invariant () const; }; @@ -29,54 +38,62 @@ bool operator!= (region_t const & a, region_t const & b) -inline size_t memoryof (region_t const & reg) -{ - return - memoryof (reg.extent) + - memoryof (reg.outer_boundaries) + - memoryof (reg.map) + - memoryof (reg.processor); -} +void +combine_regions (vector<region_t> const & oldregs, + vector<region_t> & newregs); + + + +size_t memoryof (region_t const & reg); istream & operator>> (istream & is, region_t & reg); ostream & operator<< (ostream & os, region_t const & reg); -// A pseudoregion is almost a region; it is a bbox that lives on a -// certain processor. Pseudoregions are a compact way to store -// information about what processors needs to send data to what other -// processors during synchronisation or regridding. +// A pseudoregion is almost a region; it is a bbox that belongs to a +// certain component. Pseudoregions are a compact way to store +// information about what components needs to send data to what other +// components during synchronisation or regridding. struct pseudoregion_t { ibbox extent; - int processor; + int component; pseudoregion_t () { } - pseudoregion_t (ibbox const & extent_, int const processor_) - : extent (extent_), processor (processor_) + pseudoregion_t (ibbox const & extent_, int const component_) + : extent (extent_), component (component_) { } }; +bool operator== (pseudoregion_t const & a, pseudoregion_t const & b); +inline +bool operator!= (pseudoregion_t const & a, pseudoregion_t const & b) +{ + return not (a == b); +} + inline size_t memoryof (pseudoregion_t const & p) { return memoryof (p.extent) + - memoryof (p.processor); + memoryof (p.component); } ostream & operator<< (ostream & os, pseudoregion_t const & p); + + struct sendrecv_pseudoregion_t { pseudoregion_t send, recv; sendrecv_pseudoregion_t () { } - sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_processor, - ibbox const & recv_extent, int const recv_processor) - : send (pseudoregion_t (send_extent, send_processor)), - recv (pseudoregion_t (recv_extent, recv_processor)) + sendrecv_pseudoregion_t (ibbox const & send_extent, int const send_component, + ibbox const & recv_extent, int const recv_component) + : send (pseudoregion_t (send_extent, send_component)), + recv (pseudoregion_t (recv_extent, recv_component)) { } }; diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl index 1cece9380..ac9d5b112 100644 --- a/Carpet/CarpetRegrid/interface.ccl +++ b/Carpet/CarpetRegrid/interface.ccl @@ -52,9 +52,11 @@ USES FUNCTION ConvertFromPhysicalBoundary # The true prototype of the routine below: # int Carpet_Regrid (cGH const * cctkGH, +# gh::rregs * superregss, # gh::mregs * regsss, # int force); CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_POINTER IN superregss, \ CCTK_POINTER IN regsss, \ CCTK_INT IN force) PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid_Regrid LANGUAGE C diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc index 16c12063c..442224f2c 100644 --- a/Carpet/CarpetRegrid/src/automatic.cc +++ b/Carpet/CarpetRegrid/src/automatic.cc @@ -26,15 +26,12 @@ namespace CarpetRegrid { int Automatic (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; assert (refinement_levels >= 1); - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - const int vi = CCTK_VarIndex (errorvar); @@ -60,9 +57,6 @@ namespace CarpetRegrid { (cctkGH, hh, reflevel, min(reflevels+1, (int)refinement_levels), errorgf, regs); - // make multiprocessor aware - SplitRegions (cctkGH, regs); - if (regs.size() == 0) { // remove all finer levels regss.resize(reflevel+1); @@ -77,9 +71,6 @@ namespace CarpetRegrid { } } - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc index 105b9f72c..a7229fdad 100644 --- a/Carpet/CarpetRegrid/src/baselevel.cc +++ b/Carpet/CarpetRegrid/src/baselevel.cc @@ -19,14 +19,12 @@ namespace CarpetRegrid { int BaseLevel (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; assert (refinement_levels == 1); - assert (regsss.size() == 1); - return 0; } diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc index 5bda85c85..078dd0181 100644 --- a/Carpet/CarpetRegrid/src/centre.cc +++ b/Carpet/CarpetRegrid/src/centre.cc @@ -19,7 +19,7 @@ namespace CarpetRegrid { int Centre (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; @@ -28,9 +28,6 @@ namespace CarpetRegrid { // do nothing if the levels already exist if (reflevel == refinement_levels) return 0; - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - regss.resize (refinement_levels); bvect const symmetric (symmetry_x, symmetry_y, symmetry_z); @@ -40,7 +37,7 @@ namespace CarpetRegrid { ivect rlb = hh.baseextents.at(0).at(0).lower(); ivect rub = hh.baseextents.at(0).at(0).upper(); - assert (! smart_outer_boundaries); + assert (not smart_outer_boundaries); for (size_t rl=1; rl<regss.size(); ++rl) { @@ -69,16 +66,10 @@ namespace CarpetRegrid { regs.at(0).map = Carpet::map; - // make multiprocessor aware - SplitRegions (cctkGH, regs); - regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc index 7becee434..f32ec597c 100644 --- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc +++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc @@ -23,7 +23,7 @@ namespace CarpetRegrid { int ManualCoordinateList (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; int ierr; @@ -33,9 +33,6 @@ namespace CarpetRegrid { // do nothing if the levels already exist if (reflevel == refinement_levels && !tracking) return 0; - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - jjvect nboundaryzones, is_internal, is_staggered, shiftout; ierr = GetBoundarySpecification (2*dim, &nboundaryzones[0][0], &is_internal[0][0], @@ -116,7 +113,7 @@ namespace CarpetRegrid { } // for c } // for rl - } else { // if ! smart_outer_boundaries + } else { // if not smart_outer_boundaries vector<vector<bbvect> > newobss1; if (strcmp(outerbounds, "") != 0) { @@ -159,7 +156,7 @@ namespace CarpetRegrid { } } - } // if ! smart_outer_boundaries + } // if not smart_outer_boundaries for (int rl=1; rl<refinement_levels; ++rl) { @@ -224,17 +221,11 @@ namespace CarpetRegrid { } } // if merge_overlapping_components - - // make multiprocessor aware - SplitRegions (cctkGH, regs); regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc index e6d2d2ae0..1c4cea5b0 100644 --- a/Carpet/CarpetRegrid/src/manualcoordinates.cc +++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc @@ -20,7 +20,7 @@ namespace CarpetRegrid { int ManualCoordinates (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; @@ -32,9 +32,6 @@ namespace CarpetRegrid { // do nothing if the levels already exist if (reflevel == refinement_levels) return 0; - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - regss.resize (refinement_levels); vector<rvect> lower(3), upper(3); @@ -59,16 +56,10 @@ namespace CarpetRegrid { (cctkGH, hh, rl, refinement_levels, lower.at(rl-1), upper.at(rl-1), reg, regs); - // make multiprocessor aware - SplitRegions (cctkGH, regs); - regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc index 60968e59d..fe5bdb4ab 100644 --- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc +++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc @@ -23,7 +23,7 @@ namespace CarpetRegrid { int ManualGridpointList (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; @@ -32,9 +32,6 @@ namespace CarpetRegrid { // do nothing if the levels already exist if (reflevel == refinement_levels) return 0; - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - regss.resize (refinement_levels); vector<vector<ibbox> > newbbss; @@ -116,16 +113,10 @@ namespace CarpetRegrid { ext.lower(), ext.upper(), reg, regs); } - // make multiprocessor aware - SplitRegions (cctkGH, regs); - regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc index 6481247c9..3df1f2119 100644 --- a/Carpet/CarpetRegrid/src/manualgridpoints.cc +++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc @@ -22,7 +22,7 @@ namespace CarpetRegrid { int ManualGridpoints (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_PARAMETERS; @@ -34,9 +34,6 @@ namespace CarpetRegrid { // do nothing if the levels already exist if (reflevel == refinement_levels) return 0; - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - regss.resize (refinement_levels); vector<ivect> ilower(3), iupper(3); @@ -47,7 +44,7 @@ namespace CarpetRegrid { ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin); iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax); - assert (! smart_outer_boundaries); + assert (not smart_outer_boundaries); for (size_t rl=1; rl<regss.size(); ++rl) { @@ -61,16 +58,10 @@ namespace CarpetRegrid { (cctkGH, hh, rl,refinement_levels, ilower.at(rl-1), iupper.at(rl-1), reg, regs); - // make multiprocessor aware - SplitRegions (cctkGH, regs); - regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc index 41980ea45..075fdbb93 100644 --- a/Carpet/CarpetRegrid/src/moving.cc +++ b/Carpet/CarpetRegrid/src/moving.cc @@ -20,16 +20,13 @@ namespace CarpetRegrid { int Moving (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss) + gh::rregs & regss) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; assert (refinement_levels >= 1); - assert (regsss.size() >= 1); - vector<vector<region_t> > regss = regsss.at(0); - regss.resize (refinement_levels); bvect const symmetric (symmetry_x, symmetry_y, symmetry_z); @@ -60,16 +57,10 @@ namespace CarpetRegrid { ManualCoordinates_OneLevel (cctkGH, hh, rl, refinement_levels, rlb, rub, reg, regs); - // make multiprocessor aware - SplitRegions (cctkGH, regs); - regss.at(rl) = regs; } // for rl - // make multigrid aware - MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); - return 1; } diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc index 2dd2a33a9..0a5567e2b 100644 --- a/Carpet/CarpetRegrid/src/regrid.cc +++ b/Carpet/CarpetRegrid/src/regrid.cc @@ -19,6 +19,7 @@ namespace CarpetRegrid { using namespace Carpet; CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregss_, CCTK_POINTER const regsss_, CCTK_INT force) { @@ -26,6 +27,7 @@ namespace CarpetRegrid { const cGH * const cctkGH = (const cGH *) cctkGH_; + gh::rregs & superregss = * (gh::rregs *) superregss_; gh::mregs & regsss = * (gh::mregs *) regsss_; gh const & hh = *vhh.at(Carpet::map); @@ -151,44 +153,54 @@ namespace CarpetRegrid { if (CCTK_EQUALS(refined_regions, "none")) { - do_recompose = BaseLevel (cctkGH, hh, regsss); + do_recompose = BaseLevel (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "centre")) { - do_recompose = Centre (cctkGH, hh, regsss); + do_recompose = Centre (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "manual-gridpoints")) { do_recompose - = ManualGridpoints (cctkGH, hh, regsss); + = ManualGridpoints (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) { do_recompose - = ManualCoordinates (cctkGH, hh, regsss); + = ManualCoordinates (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) { do_recompose - = ManualGridpointList (cctkGH, hh, regsss); + = ManualGridpointList (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) { do_recompose - = ManualCoordinateList (cctkGH, hh, regsss); + = ManualCoordinateList (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "moving")) { - do_recompose = Moving (cctkGH, hh, regsss); + do_recompose = Moving (cctkGH, hh, superregss); } else if (CCTK_EQUALS(refined_regions, "automatic")) { - do_recompose = Automatic (cctkGH, hh, regsss); + do_recompose = Automatic (cctkGH, hh, superregss); } else { assert (0); } + // make multiprocessor aware + vector<vector<region_t> > regss(superregss.size()); + for (size_t rl=0; rl<superregss.size(); ++rl) { +#warning "TODO: delete .processors" + SplitRegions (cctkGH, superregss.at(rl), regss.at(rl)); + } + + // make multigrid aware + MakeMultigridBoxes (cctkGH, Carpet::map, regss, regsss); + assert (regsss.size() > 0); return do_recompose; diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh index b7d727317..f37c9cb1d 100644 --- a/Carpet/CarpetRegrid/src/regrid.hh +++ b/Carpet/CarpetRegrid/src/regrid.hh @@ -31,8 +31,10 @@ namespace CarpetRegrid { // Aliased functions // CCTK_INT CarpetRegrid_Regrid (const cGH * const cctkGH, +// gh<dim>::rregs * superregss, // gh<dim>::mregs * regsss); CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregss_, CCTK_POINTER const regsss_, CCTK_INT force); } @@ -41,15 +43,15 @@ namespace CarpetRegrid { int BaseLevel (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); int Centre (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); int ManualGridpoints (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); void ManualGridpoints_OneLevel (const cGH * const cctkGH, const gh & hh, @@ -62,7 +64,7 @@ namespace CarpetRegrid { int ManualCoordinates (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); void ManualCoordinates_OneLevel (const cGH * const cctkGH, const gh & hh, @@ -84,19 +86,19 @@ namespace CarpetRegrid { int ManualGridpointList (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); int ManualCoordinateList (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); int Moving (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); int Automatic (cGH const * const cctkGH, gh const & hh, - gh::mregs & regsss); + gh::rregs & regss); void Automatic_OneLevel (const cGH * const cctkGH, const gh & hh, diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl index aa6f799da..6eb913280 100644 --- a/Carpet/CarpetRegrid2/interface.ccl +++ b/Carpet/CarpetRegrid2/interface.ccl @@ -87,21 +87,25 @@ USES FUNCTION MultiPatch_ConvertFromPhysicalBoundary # The true prototype of the routine below: # int Carpet_Regrid (cGH const * cctkGH, +# gh::rregs * superregss, # gh::mregs * regsss, # int force); -CCTK_INT FUNCTION Carpet_Regrid \ - (CCTK_POINTER_TO_CONST IN cctkGH, \ - CCTK_POINTER IN regsss, \ +CCTK_INT FUNCTION Carpet_Regrid \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_POINTER IN superregss, \ + CCTK_POINTER IN regsss, \ CCTK_INT IN force) PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid2_Regrid LANGUAGE C # The true prototype of the routine below: # int Carpet_Regrid (cGH const * cctkGH, +# vector <gh::rregs> * superregsss, # vector <gh::mregs> * regssss, # int force) -CCTK_INT FUNCTION Carpet_RegridMaps \ - (CCTK_POINTER_TO_CONST IN cctkGH, \ - CCTK_POINTER IN regssss, \ +CCTK_INT FUNCTION Carpet_RegridMaps \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_POINTER IN superregsss, \ + CCTK_POINTER IN regssss, \ CCTK_INT IN force) PROVIDES FUNCTION Carpet_RegridMaps WITH CarpetRegrid2_RegridMaps LANGUAGE C diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc index c08705d29..eaf993bc0 100644 --- a/Carpet/CarpetRegrid2/src/regrid.cc +++ b/Carpet/CarpetRegrid2/src/regrid.cc @@ -239,11 +239,13 @@ namespace CarpetRegrid2 { extern "C" { CCTK_INT CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregss_, CCTK_POINTER const regsss_, CCTK_INT const force); CCTK_INT CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregsss_, CCTK_POINTER const regssss_, CCTK_INT const force); } @@ -759,6 +761,7 @@ namespace CarpetRegrid2 { CCTK_INT CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregss_, CCTK_POINTER const regsss_, CCTK_INT const force) { @@ -882,16 +885,15 @@ namespace CarpetRegrid2 { if (do_recompose) { - gh::mregs & regsss = * static_cast <gh::mregs *> (regsss_); + gh::rregs & superregss = * static_cast <gh::rregs *> (superregss_); + gh::mregs & regsss = * static_cast <gh::mregs *> (regsss_); - // Make multigrid unaware - vector <vector <region_t> > regss = regsss.at(0); - - Regrid (cctkGH, regss); + Regrid (cctkGH, superregss); // Make multiprocessor aware + vector <vector <region_t> > regss; for (size_t rl = 0; rl < regss.size(); ++ rl) { - Carpet::SplitRegions (cctkGH, regss.at(rl)); + Carpet::SplitRegions (cctkGH, superregss.at(rl), regss.at(rl)); } // for rl // Make multigrid aware @@ -919,6 +921,7 @@ namespace CarpetRegrid2 { CCTK_INT CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_POINTER const superregsss_, CCTK_POINTER const regssss_, CCTK_INT const force) { @@ -1040,23 +1043,21 @@ namespace CarpetRegrid2 { if (do_recompose) { + vector <gh::rregs> & superregsss = + * static_cast <vector <gh::rregs> *> (superregsss_); vector <gh::mregs> & regssss = * static_cast <vector <gh::mregs> *> (regssss_); - // Make multigrid unaware - vector< vector <vector <region_t> > > regsss (Carpet::maps); - for (int m = 0; m < maps; ++ m) { - regsss.at(m) = regssss.at(m).at(0); - } - BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { - Regrid (cctkGH, regsss.at(Carpet::map)); + Regrid (cctkGH, superregsss.at(Carpet::map)); } END_MAP_LOOP; + vector< vector <vector <region_t> > > regsss (Carpet::maps); + // Count levels vector <int> rls (maps); for (int m = 0; m < maps; ++ m) { - rls.at(m) = regsss.at(m).size(); + rls.at(m) = superregsss.at(m).size(); } int maxrl = 0; for (int m = 0; m < maps; ++ m) { @@ -1064,17 +1065,21 @@ namespace CarpetRegrid2 { } // All maps must have the same number of levels for (int m = 0; m < maps; ++ m) { + superregsss.at(m).resize (maxrl); regsss.at(m).resize (maxrl); } // Make multiprocessor aware for (int rl = 0; rl < maxrl; ++ rl) { - vector <vector <region_t> > regss (maps); + vector <vector <region_t> > superregss (maps); for (int m = 0; m < maps; ++ m) { - regss.at(m) = regsss.at(m).at(rl); + superregss.at(m) = superregsss.at(m).at(rl); } - Carpet::SplitRegionsMaps (cctkGH, regss); + vector <vector <region_t> > regss (maps); + Carpet::SplitRegionsMaps (cctkGH, superregss, regss); + for (int m = 0; m < maps; ++ m) { + superregsss.at(m).at(rl) = superregss.at(m); regsss.at(m).at(rl) = regss.at(m); } } // for rl |