From 4fd01aecd88bd63bf3f6a79a63657915992a42b6 Mon Sep 17 00:00:00 2001 From: schnetter <> Date: Mon, 11 Mar 2002 12:17:00 +0000 Subject: Added stream input routines for some CarpetLib containers. Added stream input routines for some CarpetLib containers. The regridder now has to explicitly say which boundaries are outer, and which are internal. This will make outer boundaries on fine grid possible, and is also necessary when there are multiple grid patches. Started to add support for arbitrariliy many user-specified refinement regions. Not yet finished. The Carpet driver can now handle multiple grid patches. Added example files for multiple grid patches. They use initial data that does not "fit" the boundary conditions, and they don't use multiple refinement levels so far. Removed old and unused example files in CarpetLib. darcs-hash:20020311121709-07bb3-18594c42bd7a958ee0840d29e158a343208f5711.gz --- Carpet/CarpetRegrid/param.ccl | 22 +++- Carpet/CarpetRegrid/src/regrid.cc | 236 ++++++++++++++++++++++---------------- Carpet/CarpetRegrid/src/regrid.hh | 3 +- 3 files changed, 159 insertions(+), 102 deletions(-) (limited to 'Carpet/CarpetRegrid') diff --git a/Carpet/CarpetRegrid/param.ccl b/Carpet/CarpetRegrid/param.ccl index e6902bc80..951e7949b 100644 --- a/Carpet/CarpetRegrid/param.ccl +++ b/Carpet/CarpetRegrid/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn CarpetRegrid -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.3 2002/01/11 17:37:13 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.4 2002/03/11 13:17:15 schnetter Exp $ @@ -25,6 +25,8 @@ KEYWORD refined_regions "Regions where the grid is refined" "centre" :: "Refine around the centre of the grid only" "manual-gridpoints" :: "Refine the regions specified by integer grid points l[123]i[xyz]{min,max}" "manual-coordinates" :: "Refine the regions specified by coordinates l[123][xyz]{min,max}" + "manual-gridpoint-list" :: "" + "manual-coordinate-list" :: "" "automatic" :: "Refine automatically" } "centre" @@ -202,6 +204,24 @@ CCTK_REAL l3zmax "Upper boundary of level 3 box in z-direction" +# Refinement criteria for manual-gridpoint-list + +CCTK_STRING gridpoints "List of bounding box gridpoints" +{ + .* :: "[ [ ([,,]:[,,]:[,,]), ... ], ... ]" +} "" + + + +# Refinement criteria for manual-coordinate-list + +CCTK_STRING coordinates "List of bounding box coordinates" +{ + .* :: "[ [ ([,,]:[,,]:[,,]), ... ], ... ]" +} "" + + + # Refinement criteria for automatic refining CCTK_INT minwidth "Minimum width of refined region" diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc index 1b0b95d87..ad17ca833 100644 --- a/Carpet/CarpetRegrid/src/regrid.cc +++ b/Carpet/CarpetRegrid/src/regrid.cc @@ -19,7 +19,7 @@ #include "carpet.hh" #include "regrid.hh" -static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.8 2002/03/06 17:47:58 schnetter Exp $"; +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.9 2002/03/11 13:17:15 schnetter Exp $"; @@ -30,6 +30,13 @@ namespace CarpetRegrid { + typedef vect ivect; + typedef bbox ibbox; + + typedef vect,dim> bvect; + + + int CarpetRegridStartup () { RegisterRegridRoutine (CarpetRegridRegrid); @@ -40,6 +47,7 @@ namespace CarpetRegrid { int CarpetRegridRegrid (const cGH * const cctkGH, gh::rexts& bbsss, + gh::rbnds& obss, gh::rprocs& pss) { DECLARE_CCTK_PARAMETERS; @@ -56,9 +64,6 @@ namespace CarpetRegrid { #endif assert (mglevel == -1); - // Return if this is the finest possible level - if (reflevel == maxreflevels-1) return 0; - // Return if we want to regrid during initial data only, and this // is not the time for initial data if (regrid_every == 0 && cctkGH->cctk_iteration != 0) return 0; @@ -69,7 +74,7 @@ namespace CarpetRegrid { return 0; } - list > bbl; + list bbl; if (CCTK_EQUALS(refined_regions, "none")) { @@ -85,13 +90,13 @@ namespace CarpetRegrid { CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels"); } assert (refinement_levels<=4); - vector > lower(3), upper(3); - lower[0] = vect (l1ixmin, l1iymin, l1izmin); - upper[0] = vect (l1ixmax, l1iymax, l1izmax); - lower[1] = vect (l2ixmin, l2iymin, l2izmin); - upper[1] = vect (l2ixmax, l2iymax, l2izmax); - lower[2] = vect (l3ixmin, l3iymin, l3izmin); - upper[2] = vect (l3ixmax, l3iymax, l3izmax); + vector lower(3), upper(3); + lower[0] = ivect (l1ixmin, l1iymin, l1izmin); + upper[0] = ivect (l1ixmax, l1iymax, l1izmax); + lower[1] = ivect (l2ixmin, l2iymin, l2izmin); + upper[1] = ivect (l2ixmax, l2iymax, l2izmax); + lower[2] = ivect (l3ixmin, l3iymin, l3izmin); + upper[2] = ivect (l3ixmax, l3iymax, l3izmax); MakeRegions_AsSpecified (cctkGH, refinement_levels, lower, upper, bbl); } else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) { @@ -109,6 +114,14 @@ namespace CarpetRegrid { upper[2] = vect (l3xmax, l3ymax, l3zmax); MakeRegions_AsSpecified (cctkGH, refinement_levels, lower, upper, bbl); + } else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) { + + abort (); + + } else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) { + + abort (); + } else if (CCTK_EQUALS(refined_regions, "automatic")) { const int vi = CCTK_VarIndex (errorvar); @@ -135,9 +148,10 @@ namespace CarpetRegrid { abort(); } - vector > bbs; + // transform bbox list into bbox vector + vector bbs; bbs.reserve (bbl.size()); - for (list >::const_iterator it = bbl.begin(); + for (list::const_iterator it = bbl.begin(); it != bbl.end(); ++it) { bbs.push_back (*it); @@ -145,21 +159,40 @@ namespace CarpetRegrid { // TODO: ensure nesting properties - SplitRegions (cctkGH, bbs); + // TODO: set outer boundaries correctly + vector obs; + obs.resize (bbs.size()); + for (vector::iterator it = obs.begin(); + it != obs.end(); + ++it) { + *it = bvect(vect(false)); + } + + // make multiprocessor aware + SplitRegions (cctkGH, bbs, obs); + // make multigrid aware gh::cexts bbss; bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels); + // insert into grid hierarchy bbsss = hh->extents; + obss = hh->outer_boundaries; if (bbss.size() == 0) { + // remove all finer levels // TODO: this might not work bbsss.resize(reflevel+1); + obss.resize(reflevel+1); } else { assert (reflevel < (int)bbsss.size()); if (reflevel+1 == (int)bbsss.size()) { + // add a finer level bbsss.push_back (bbss); + obss.push_back (obs); } else { + // change a finer level bbsss[reflevel+1] = bbss; + obss[reflevel+1] = obs; } } @@ -170,7 +203,7 @@ namespace CarpetRegrid { - void MakeRegions_BaseLevel (const cGH* cctkGH, list >& bbl) + void MakeRegions_BaseLevel (const cGH* cctkGH, list& bbl) { assert (bbl.empty()); } @@ -181,27 +214,27 @@ namespace CarpetRegrid { // how the hierarchy should be refined. But the result of this // routine is rather arbitrary. void MakeRegions_RefineCentre (const cGH* cctkGH, const int reflevels, - list >& bbl) + list& bbl) { assert (bbl.empty()); if (reflevel+1 >= reflevels) return; // note: what this routine calls "ub" is "ub+str" elsewhere - vect rstr = hh->baseextent.stride(); - vect rlb = hh->baseextent.lower(); - vect rub = hh->baseextent.upper() + rstr; + ivect rstr = hh->baseextent.stride(); + ivect rlb = hh->baseextent.lower(); + ivect rub = hh->baseextent.upper() + rstr; for (int rl=0; rl oldrlb = rlb; - const vect oldrub = rub; + const ivect oldrlb = rlb; + const ivect oldrub = rub; // calculate extent and centre - const vect rextent = rub - rlb; - const vect rcentre = rlb + (rextent / 2 / rstr) * rstr; + const ivect rextent = rub - rlb; + const ivect rcentre = rlb + (rextent / 2 / rstr) * rstr; // calculate new extent assert (all(rextent % hh->reffact == 0)); - const vect newrextent = rextent / hh->reffact; + const ivect newrextent = rextent / hh->reffact; // refined boxes have smaller stride assert (all(rstr%hh->reffact == 0)); rstr /= hh->reffact; @@ -219,32 +252,32 @@ namespace CarpetRegrid { assert (all(rlb >= oldrlb && rub < oldrub)); } - bbl.push_back (bbox(rlb, rub-rstr, rstr)); + bbl.push_back (ibbox(rlb, rub-rstr, rstr)); } static void MakeRegions_AsSpecified_OneLevel (const cGH* cctkGH, const int reflevels, - const vect ilower, - const vect iupper, - list >& bbl) + const ivect ilower, + const ivect iupper, + list& bbl) { assert (bbl.empty()); if (reflevel+1 >= reflevels) return; - const vect rstr = hh->baseextent.stride(); - const vect rlb = hh->baseextent.lower(); - const vect rub = hh->baseextent.upper(); + const ivect rstr = hh->baseextent.stride(); + const ivect rlb = hh->baseextent.lower(); + const ivect rub = hh->baseextent.upper(); const int rl = reflevel+1; const int levfac = ipow(hh->reffact, rl); assert (all (rstr % levfac == 0)); - const vect str (rstr / levfac); - const vect lb (ilower); - const vect ub (iupper); + const ivect str (rstr / levfac); + const ivect lb (ilower); + const ivect ub (iupper); if (! all(lb>=rlb && ub<=rub)) { CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, "The refinement region boundaries for refinement level #%d are not within the main grid", rl); @@ -261,22 +294,22 @@ namespace CarpetRegrid { assert (all(lb<=ub)); assert (all(lb%str==0 && ub%str==0)); - bbl.push_back (bbox(lb, ub, str)); + bbl.push_back (ibbox(lb, ub, str)); } void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels, - const vector > lower, - const vector > upper, - list >& bbl) + const vector lower, + const vector upper, + list& bbl) { if (reflevel+1 >= reflevels) return; const int rl = reflevel+1; - const vect ilower = lower[rl-1]; - const vect iupper = upper[rl-1]; + const ivect ilower = lower[rl-1]; + const ivect iupper = upper[rl-1]; MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels, ilower, iupper, bbl); } @@ -286,7 +319,7 @@ namespace CarpetRegrid { void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels, const vector > lower, const vector > upper, - list >& bbl) + list& bbl) { if (reflevel+1 >= reflevels) return; @@ -299,12 +332,15 @@ namespace CarpetRegrid { global_upper[d] = 1; } } - const vect global_extent (hh->baseextent.upper() - hh->baseextent.lower() + hh->baseextent.stride() * (dd->lghosts + dd->ughosts)); + const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower() + hh->baseextent.stride() * (dd->lghosts + dd->ughosts)); const int rl = reflevel+1; - const vect ilower = vect(map((CCTK_REAL(*)(CCTK_REAL))(floor), (lower[rl-1] - global_lower) / (global_upper - global_lower) * vect(global_extent) + 0.5)); - const vect iupper = vect(map((CCTK_REAL(*)(CCTK_REAL))(floor), (upper[rl-1] - global_lower) / (global_upper - global_lower) * vect(global_extent) + 0.5)); + const vect scale = vect(global_extent) / (global_upper - global_lower); + const vect rlower = (lower[rl-1] - global_lower) * scale; + const vect rupper = (upper[rl-1] - global_lower) * scale; + const ivect ilower = ivect(map(floor, vect(rlower) + 0.5)); + const ivect iupper = ivect(map(floor, vect(rupper) + 0.5)); MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels, ilower, iupper, bbl); } @@ -312,10 +348,10 @@ namespace CarpetRegrid { static void - MakeRegions_Adaptively_Recombine (list >& bbl1, - list >& bbl2, - list >& bbl, - const bbox& iface, + MakeRegions_Adaptively_Recombine (list& bbl1, + list& bbl2, + list& bbl, + const ibbox& iface, const int dir) { assert (!iface.empty()); @@ -325,23 +361,23 @@ namespace CarpetRegrid { int numcombinedboxes = 0; int oldnumpoints = 0; - for (list >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { + for (list::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { oldnumpoints += ibb->num_points(); } - for (list >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) { + for (list::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) { oldnumpoints += ibb1->num_points(); } - for (list >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) { + for (list::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) { oldnumpoints += ibb2->num_points(); } #if 0 // remember old bounding boxes bboxset oldboxes; - for (list >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) { + for (list::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) { oldboxes += *ibb1; } - for (list >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) { + for (list::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) { oldboxes += *ibb2; } #endif @@ -355,18 +391,18 @@ namespace CarpetRegrid { cout << bbl2 << endl; #endif - const vect lo = iface.lower(); - const vect up = iface.upper(); - const vect str = iface.stride(); + const ivect lo = iface.lower(); + const ivect up = iface.upper(); + const ivect str = iface.stride(); { // prune boxes on the left - list >::iterator ibb1 = bbl1.begin(); + list::iterator ibb1 = bbl1.begin(); while (ibb1 != bbl1.end()) { // is this bbox just to the left of the interface? -// const vect lo1 = ibb1->lower(); - const vect up1 = ibb1->upper(); - const vect str1 = ibb1->stride(); +// const ivect lo1 = ibb1->lower(); + const ivect up1 = ibb1->upper(); + const ivect str1 = ibb1->stride(); assert (up1[dir]+str1[dir] <= lo[dir]); assert (all(str1 == str)); if (up1[dir]+str1[dir] < lo[dir]) { @@ -381,12 +417,12 @@ namespace CarpetRegrid { { // prune boxes on the right - list >::iterator ibb2 = bbl2.begin(); + list::iterator ibb2 = bbl2.begin(); while (ibb2 != bbl2.end()) { // is this bbox just to the right of the interface? - const vect lo2 = ibb2->lower(); -// const vect up2 = ibb2->upper(); - const vect str2 = ibb2->stride(); + const ivect lo2 = ibb2->lower(); +// const ivect up2 = ibb2->upper(); + const ivect str2 = ibb2->stride(); assert (up[dir] <= lo2[dir]); assert (all(str2 == str)); if (up[dir] < lo2[dir]) { @@ -401,31 +437,31 @@ namespace CarpetRegrid { { // walk all boxes on the left - list >::iterator ibb1 = bbl1.begin(); + list::iterator ibb1 = bbl1.begin(); while (ibb1 != bbl1.end()) { - vect lo1 = ibb1->lower(); - vect up1 = ibb1->upper(); - vect str1 = ibb1->stride(); + ivect lo1 = ibb1->lower(); + ivect up1 = ibb1->upper(); + ivect str1 = ibb1->stride(); assert (up1[dir]+str1[dir] == lo[dir]); lo1[dir] = lo[dir]; up1[dir] = up[dir]; - const bbox iface1 (lo1,up1,str1); + const ibbox iface1 (lo1,up1,str1); { // walk all boxes on the right - list >::iterator ibb2 = bbl2.begin(); + list::iterator ibb2 = bbl2.begin(); while (ibb2 != bbl2.end()) { - vect lo2 = ibb2->lower(); - vect up2 = ibb2->upper(); - vect str2 = ibb2->stride(); + ivect lo2 = ibb2->lower(); + ivect up2 = ibb2->upper(); + ivect str2 = ibb2->stride(); assert (lo2[dir] == up[dir]); lo2[dir] = lo[dir]; up2[dir] = up[dir]; - const bbox iface2 (lo2,up2,str2); + const ibbox iface2 (lo2,up2,str2); // check for a match if (iface1 == iface2) { - const bbox combined (ibb1->lower(), ibb2->upper(), str); + const ibbox combined (ibb1->lower(), ibb2->upper(), str); bbl.push_back (combined); ibb1 = bbl1.erase(ibb1); ibb2 = bbl2.erase(ibb2); @@ -452,7 +488,7 @@ namespace CarpetRegrid { assert (newnumboxes + numcombinedboxes == oldnumboxes); int newnumpoints = 0; - for (list >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { + for (list::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { newnumpoints += ibb->num_points(); } assert (newnumpoints == oldnumpoints); @@ -460,7 +496,7 @@ namespace CarpetRegrid { #if 0 // find new bounding boxes bboxset newboxes; - for (list >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { + for (list::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { newboxes += *ibb; } @@ -482,8 +518,8 @@ namespace CarpetRegrid { const CCTK_REAL minfraction, const CCTK_REAL maxerror, const data& error, - list >& bbl, - const bbox& region) + list& bbl, + const ibbox& region) { // Just to be sure assert (! region.empty()); @@ -492,7 +528,7 @@ namespace CarpetRegrid { // (this doesn't work yet on multiple processors) assert (CCTK_nProcs(cctkGH)==1); int cnt = 0; - for (bbox::iterator it=region.begin(); it!=region.end(); ++it) { + for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) { if (error[*it] > maxerror) ++cnt; } const CCTK_REAL fraction = (CCTK_REAL)cnt / region.num_points(); @@ -502,30 +538,30 @@ namespace CarpetRegrid { // Don't refine } else if (width < 2*minwidth || fraction >= minfraction) { // Refine the whole region - const vect lo(region.lower()); - const vect up(region.upper()); - const vect str(region.stride()); - bbl.push_back (bbox(lo,up+str-str/reffact,str/reffact)); + const ivect lo(region.lower()); + const ivect up(region.upper()); + const ivect str(region.stride()); + bbl.push_back (ibbox(lo,up+str-str/reffact,str/reffact)); // cout << "MRA: Refining to " << bbl.back() << " size " << bbl.back().num_points() << " fraction " << fraction << endl; } else { // Split the region and check recursively const int dir = maxloc(region.shape()); - const vect lo(region.lower()); - const vect up(region.upper()); - const vect str(region.stride()); - vect lo1(lo), lo2(lo); - vect up1(up), up2(up); + const ivect lo(region.lower()); + const ivect up(region.upper()); + const ivect str(region.stride()); + ivect lo1(lo), lo2(lo); + ivect up1(up), up2(up); const int mgstr = ipow(hh->mgfact, mglevels); // honour multigrid factors const int step = str[dir]*mgstr; lo2[dir] = ((up[dir]+lo[dir])/2/step)*step; up1[dir] = lo2[dir]-str[dir]; - const bbox region1(lo1,up1,str); - const bbox region2(lo2,up2,str); + const ibbox region1(lo1,up1,str); + const ibbox region2(lo2,up2,str); assert (region1.is_contained_in(region)); assert (region2.is_contained_in(region)); assert ((region1 & region2).empty()); assert (region1 + region2 == region); - list > bbl1, bbl2; + list bbl1, bbl2; MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror, error, bbl1, region1); MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror, @@ -533,7 +569,7 @@ namespace CarpetRegrid { // Combine regions if possible up2 += str-str/reffact; up2[dir] = lo2[dir]; - const bbox iface(lo2,up2,str/reffact); + const ibbox iface(lo2,up2,str/reffact); MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface, dir); } @@ -545,7 +581,7 @@ namespace CarpetRegrid { const CCTK_REAL minfraction, const CCTK_REAL maxerror, const gf& error, - list >& bbl) + list& bbl) { assert (bbl.empty()); @@ -562,7 +598,7 @@ namespace CarpetRegrid { // cout << endl << "MRA: Choosing regions to refine in " << hh->components(rl) << " components" << endl; for (int c=0; ccomponents(rl); ++c) { - const bbox region = hh->extents[rl][c][ml]; + const ibbox region = hh->extents[rl][c][ml]; assert (! region.empty()); const data& errdata = *error(tl,rl,c,ml); @@ -571,10 +607,10 @@ namespace CarpetRegrid { errdata, bbl, region); } - int numpoints = 0; - for (list >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { - numpoints += ibb->num_points(); - } +// int numpoints = 0; +// for (list::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) { +// numpoints += ibb->num_points(); +// } // cout << "MRA: Chose " << bbl.size() << " regions with a total size of " << numpoints << " to refine." << endl << endl; } diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh index cef62a2fb..402e4f8df 100644 --- a/Carpet/CarpetRegrid/src/regrid.hh +++ b/Carpet/CarpetRegrid/src/regrid.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.3 2002/01/11 17:37:14 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.4 2002/03/11 13:17:16 schnetter Exp $ #ifndef REGRID_HH #define REGRID_HH @@ -28,6 +28,7 @@ namespace CarpetRegrid { int CarpetRegridRegrid (const cGH * const cctkGH, gh::rexts& bbsss, + gh::rbnds& obss, gh::rprocs& pss); -- cgit v1.2.3