diff options
Diffstat (limited to 'Carpet')
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 |