aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Carpet/CarpetLib/src/gh.cc396
-rw-r--r--Carpet/CarpetLib/src/gh.hh87
2 files changed, 250 insertions, 233 deletions
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 1fa06ee2d..e85d010d8 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -18,46 +18,163 @@ using namespace std;
// Constructors
-gh::gh (const vector<ivect> & reffacts_, const centering refcent_,
- const int mgfact_, const centering mgcent_,
- const ibbox baseextent_)
+gh::
+gh (vector<ivect> const & reffacts_, centering const refcent_,
+ int const mgfact_, centering const mgcent_,
+ vector<vector<ibbox> > const & baseextents_,
+ i2vect const & boundary_width_)
: reffacts(reffacts_), refcent(refcent_),
mgfact(mgfact_), mgcent(mgcent_),
- baseextent(baseextent_)
+ baseextents(baseextents_),
+ boundary_width(boundary_width_)
{
assert (reffacts.size() >= 1);
assert (all (reffacts.front() == 1));
- for (size_t n = 1; n < reffacts.size(); ++ n) {
- assert (all (reffacts.AT(n) >= reffacts.AT(n-1)));
- assert (all (reffacts.AT(n) % reffacts.AT(n-1) == 0));
+ for (int rl = 1; rl < (int)reffacts.size(); ++ rl) {
+ assert (all (reffacts.AT(rl) >= reffacts.AT(rl-1)));
+ assert (all (reffacts.AT(rl) % reffacts.AT(rl-1) == 0));
+ }
+ assert (refcent == vertex_centered or refcent == cell_centered);
+
+ assert (mgfact >= 1);
+ assert (mgcent == vertex_centered or mgcent == cell_centered);
+
+ assert (baseextents.size() >= 1);
+ assert (baseextents.AT(0).size() >= 1);
+ assert (baseextents.AT(0).size() == reffacts.size());
+ for (int ml = 1; ml < (int)baseextents.size(); ++ ml) {
+ assert (baseextents.AT(ml).size() == baseextents.AT(ml-1).size());
+ }
+ for (int ml = 0; ml < (int)baseextents.size(); ++ ml) {
+ for (int rl = 1; rl < (int)baseextents.AT(ml).size(); ++ rl) {
+ ibbox const & cbox = baseextents.AT(ml).AT(rl-1);
+ ibbox const & fbox = baseextents.AT(ml).AT(rl);
+ assert (all (cbox.stride() * reffacts.AT(rl-1) ==
+ fbox.stride() * reffacts.AT(rl)));
+ }
+ }
+
+ assert (all (all (boundary_width >= 0)));
+ for (int ml = 0; ml < (int)baseextents.size(); ++ ml) {
+ for (int rl = 0; rl < (int)baseextents.AT(ml).size(); ++ rl) {
+ ibbox const & box = baseextents.AT(ml).AT(rl);
+ assert (all (box.shape() / box.stride() >
+ boundary_width[0] + boundary_width[1]));
+ }
}
- assert (all (baseextent.stride() % reffacts.AT(reffacts.size()-1) == 0));
}
// Destructors
-gh::~gh () { }
+gh::
+~gh ()
+{
+}
+
+
// Modifiers
-void gh::regrid (mregs const & regs)
+void
+gh::
+regrid (mregs const & regs)
{
DECLARE_CCTK_PARAMETERS;
- // Save the old grid hierarchy
- _oldregions.clear ();
- swap (_oldregions, _regions);
- _regions = regs;
+ // Save the grid hierarchy
+ oldregions.clear ();
+ swap (oldregions, regions);
+ regions = regs;
+
+
// Consistency checks
- // nota bene: there might be 0 refinement levels.
+ // Note: there may be zero refinement levels
- check_multigrid_consistency ();
- check_component_consistency ();
- check_base_grid_extent ();
- check_refinement_levels ();
+ // Check multigrid consistency
+ assert (mglevels()>0);
+ for (int ml=1; ml<mglevels(); ++ml) {
+ for (int rl=0; rl<reflevels(); ++rl) {
+ for (int c=0; c<components(rl); ++c) {
+ assert (all(extent(ml,rl,c).stride() ==
+ mgfact * extent(ml-1,rl,c).stride()));
+ assert (extent(ml,rl,c)
+ .contracted_for(extent(ml-1,rl,c))
+ .is_contained_in(extent(ml-1,rl,c)));
+ }
+ }
+ }
- calculate_base_extents_of_all_levels ();
+ // Check component consistency
+ for (int ml=0; ml<mglevels(); ++ml) {
+ for (int rl=0; rl<reflevels(); ++rl) {
+ assert (components(rl)>0);
+ for (int c=0; c<components(rl); ++c) {
+ ibbox const & b = extent(ml,rl,c);
+ ibbox const & b0 = extent(ml,rl,0);
+ assert (all(b.stride() == b0.stride()));
+ assert (b.is_aligned_with(b0));
+ for (int cc=c+1; cc<components(rl); ++cc) {
+ assert ((b & extent(ml,rl,cc)).empty());
+ }
+ }
+ }
+ }
+
+ // Check base grid extent
+ for (int ml=0; ml<mglevels(); ++ml) {
+ for (int rl=0; rl<reflevels(); ++rl) {
+ for (int c=0; c<components(rl); ++c) {
+ assert (extent(ml,rl,c).is_contained_in(baseextents.at(ml).at(rl)));
+ }
+ }
+ }
+ // Check proper nesting, i.e., whether the fine grids are contained
+ // in the coarser grids
+ {
+ bool have_error = false;
+ for (int ml=0; ml<mglevels(); ++ml) {
+ for (int rl=1; rl<reflevels(); ++rl) {
+ assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
+ extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
+ // Check contained-ness:
+ // first take all coarse grids
+ ibset coarse;
+ for (int c=0; c<components(rl-1); ++c) {
+ coarse += extent(ml,rl-1,c);
+ }
+ coarse.normalize();
+ // then check all fine grids
+ for (int c=0; c<components(rl); ++c) {
+ ibbox const & fine =
+ extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
+ if (not (fine <= coarse)) {
+ if (not have_error) {
+ cout << "The following components are not properly nested, i.e.," << endl
+ << "they are not contained within the next coarser level's components:" << endl;
+ have_error = true;
+ }
+ cout << " ml " << ml << " rl " << rl << " c " << c << ": "
+ << fine << endl;
+ }
+ } // for c
+ } // for rl
+ } // for ml
+ if (have_error) {
+ cout << "The grid hierarchy is:" << endl;
+ for (int ml=0; ml<mglevels(); ++ml) {
+ for (int rl=0; rl<reflevels(); ++rl) {
+ for (int c=0; c<components(rl); ++c) {
+ cout << " ml " << ml << " rl " << rl << " c " << c << ": "
+ << extent(ml,rl,c) << endl;
+ }
+ }
+ }
+ CCTK_WARN (0, "The refinement hierarchy is not properly nested.");
+ }
+ }
+
+ // Output
if (output_bboxes) {
do_output_bboxes (cout);
do_output_bases (cout);
@@ -73,12 +190,16 @@ void gh::regrid (mregs const & regs)
}
}
-bool gh::recompose (const int rl,
- const bool do_prolongate)
+
+
+bool
+gh::
+recompose (int const rl,
+ bool const do_prolongate)
{
// Handle changes in number of mglevels
- if (_oldregions.size() != _regions.size()) {
- _oldregions.resize (_regions.size());
+ if (oldregions.size() != regions.size()) {
+ oldregions.resize (regions.size());
}
bool const do_recompose = level_did_change(rl);
@@ -92,8 +213,8 @@ bool gh::recompose (const int rl,
// Overwrite old with new grid hierarchy
for (int ml=0; ml<mglevels(); ++ml) {
- _oldregions.AT(ml).resize (_regions.AT(ml).size());
- _oldregions.AT(ml).AT(rl) = _regions.AT(ml).AT(rl);
+ oldregions.AT(ml).resize (regions.AT(ml).size());
+ oldregions.AT(ml).AT(rl) = regions.AT(ml).AT(rl);
}
}
@@ -101,151 +222,42 @@ bool gh::recompose (const int rl,
return do_recompose;
}
-bool gh::level_did_change (const int rl) const
+
+
+bool
+gh::
+level_did_change (int const rl)
+ const
{
// Find out whether this level changed
- if (_regions.size() != _oldregions.size()) return true;
+ if (regions.size() != oldregions.size()) return true;
for (int ml=0; ml<mglevels(); ++ml) {
assert (rl>=0 and rl<reflevels());
- if (rl >= (int)_oldregions.AT(ml).size()) return true;
- if (_regions.AT(ml).AT(rl).size() != _oldregions.AT(ml).AT(rl).size()) {
+ if (rl >= (int)oldregions.AT(ml).size()) return true;
+ if (regions.AT(ml).AT(rl).size() != oldregions.AT(ml).AT(rl).size()) {
return true;
}
for (int c=0; c<components(rl); ++c) {
- if (_regions.AT(ml).AT(rl).AT(c) != _oldregions.AT(ml).AT(rl).AT(c)) {
+ if (regions.AT(ml).AT(rl).AT(c) != oldregions.AT(ml).AT(rl).AT(c)) {
return true;
}
} // for c
} // for ml
return false;
}
-
-void gh::check_multigrid_consistency ()
-{
- assert (mglevels()>0);
- for (int ml=1; ml<mglevels(); ++ml) {
- for (int rl=0; rl<reflevels(); ++rl) {
- for (int c=0; c<components(rl); ++c) {
- assert (all(extent(ml,rl,c).stride()
- == ivect(mgfact) * extent(ml-1,rl,c).stride()));
- // TODO: put the check back in, taking outer boundaries into
- // account
-#if 0
- assert (extent(ml,rl,c)
- .contracted_for(extent(ml-1,rl,c))
- .is_contained_in(extent(ml-1,rl,c)));
-#endif
- }
- }
- }
-}
-
-void gh::check_component_consistency ()
-{
- for (int ml=0; ml<mglevels(); ++ml) {
- for (int rl=0; rl<reflevels(); ++rl) {
- assert (components(rl)>0);
- for (int c=0; c<components(rl); ++c) {
- const ibbox &b = extent(ml,rl,c);
- const ibbox &b0 = extent(ml,rl,0);
- assert (all(b.stride() == b0.stride()));
- assert (b.is_aligned_with(b0));
- for (int cc=c+1; cc<components(rl); ++cc) {
- assert ((b & extent(ml,rl,cc)).empty());
- }
- }
- }
- }
-}
-void gh::check_base_grid_extent ()
-{
- if (reflevels()>0) {
- for (int c=0; c<components(0); ++c) {
- // TODO: put the check back in, taking outer boundaries into
- // account
-#if 0
- assert (extent(0,c,0).is_contained_in(baseextent));
-#endif
- }
- }
-}
-
-// Check proper nesting, i.e., whether the fine grids are contained in
-// the coarser grids
-void gh::check_refinement_levels ()
-{
- bool have_error = false;
- for (int ml=0; ml<mglevels(); ++ml) {
- for (int rl=1; rl<reflevels(); ++rl) {
- assert (all(extent(ml,rl-1,0).stride()
- == ((reffacts.AT(rl) / reffacts.AT(rl-1))
- * extent(ml,rl,0).stride())));
- // Check contained-ness:
- // first take all coarse grids ...
- ibset all;
- for (int c=0; c<components(rl-1); ++c) {
- all |= extent(ml,rl-1,c);
- }
-#if 0
- // ... remember their size ...
- const int sz = all.size();
- // ... then add the coarsified fine grids ...
- for (int c=0; c<components(rl); ++c) {
- all |= extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
- }
- // ... and then check the sizes:
- assert (all.size() == sz);
-#endif
- for (int c=0; c<components(rl); ++c) {
- ibset const finebox
- = (extent(ml,rl,c).contracted_for (extent(ml,rl-1,0)));
- if (! (finebox <= all)) {
- if (! have_error) {
- cout << "The following components are not properly nested, i.e., they are not" << endl
- << "contained within their next coarser level's components:" << endl;
- }
- cout << " ml " << ml << " rl " << rl << " c " << c << endl;
- have_error = true;
- }
- }
- }
- }
- if (have_error) {
- cout << "The grid hierarchy looks like:" << endl;
- for (int ml=0; ml<mglevels(); ++ml) {
- for (int rl=0; rl<reflevels(); ++rl) {
- for (int c=0; c<components(rl); ++c) {
- cout << " ml " << ml << " rl " << rl << " c " << c << ": "
- << extent(ml,rl,c) << endl;
- }
- }
- }
- CCTK_WARN (0, "The refinement hierarchy is not properly nested.");
- }
-}
-void gh::calculate_base_extents_of_all_levels ()
-{
- _bases.resize(mglevels());
- for (int ml=0; ml<mglevels(); ++ml) {
- _bases.AT(ml).resize(reflevels());
- for (int rl=0; rl<reflevels(); ++rl) {
- _bases.AT(ml).AT(rl) = ibbox();
- ibbox &bb = _bases.AT(ml).AT(rl);
- for (int c=0; c<components(rl); ++c) {
- bb = bb.expanded_containing(extent(ml,rl,c));
- }
- }
- }
-}
// Accessors
-int gh::local_components (const int rl) const
+
+int
+gh::
+local_components (int const rl)
+ const
{
int lc = 0;
for (int c=0; c<components(rl); ++c) {
- if (is_local(rl,c)) ++lc;
+ lc += is_local(rl,c);
}
return lc;
}
@@ -253,75 +265,99 @@ int gh::local_components (const int rl) const
// Time hierarchy management
-void gh::add (th* t)
+
+void
+gh::
+add (th * const t)
{
- ths.push_back(t);
+ ths.push_back (t);
}
-void gh::remove (th* t)
+void
+gh::
+remove (th * const t)
{
- ths.remove(t);
+ ths.remove (t);
}
// Data hierarchy management
-void gh::add (dh* d)
+
+void
+gh::
+add (dh * const d)
{
- dhs.push_back(d);
+ dhs.push_back (d);
}
-void gh::remove (dh* d)
+void
+gh::
+remove (dh * const d)
{
- dhs.remove(d);
+ dhs.remove (d);
}
-void gh::do_output_bboxes (ostream& os) const
+
+// Output
+
+void
+gh::
+do_output_bboxes (ostream & os)
+ const
{
for (int ml=0; ml<mglevels(); ++ml) {
for (int rl=0; rl<reflevels(); ++rl) {
for (int c=0; c<components(rl); ++c) {
- os << endl;
- os << "gh bboxes:" << endl;
- os << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
- os << "extent=" << extent(ml,rl,c) << endl;
- os << "outer_boundaries=" << outer_boundaries(rl,c) << endl;
- os << "refinement_boundaries=" << refinement_boundaries(rl,c) << endl;
- os << "processor=" << processor(rl,c) << endl;
+ os << eol;
+ os << "gh bboxes:" << eol;
+ os << "ml=" << ml << " rl=" << rl << " c=" << c << eol;
+ os << "extent=" << extent(ml,rl,c) << eol;
+ os << "outer_boundaries=" << outer_boundaries(rl,c) << eol;
+ os << "processor=" << processor(rl,c) << eol;
}
}
}
}
-void gh::do_output_bases (ostream& os) const
+void
+gh::
+do_output_bases (ostream & os)
+ const
{
for (int ml=0; ml<mglevels(); ++ml) {
for (int rl=0; rl<reflevels(); ++rl) {
if (components(rl)>0) {
- os << endl;
- os << "gh bases:" << endl;
- os << "ml=" << ml << " rl=" << rl << endl;
- os << "base=" << bases().AT(ml).AT(rl) << endl;
+ os << eol;
+ os << "gh bases:" << eol;
+ os << "ml=" << ml << " rl=" << rl << eol;
+ os << "base=" << baseextents.AT(ml).AT(rl) << eol;
}
}
}
}
-ostream& gh::output (ostream& os) const
+ostream &
+gh::
+output (ostream & os)
+ const
{
os << "gh:"
<< "reffacts=" << reffacts << ",refcentering=" << refcent << ","
<< "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
- << "regions=" << regions() << ","
+ << "baseextents=" << baseextents << ","
+ << "boundary_width=" << boundary_width << ","
+ << "regions=" << regions << ","
<< "dhs={";
- const char * sep = "";
- for (list<dh*>::const_iterator d = dhs.begin();
- d != dhs.end(); ++d)
{
- os << sep;
- (*d)->output(os);
- sep = ",";
+ bool isfirst = true;
+ for (list<dh*>::const_iterator
+ d = dhs.begin(); d != dhs.end(); ++ d, isfirst = false)
+ {
+ if (not isfirst) os << ",";
+ os << *d;
+ }
}
os << "}";
return os;
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 619aeebf0..6346e1509 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -22,9 +22,6 @@ class dh;
class th;
class gh;
-// Output
-ostream& operator<< (ostream& os, const gh& h);
-
// A refinement hierarchy, where higher levels are finer than the base
@@ -47,15 +44,11 @@ public: // should be readonly
const int mgfact; // default multigrid factor
const centering mgcent; // default (vertex or cell centered)
- const ibbox baseextent;
-
-private:
- vector<vector<ibbox> > _bases; // [ml][rl]
+ vector<vector<ibbox> > baseextents; // [ml][rl]
+ const i2vect boundary_width;
- // Extents and properties of all grids
- mregs _regions;
- // A copy, used during regridding
- mregs _oldregions;
+ mregs regions; // extents and properties of all grids
+ mregs oldregions; // a copy, used during regridding
list<th*> ths; // list of all time hierarchies
list<dh*> dhs; // list of all data hierarchies
@@ -63,67 +56,60 @@ private:
public:
// Constructors
- gh (const vector<ivect> & reffacts, const centering refcent,
- const int mgfact, const centering mgcent,
- const ibbox baseextent);
+ gh (vector<ivect> const & reffacts, centering refcent,
+ int mgfact, centering mgcent,
+ vector<vector<ibbox> > const & baseextents,
+ i2vect const & boundary_width);
// Destructors
~gh ();
// Modifiers
void regrid (mregs const & regs);
- bool recompose (const int rl,
- const bool do_prolongate);
+ bool recompose (int rl, bool do_prolongate);
+
private:
- bool level_did_change (const int rl) const;
+
+ bool level_did_change (int rl) const;
// Accessors
public:
- mregs const & regions () const
- {
- return _regions;
- }
-
- const vector<vector<ibbox> > & bases() const
- {
- return _bases;
- }
- ibbox extent (const int m, const int rl, const int c) const
+ ibbox const & extent (const int ml, const int rl, const int c) const
{
- return _regions.AT(m).AT(rl).AT(c).extent;
+ return regions.AT(ml).AT(rl).AT(c).extent;
}
- b2vect outer_boundaries (const int rl, const int c) const
+ ibbox const & baseextent (const int ml, const int rl) const
{
- return _regions.AT(0).AT(rl).AT(c).outer_boundaries;
+ return baseextents.AT(ml).AT(rl);
}
-
- b2vect refinement_boundaries (const int rl, const int c) const
+
+ b2vect const & outer_boundaries (const int rl, const int c) const
{
- return _regions.AT(0).AT(rl).AT(c).refinement_boundaries;
+ return regions.AT(0).AT(rl).AT(c).outer_boundaries;
}
-
+
int processor (const int rl, const int c) const
{
- return _regions.AT(0).AT(rl).AT(c).processor;
+ return regions.AT(0).AT(rl).AT(c).processor;
}
int mglevels () const
{
- return (int)_regions.size();
+ return (int)regions.size();
}
int reflevels () const
{
if (mglevels() == 0) return 0;
- return (int)_regions.AT(0).size();
+ return (int)regions.AT(0).size();
}
int components (const int rl) const
{
- return (int)_regions.AT(0).AT(rl).size();
+ return (int)regions.AT(0).AT(rl).size();
}
bool is_local (const int rl, const int c) const
@@ -134,32 +120,27 @@ public:
int local_components (const int rl) const;
// Time hierarchy management
- void add (th* t);
- void remove (th* t);
+ void add (th * t);
+ void remove (th * t);
// Data hierarchy management
- void add (dh* d);
- void remove (dh* d);
+ void add (dh * d);
+ void remove (dh * d);
// Output
- ostream& output (ostream& os) const;
+ ostream & output (ostream & os) const;
private:
- void check_multigrid_consistency ();
- void check_component_consistency ();
- void check_base_grid_extent ();
- void check_refinement_levels ();
- void calculate_base_extents_of_all_levels ();
- void do_output_bboxes (ostream& os) const;
- void do_output_bases (ostream& os) const;
+ void do_output_bboxes (ostream & os) const;
+ void do_output_bases (ostream & os) const;
};
-inline ostream& operator<< (ostream& os, const gh& h) {
- h.output(os);
- return os;
+inline ostream & operator<< (ostream & os, gh const & h)
+{
+ return h.output(os);
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-06 22:07:07 +0000