From 08ebc90094302d03573a613684b1b051d8180316 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 27 Apr 2010 12:02:22 -0500 Subject: CarpetLib: Begin to add new load balancing mechanism Add new class bintree. Add new routine balance. --- Carpet/CarpetLib/src/balance.cc | 495 ++++++++++++++++++++++++++++++++++++ Carpet/CarpetLib/src/balance.hh | 29 +++ Carpet/CarpetLib/src/bintree.cc | 489 +++++++++++++++++++++++++++++++++++ Carpet/CarpetLib/src/bintree.hh | 172 +++++++++++++ Carpet/CarpetLib/src/make.code.defn | 5 +- Carpet/CarpetLib/src/region.cc | 44 ++++ Carpet/CarpetLib/src/region.hh | 4 + 7 files changed, 1237 insertions(+), 1 deletion(-) create mode 100644 Carpet/CarpetLib/src/balance.cc create mode 100644 Carpet/CarpetLib/src/balance.hh create mode 100644 Carpet/CarpetLib/src/bintree.cc create mode 100644 Carpet/CarpetLib/src/bintree.hh (limited to 'Carpet/CarpetLib') diff --git a/Carpet/CarpetLib/src/balance.cc b/Carpet/CarpetLib/src/balance.cc new file mode 100644 index 000000000..e9f5cad6a --- /dev/null +++ b/Carpet/CarpetLib/src/balance.cc @@ -0,0 +1,495 @@ +#include + +#include +#include +#include +#include +#include +#include + +#include "region.hh" + +#include "balance.hh" + +using namespace std; + + + +namespace CarpetLib { + + // Interface for one item + struct item_ifc { + CCTK_REAL load() const; + item_ifc split (CCTK_REAL ratio_new_over_old); + }; + + // Declaration of template functions implementing an equivalen + // interface + template + CCTK_REAL + item_load (item_t const& item); + + template + item_t item_split (item_t& item, CCTK_REAL ratio_new_over_old); + + // Default definitions of these template functions using the + // interface above + template + CCTK_REAL + item_load (item_t const& item) + { + return item.load(); + } + + template + item_t item_split (item_t& item, CCTK_REAL ratio_new_over_old) + { + return item.split (ratio_new_over_old); + } + + + + // A collection of items + template + class items_t { + typedef list coll_t; + coll_t items; + typename coll_t::iterator find_largest_item (); + public: + items_t (); + items_t (vector const& items_); + void add_item (item_t const& item); + bool empty () const; + size_t size () const; + CCTK_REAL load () const; + item_t& get_one_item (); + item_t& get_largest_item (); + item_t get_and_remove_largest_item (); + void copy_out (vector& dst) const; + }; + + // A worker + template + class worker_t: public items_t { + public: + CCTK_REAL strength () const; + CCTK_REAL business () const; + }; + + // A collection of workers + template + class workers_t { + typedef vector > coll_t; + coll_t workers; + typename coll_t::iterator find_least_busy_worker (); + typename coll_t::iterator find_most_busy_worker (); + public: + workers_t (int nworkers); + bool empty () const; + size_t size () const; + CCTK_REAL load () const; + CCTK_REAL strength () const; + CCTK_REAL ideal_business () const; + CCTK_REAL imbalance () const; + worker_t& get_least_busy_worker (); + worker_t& get_most_busy_worker (); + void ensure_same_size (); + void copy_out (vector >& dst) const; + }; + + + + template + items_t::items_t () + { + } + + template + items_t::items_t (vector const& items_) + { + for (typename vector::const_iterator + p = items_.begin(); p != items_.end(); ++p) + { + add_item (*p); + } + } + + template + void + items_t::add_item (item_t const& item) + { + items.push_back (item); + } + + template + bool + items_t::empty () const + { + return items.empty(); + } + + template + size_t + items_t::size () const + { + return items.size(); + } + + template + CCTK_REAL + items_t::load () const + { + CCTK_REAL total_load = 0.0; + for (typename coll_t::const_iterator + p = items.begin(); p != items.end(); ++p) + { + total_load += item_load (*p); + } + return total_load; + } + + template + typename items_t::coll_t::iterator + items_t::find_largest_item () + { + typename coll_t::iterator max_item = items.end(); + CCTK_REAL max_load = -1.0; + for (typename coll_t::iterator p = items.begin(); p != items.end(); ++p) { + if (item_load (*p) > max_load) { + max_item = p; + max_load = item_load (*max_item); + } + } + return max_item; + } + + template + item_t& + items_t::get_one_item () + { + assert (not empty()); + return items.front(); + } + + template + item_t& + items_t::get_largest_item () + { + typename coll_t::iterator const max_item = find_largest_item(); + assert (max_item != items.end()); + return *max_item; + } + + template + item_t + items_t::get_and_remove_largest_item () + { + typename coll_t::iterator const max_item = find_largest_item(); + assert (max_item != items.end()); + item_t const item = *max_item; + items.erase (max_item); + return item; + } + + template + void + items_t::copy_out (vector& dst) const + { + dst.resize (items.size()); + copy (items.begin(), items.end(), dst.begin()); + } + + + + template + CCTK_REAL + worker_t::strength () const + { + return 1.0; // All workers have the same strength + } + + template + CCTK_REAL + worker_t::business () const + { + return this->load() / strength(); + } + + + + template + workers_t::workers_t (int const nworkers) + : workers (nworkers) + { + } + + template + bool + workers_t::empty () const + { + return workers.empty(); + } + + template + size_t + workers_t::size () const + { + return workers.size(); + } + + template + CCTK_REAL + workers_t::load () const + { + CCTK_REAL total_load = 0.0; + for (typename coll_t::const_iterator + w = workers.begin(); w != workers.end(); ++w) + { + total_load += w->load(); + } + return total_load; + } + + template + CCTK_REAL + workers_t::strength () const + { + CCTK_REAL total_strength = 0.0; + for (typename coll_t::const_iterator + w = workers.begin(); w != workers.end(); ++w) + { + total_strength += w->strength(); + } + return total_strength; + } + + template + CCTK_REAL + workers_t::ideal_business () const + { + return load() / strength(); + } + + template + CCTK_REAL + workers_t::imbalance () const + { + assert (not empty()); + CCTK_REAL max_load = 0.0; + CCTK_REAL avg_load = 0.0; + for (typename coll_t::const_iterator + w = workers.begin(); w != workers.end(); ++w) + { + max_load = max (max_load, w->business()); + avg_load += w->business(); + } + avg_load /= size(); + return max_load - avg_load; + } + + template + typename workers_t::coll_t::iterator + workers_t::find_least_busy_worker () + { + typename coll_t::iterator min_worker = workers.end(); + CCTK_REAL min_business = numeric_limits::max(); + for (typename coll_t::iterator + w = workers.begin(); w != workers.end(); ++w) + { + if (w->business() < min_business) { + min_worker = w; + min_business = min_worker->business(); + } + } + return min_worker; + } + + template + typename workers_t::coll_t::iterator + workers_t::find_most_busy_worker () + { + typename coll_t::iterator max_worker = workers.end(); + CCTK_REAL max_business = 0.0; + for (typename coll_t::iterator + w = workers.begin(); w != workers.end(); ++w) + { + if (w->business() > max_business) { + max_worker = w; + max_business = max_worker->business(); + } + } + return max_worker; + } + + template + worker_t& + workers_t::get_least_busy_worker () + { + typename coll_t::iterator const min_worker = find_least_busy_worker(); + assert (min_worker != workers.end()); + return *min_worker; + } + + template + worker_t& + workers_t::get_most_busy_worker () + { + typename coll_t::iterator const max_worker = find_most_busy_worker(); + assert (max_worker != workers.end()); + return *max_worker; + } + + template + void + workers_t::ensure_same_size () + { + if (empty()) return; // nothing to do + + size_t max_items = 0; + typename coll_t::iterator nonempty_worker = workers.end(); + for (typename coll_t::iterator + w = workers.begin(); w != workers.end(); ++w) + { + if (nonempty_worker == workers.end() and not w->empty()) { + nonempty_worker = w; + } + max_items = max (max_items, w->size()); + } + if (max_items == 0) return; // load is already equal + + for (typename coll_t::iterator + w = workers.begin(); w != workers.end(); ++w) + { + // find a worker who has an item + typename coll_t::iterator const worker = w->empty() ? nonempty_worker : w; + while (w->size() < max_items) { + CCTK_REAL const ratio = 0.0; // create empty fill-up items + w->add_item (item_split (w->get_one_item(), ratio)); + } + } + + for (typename coll_t::const_iterator + w = workers.begin(); w != workers.end(); ++w) + { + assert (w->size() == max_items); + } + } + + template + void + workers_t::copy_out (vector >& dst) const + { + dst.resize (workers.size()); + for (size_t w = 0; w < workers.size(); ++w) { + workers.at(w).copy_out (dst.at(w)); + } + } + + + + ////////////////////////////////////////////////////////////////////////////// + + + + template + void + assign_item (items_t& items, workers_t& workers) + { + // Assign the largest item to the least busy worker + item_t const item = items.get_and_remove_largest_item (); + worker_t& worker = workers.get_least_busy_worker (); + worker.add_item (item); + } + + template + void + split_and_distribute (workers_t& workers) + { + // Split the largest item of the most busy worker and give the + // remainder to another worker + worker_t& worker = workers.get_most_busy_worker (); + item_t& item = worker.get_largest_item (); + + // Determine how to split the item + CCTK_REAL const imbalance = worker.business() - workers.ideal_business(); + // Should we even be here? + assert (imbalance > 0.0); + CCTK_REAL const item_business = item_load (item) / worker.strength(); + // This should be the largest item! + assert (item_business > 0.0); + // Determine how much of the item to give away + CCTK_REAL const ratio = imbalance / item_business; + // A ratio of one or more indicates that the item should be given + // away in its entirety -- which would mean that something went + // wrong before + assert (ratio < 1.0); + + // Split the item... + item_t const new_item = item_split (item, ratio); + // ...and give the remainder to someone else + worker_t& new_worker = workers.get_least_busy_worker (); + // This should be someone else! + assert (&new_worker != &worker); + new_worker.add_item (new_item); + } + + + + template + void + balance (vector const& items_, + vector >& split_items, + int const nworkers, + CCTK_REAL const max_imbalance, + bool const ensure_same_size) + { + assert (max_imbalance >= 0.0); + + if (items_.empty()) return; // nothing to do + items_t items (items_); + assert (nworkers > 0); + workers_t workers (nworkers); + + // Assign items + while (not items.empty()) { + assign_item (items, workers); + } + + // TODO: To parallelise this: Group workers into groups, and + // assign work to these groups. Then balance the load + // recursively. + + // Balance items + for (;;) { + // Measure imbalance + CCTK_REAL const imbalance = workers.imbalance(); + // Are we done? + if (imbalance <= max_imbalance) break; + // Should we even be here? + assert (workers.size() > 1); + + // Do something + split_and_distribute (workers); + + // Ensure progress + assert (workers.imbalance() < imbalance); + } + + if (ensure_same_size) { + workers.ensure_same_size (); + } + + workers.copy_out (split_items); + } + + + + template + void + balance (vector const& items_, + vector >& split_items, + int const nworkers, + CCTK_REAL const max_imbalance, + bool const ensure_same_size); + +} // namespace CarpetLib diff --git a/Carpet/CarpetLib/src/balance.hh b/Carpet/CarpetLib/src/balance.hh new file mode 100644 index 000000000..ccf7dad91 --- /dev/null +++ b/Carpet/CarpetLib/src/balance.hh @@ -0,0 +1,29 @@ +#ifndef BALANCE_HH +#define BALANCE_HH + +#include + +#include + +using namespace std; + + + +namespace CarpetLib { + + // This routine splits N items over P workers. It can split the + // items if necessary to ensure a maximum imbalance, and it can + // ensure (artificially) that each worker receives the same number + // of items. + + template + void + balance (vector const& items_, + vector >& split_items_, + int nworkers, + CCTK_REAL max_imbalance, + bool ensure_same_size); + +} // namespace CarpetLib + +#endif // #ifndef BALANCE_HH diff --git a/Carpet/CarpetLib/src/bintree.cc b/Carpet/CarpetLib/src/bintree.cc new file mode 100644 index 000000000..db3b27e1c --- /dev/null +++ b/Carpet/CarpetLib/src/bintree.cc @@ -0,0 +1,489 @@ +#include +#include +#include + +#include "defs.hh" +#include "region.hh" + +#include "bintree.hh" + + + +// Create an empty tree +template +bintree::bintree () + : type (type_empty) +{ + assert (invariant()); + // This is unused + assert (0); +} + + + +// Create a tree branch from a list of bounds and subtrees +template +bintree::bintree (int const dir_, vect const & bounds_, + vect const & subtrees_) + : type (type_branch), dir (dir_), bounds (bounds_), subtrees (subtrees_) +{ + assert (dir>=0 and dir +bintree::bintree (P const & p_) + : type (type_leaf), p (p_) +{ + assert (invariant()); +} + + + +// Create a tree as copy from another tree +template +bintree::bintree (bintree const & t) + : type (t.type) +{ + switch (type) { + case type_empty: + // do nothing + break; + case type_branch: + dir = t.dir; + bounds = t.bounds; + for (int i=0; i<2; ++i) { + subtrees[i] = new bintree (*t.subtrees[i]); + } + break; + case type_leaf: + p = t.p; + break; + default: + assert (0); + } + assert (invariant()); +} + + + +// Copy a tree +template +bintree & +bintree::operator= (bintree const & t) +{ + assert (invariant()); + if (is_branch()) { + for (int i=0; i<2; ++i) { + delete subtrees[i]; + } + } + type = t.type; + switch (type) { + case type_empty: + // do nothing + break; + case type_branch: + dir = t.dir; + bounds = t.bounds; + for (int i=0; i<2; ++i) { + subtrees[i] = new bintree (*t.subtrees[i]); + } + break; + case type_leaf: + p = t.p; + break; + default: + assert (0); + } + assert (invariant()); + return *this; +} + + + +// Delete a tree (and its subtrees) +template +bintree::~bintree () +{ + assert (invariant()); + if (is_branch()) { + for (int i=0; i<2; ++i) { + delete subtrees[i]; + } + } +} + + + +// Compare trees +template +bool +bintree::operator== (bintree 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 (any (bounds != t.bounds)) return false; + if (any (subtrees != t.subtrees)) return false; + break; + case type_leaf: + return p == t.p; + default: + assert (0); + } + return true; +} + + + +// Invariant +template +bool +bintree::invariant () const +{ + return empty() + is_branch() + is_leaf() == 1; +} + + + +// Find the leaf payload corresponding to a position +template +P const * +bintree::search (tvect const & ipos) const +{ + assert (not empty()); + if (is_leaf()) return & p; + int const i = ::asearch (ipos[dir], bounds); + if (i<0 or i>=2) return NULL; // not found + return subtrees[i]->search(ipos); +} + +template +P * +bintree::search (tvect const & ipos) +{ + assert (not empty()); + if (is_leaf()) return & p; + int const i = ::asearch (ipos[dir], bounds); + if (i<0 or i>=2) return NULL; // not found + return subtrees[i]->search(ipos); +} + + + +// Const iterator +template +bintree::const_iterator::const_iterator (bintree const & f_) + : f(f_), i(0), it(0) +{ + if (f.is_branch()) { + it = new const_iterator (* f.subtrees[i]); + while ((*it).done()) { + delete it; + it = 0; + ++ i; + if (done()) break; + // to do: use a new function "reset iterator" instead + it = new const_iterator (* f.subtrees[i]); + } + assert (done() or not (*it).done()); + } +} + +template +bintree::const_iterator::const_iterator (bintree const & f_, int) + : f(f_), it(0) +{ + if (f.empty()) { + i = 0; + } else if (f.is_leaf()) { + i = 1; + } else { + i = 2; + } +} + +template +bintree::const_iterator::~const_iterator () +{ + if (it) { + delete it; + } +} + +template +bintree const & +bintree::const_iterator::operator* () const +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + return f; + } else { + assert (it); + return **it; + } +} + +template +bool +bintree::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 bintree::const_iterator & +bintree::const_iterator::operator++ () +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + ++ i; + } else { + ++ *it; + while ((*it).done()) { + delete it; + it = 0; + ++ i; + if (done()) break; + // to do: use a new function "reset iterator" instead + it = new const_iterator (* f.subtrees[i]); + } + assert (done() or not (*it).done()); + } + return *this; +} + +template +bool +bintree::const_iterator::done () const +{ + if (f.empty()) { + return true; + } else if (f.is_leaf()) { + return i > 0; + } else { + return i == 2; + } +} + + + +// Non-const iterator +template +bintree::iterator::iterator (bintree & f_) + : f(f_), i(0), it(0) +{ + if (f.is_branch()) { + it = new iterator (* f.subtrees[i]); + while ((*it).done()) { + delete it; + it = 0; + ++ i; + if (done()) break; + // to do: use a new function "reset iterator" instead + it = new iterator (* f.subtrees[i]); + } + assert (done() or not (*it).done()); + } +} + +template +bintree::iterator::iterator (bintree & f_, int) + : f(f_), it(0) +{ + if (f.empty()) { + i = 0; + } else if (f.is_leaf()) { + i = 1; + } else { + i = 2; + } +} + +template +bintree::iterator::~iterator () +{ + if (it) { + delete it; + } +} + +template +bintree & +bintree::iterator::operator* () const +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + return f; + } else { + assert (it); + return **it; + } +} + +template +bool +bintree::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 bintree::iterator & +bintree::iterator::operator++ () +{ + assert (not done()); + assert (not f.empty()); + if (f.is_leaf()) { + ++ i; + } else { + ++ *it; + while ((*it).done()) { + delete it; + it = 0; + ++ i; + if (done()) break; + // to do: use a new function "reset iterator" instead + it = new iterator (* f.subtrees[i]); + } + assert (done() or not (*it).done()); + } + return *this; +} + +template +bool +bintree::iterator::done () const +{ + if (f.empty()) { + return true; + } else if (f.is_leaf()) { + return i > 0; + } else { + return i == 2; + } +} + + + +// Memory usage +template +size_t +bintree::memory () const +{ + size_t size = sizeof *this; + if (is_branch()) { + for (int i=0; i<2; ++i) { + size += memoryof(*subtrees[i]); + } + } + return size; +} + + + +// Output helper +template +void +bintree::output (ostream & os) const +{ + os << "bintree{"; + if (empty()) { + os << "empty"; + } else if (is_branch()) { + os << "branch:" + << "dir=" << dir << "," + << "subtrees=["; + for (int i=0; i<2; ++i) { + os << bounds[i] << ":[" << i << "]=" << *subtrees[i] << ":"; + } + os << bounds[2] << "]"; + } else { + os << "leaf:" + << "payload=" << p; + } + os << "}"; +} + + + +// Generic arithmetic search +template +static +int asearch (T const t, vect const & ts) +{ + int imin = 0; + int imax = D - 1; + if (imax <= imin) { + return imin; + } + T tmin = ts[imin]; + if (t < tmin) { + return -1; + } + T tmax = ts[imax]; + if (t >= tmax) { + return imax; + } + int isize = imax - imin; + for (;;) { + // check loop invariants + assert (imin < imax); + assert (t >= tmin); + assert (t < tmax); + if (imax == imin+1) { + 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) { + return imax - 1; + } + assert (iguess > imin and iguess < imax); + T const tguess = ts[iguess]; + if (t < tguess) { + imax = iguess; + tmax = tguess; + } else { + imin = iguess; + tmin = tguess; + } + // check loop variant + int const newisize = imax - imin; + assert (newisize < isize); + isize = newisize; + } +} + + + +template class bintree ; + +template size_t memoryof (bintree const & f); + +template ostream & operator<< (ostream & os, bintree const & f); diff --git a/Carpet/CarpetLib/src/bintree.hh b/Carpet/CarpetLib/src/bintree.hh new file mode 100644 index 000000000..7a7927d9a --- /dev/null +++ b/Carpet/CarpetLib/src/bintree.hh @@ -0,0 +1,172 @@ +#ifndef BINTREE_HH +#define BINTREE_HH + +#include +#include +#include +#include +#include + +#include + +using namespace std; + + + +// This is a binary tree data structure, i.e. a tree data structure +// which decomposes a cuboid domain into two non-overlapping cuboid +// subdomains. 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 +static +int asearch (T t, vect const & ts); + + +// T: index space (usually integer, or maybe real) +// D: number of dimensions (rank) +// P: payload (e.g. region_t) +template +class bintree { + +public: + + // Short name for a small vector + typedef vect 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: + // 3 bounds, splitting the domain + vect bounds; + // 2 pointers to subtrees + vect subtrees; + + // If this is a leaf: + // just the payload + P p; + +public: + + // Create an empty tree + bintree (); + + // Create a tree branch from a list of bounds and subtrees + bintree (int dir_, vect const & bounds_, + vect const & subtrees_); + + // Create a tree leaf from a payload + bintree (P const & p_); + + // Create a tree as copy from another tree + bintree (bintree const & t); + + // Copy a tree + bintree & operator= (bintree const & t); + + // Delete a tree (and its subtrees) + ~bintree (); + + // 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== (bintree const & t) const; + bool operator!= (bintree 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 { + bintree const & f; + int i; + const_iterator * it; + public: + const_iterator (bintree const & f_); + const_iterator (bintree const & f_, int); + ~const_iterator (); + private: + const_iterator (const_iterator const &); + const_iterator & operator= (const_iterator const &); + public: + bintree 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 { + bintree & f; + int i; + iterator * it; + public: + iterator (bintree & f_); + iterator (bintree & f_, int); + ~iterator (); + private: + iterator (iterator const &); + iterator & operator= (iterator const &); + public: + bintree & operator* () const; + bool operator== (iterator const & it2) const; + bool operator!= (iterator const & it2) const + { return not (*this == it2); } + iterator & operator++ (); + bool done () const; + }; + + // Memory usage + size_t memory () const CCTK_ATTRIBUTE_PURE; + + // Output helper + void output (ostream & os) const; +}; + + + +// Memory usage +template +inline size_t memoryof (bintree const & f) CCTK_ATTRIBUTE_PURE; +template +inline size_t memoryof (bintree const & f) { return f.memory(); } + + + +// Output +template +ostream & +operator<< (ostream & os, bintree const & f) +{ + f.output (os); + return os; +} + + + +#endif // #ifndef BINTREE_HH diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index ccf896ae5..96d621172 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -1,8 +1,11 @@ # Main make.code.defn file for thorn CarpetLib -*-Makefile-*- # Source files in this directory -SRCS = bbox.cc \ +SRCS = balance.cc \ + bbox.cc \ bboxset.cc \ + bboxtree.cc \ + bintree.cc \ commstate.cc \ data.cc \ defs.cc \ diff --git a/Carpet/CarpetLib/src/region.cc b/Carpet/CarpetLib/src/region.cc index 00dabc410..eaf04d490 100644 --- a/Carpet/CarpetLib/src/region.cc +++ b/Carpet/CarpetLib/src/region.cc @@ -88,6 +88,50 @@ operator== (region_t const & a, region_t const & b) +// Assign a load to a region +CCTK_REAL +region_t::load () const +{ + return extent.size(); +} + + + +// Split a region into two +region_t +region_t::split (CCTK_REAL const ratio_new_over_old) +{ + assert (ratio_new_over_old >= 0 and ratio_new_over_old <= 1); + if (extent.empty()) { + // Don't do anything for empty regions + return *this; + } + // Choose a direction (prefer the z direction) + int const idir = maxloc1 (extent.shape()); + int const np = extent.shape()[idir]; + // Keep the lower part, and split off the upper part + int const new_np = floor (np * ratio_new_over_old + 0.5); + int const keep_np = np - new_np; + // Calculate new region extents + ivect const lo = extent.lower(); + ivect const up = extent.upper(); + ivect const str = extent.stride(); + ivect const cut = lo + str * ivect::dir(idir) * keep_np; + + // Copy the region + region_t newreg = *this; + // Set new extents + extent = ibbox (lo, cut-str, str); + newreg.extent = ibbox (cut, up, str); + // Mark cutting boundary as not outer boundary + outer_boundaries[idir][1] = false; + newreg.outer_boundaries[idir][0] = false; + + return newreg; +} + + + // 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 diff --git a/Carpet/CarpetLib/src/region.hh b/Carpet/CarpetLib/src/region.hh index ed14b024d..f47c38617 100644 --- a/Carpet/CarpetLib/src/region.hh +++ b/Carpet/CarpetLib/src/region.hh @@ -26,6 +26,10 @@ struct region_t { ~region_t (); bool invariant () const CCTK_ATTRIBUTE_PURE; + + // For regridding + CCTK_REAL load () const; + region_t split (CCTK_REAL ratio_new_over_old); }; -- cgit v1.2.3