From ece47455bad491c45b3136362e9239f505de23b9 Mon Sep 17 00:00:00 2001 From: eschnett <> Date: Thu, 1 Mar 2001 12:40:00 +0000 Subject: Initial revision darcs-hash:20010301124010-f6438-fca5ed1e25f84efd816aa0d13fc23b58add7195d.gz --- Carpet/CarpetLib/README | 5 +- Carpet/CarpetLib/interface.ccl | 23 +- Carpet/CarpetLib/param.ccl | 24 +- Carpet/CarpetLib/schedule.ccl | 2 +- Carpet/CarpetLib/src/amr.cc | 62 ++ Carpet/CarpetLib/src/bbox.cc | 150 ++-- Carpet/CarpetLib/src/bbox.hh | 142 ++-- Carpet/CarpetLib/src/bboxset.cc | 182 ++--- Carpet/CarpetLib/src/bboxset.hh | 98 +-- Carpet/CarpetLib/src/data.cc | 1529 +++++++---------------------------- Carpet/CarpetLib/src/data.hh | 160 ++-- Carpet/CarpetLib/src/defs.cc | 184 ++--- Carpet/CarpetLib/src/defs.hh | 199 +---- Carpet/CarpetLib/src/dh.cc | 684 ++++------------ Carpet/CarpetLib/src/dh.hh | 99 +-- Carpet/CarpetLib/src/dist.cc | 114 +-- Carpet/CarpetLib/src/dist.hh | 188 +++-- Carpet/CarpetLib/src/gdata.cc | 446 +--------- Carpet/CarpetLib/src/gdata.hh | 256 +++--- Carpet/CarpetLib/src/gf.cc | 95 ++- Carpet/CarpetLib/src/gf.hh | 51 +- Carpet/CarpetLib/src/ggf.cc | 663 ++++++--------- Carpet/CarpetLib/src/ggf.hh | 224 +++-- Carpet/CarpetLib/src/gh.cc | 231 +++--- Carpet/CarpetLib/src/gh.hh | 102 +-- Carpet/CarpetLib/src/instantiate | 163 +--- Carpet/CarpetLib/src/io.cc | 44 + Carpet/CarpetLib/src/io.hh | 61 ++ Carpet/CarpetLib/src/make.code.defn | 42 +- Carpet/CarpetLib/src/th.cc | 75 +- Carpet/CarpetLib/src/th.hh | 70 +- Carpet/CarpetLib/src/vect.cc | 80 +- Carpet/CarpetLib/src/vect.hh | 443 ++-------- Carpet/CarpetLib/src/wave.cc | 468 +++++++++++ Carpet/CarpetLib/src/wave.hh | 122 +++ 35 files changed, 2664 insertions(+), 4817 deletions(-) create mode 100644 Carpet/CarpetLib/src/amr.cc create mode 100644 Carpet/CarpetLib/src/io.cc create mode 100644 Carpet/CarpetLib/src/io.hh create mode 100644 Carpet/CarpetLib/src/wave.cc create mode 100644 Carpet/CarpetLib/src/wave.hh (limited to 'Carpet/CarpetLib') diff --git a/Carpet/CarpetLib/README b/Carpet/CarpetLib/README index 050a838b0..d71eab1f1 100644 --- a/Carpet/CarpetLib/README +++ b/Carpet/CarpetLib/README @@ -1,8 +1,7 @@ Cactus Code Thorn CarpetLib -Authors : Erik Schnetter -CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/README,v 1.3 2004/01/25 14:57:29 schnetter Exp $ +Authors : ... +CVS info : $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/README,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -------------------------------------------------------------------------- Purpose of the thorn: -This thorn contains the backend library that provides mesh refinement. diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl index 593228e75..74f358587 100644 --- a/Carpet/CarpetLib/interface.ccl +++ b/Carpet/CarpetLib/interface.ccl @@ -1,22 +1,5 @@ # Interface definition for thorn CarpetLib -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/interface.ccl,v 1.5 2004/05/04 22:09:54 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/interface.ccl,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -IMPLEMENTS: CarpetLib - -includes header: defs.hh in defs.hh -includes header: dist.hh in dist.hh - -includes header: bbox.hh in bbox.hh -includes header: bboxset.hh in bboxset.hh -includes header: vect.hh in vect.hh - -includes header: data.hh in data.hh -includes header: gdata.hh in gdata.hh - -includes header: dh.hh in dh.hh -includes header: gf.hh in gf.hh -includes header: ggf.hh in ggf.hh -includes header: gh.hh in gh.hh -includes header: th.hh in th.hh - -includes header: operators.hh in operators.hh +implements: CarpetLib +inherits: FlexIO diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index 1b4e14a87..3f7d08893 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -1,24 +1,2 @@ # Parameter definitions for thorn CarpetLib -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/param.ccl,v 1.8 2004/05/21 18:13:41 schnetter Exp $ - -private: - -BOOLEAN verbose "Print info to the screen" STEERABLE=always -{ -} "no" - -BOOLEAN check_array_accesses "Check all array accesses in Fortran" STEERABLE=always -{ -} "no" - -BOOLEAN barriers "Insert barriers at strategic places for debugging purposes (slows down execution)" STEERABLE=always -{ -} "no" - -BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=always -{ -} "no" - -BOOLEAN save_memory_during_regridding "Save some memory during regridding at the expense of speed" -{ -} "no" +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/param.ccl,v 1.1 2001/03/01 13:40:10 eschnett Exp $ diff --git a/Carpet/CarpetLib/schedule.ccl b/Carpet/CarpetLib/schedule.ccl index ef0ff81a6..f3a51901c 100644 --- a/Carpet/CarpetLib/schedule.ccl +++ b/Carpet/CarpetLib/schedule.ccl @@ -1,2 +1,2 @@ # Schedule definitions for thorn CarpetLib -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/schedule.ccl,v 1.2 2003/09/19 16:06:41 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/schedule.ccl,v 1.1 2001/03/01 13:40:10 eschnett Exp $ diff --git a/Carpet/CarpetLib/src/amr.cc b/Carpet/CarpetLib/src/amr.cc new file mode 100644 index 000000000..71cd9b723 --- /dev/null +++ b/Carpet/CarpetLib/src/amr.cc @@ -0,0 +1,62 @@ +/*************************************************************************** + amr.cc - Wavetoy example driver + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/amr.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include +#include +#include +#include + +#include "dist.hh" +#include "vect.hh" + +#include "wave.hh" + +int main (int argc, char *argv[]) { + + dist::init (argc, argv); + + typedef vect ivect; + + if (argc != 2) { + cerr << "Error: wrong arguments" << endl + << "Synopsis: " << argv[0] << " " + << endl; + exit(1); + } + + const int ref_fact = 2; + const int ref_levs = atoi(argv[1]); + assert (ref_levs>0); + const int mg_fact = 2; + const int mg_levs = 1; + const ivect lb(-24), ub(24), str((int)rint(pow(ref_fact, ref_levs-1))); + const int tstr = (int)rint(pow(ref_fact, ref_levs-1)); + const int nsteps = 16 / tstr; + + wave w(ref_fact, ref_levs, mg_fact, mg_levs, lb, ub, str, tstr); + w.init(); + for (int i=0; i - +/*************************************************************************** + bbox.cc - Bounding boxes + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include -#include #include "defs.hh" #include "vect.hh" -#include "bbox.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) || !defined(BBOX_HH) +# include "bbox.hh" +#endif // Constructors template -bbox::bbox (): _lower(T(1)), _upper(T(0)), _stride(T(1)) { } +bbox::bbox (): _lower(1), _upper(0), _stride(1) { } template bbox::bbox (const bbox& b) @@ -34,36 +50,8 @@ bbox::bbox (const vect& lower, const vect& upper, const vect& stride) : _lower(lower), _upper(upper), _stride(stride) { - assert (all(_stride>T(0))); - assert (all((_upper-_lower)%_stride == T(0))); - if (numeric_limits::is_integer && numeric_limits::is_signed) { - // prevent accidental wrap-around - assert (all(_lower < numeric_limits::max() / 2)); - assert (all(_lower > numeric_limits::min() / 2)); - assert (all(_upper < numeric_limits::max() / 2)); - assert (all(_upper > numeric_limits::min() / 2)); - } -} - -// Accessors -template -T bbox::size () const { - if (empty()) return 0; -// return prod((shape()+stride()-1)/stride()); - const vect sh((shape()+stride()-1)/stride()); - T sz = 1, max = numeric_limits::max(); - for (int d=0; d -bool bbox::contains (const vect& x) const { - return all(x>=lower() && x<=upper()); + assert (all(stride>=1)); + assert (all((upper-lower)%stride==0)); } // Operators @@ -126,23 +114,20 @@ bbox bbox::operator& (const bbox& b) const { // Containment template -bool bbox::is_contained_in (const bbox& b) const { - if (empty()) return true; +bool bbox::contained_in (const bbox& b) const { // no alignment check return all(lower()>=b.lower() && upper()<=b.upper()); } // Alignment check template -bool bbox::is_aligned_with (const bbox& b) const { - return all(stride()==b.stride() && (lower()-b.lower()) % stride() == T(0)); +bool bbox::aligned_with (const bbox& b) const { + return all(stride()==b.stride() && (lower()-b.lower()) % stride() == 0); } // Expand the bbox a little by multiples of the stride template bbox bbox::expand (const vect& lo, const vect& hi) const { - // Allow expansion only into directions where the extent is not negative - assert (all(lower()<=upper() || (lo==T(0) && hi==T(0)))); const vect str = stride(); const vect lb = lower() - lo * str; const vect ub = upper() + hi * str; @@ -152,7 +137,6 @@ bbox bbox::expand (const vect& lo, const vect& hi) const { // Find the smallest b-compatible box around *this template bbox bbox::expanded_for (const bbox& b) const { - if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect str = b.stride(); const vect loff = ((lower() - b.lower()) % str + str) % str; const vect uoff = ((upper() - b.lower()) % str + str) % str; @@ -164,7 +148,6 @@ bbox bbox::expanded_for (const bbox& b) const { // Find the largest b-compatible box inside *this template bbox bbox::contracted_for (const bbox& b) const { - if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect str = b.stride(); const vect loff = ((lower() - b.lower()) % str + str) % str; const vect uoff = ((upper() - b.lower()) % str + str) % str; @@ -173,23 +156,11 @@ bbox bbox::contracted_for (const bbox& b) const { return bbox(lo,up,str); } -// Smallest bbox containing both boxes -template -bbox bbox::expanded_containing (const bbox& b) const { - if (empty()) return b; - if (b.empty()) return *this; - assert (is_aligned_with(b)); - const vect lo = min(lower(), b.lower()); - const vect up = max(upper(), b.upper()); - const vect str = min(stride(), b.stride()); - return bbox(lo,up,str); -} - // Iterators template bbox::iterator::iterator (const bbox& box, const vect& pos) : box(box), pos(pos) { - if (box.empty()) this->pos=box.upper(); + if (box.empty()) this->pos=box.upper()+box.stride(); } template @@ -198,73 +169,42 @@ bool bbox::iterator::operator!= (const iterator& i) const { } template -typename bbox::iterator& bbox::iterator::operator++ () { +bbox::iterator& bbox::iterator::operator++ () { for (int d=0; d -typename bbox::iterator bbox::begin () const { +bbox::iterator bbox::begin () const { return iterator(*this, lower()); } template -typename bbox::iterator bbox::end () const { - return iterator(*this, lower()); -} - - - -// Input -template -void bbox::input (istream& is) { - try { - skipws (is); - consume (is, '('); - is >> _lower; - skipws (is); - consume (is, ':'); - is >> _upper; - skipws (is); - consume (is, ':'); - is >> _stride; - skipws (is); - consume (is, ')'); - } catch (input_error &err) { - cout << "Input error while reading a bbox" << endl; - throw err; - } - if (any(_stride<=T(0))) { - cout << "While reading the bbox " << *this << ":" << endl - << " The stride is not positive." << endl; - throw input_error(); - } - if (any((_upper-_lower)%_stride != T(0))) { - cout << "While reading the bbox " << *this << ":" << endl - << " The stride does not evenly divide the extent." << endl; - throw input_error(); - } - assert (all(_stride>T(0))); - assert (all((_upper-_lower)%_stride == T(0))); +bbox::iterator bbox::end () const { + return iterator(*this, upper()+stride()); } // Output template -void bbox::output (ostream& os) const { - os << "(" << lower() << ":" << upper() << ":" << stride() << ")"; +ostream& operator<< (ostream& os, const bbox& b) { + os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")"; + return os; } -// Note: We need all dimensions all the time. -template class bbox; +#if defined(TMPL_EXPLICIT) template class bbox; +template ostream& operator<< (ostream& os, const bbox& b); template class bbox; +template ostream& operator<< (ostream& os, const bbox& b); template class bbox; -template class bbox; +template ostream& operator<< (ostream& os, const bbox& b); +#endif diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index e23d000ac..3529f944e 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.hh @@ -1,4 +1,22 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.15 2004/04/18 13:03:44 schnetter Exp $ +/*************************************************************************** + bbox.hh - Bounding boxes + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef BBOX_HH #define BBOX_HH @@ -8,77 +26,52 @@ #include "defs.hh" #include "vect.hh" -using namespace std; - -// Forward declaration +// Forward definition template class bbox; -// Input/Output -template -istream& operator>> (istream& is, bbox& b); -template +// Output +template ostream& operator<< (ostream& os, const bbox& b); -/** - * A bounding box, i.e. a rectangle with lower and upper bound and a - * stride. - */ +// Bounding box class template class bbox { // Fields - - /** Bounding box bounds and stride. The bounds are inclusive. */ - vect _lower, _upper, _stride; + vect _lower, _upper, _stride;// bounds are inclusive public: // Constructors - - /** Construct an empty bbox. */ bbox (); - - /** Copy constructor. */ bbox (const bbox& b); - - /** Assignment operator. */ bbox& operator= (const bbox& b); - - /** Create a bbox from bounds and stride. */ bbox (const vect& lower, const vect& upper, const vect& stride); // Accessors - // (Don't return references; *this might be a temporary) - - /** Get lower bound. */ - vect lower () const { return _lower; } - - /** Get upper bound. */ - vect upper () const { return _upper; } - - /** Get stride. */ - vect stride () const { return _stride; } - - /** Get the shape (or extent). */ + const vect& lower () const { return _lower; } + const vect& upper () const { return _upper; } + const vect& stride () const { return _stride; } vect shape () const { return _upper - _lower + _stride; } - /** Determine whether the bbox is empty. */ bool empty() const { return any(lower()>upper()); } - /** Return the size, which is the product of the shape. */ - T size () const; - - // Queries + T size () const { + if (empty()) return 0; + return prod(shape()); + } - /** Find out whether the bbox contains the point x. */ - bool contains (const vect& x) const; + T num_points () const { + if (empty()) return 0; + return prod((shape()+stride()-1)/stride()); + } // Operators bool operator== (const bbox& b) const; @@ -88,84 +81,47 @@ public: bool operator<= (const bbox& b) const; bool operator>= (const bbox& b) const; - /** Calculate the intersection (the set of common points) with the - bbox b. */ + // Intersection bbox operator& (const bbox& b) const; - /** Find out whether this bbox is contained in the bbox b. */ - bool is_contained_in (const bbox& b) const; + // Containment + bool contained_in (const bbox& b) const; - /** Find out whether this bbox is aligned with the bbox b. - ("aligned" means that both bboxes have the same stride and that - their boundaries are commesurate.) */ - bool is_aligned_with (const bbox& b) const; + // Alignment check + bool aligned_with (const bbox& b) const; - /** Expand (enlarge) the bbox by multiples of the stride. */ + // Expand the bbox a little by multiples of the stride bbox expand (const vect& lo, const vect& hi) const; - /** Find the smallest b-compatible box around this bbox. - ("compatible" means having the same stride.) */ + // Find the smallest b-compatible box around *this bbox expanded_for (const bbox& b) const; - /** Find the largest b-compatible box inside this bbox. */ + // Find the largest b-compatible box inside *this bbox contracted_for (const bbox& b) const; - /** Find the smallest bbox containing both boxes. */ - bbox expanded_containing (const bbox& b) const; - // Iterators - - /** An iterator over all points in a bbox. */ class iterator { protected: - /** The bbox over which we iterate. */ - const bbox& box; - /** Current position. */ + bbox box; vect pos; public: - /** Constructor. */ iterator (const bbox& box, const vect& pos); - /** Accessor. */ const vect& operator* () const { return pos; } - /** Check whether the position is the same. */ bool operator!= (const iterator& i) const; - /** Advance. */ iterator& operator++ (); }; - /** Create an iterator that points to the first point in a bbox. */ iterator begin () const; - /** Create an iterator that points "after the last point" in a bbox, - which means that it also points to the first point. */ iterator end () const; - // Input/Output helpers - void input (istream& is); - void output (ostream& os) const; + // Output + friend ostream& operator<< <>(ostream& os, const bbox& b); }; -// Input - -/** Read a formatted bbox from a stream. */ -template -inline istream& operator>> (istream& is, bbox& b) { - b.input(is); - return is; -} - - - -// Output - -/** Write a bbox formatted to a stream. */ -template -inline ostream& operator<< (ostream& os, const bbox& b) { - b.output(os); - return os; -} - - +#if defined(TMPL_IMPLICIT) +# include "bbox.cc" +#endif #endif // BBOX_HH diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 94cf1997e..1dd685178 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -1,17 +1,32 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.16 2004/06/13 22:46:48 schnetter Exp $ - -#include - +/*************************************************************************** + bboxset.cc - Sets of bounding boxes + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include -#include #include -#include #include "defs.hh" -#include "bboxset.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) && !defined(BBOXSET_HH) +# include "bboxset.hh" +#endif @@ -23,7 +38,7 @@ bboxset::bboxset () { template bboxset::bboxset (const box& b) { - if (!b.empty()) bs.insert(b); + bs.insert(b); assert (invariant()); } @@ -44,9 +59,11 @@ template bool bboxset::invariant () const { for (const_iterator bi=begin(); bi!=end(); ++bi) { if ((*bi).empty()) return false; - if (! (*bi).is_aligned_with(*bs.begin())) return false; + if (! (*bi).aligned_with(*bs.begin())) return false; // check for overlap (quadratic -- expensive) - for (const_iterator bi2=begin(); bi2!=bi; ++bi2) { + int cnt=0; + for (const_iterator bi2=bi; bi2!=end(); ++bi2) { + if (!cnt++) continue; if (! ((*bi2) & (*bi)).empty()) return false; } } @@ -58,67 +75,7 @@ bool bboxset::invariant () const { // Normalisation template void bboxset::normalize () { - assert (invariant()); - const int num_initial_boxes = bs.size(); - int num_combined_boxes = 0; - stack todo, done; - for (typename set::const_iterator elt = bs.begin(); elt != bs.end(); ++elt) { - done.push (*elt); - } - // TODO: This will not catch all cases where bboxes can be combined. - for (int d=0; d(); - while (! todo.empty()) { - restart:; - box item = todo.top(); - todo.pop(); - stack work = done; - done = stack(); - while (! work.empty()) { - box comp = work.top(); - work.pop(); - { - assert (all(comp.stride() == item.stride())); - if (comp.upper()[d] + item.stride()[d] == item.lower()[d]) { - if (all((comp.lower() == item.lower() - && comp.upper() == item.upper()).replace (d, true))) { - box newbox = box(comp.lower(), item.upper(), item.stride()); - todo.push (newbox); - while (! work.empty()) { - done.push (work.top()); - work.pop(); - } - ++num_combined_boxes; - goto restart; - } - } - if (item.upper()[d] + item.stride()[d] == comp.lower()[d]) { - if (all((comp.lower() == item.lower() - && comp.upper() == item.upper()).replace (d, true))) { - box newbox = box(item.lower(), comp.upper(), item.stride()); - todo.push (newbox); - while (! work.empty()) { - done.push (work.top()); - work.pop(); - } - ++num_combined_boxes; - goto restart; - } - } - } - done.push (comp); - } // while work - done.push (item); - } // while todo - } // for d - bs.clear(); - while (! done.empty()) { - bs.insert (done.top()); - done.pop(); - } - const int num_final_boxes = bs.size(); - assert (num_initial_boxes - num_combined_boxes == num_final_boxes); + // TODO assert (invariant()); } @@ -129,9 +86,7 @@ template T bboxset::size () const { T s=0; for (const_iterator bi=begin(); bi!=end(); ++bi) { - const T bs = (*bi).size(); - assert (numeric_limits::max() - bs >= s); - s += bs; + s += (*bi).size(); } return s; } @@ -176,16 +131,6 @@ bboxset bboxset::operator+ (const bboxset& s) const { return r; } -template -bboxset bboxset::plus (const bbox& b1, const bbox& b2) { - return bboxset(b1) + b2; -} - -template -bboxset bboxset::plus (const bbox& b, const bboxset& s) { - return s + b; -} - // Union @@ -266,8 +211,8 @@ bboxset& bboxset::operator&= (const bboxset& s) { // Difference template -bboxset bboxset::minus (const bbox& b1, const bbox& b2) { - assert (b1.is_aligned_with(b2)); +bboxset operator- (const bbox& b1, const bbox& b2) { + assert (b1.aligned_with(b2)); if (b1.empty()) return bboxset(); if (b2.empty()) return bboxset(b1); const vect str = b1.stride(); @@ -337,7 +282,7 @@ bboxset bboxset::operator- (const bboxset& s) const { } template -bboxset bboxset::minus (const bbox& b, const bboxset& s) { +bboxset operator- (const bbox& b, const bboxset& s) { bboxset r = bboxset(b) - s; assert (r.invariant()); return r; @@ -345,49 +290,30 @@ bboxset bboxset::minus (const bbox& b, const bboxset& s) { -// Equality -template -bool bboxset::operator<= (const bboxset& s) const { - return (*this - s).empty(); -} - -template -bool bboxset::operator< (const bboxset& s) const { - return (*this - s).empty() && ! (s - *this).empty(); -} - -template -bool bboxset::operator>= (const bboxset& s) const { - return s <= *this; -} - -template -bool bboxset::operator> (const bboxset& s) const { - return s < *this; -} - -template -bool bboxset::operator== (const bboxset& s) const { - return (*this <= s) && (*this >= s); -} - -template -bool bboxset::operator!= (const bboxset& s) const { - return ! (*this == s); -} - - - // Output template -void bboxset::output (ostream& os) const { - T Tdummy; - os << "bboxset<" << typestring(Tdummy) << "," << D << ">:" - << "size=" << size() << "," - << "setsize=" << setsize() << "," - << "set=" << bs; +ostream& operator<< (ostream& os, const bboxset& s) { +// os << "bboxset<" STR(T) "," << D << ">:size=" << s.size() << "," +// << "set=" << s.bs; + os << s.bs; + return os; } +#if defined(TMPL_EXPLICIT) +template class bboxset; +template bboxset operator- (const bbox& b1, const bbox& b2); +template bboxset operator- (const bbox& b, const bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& b); + +template class bboxset; +template bboxset operator- (const bbox& b1, const bbox& b2); +template bboxset operator- (const bbox& b, const bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& b); + template class bboxset; +template bboxset operator- (const bbox& b1, const bbox& b3); +template bboxset operator- (const bbox& b, const bboxset& s); +template ostream& operator<< (ostream& os, const bboxset& b); +#endif diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index a94c8940a..009d9afdd 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.hh @@ -1,10 +1,27 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.11 2003/09/19 16:06:41 schnetter Exp $ +/*************************************************************************** + bboxset.hh - Sets of bounding boxes + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef BBOXSET_HH #define BBOXSET_HH -#include - +#include #include #include @@ -12,22 +29,15 @@ #include "defs.hh" #include "vect.hh" -using namespace std; - -// Forward declaration +// Forward definition template class bboxset; -// template -// bboxset operator+ (const bbox& b1, const bbox& b2); -// template -// bboxset operator+ (const bbox& b, const bboxset& s); - -// template -// bboxset operator- (const bbox& b1, const bbox& b2); -// template -// bboxset operator- (const bbox& b, const bboxset& s); +template +bboxset operator- (const bbox& b1, const bbox& b2); +template +bboxset operator- (const bbox& b, const bboxset& s); // Output template @@ -67,15 +77,12 @@ public: // Accessors bool empty () const { return bs.empty(); } T size () const; - int setsize () const { return bs.size(); } // Add (bboxes that don't overlap) bboxset& operator+= (const box& b); bboxset& operator+= (const bboxset& s); bboxset operator+ (const box& b) const; bboxset operator+ (const bboxset& s) const; - static bboxset plus (const box& b1, const box& b2); - static bboxset plus (const box& b, const bboxset& s); // Union bboxset& operator|= (const box& b); @@ -91,66 +98,27 @@ public: // Difference // friend bboxset operator- (const box& b1, const box& b2); - static bboxset minus (const box& b1, const box& b2); bboxset operator- (const box& b) const; bboxset& operator-= (const box& b); bboxset& operator-= (const bboxset& s); bboxset operator- (const bboxset& s) const; // friend bboxset operator- (const box& b, const bboxset& s); - static bboxset minus (const box& b, const bboxset& s); - - // Equality - bool operator== (const bboxset& s) const; - bool operator!= (const bboxset& s) const; - bool operator< (const bboxset& s) const; - bool operator<= (const bboxset& s) const; - bool operator> (const bboxset& s) const; - bool operator>= (const bboxset& s) const; // Iterators - typedef typename bset::const_iterator const_iterator; - typedef typename bset::iterator iterator; + typedef bset::const_iterator const_iterator; + typedef bset::iterator iterator; - const_iterator begin () const { return bs.begin(); } - const_iterator end () const { return bs.end(); } -// iterator begin () const { return bs.begin(); } -// iterator end () const { return bs.end(); } + iterator begin () const { return bs.begin(); } + iterator end () const { return bs.end(); } // Output - void output (ostream& os) const; + friend ostream& operator<< <>(ostream& os, const bboxset& s); }; -template -inline bboxset operator+ (const bbox& b1, const bbox& b2) { - return bboxset::plus(b1,b2); -} - -template -inline bboxset operator+ (const bbox& b, const bboxset& s) { - return bboxset::plus(b,s); -} - -template -inline bboxset operator- (const bbox& b1, const bbox& b2) { - return bboxset::minus(b1,b2); -} - -template -inline bboxset operator- (const bbox& b, const bboxset& s) { - return bboxset::minus(b,s); -} - - - -// Output -template -inline ostream& operator<< (ostream& os, const bboxset& s) { - s.output(os); - return os; -} - - +#if defined(TMPL_IMPLICIT) +# include "bboxset.cc" +#endif #endif // BBOXSET_HH diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 5c9d219e7..18fb2ac98 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -1,499 +1,145 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.55 2004/05/21 18:13:41 schnetter Exp $ - -#include -#include -#include -#include - -#include -#include -#include -#include +/*************************************************************************** + data.cc - Data storage + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include +#include #include -#include #include -#include "cctk.h" - #include "bbox.hh" #include "defs.hh" #include "dist.hh" #include "vect.hh" -#include "data.hh" - -using namespace std; - - - -static size_t total_allocated_bytes; // total number of allocated bytes +#if !defined(TMPL_IMPLICIT) || !defined(DATA_HH) +# include "data.hh" +#endif // Constructors template -data::data (const int varindex_, const operator_type transport_operator_, - const int vectorlength, const int vectorindex, - data* const vectorleader) - : gdata(varindex_, transport_operator_), - _storage(NULL), _allocated_bytes(0), - vectorlength(vectorlength), vectorindex(vectorindex), - vectorleader(vectorleader) -{ - assert (vectorlength>=1); - assert (vectorindex>=0 && vectorindexregister_client (vectorindex); -} +data::data () + : _storage(0) +{ } template -data::data (const int varindex_, const operator_type transport_operator_, - const int vectorlength, const int vectorindex, - data* const vectorleader, - const ibbox& extent_, const int proc_) - : gdata(varindex_, transport_operator_), - _storage(NULL), _allocated_bytes(0), - vectorlength(vectorlength), vectorindex(vectorindex), - vectorleader(vectorleader) +data::data (const ibbox& extent, const int proc) + : _storage(0) { - assert (vectorlength>=1); - assert (vectorindex>=0 && vectorindexregister_client (vectorindex); - allocate(extent_, proc_); + allocate(extent, proc); } // Destructors template -data::~data () -{ - if (vectorleader) vectorleader->unregister_client (vectorindex); - if (vectorindex==0) assert (! has_clients()); +data::~data () { free(); } // Pseudo constructors template -data* data::make_typed (const int varindex_, - const operator_type transport_operator_) - const -{ - return new data(varindex_, transport_operator_); -} - - - -// Vector mamagement -template -void data::register_client (const int index) -{ - assert (! vectorclients.at(index)); - vectorclients.at(index) = true; +data* data::make_typed (const ibbox& extent, const int proc) const { + return new data(extent, proc); } -template -void data::unregister_client (const int index) -{ - assert (vectorclients.at(index)); - vectorclients.at(index) = false; -} - -template -bool data::has_clients () -{ - bool retval = false; - for (size_t n=0; n -void data::getmem (const size_t nelems) -{ - const size_t nbytes = nelems * sizeof(T); - try { - assert (this->_allocated_bytes == 0); - this->_storage = new T[nelems]; - this->_allocated_bytes = nbytes; - } catch (...) { - T Tdummy; - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Failed to allocate %.0f bytes (%.3f MB) of memory for type %s. %.0f bytes (%.3f MB) are currently allocated.", - (double)nbytes, nbytes/1.0e6, - typestring(Tdummy), - (double)total_allocated_bytes, total_allocated_bytes/1.0e6); - } - total_allocated_bytes += nbytes; -} - - - -template -void data::freemem () -{ - delete [] _storage; - total_allocated_bytes -= this->_allocated_bytes; - this->_allocated_bytes = 0; -} - - - -template -void data::allocate (const ibbox& extent_, - const int proc_, - void* const mem) -{ - assert (!this->_has_storage); - this->_has_storage = true; - // prevent accidental wrap-around - assert (all(extent_.lower() < numeric_limits::max() / 2)); - assert (all(extent_.lower() > numeric_limits::min() / 2)); - assert (all(extent_.upper() < numeric_limits::max() / 2)); - assert (all(extent_.upper() > numeric_limits::min() / 2)); +void data::allocate (const ibbox& extent, const int proc) { + assert (!_has_storage); + _has_storage = true; // data - this->_extent = extent_; - this->_shape = max(ivect(0), this->_extent.shape() / this->_extent.stride()); - this->_size = 1; + _extent = extent; + _shape = _extent.shape() / _extent.stride(); + _size = 1; for (int d=0; d_stride[d] = this->_size; - assert (this->_shape[d]==0 || this->_size <= INT_MAX / this->_shape[d]); - this->_size *= this->_shape[d]; + _stride[d] = _size; + _size *= _shape[d]; } - this->_proc = proc_; + _proc = proc; int rank; MPI_Comm_rank (dist::comm, &rank); - if (rank==this->_proc) { - this->_owns_storage = !mem; - if (this->_owns_storage) { - if (this->vectorindex == 0) { - assert (! this->vectorleader); - getmem (this->vectorlength * this->_size); - } else { - assert (this->vectorleader); - this->_storage = this->vectorleader->vectordata (this->vectorindex); - } - } else { - this->_storage = (T*)mem; - } - } else { - assert (!mem); + if (rank==_proc) { + _storage = new T[_size]; } } template -void data::free () -{ - if (this->_storage && this->_owns_storage && this->vectorindex==0) { - freemem (); - } +void data::free () { + if (_storage) delete [] _storage; _storage = 0; - this->_has_storage = false; + _has_storage = false; } template -void data::transfer_from (gdata* gsrc) -{ - assert (this->vectorlength==1); +void data::transfer_from (generic_data* gsrc) { data* src = (data*)gsrc; - assert (src->vectorlength==1); assert (!_storage); *this = *src; - *src = data(this->varindex, this->transport_operator); + *src = data(); } -template -T* data::vectordata (const int vectorindex) const -{ - assert (this->vectorindex==0); - assert (! this->vectorleader); - assert (vectorindex>=0 && vectorindexvectorlength); - assert (this->_storage && this->_owns_storage); - return this->_storage + vectorindex * this->_size; -} - - - // Processor management -template -void data::change_processor (comm_state& state, - const int newproc, - void* const mem) -{ - switch (state.thestate) { - case state_recv: - change_processor_recv (newproc, mem); - break; - case state_send: - change_processor_send (newproc, mem); - break; - case state_wait: - change_processor_wait (newproc, mem); - break; - default: - assert(0); - } -} - - template -void data::change_processor_recv (const int newproc, void* const mem) -{ - assert (!this->comm_active); - this->comm_active = true; - - if (newproc == this->_proc) { - assert (!mem); - return; - } - - if (this->_has_storage) { +void data::change_processor (const int newproc) { + if (newproc == _proc) return; + + if (_has_storage) { int rank; MPI_Comm_rank (dist::comm, &rank); if (rank == newproc) { // copy from other processor assert (!_storage); - this->_owns_storage = !mem; - if (this->_owns_storage) { - getmem (this->_size); - } else { - _storage = (T*)mem; - } - - const double wtime1 = MPI_Wtime(); - T dummy; - MPI_Irecv (_storage, this->_size, dist::datatype(dummy), this->_proc, - this->tag, dist::comm, &this->request); - const double wtime2 = MPI_Wtime(); - this->wtime_irecv += wtime2 - wtime1; - - } else if (rank == this->_proc) { - // copy to other processor - - } else { - assert (!mem); - assert (!_storage); - } - } -} - + _storage = new T[_size]; + T dummy; + MPI_Status status; + MPI_Recv (_storage, _size, dist::datatype(dummy), _proc, + dist::tag, dist::comm, &status); -template -void data::change_processor_send (const int newproc, void* const mem) -{ - assert (this->comm_active); - - if (newproc == this->_proc) { - assert (!mem); - return; - } - - if (this->_has_storage) { - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == newproc) { - // copy from other processor - - } else if (rank == this->_proc) { + } else if (rank == _proc) { // copy to other processor - - assert (!mem); + assert (_storage); - - const double wtime1 = MPI_Wtime(); T dummy; - MPI_Isend (_storage, this->_size, dist::datatype(dummy), newproc, - this->tag, dist::comm, &this->request); - const double wtime2 = MPI_Wtime(); - this->wtime_isend += wtime2 - wtime1; - - } else { - assert (!mem); - assert (!_storage); - } - } -} - + MPI_Send (_storage, _size, dist::datatype(dummy), newproc, + dist::tag, dist::comm); - -template -void data::change_processor_wait (const int newproc, void* const mem) -{ - assert (this->comm_active); - this->comm_active = false; - - if (newproc == this->_proc) { - assert (!mem); - return; - } - - if (this->_has_storage) { - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == newproc) { - // copy from other processor - - const double wtime1 = MPI_Wtime(); - MPI_Status status; - MPI_Wait (&this->request, &status); - const double wtime2 = MPI_Wtime(); - this->wtime_irecvwait += wtime2 - wtime1; - - } else if (rank == this->_proc) { - // copy to other processor - - assert (!mem); - assert (_storage); - - const double wtime1 = MPI_Wtime(); - MPI_Status status; - MPI_Wait (&this->request, &status); - const double wtime2 = MPI_Wtime(); - this->wtime_isendwait += wtime2 - wtime1; - - if (this->_owns_storage) { - freemem (); - } + delete [] _storage; _storage = 0; - + } else { - assert (!mem); assert (!_storage); } } - - this->_proc = newproc; -} - - -// Data manipulators -template -void data -::copy_from_innerloop (const gdata* gsrc, const ibbox& box) -{ - const data* src = (const data*)gsrc; - assert (this->has_storage() && src->has_storage()); - assert (all(box.lower()>=this->extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=this->extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==this->extent().stride() - && box.stride()==src->extent().stride())); - assert (all((box.lower()-this->extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - assert (this->proc() == src->proc()); - - const int groupindex = CCTK_GroupIndexFromVarI(this->varindex); - const int group_tags_table = CCTK_GroupTagsTableI(groupindex); - assert (group_tags_table >= 0); - - // Disallow this. - T Tdummy; - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "There is no copy operator available for the variable type %s, dimension %d.", - typestring(Tdummy), D); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == this->proc()); - - for (typename ibbox::iterator it=box.begin(); it!=box.end(); ++it) { - const ivect index = *it; - (*this)[index] = (*src)[index]; - } - + _proc = newproc; } - - +// Data manipulators template -void data -::interpolate_from_innerloop (const vector*> gsrcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (this->has_storage()); - assert (all(box.lower()>=this->extent().lower())); - assert (all(box.upper()<=this->extent().upper())); - assert (all(box.stride()==this->extent().stride())); - assert (all((box.lower()-this->extent().lower())%box.stride() == 0)); - vector srcs(gsrcs.size()); - for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t]; - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs[t]->has_storage()); - assert (all(box.lower()>=srcs[t]->extent().lower())); - assert (all(box.upper()<=srcs[t]->extent().upper())); - assert (this->proc() == srcs[t]->proc()); - } - assert (order_space >= 0); - assert (order_time >= 0); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == this->proc()); - - assert (this->varindex >= 0); - const int groupindex = CCTK_GroupIndexFromVarI (this->varindex); - assert (groupindex >= 0); - char* groupname = CCTK_GroupName(groupindex); - T Tdummy; - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "There is no interpolator available for the group \"%s\" with variable type %s, dimension %d, spatial interpolation order %d, temporal interpolation order %d.", - groupname, typestring(Tdummy), D, order_space, order_time); - ::free (groupname); -} - - - -extern "C" { - void CCTK_FCALL CCTK_FNAME(copy_3d_int4) - (const CCTK_INT4* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_INT4* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(copy_3d_real8) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(copy_3d_complex16) - (const CCTK_COMPLEX16* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_COMPLEX16* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); -} - -template<> -void data -::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box) -{ +void data::copy_from (const generic_data* gsrc, const ibbox& box) { const data* src = (const data*)gsrc; - assert (has_storage() && src->has_storage()); + assert (_has_storage && src->_has_storage); assert (all(box.lower()>=extent().lower() && box.lower()>=src->extent().lower())); assert (all(box.upper()<=extent().upper() @@ -502,887 +148,278 @@ void data && box.stride()==src->extent().stride())); assert (all((box.lower()-extent().lower())%box.stride() == 0 && (box.lower()-src->extent().lower())%box.stride() == 0)); - - assert (proc() == src->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const ibbox& sext = src->extent(); - const ibbox& dext = extent(); - - int srcshp[3], dstshp[3]; - int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; - - for (int d=0; d<3; ++d) { - srcshp[d] = (sext.shape() / sext.stride())[d]; - dstshp[d] = (dext.shape() / dext.stride())[d]; - - srcbbox[0][d] = sext.lower()[d]; - srcbbox[1][d] = sext.upper()[d]; - srcbbox[2][d] = sext.stride()[d]; + + if (_proc == src->_proc) { + // copy on same processor + + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == _proc) { + + for (ibbox::iterator it=box.begin(); it!=box.end(); ++it) { + const ivect index = *it; + (*this)[index] = (*src)[index]; + } + + } - dstbbox[0][d] = dext.lower()[d]; - dstbbox[1][d] = dext.upper()[d]; - dstbbox[2][d] = dext.stride()[d]; + } else { - regbbox[0][d] = box.lower()[d]; - regbbox[1][d] = box.upper()[d]; - regbbox[2][d] = box.stride()[d]; - } - - assert (all(dext.stride() == box.stride())); - if (all(sext.stride() == dext.stride())) { - CCTK_FNAME(copy_3d_int4) ((const CCTK_INT4*)src->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_INT4*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); + // copy to different processor + data* tmp = new data(box, src->_proc); + tmp->copy_from (src, box); + tmp->change_processor (_proc); + copy_from (tmp, box); + delete tmp; - } else { - assert (0); } } -template<> -void data -::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box) +template +void data::interpolate_from (const generic_data* gsrc, + const ibbox& box) { const data* src = (const data*)gsrc; - assert (has_storage() && src->has_storage()); + assert (_has_storage && src->_has_storage); + assert (all(box.lower()>=extent().lower() + && box.upper()<=extent().upper())); assert (all(box.lower()>=extent().lower() && box.lower()>=src->extent().lower())); assert (all(box.upper()<=extent().upper() && box.upper()<=src->extent().upper())); assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); + /* && box.stride()<=src->extent().stride() */ )); assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - assert (proc() == src->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const ibbox& sext = src->extent(); - const ibbox& dext = extent(); - - int srcshp[3], dstshp[3]; - int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; - - for (int d=0; d<3; ++d) { - srcshp[d] = (sext.shape() / sext.stride())[d]; - dstshp[d] = (dext.shape() / dext.stride())[d]; + /* && (box.lower()-src->extent().lower())%box.stride() == 0 */ )); + + if (_proc == src->_proc) { + // interpolate on same processor - srcbbox[0][d] = sext.lower()[d]; - srcbbox[1][d] = sext.upper()[d]; - srcbbox[2][d] = sext.stride()[d]; + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == _proc) { + + for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { + const ivect& pos = *posi; + + // get box around current position + const ibbox frombox + = ibbox(pos,pos,extent().stride()).expanded_for(src->extent()); + + // interpolate from box to position + T sum = 0; + for (ibbox::iterator fromposi=frombox.begin(); + fromposi!=frombox.end(); ++fromposi) + { + const ivect& frompos = *fromposi; + + // interpolation weight + const ivect str = src->extent().stride(); + const T f = prod(vect(str - abs(pos - frompos)) + / vect(str)); + sum += f * (*src)[frompos]; + } + (*this)[pos] = sum; + + } // for pos + + } - dstbbox[0][d] = dext.lower()[d]; - dstbbox[1][d] = dext.upper()[d]; - dstbbox[2][d] = dext.stride()[d]; + } else { + // interpolate from other processor - regbbox[0][d] = box.lower()[d]; - regbbox[1][d] = box.upper()[d]; - regbbox[2][d] = box.stride()[d]; - } - - assert (all(dext.stride() == box.stride())); - if (all(sext.stride() == dext.stride())) { - CCTK_FNAME(copy_3d_real8) ((const CCTK_REAL8*)src->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); + data* tmp = new data(box, src->_proc); + tmp->interpolate_from (src, box); + tmp->change_processor (_proc); + copy_from (tmp, box); + delete tmp; - } else { - assert (0); } } -template<> -void data -::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box) +template +void data::interpolate_from (const generic_data* gsrc, + const double sfact, + const generic_data* gtrc, + const double tfact, + const ibbox& box) { const data* src = (const data*)gsrc; - assert (has_storage() && src->has_storage()); + const data* trc = (const data*)gtrc; + assert (_has_storage && src->_has_storage && trc->_has_storage); assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); + && box.upper()<=extent().upper())); + assert (all(box.lower()>=extent().lower() + && box.lower()>=src->extent().lower() + && box.lower()>=trc->extent().lower())); assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); + && box.upper()<=src->extent().upper() + && box.upper()<=trc->extent().upper())); assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); + /* && box.stride()<=src->extent().stride() */ + /* && trc->extent().stride()==src->extent().stride() */ )); assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - assert (proc() == src->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const ibbox& sext = src->extent(); - const ibbox& dext = extent(); - - int srcshp[3], dstshp[3]; - int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; - - for (int d=0; d<3; ++d) { - srcshp[d] = (sext.shape() / sext.stride())[d]; - dstshp[d] = (dext.shape() / dext.stride())[d]; - - srcbbox[0][d] = sext.lower()[d]; - srcbbox[1][d] = sext.upper()[d]; - srcbbox[2][d] = sext.stride()[d]; - - dstbbox[0][d] = dext.lower()[d]; - dstbbox[1][d] = dext.upper()[d]; - dstbbox[2][d] = dext.stride()[d]; + && (box.lower()-src->extent().lower())%box.stride() == 0 + && (box.lower()-trc->extent().lower())%box.stride() == 0)); + + if (_proc == src->_proc && _proc == trc->_proc) { + // interpolate on same processor - regbbox[0][d] = box.lower()[d]; - regbbox[1][d] = box.upper()[d]; - regbbox[2][d] = box.stride()[d]; - } - - assert (all(dext.stride() == box.stride())); - if (all(sext.stride() == dext.stride())) { - CCTK_FNAME(copy_3d_complex16) ((const CCTK_COMPLEX16*)src->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_COMPLEX16*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, - dstbbox, - regbbox); + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == _proc) { + + for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) { + const ivect& pos = *posi; + + // get box around current position + const ibbox frombox + = ibbox(pos,pos,extent().stride()).expanded_for(src->extent()); + + // interpolate from box to position + T src_sum = 0; + T trc_sum = 0; + for (ibbox::iterator fromposi=frombox.begin(); + fromposi!=frombox.end(); ++fromposi) + { + const ivect& frompos = *fromposi; + + // interpolation weight + const ivect str = src->extent().stride(); + const T f = prod(vect(str - abs(pos - frompos)) + / vect(str)); + src_sum += f * (*src)[frompos]; + trc_sum += f * (*trc)[frompos]; + } + (*this)[pos] = (T)sfact * src_sum + (T)tfact * trc_sum; + + } // for pos + + } } else { - assert (0); + // interpolate from other processors + + data* smp = new data(box, src->_proc); + smp->copy_from (src, box); + data* tmp = new data(box, trc->_proc); + tmp->copy_from (trc, box); + interpolate_from (smp, sfact, tmp, tfact, box); + delete smp; + delete tmp; + } } - - -extern "C" { - - void CCTK_FCALL CCTK_FNAME(restrict_3d_real8) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(restrict_3d_real8_rf2) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - - - - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_rf2) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3_rf2) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_minmod) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_eno) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o5) - (const CCTK_REAL8* src, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_rf2) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_minmod) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_eno) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o5) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_rf2) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3_rf2) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_minmod) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_eno) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); - void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o5) - (const CCTK_REAL8* src1, const CCTK_REAL8& t1, - const CCTK_REAL8* src2, const CCTK_REAL8& t2, - const CCTK_REAL8* src3, const CCTK_REAL8& t3, - const int& srciext, const int& srcjext, const int& srckext, - CCTK_REAL8* dst, const CCTK_REAL8& t, - const int& dstiext, const int& dstjext, const int& dstkext, - const int srcbbox[3][3], - const int dstbbox[3][3], - const int regbbox[3][3]); -} - -template<> -void data -::interpolate_from_innerloop (const vector*> gsrcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) +// Output +template +template +void data::write_ascii (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const { - const CCTK_REAL eps = 1.0e-10; + assert (_has_storage); - assert (has_storage()); - assert (all(box.lower()>=extent().lower())); - assert (all(box.upper()<=extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - vector srcs(gsrcs.size()); - for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t]; - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs[t]->has_storage()); - assert (all(box.lower()>=srcs[t]->extent().lower())); - assert (all(box.upper()<=srcs[t]->extent().upper())); - } - - assert (proc() == srcs[0]->proc()); - - int rank; - MPI_Comm_rank (dist::comm, &rank); - assert (rank == proc()); - - const ibbox& sext = srcs[0]->extent(); - const ibbox& dext = extent(); - - int srcshp[3], dstshp[3]; - int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; - - for (int d=0; d<3; ++d) { - srcshp[d] = (sext.shape() / sext.stride())[d]; - dstshp[d] = (dext.shape() / dext.stride())[d]; - - srcbbox[0][d] = sext.lower()[d]; - srcbbox[1][d] = sext.upper()[d]; - srcbbox[2][d] = sext.stride()[d]; + if (_proc==0) { + // output on processor 0 - dstbbox[0][d] = dext.lower()[d]; - dstbbox[1][d] = dext.upper()[d]; - dstbbox[2][d] = dext.stride()[d]; - - regbbox[0][d] = box.lower()[d]; - regbbox[1][d] = box.upper()[d]; - regbbox[2][d] = box.stride()[d]; - } - - // Check that the times are consistent - assert (times.size() > 0); - CCTK_REAL min_time = times[0]; - CCTK_REAL max_time = times[0]; - for (size_t tl=1; tl 1.2); - min_time = min(min_time, times[tl]); - max_time = max(max_time, times[tl]); - } - if (time < min_time - eps || time > max_time + eps) { - ostringstream buf; - buf << "Internal error: extrapolation in time." - << " time=" << time - << " times=" << times; - CCTK_WARN (0, buf.str().c_str()); - } - - // Is it necessary to interpolate in time? - if (times.size() > 1) { - for (size_t tl=0; tl 1.4); - if (abs(times[tl] - time) < eps) { - // It is not. - vector*> my_gsrcs(1); - vector my_times(1); - my_gsrcs[0] = gsrcs[tl]; - my_times[0] = times[tl]; - const int my_order_time = 0; - this->interpolate_from_innerloop - (my_gsrcs, my_times, box, time, order_space, my_order_time); - return; - } - } - } - - assert (all(dext.stride() == box.stride())); - if (all(sext.stride() < dext.stride())) { - // Restrict - - assert (times.size() == 1); - assert (abs(times[0] - time) < eps); - - switch (transport_operator) { - - case op_Lagrange: - case op_TVD: - case op_ENO: - assert (srcs.size() == 1); - if (all (dext.stride() == sext.stride() * 2)) { - CCTK_FNAME(restrict_3d_real8_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(restrict_3d_real8) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == 0) { - default: - assert (0); + const vect lo = extent().lower()[dirs]; + const vect up = extent().upper()[dirs]; + const vect str = extent().stride()[dirs]; + const bbox ext(lo,up,str); - } // switch (transport_operator) - - } else if (all(sext.stride() > dext.stride())) { - // Prolongate - - switch (transport_operator) { + ofstream file(name.c_str(), ios::app); + assert (file.good()); - case op_Lagrange: - switch (order_time) { - - case 0: - assert (times.size() == 1); - assert (abs(times[0] - time) < eps); - assert (srcs.size()>=1); - switch (order_space) { - case 0: - case 1: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 2: - case 3: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_o3_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - case 1: - assert (srcs.size()>=2); - switch (order_space) { - case 0: - case 1: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_2tl_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8_2tl) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 2: - case 3: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_2tl_o3_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8_2tl_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_2tl_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - case 2: - assert (srcs.size()>=3); - switch (order_space) { - case 0: - case 1: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_3tl_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8_3tl) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 2: - case 3: - if (all (sext.stride() == dext.stride() * 2)) { - CCTK_FNAME(prolongate_3d_real8_3tl_o3_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } else { - CCTK_FNAME(prolongate_3d_real8_3tl_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - } - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_3tl_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - default: - assert (0); - } // switch (order_time) - break; + file << "# time=" << time << " tl=" << tl << " rl=" << rl + << " c=" << c << " ml=" << ml << " "; + assert (DD<=3); + for (int d=0; dstorage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); -// CCTK_FNAME(prolongate_3d_real8_minmod) -// ((const CCTK_REAL8*)srcs[0]->storage(), -// srcshp[0], srcshp[1], srcshp[2], -// (CCTK_REAL8*)storage(), -// dstshp[0], dstshp[1], dstshp[2], -// srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - case 1: - switch (order_space) { - case 0: - case 1: - CCTK_WARN (0, "There is no stencil for op=\"TVD\" with order_space=1"); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_2tl_rf2) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); -// CCTK_FNAME(prolongate_3d_real8_2tl_minmod) -// ((const CCTK_REAL8*)srcs[0]->storage(), times[0], -// (const CCTK_REAL8*)srcs[1]->storage(), times[1], -// srcshp[0], srcshp[1], srcshp[2], -// (CCTK_REAL8*)storage(), time, -// dstshp[0], dstshp[1], dstshp[2], -// srcbbox, dstbbox, regbbox); - break; - default: - assert (0); + for (bbox::iterator it=ext.begin(); it!=ext.end(); ++it) { + ivect index(0); + for (int d=0; dstorage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); -// CCTK_FNAME(prolongate_3d_real8_3tl_minmod) -// ((const CCTK_REAL8*)srcs[0]->storage(), times[0], -// (const CCTK_REAL8*)srcs[1]->storage(), times[1], -// (const CCTK_REAL8*)srcs[2]->storage(), times[2], -// srcshp[0], srcshp[1], srcshp[2], -// (CCTK_REAL8*)storage(), time, -// dstshp[0], dstshp[1], dstshp[2], -// srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - default: - assert (0); - } // switch (order_time) - break; + } - case op_ENO: - switch (order_time) { - case 0: - assert (times.size() == 1); - assert (abs(times[0] - time) < eps); - switch (order_space) { - case 0: - case 1: - CCTK_WARN (0, "There is no stencil for op=\"ENO\" with order_space=1"); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_eno) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - case 1: - switch (order_space) { - case 0: - case 1: - CCTK_WARN (0, "There is no stencil for op=\"ENO\" with order_space=1"); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_2tl_eno) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - case 2: - switch (order_space) { - case 0: - case 1: - CCTK_WARN (0, "There is no stencil for op=\"ENO\" with order_space=1"); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_3tl_eno) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - default: - assert (0); - } // switch (order_time) - break; + file.close(); + assert (file.good()); - default: - assert(0); - } // switch (transport_operator) + } } else { - assert (0); + // copy to processor 0 and output there + + data* tmp = new data(_extent, 0); + tmp->copy_from (this, _extent); + tmp->write_ascii (name, time, dirs, tl, rl, c, ml); + delete tmp; + } } - - // Output template -ostream& data::output (ostream& os) const -{ - T Tdummy; - os << "data<" << typestring(Tdummy) << "," << D << ">:" - << "extent=" << this->extent() << "," - << "stride=" << this->stride() << ",size=" << this->size(); +ostream& data::out (ostream& os) const { + os << "data<" STR(T) "," << D << ">:" + << "extent=" << extent() << "," + << "stride=" << stride() << ",size=" << size(); return os; } -#define INSTANTIATE(T) \ -template class data; +#if defined(TMPL_EXPLICIT) + +#define INSTANTIATE(T) \ + \ +template data; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; \ + \ +template data; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; \ + \ +template data; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; \ +template void data::write_ascii (const string name, const double time, \ + const vect& dirs, \ + const int tl, const int rl, \ + const int c, const int ml) const; #include "instantiate" #undef INSTANTIATE + +#endif diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index ba44556ff..c528833e1 100644 --- a/Carpet/CarpetLib/src/data.hh +++ b/Carpet/CarpetLib/src/data.hh @@ -1,130 +1,156 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.23 2004/06/08 22:57:22 schnetter Exp $ +/*************************************************************************** + data.hh - Data storage + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef DATA_HH #define DATA_HH #include -#include #include -#include - -#include "cctk.h" #include "defs.hh" #include "dist.hh" #include "bbox.hh" +#include "bboxset.hh" #include "gdata.hh" #include "vect.hh" -using namespace std; + + +// Forward definition +template class data; + +// Output +template +ostream& operator<< (ostream& os, const data& d); -// A distributed multi-dimensional array +// A real data storage template -class data: public gdata -{ +class data: public generic_data { // Types typedef vect ivect; typedef bbox ibbox; + typedef bboxset ibset; // Fields - T* _storage; // the data (if located on this processor) - size_t _allocated_bytes; // number of allocated bytes - - // For vector groups with contiguous storage - int vectorlength; // number of vector elements - int vectorindex; // index of this vector element - data* vectorleader; // if index!=0: first vector element - vector vectorclients; // if index==0: registered elements - - void register_client (int index); - void unregister_client (int index); - bool has_clients (); - + T* restrict _storage; // the data (if located on this processor) + public: // Constructors - data (const int varindex = -1, - const operator_type transport_operator = op_error, - const int vectorlength = 1, const int vectorindex = 0, - data* const vectorleader = NULL); - data (const int varindex, const operator_type transport_operator, - const int vectorlength, const int vectorindex, - data* const vectorleader, - const ibbox& extent, const int proc); + data (); + data (const ibbox& extent, const int proc); // Destructors virtual ~data (); // Pseudo constructors - virtual data* make_typed (const int varindex, - const operator_type transport_operator) const; + virtual data* make_typed (const ibbox& extent, const int proc) const; // Storage management -private: - void getmem (const size_t nelems); - void freemem (); -public: - virtual void allocate (const ibbox& extent, const int proc, - void* const mem=0); + virtual void allocate (const ibbox& extent, const int proc); virtual void free (); - virtual void transfer_from (gdata* gsrc); - -private: - T* vectordata (const int vectorindex) const; -public: + virtual void transfer_from (generic_data* gsrc); // Processor management - virtual void change_processor (comm_state& state, - const int newproc, void* const mem=0); -private: - virtual void change_processor_recv (const int newproc, void* const mem=0); - virtual void change_processor_send (const int newproc, void* const mem=0); - virtual void change_processor_wait (const int newproc, void* const mem=0); -public: + virtual void change_processor (const int newproc); // Accessors - virtual const void* storage () const - { - assert (this->_has_storage); + virtual const T* storage () const { + assert (_has_storage); return _storage; } - virtual void* storage () { - assert (this->_has_storage); + virtual T* storage () { + assert (_has_storage); return _storage; } // Data accessors - const T& operator[] (const ivect& index) const - { + const T& operator[] (const ivect& index) const { assert (_storage); return _storage[offset(index)]; } - T& operator[] (const ivect& index) - { + T& operator[] (const ivect& index) { assert (_storage); return _storage[offset(index)]; } // Data manipulators - void copy_from_innerloop (const gdata* gsrc, - const ibbox& box); - void interpolate_from_innerloop (const vector*> gsrcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); + virtual void copy_from (const generic_data* gsrc, const ibbox& b); + virtual void interpolate_from (const generic_data* gsrc, + const ibbox& box); + virtual void interpolate_from (const generic_data* gsrc, + const double sfact, + const generic_data* gtrc, + const double tfact, + const ibbox& box); + // Output + template + void write_ascii (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const; +protected: + virtual void write_ascii_1 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) const { + write_ascii (name, time, dirs, tl, rl, c, ml); + } + virtual void write_ascii_2 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) const { + write_ascii (name, time, dirs, tl, rl, c, ml); + } + virtual void write_ascii_3 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) const { + write_ascii (name, time, dirs, tl, rl, c, ml); + } +// void write_ieee (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; +// void write_hdf (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; +// void write_h5 (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; public: // Output - ostream& output (ostream& os) const; + ostream& out (ostream& os) const; }; +#if defined(TMPL_IMPLICIT) +# include "data.cc" +#endif + #endif // DATA_HH diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index 69cd22fca..94315703d 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.cc @@ -1,105 +1,40 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.19 2004/05/21 18:13:41 schnetter Exp $ - -#include -#include - +/*************************************************************************** + defs.cc - Commonly used definitions + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include #include #include -#include #include -#include "cctk.h" - -#include "defs.hh" - -using namespace std; - - - -template -T ipow (T x, int y) { - if (y<0) { - y = -y; - x = T(1)/x; - } - T res = T(1); - for (;;) { - if (y%2) res *= x; - y /= 2; - if (y==0) break; - x *= x; - } - return res; -} - - - -void skipws (istream& is) { - while (is.good() && isspace(is.peek())) { - is.get(); - } -} - - - -void expect (istream& is, const char c) { - if (is.peek() == c) return; - cout << "While reading characters from a stream:" << endl - << " Character '" << c << "' expected, but not found." << endl - << " The next up to 100 available characters are \""; - for (int i=0; i<100; ++i) { - const int uc = is.get(); - if (uc<0) break; - cout << (unsigned char)uc; - } - cout << "\"." << endl; - throw input_error(); -} - - - -void consume (istream& is, const char c) { - expect (is, c); - is.get(); -} - - - -// Vector input -template -istream& input (istream& is, vector& v) { - v.clear(); - try { - skipws (is); - consume (is, '['); - skipws (is); - while (is.good() && is.peek() != ']') { - T elem; - is >> elem; - v.push_back (elem); - skipws (is); - if (is.peek() != ',') break; - is.get(); - skipws (is); - } - skipws (is); - consume (is, ']'); - } catch (input_error &err) { - cout << "Input error while reading a vector<>" << endl - << " The following elements have been read so far: " << v << endl; - throw err; - } - return is; -} +#if !defined(TMPL_IMPLICIT) || !defined(DEFS_HH) +# include "defs.hh" +#endif // List output template -ostream& output (ostream& os, const list& l) { +ostream& operator<< (ostream& os, const list& l) { os << "["; - for (typename list::const_iterator ti=l.begin(); ti!=l.end(); ++ti) { + for (list::const_iterator ti=l.begin(); ti!=l.end(); ++ti) { if (ti!=l.begin()) os << ","; os << *ti; } @@ -109,9 +44,9 @@ ostream& output (ostream& os, const list& l) { // Set output template -ostream& output (ostream& os, const set& s) { +ostream& operator<< (ostream& os, const set& s) { os << "{"; - for (typename set::const_iterator ti=s.begin(); ti!=s.end(); ++ti) { + for (set::const_iterator ti=s.begin(); ti!=s.end(); ++ti) { if (ti!=s.begin()) os << ","; os << *ti; } @@ -119,24 +54,12 @@ ostream& output (ostream& os, const set& s) { return os; } -// Stack output -template -ostream& output (ostream& os, const stack& s) { - stack s2 (s); - list l; - while (! s2.empty()) { - l.insert (l.begin(), s2.top()); - s2.pop(); - } - return output (os, l); -} - // Vector output template -ostream& output (ostream& os, const vector& v) { +ostream& operator<< (ostream& os, const vector& v) { os << "["; int cnt=0; - for (typename vector::const_iterator ti=v.begin(); ti!=v.end(); ++ti) { + for (vector::const_iterator ti=v.begin(); ti!=v.end(); ++ti) { if (ti!=v.begin()) os << ","; os << cnt++ << ":" << *ti; } @@ -146,31 +69,28 @@ ostream& output (ostream& os, const vector& v) { +#if defined(TMPL_EXPLICIT) #include "bbox.hh" #include "bboxset.hh" -template int ipow (int x, int y); - -template istream& input (istream& os, vector >& v); -template istream& input (istream& os, vector >& v); -template istream& input (istream& os, vector > >& v); -template istream& input (istream& os, vector > >& v); -template istream& input (istream& os, vector,3> >& v); -template istream& input (istream& os, vector,3> > >& v); - -template ostream& output (ostream& os, const list >& l); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const stack >& s); -template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector& v); -template ostream& output (ostream& os, const vector >& v); -template ostream& output (ostream& os, const vector >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector >& v); -template ostream& output (ostream& os, const vector >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector,3> > >& v); -template ostream& output (ostream& os, const vector,3> >& v); -template ostream& output (ostream& os, const vector > > >& v); +template ostream& operator<< (ostream& os, const list >& l); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > > >& v); + +template ostream& operator<< (ostream& os, const list >& l); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > > >& v); + +template ostream& operator<< (ostream& os, const list >& l); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const set >& s); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > >& v); +template ostream& operator<< (ostream& os, const vector > > >& v); +#endif diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh index 4297ba8c6..741678370 100644 --- a/Carpet/CarpetLib/src/defs.hh +++ b/Carpet/CarpetLib/src/defs.hh @@ -1,4 +1,22 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.13 2004/03/01 19:43:08 schnetter Exp $ +/*************************************************************************** + defs.hh - Commonly used definitions + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef DEFS_HH #define DEFS_HH @@ -7,188 +25,35 @@ #include #endif -#include - #include -#include -#include #include #include -#include #include -#include "cctk.h" +// Stringification +#define STR(s) #s + +// Fortranification +#define FORTRAN_NAME(x) x##_ -using namespace std; +// Fortran style function arguments +#define restrict __restrict__ // A general type enum centering { vertex_centered, cell_centered }; -// Useful helper template inline T square (const T& x) { return x*x; } -// Another useful helper -template -T ipow (T x, int y); - - - -// Input streams -struct input_error { }; -void skipws (istream& is); -void expect (istream& is, const char c); -void consume (istream& is, const char c); - - - -// Names for types - -#if 0 - -inline const char * typestring (const char& dummy) -{ return "char"; } - -inline const char * typestring (const signed char& dummy) -{ return "signed char"; } - -inline const char * typestring (const unsigned char& dummy) -{ return "unsigned char"; } - -inline const char * typestring (const short& dummy) -{ return "short"; } - -inline const char * typestring (const unsigned short& dummy) -{ return "unsigned short"; } - -inline const char * typestring (const int& dummy) -{ return "int"; } - -inline const char * typestring (const unsigned int& dummy) -{ return "unsigned int"; } - -inline const char * typestring (const long& dummy) -{ return "long"; } - -inline const char * typestring (const unsigned long& dummy) -{ return "unsigned long"; } - -inline const char * typestring (const long long& dummy) -{ return "long long"; } - -inline const char * typestring (const unsigned long long& dummy) -{ return "unsigned long long"; } - -inline const char * typestring (const float& dummy) -{ return "float"; } - -inline const char * typestring (const double& dummy) -{ return "double"; } - -inline const char * typestring (const long double& dummy) -{ return "long double"; } - -inline const char * typestring (const complex& dummy) -{ return "complex"; } - -inline const char * typestring (const complex& dummy) -{ return "complex"; } - -inline const char * typestring (const complex& dummy) -{ return "complex"; } - -#else - -# ifdef CCTK_INT1 -inline const char * typestring (const CCTK_INT1& dummy) -{ return "CCTK_INT1"; } -# endif - -# ifdef CCTK_INT2 -inline const char * typestring (const CCTK_INT2& dummy) -{ return "CCTK_INT2"; } -# endif - -# ifdef CCTK_INT4 -inline const char * typestring (const CCTK_INT4& dummy) -{ return "CCTK_INT4"; } -# endif - -# ifdef CCTK_INT8 -inline const char * typestring (const CCTK_INT8& dummy) -{ return "CCTK_INT8"; } -# endif - -# ifdef CCTK_REAL4 -inline const char * typestring (const CCTK_REAL4& dummy) -{ return "CCTK_REAL4"; } -# endif - -# ifdef CCTK_REAL8 -inline const char * typestring (const CCTK_REAL8& dummy) -{ return "CCTK_REAL8"; } -# endif - -# ifdef CCTK_REAL16 -inline const char * typestring (const CCTK_REAL16& dummy) -{ return "CCTK_REAL16"; } -# endif - -# ifdef CCTK_REAL4 -inline const char * typestring (const CCTK_COMPLEX8& dummy) -{ return "CCTK_COMPLEX8"; } -# endif - -# ifdef CCTK_REAL8 -inline const char * typestring (const CCTK_COMPLEX16& dummy) -{ return "CCTK_COMPLEX16"; } -# endif - -# ifdef CCTK_REAL16 -inline const char * typestring (const CCTK_COMPLEX32& dummy) -{ return "CCTK_COMPLEX32"; } -# endif - -#endif - - - -// Container input -template istream& input (istream& is, vector& v); - -template -inline istream& operator>> (istream& is, vector& v) { - return input(is,v); -} - - - // Container output -template ostream& output (ostream& os, const list& l); -template ostream& output (ostream& os, const set& s); -template ostream& output (ostream& os, const stack& s); -template ostream& output (ostream& os, const vector& v); - -template -inline ostream& operator<< (ostream& os, const list& l) { - return output(os,l); -} +template ostream& operator<< (ostream& os, const list& l); +template ostream& operator<< (ostream& os, const set& s); +template ostream& operator<< (ostream& os, const vector& v); -template -inline ostream& operator<< (ostream& os, const set& s) { - return output(os,s); -} - -template -inline ostream& operator<< (ostream& os, const stack& s) { - return output(os,s); -} - -template -inline ostream& operator<< (ostream& os, const vector& v) { - return output(os,v); -} +#if defined(TMPL_IMPLICIT) +# include "defs.cc" +#endif #endif // DEFS_HH diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index e7e8778ff..804d52f4c 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1,91 +1,80 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.57 2004/09/17 16:37:26 schnetter Exp $ +/*************************************************************************** + dh.cc - Data Hierarchy + a grid hierarchy plus ghost zones + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu -#include + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -#include "cctk.h" -#include "cctk_Parameters.h" + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include "defs.hh" -#include "dist.hh" #include "ggf.hh" #include "vect.hh" -#include "dh.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) || !defined(DH_HH) +# include "dh.hh" +#endif // Constructors template -dh::dh (gh& h, - const ivect& lghosts, const ivect& ughosts, - const int prolongation_order_space, const int buffer_width) - : h(h), - lghosts(lghosts), ughosts(ughosts), - prolongation_order_space(prolongation_order_space), - buffer_width(buffer_width) +dh::dh (gh& h, const ivect& lghosts, const ivect& ughosts) + : h(h), lghosts(lghosts), ughosts(ughosts) { assert (all(lghosts>=0 && ughosts>=0)); - assert (prolongation_order_space>=0); - assert (buffer_width>=0); h.add(this); - CHECKPOINT; - recompose (false); + recompose(); } // Destructors template -dh::~dh () -{ - CHECKPOINT; +dh::~dh () { h.remove(this); } -// Helpers -template -int dh::prolongation_stencil_size () const { - assert (prolongation_order_space>=0); - return prolongation_order_space/2; -} - // Modifiers template -void dh::recompose (const bool do_prolongate) { - DECLARE_CCTK_PARAMETERS; - - CHECKPOINT; +void dh::recompose () { boxes.clear(); boxes.resize(h.reflevels()); for (int rl=0; rl::recompose (const bool do_prolongate) { // Sync boxes const int cs = h.components(rl); - boxes.at(rl).at(c).at(ml).send_sync.resize(cs); - boxes.at(rl).at(c).at(ml).recv_sync.resize(cs); + boxes[rl][c][ml].send_sync.resize(cs); + boxes[rl][c][ml].recv_sync.resize(cs); // Refinement boxes if (rl>0) { const int csm1 = h.components(rl-1); - boxes.at(rl).at(c).at(ml).send_ref_coarse.resize(csm1); - boxes.at(rl).at(c).at(ml).recv_ref_coarse.resize(csm1); - boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse.resize(csm1); + boxes[rl][c][ml].send_ref_coarse.resize(csm1); + boxes[rl][c][ml].recv_ref_coarse.resize(csm1); + boxes[rl][c][ml].recv_ref_bnd_coarse.resize(csm1); } if (rl::recompose (const bool do_prolongate) { for (int rl=0; rl0) { - const ibbox intrf = boxes.at(rl).at(c).at(ml-1).interior; - const ibbox extrf = boxes.at(rl).at(c).at(ml-1).exterior; + const ibbox intrf = boxes[rl][c][ml-1].interior; + const ibbox extrf = boxes[rl][c][ml-1].exterior; // Restriction (interior) { // (the restriction must fill all of the interior of the // coarse grid, and may use the exterior of the fine grid) const ibbox recv = intr; - assert (intr.empty() || ! recv.empty()); const ibbox send = recv.expanded_for(extrf); - assert (intr.empty() || ! send.empty()); - // TODO: put the check back in, taking outer boundaries - // into account -#if 0 - assert (send.is_contained_in(extrf)); -#endif - boxes.at(rl).at(c).at(ml-1).send_mg_coarse.push_back(send); - boxes.at(rl).at(c).at(ml ).recv_mg_fine .push_back(recv); + assert (send.contained_in(extrf)); + boxes[rl][c][ml-1].send_mg_coarse.push_back(send); + boxes[rl][c][ml ].recv_mg_fine .push_back(recv); } // Prolongation (interior) { @@ -172,195 +145,58 @@ void dh::recompose (const bool do_prolongate) { // grid, and may fill only the interior of the fine grid, // and the bbox must be as large as possible) const ibbox recv = extr.contracted_for(intrf) & intrf; - assert (intr.empty() || ! recv.empty()); const ibbox send = recv.expanded_for(extr); - assert (intr.empty() || ! send.empty()); - boxes.at(rl).at(c).at(ml-1).recv_mg_coarse.push_back(recv); - boxes.at(rl).at(c).at(ml ).send_mg_fine .push_back(send); + boxes[rl][c][ml-1].recv_mg_coarse.push_back(recv); + boxes[rl][c][ml ].send_mg_fine .push_back(send); } - } // if not finest multigrid level - - } // for ml - } // for c - } // for rl - - for (int rl=0; rlbegin(); - li!=lvi->end(); ++li) { - recvs -= *li; - } - } - recvs.normalize(); - assert (recvs.setsize() <= 1); - if (recvs.setsize() == 1) { - const ibbox recv = *recvs.begin(); - const ibbox send = recv.expanded_for(extr); - assert (! send.empty()); - assert (send.is_contained_in(extr)); - boxes.at(rl+1).at(cc).at(ml).recv_ref_coarse.at(c ) - .push_back(recv); - boxes.at(rl ).at(c ).at(ml).send_ref_fine .at(cc) - .push_back(send); - } + const ibbox recv = extr.contracted_for(intrf) & intrf; + const ibbox send = recv.expanded_for(extr); + assert (send.contained_in(extr)); + boxes[rl+1][cc][ml].recv_ref_coarse[c ].push_back(recv); + boxes[rl ][c ][ml].send_ref_fine [cc].push_back(send); } - - } // for cc - } // if not finest refinement level - - } // for ml - } // for c - } // for rl - - for (int rl=0; rlbegin(); - li!=lvi->end(); ++li) { - pbndsf -= *li; - } - } - pbndsf.normalize(); - } - // Buffer zones - ibset buffers; - { - for (typename ibset::const_iterator pbi=pbndsf.begin(); - pbi!=pbndsf.end(); ++pbi) { - buffers |= (*pbi).expand(buffer_width, buffer_width) & extrf; - } - buffers.normalize(); - } - // Add boundaries - const ibbox maxrecvs - = extr.expand(-pss,-pss).contracted_for(extrf); - ibset recvs = buffers & maxrecvs; - recvs.normalize(); - { - // Do not prolongate what is already prolongated - const iblistvect& rrbc - = boxes.at(rl+1).at(cc).at(ml).recv_ref_bnd_coarse; - for (typename iblistvect::const_iterator lvi=rrbc.begin(); - lvi!=rrbc.end(); ++lvi) { - for (typename iblist::const_iterator li=lvi->begin(); - li!=lvi->end(); ++li) { - recvs -= *li; - } - } - recvs.normalize(); - } - { - for (typename ibset::const_iterator ri = recvs.begin(); - ri != recvs.end(); ++ri) { - const ibbox & recv = *ri; - const ibbox send = recv.expanded_for(extr); - assert (! send.empty()); - assert (send.is_contained_in(extr)); - assert (send.is_contained_in(extr.expand(-pss,-pss))); - boxes.at(rl+1).at(cc).at(ml).recv_ref_bnd_coarse.at(c ) - .push_back(recv); - boxes.at(rl ).at(c ).at(ml).send_ref_bnd_fine .at(cc) - .push_back(send); - } + ibset bndsf = boxes[rl+1][cc][ml].boundaries; + // coarsify boundaries of fine component + for (ibset::const_iterator bi=bndsf.begin(); + bi!=bndsf.end(); ++bi) { + const ibbox& bndf = *bi; + // (the prolongation may use the exterior of the + // coarse grid, and must fill all of the boundary of + // the fine grid) + const ibbox recv = extr.contracted_for(bndf) & bndf; + const ibbox send = recv.expanded_for(extr); + assert (send.contained_in(extr)); + boxes[rl+1][cc][ml].recv_ref_bnd_coarse[c ].push_back(recv); + boxes[rl ][c ][ml].send_ref_bnd_fine [cc].push_back(send); } } - - } // for cc - } // if not finest refinement level - - } // for ml - } // for c - } // for rl - - for (int rl=0; rl::recompose (const bool do_prolongate) { for (int rl=0; rlbegin(); - li!=lvi->end(); ++li) { - sync_not -= *li; - recv_not -= *li; - } - } - - // Subtract boxes received during prolongation - const iblistvect& recv_ref_bnd_coarse - = boxes.at(rl).at(c).at(ml).recv_ref_bnd_coarse; - for (typename iblistvect::const_iterator - lvi=recv_ref_bnd_coarse.begin(); - lvi!=recv_ref_bnd_coarse.end(); ++lvi) { - for (typename iblist::const_iterator li=lvi->begin(); - li!=lvi->end(); ++li) { - recv_not -= *li; - } - } - - } // for ml - } // for c - } // for rl - - // Calculate bases - bases.resize(h.reflevels()); - for (int rl=0; rl:" - << "ghosts=[" << lghosts << "," << ughosts << "]," + << "ghosts=[" << d.lghosts << "," << d.ughosts << "]," << "gfs={"; int cnt=0; - for (typename list*>::const_iterator f = gfs.begin(); - f != gfs.end(); ++f) { + for (list*>::const_iterator f = d.gfs.begin(); + f != d.gfs.end(); ++f) { if (cnt++) os << ","; - (*f)->output(os); + os << **f; } os << "}"; + return os; } +#if defined(TMPL_EXPLICIT) +template class dh<1>; +template ostream& operator<< (ostream& os, const dh<1>& d); + +template class dh<2>; +template ostream& operator<< (ostream& os, const dh<2>& d); + template class dh<3>; +template ostream& operator<< (ostream& os, const dh<3>& d); +#endif diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index 22a9bd425..0d3c4c716 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -1,10 +1,28 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.20 2004/08/07 19:47:11 schnetter Exp $ +/*************************************************************************** + dh.hh - Data Hierarchy + A grid hierarchy plus ghost zones + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef DH_HH #define DH_HH -#include - +#include #include #include #include @@ -16,13 +34,11 @@ #include "gh.hh" #include "vect.hh" -using namespace std; - // Forward declaration -template class ggf; template class dh; +template class generic_gf; // Output template @@ -41,16 +57,6 @@ class dh { typedef list iblist; typedef vector iblistvect; // vector of lists -public: - - // in here, the term "boundary" means both ghost zones and - // refinement boundaries, but does not refer to outer (physical) - // boundaries. - - // ghost zones, refinement boundaries, and outer boundaries are not - // used as sources for synchronisation. this design choice might - // not be good. - struct dboxes { ibbox exterior; // whole region (including boundaries) @@ -64,83 +70,50 @@ public: iblistvect recv_ref_fine; iblistvect recv_ref_coarse; iblistvect send_sync; // send while syncing - iblistvect send_ref_bnd_fine; // sent to finer grids + iblistvect send_ref_bnd_fine; ibset boundaries; // boundaries - iblistvect recv_sync; // received while syncing - iblistvect recv_ref_bnd_coarse; // received from coarser grids - ibset sync_not; // not received while syncing (outer boundary of that level) - ibset recv_not; // not received while syncing or prolongating (globally outer boundary) - -#if 0 - // after regridding: - iblistvect prev_send; // sent from previous dh - iblistvect recv_prev; // received from previous dh - iblistvect send_prev_fine; // sent to finer - iblistvect recv_prev_coarse; // received from coarser -#endif - }; - -private: - - struct dbases { - ibbox exterior; // whole region (including boundaries) - ibbox interior; // interior (without boundaries) - ibset boundaries; // boundaries + iblistvect recv_sync; // receive while syncing + iblistvect recv_ref_bnd_coarse; }; typedef vector mboxes; // ... for each multigrid level typedef vector cboxes; // ... for each component typedef vector rboxes; // ... for each refinement level - typedef vector mbases; // ... for each multigrid level - typedef vector rbases; // ... for each refinement level - public: // should be readonly // Fields - gh& h; // hierarchy + gh &h; // hierarchy ivect lghosts, ughosts; // ghost zones - int prolongation_order_space; // order of spatial prolongation operator - int buffer_width; // buffer inside refined grids - rboxes boxes; - rbases bases; - list*> gfs; // list of all grid functions + list*> gfs; public: // Constructors - dh (gh& h, const ivect& lghosts, const ivect& ughosts, - int prolongation_order_space, int buffer_width); + dh (gh& h, const ivect& lghosts, const ivect& ughosts); // Destructors - virtual ~dh (); - - // Helpers - int prolongation_stencil_size () const; + ~dh (); // Modifiers - void recompose (const bool do_prolongate); + void recompose (); // Grid function management - void add (ggf* f); - void remove (ggf* f); + void add (generic_gf* f); + void remove (generic_gf* f); // Output - virtual void output (ostream& os) const; + friend ostream& operator<< <> (ostream& os, const dh& d); }; -template -inline ostream& operator<< (ostream& os, const dh& d) { - d.output(os); - return os; -} - - +#if defined(TMPL_IMPLICIT) +# include "dh.cc" +#endif #endif // DH_HH diff --git a/Carpet/CarpetLib/src/dist.cc b/Carpet/CarpetLib/src/dist.cc index 7e41e3cdb..fc847a523 100644 --- a/Carpet/CarpetLib/src/dist.cc +++ b/Carpet/CarpetLib/src/dist.cc @@ -1,83 +1,53 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.cc,v 1.8 2004/07/08 12:36:01 tradke Exp $ +/*************************************************************************** + dist.cc - Helpers for distributed computing + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu -#include + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -#include + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ -#include "cctk.h" -#include "cctk_Parameters.h" +#include + +#include #include "defs.hh" -#include "dist.hh" +#if !defined(TMPL_IMPLICIT) || !defined(DIST_HH) +# include "dist.hh" +#endif -using namespace std; +MPI_Comm dist::comm; -namespace dist { - - MPI_Comm comm = MPI_COMM_NULL; - -#if 0 - MPI_Datatype mpi_complex_float; - MPI_Datatype mpi_complex_double; - MPI_Datatype mpi_complex_long_double; -#else - MPI_Datatype mpi_complex8; - MPI_Datatype mpi_complex16; - MPI_Datatype mpi_complex32; -#endif - - void init (int& argc, char**& argv) { - MPI_Init (&argc, &argv); - pseudoinit(); - } - - void pseudoinit () { - comm = MPI_COMM_WORLD; - -#if 0 - MPI_Type_contiguous (2, MPI_FLOAT, &mpi_complex_float); - MPI_Type_commit (&mpi_complex_float); - MPI_Type_contiguous (2, MPI_DOUBLE, &mpi_complex_double); - MPI_Type_commit (&mpi_complex_double); - MPI_Type_contiguous (2, MPI_LONG_DOUBLE, &mpi_complex_long_double); - MPI_Type_commit (&mpi_complex_long_double); -#else -# ifdef CCTK_REAL4 - CCTK_REAL4 dummy4; - MPI_Type_contiguous (2, datatype(dummy4), &mpi_complex8); - MPI_Type_commit (&mpi_complex8); -# endif -# ifdef CCTK_REAL8 - CCTK_REAL8 dummy8; - MPI_Type_contiguous (2, datatype(dummy8), &mpi_complex16); - MPI_Type_commit (&mpi_complex16); -# endif -# ifdef CCTK_REAL16 - CCTK_REAL16 dummy16; - MPI_Type_contiguous (2, datatype(dummy16), &mpi_complex32); - MPI_Type_commit (&mpi_complex32); -# endif + + +void dist::init (int& argc, char**& argv) { + MPI_Init (&argc, &argv); + comm = MPI_COMM_WORLD; +} + +void dist::pseudoinit () { + comm = MPI_COMM_WORLD; +} + +void dist::finalize () { + MPI_Finalize (); +} + + + +#if defined(TMPL_EXPLICIT) #endif - } - - void finalize () { - MPI_Finalize (); - } - - void checkpoint (const char* file, int line) { - DECLARE_CCTK_PARAMETERS; - if (verbose) { - int rank; - MPI_Comm_rank (comm, &rank); - printf ("CHECKPOINT: processor %d, file %s, line %d\n", - rank, file, line); - } - if (barriers) { - MPI_Barrier (comm); - } - } - -} // namespace dist diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh index ab89ae580..09645a716 100644 --- a/Carpet/CarpetLib/src/dist.hh +++ b/Carpet/CarpetLib/src/dist.hh @@ -1,120 +1,112 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.10 2004/03/01 19:43:39 schnetter Exp $ +/*************************************************************************** + dist.hh - Helpers for distributed computing + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef DIST_HH #define DIST_HH -#include -#include -#include - -#include +#include +#include #include -#include "cctk.h" - #include "defs.hh" -using namespace std; - -namespace dist { - - extern MPI_Comm comm; - -#if 0 - extern MPI_Datatype mpi_complex_float; - extern MPI_Datatype mpi_complex_double; - extern MPI_Datatype mpi_complex_long_double; -#else - extern MPI_Datatype mpi_complex8; - extern MPI_Datatype mpi_complex16; - extern MPI_Datatype mpi_complex32; -#endif - - void init (int& argc, char**& argv); - void pseudoinit (); - void finalize (); - - // Debugging output -#define CHECKPOINT dist::checkpoint(__FILE__, __LINE__) - void checkpoint (const char* file, int line); - - - - // Datatype helpers - inline MPI_Datatype datatype (const char& dummy) - { return MPI_CHAR; } - - inline MPI_Datatype datatype (const signed char& dummy) - { return MPI_UNSIGNED_CHAR; } - - inline MPI_Datatype datatype (const unsigned char& dummy) - { return MPI_BYTE; } - - inline MPI_Datatype datatype (const short& dummy) - { return MPI_SHORT; } - - inline MPI_Datatype datatype (const unsigned short& dummy) - { return MPI_UNSIGNED_SHORT; } - - inline MPI_Datatype datatype (const int& dummy) - { return MPI_INT; } - - inline MPI_Datatype datatype (const unsigned int& dummy) - { return MPI_UNSIGNED; } - - inline MPI_Datatype datatype (const long& dummy) - { return MPI_LONG; } - - inline MPI_Datatype datatype (const unsigned long& dummy) - { return MPI_UNSIGNED_LONG; } - - inline MPI_Datatype datatype (const long long& dummy) - { return MPI_LONG_LONG_INT; } +struct dist { - inline MPI_Datatype datatype (const float& dummy) - { return MPI_FLOAT; } + static const int tag = 1; - inline MPI_Datatype datatype (const double& dummy) - { return MPI_DOUBLE; } + static MPI_Comm comm; - inline MPI_Datatype datatype (const long double& dummy) - { return MPI_LONG_DOUBLE; } + static void init (int& argc, char**& argv); + static void pseudoinit (); + static void finalize (); -#if 0 + template + static MPI_Datatype datatype (const T& dummy); - inline MPI_Datatype datatype (const complex& dummy) - { return mpi_complex_float; } - - inline MPI_Datatype datatype (const complex& dummy) - { return mpi_complex_double; } - - inline MPI_Datatype datatype (const complex& dummy) - { return mpi_complex_long_double; } - -#else - -# ifdef CCTK_REAL4 - inline MPI_Datatype datatype (const CCTK_COMPLEX8& dummy) - { return mpi_complex8; } -# endif - -# ifdef CCTK_REAL8 - inline MPI_Datatype datatype (const CCTK_COMPLEX16& dummy) - { return mpi_complex16; } -# endif - -# ifdef CCTK_REAL16 - inline MPI_Datatype datatype (const CCTK_COMPLEX32& dummy) - { return mpi_complex32; } -# endif +}; + + + +template +inline MPI_Datatype dist::datatype (const T& dummy) +{ abort(); return -1; } + +template<> +inline MPI_Datatype dist::datatype (const char& dummy) +{ return MPI_CHAR; } + +template<> +inline MPI_Datatype dist::datatype (const signed char& dummy) +{ return MPI_UNSIGNED_CHAR; } -#endif +template<> +inline MPI_Datatype dist::datatype (const unsigned char& dummy) +{ return MPI_BYTE; } + +template<> +inline MPI_Datatype dist::datatype (const short& dummy) +{ return MPI_SHORT; } + +template<> +inline MPI_Datatype dist::datatype (const unsigned short& dummy) +{ return MPI_UNSIGNED_SHORT; } + +template<> +inline MPI_Datatype dist::datatype (const int& dummy) +{ return MPI_INT; } + +template<> +inline MPI_Datatype dist::datatype (const unsigned int& dummy) +{ return MPI_UNSIGNED; } + +template<> +inline MPI_Datatype dist::datatype (const long& dummy) +{ return MPI_LONG; } + +template<> +inline MPI_Datatype dist::datatype (const unsigned long& dummy) +{ return MPI_UNSIGNED_LONG; } -} // namespace dist +template<> +inline MPI_Datatype dist::datatype (const long long& dummy) +{ return MPI_LONG_LONG_INT; } +template<> +inline MPI_Datatype dist::datatype (const float& dummy) +{ return MPI_FLOAT; } +template<> +inline MPI_Datatype dist::datatype (const double& dummy) +{ return MPI_DOUBLE; } + +template<> +inline MPI_Datatype dist::datatype (const long double& dummy) +{ return MPI_LONG_DOUBLE; } + + + +#if defined(TMPL_IMPLICIT) +# include "dist.cc" +#endif #endif // DIST_HH diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index f05c88122..0123b96a1 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -1,438 +1,64 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.29 2004/04/07 16:58:07 schnetter Exp $ +/*************************************************************************** + gdata.cc - description + ------------------- + begin : Wed Jul 19 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu -#include -#include + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -#include + ***************************************************************************/ -#include "cctk.h" -#include "cctk_Parameters.h" +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ -#include "util_ErrorCodes.h" -#include "util_Table.h" +#include +#include #include "bbox.hh" #include "defs.hh" #include "dist.hh" #include "vect.hh" -#include "gdata.hh" - -using namespace std; - - - -// Communication state control -template -comm_state::comm_state () - : thestate(state_recv), - current(0) -{ -} - -template -void comm_state::step () -{ - assert (thestate!=state_done); - assert (current==tmps.size()); - thestate=astate(size_t(thestate)+1); - current=0; -} - -template -bool comm_state::done () -{ - return thestate==state_done; -} - -template -comm_state::~comm_state () -{ - assert (thestate==state_recv || thestate==state_done); -} - - - -// Hand out the next MPI tag -static int nexttag () -{ - static int last = 100; - ++last; - if (last > 30000) last = 100; - return last; -} +#if !defined(TMPL_IMPLICIT) || !defined(GDATA_HH) +# include "gdata.hh" +#endif // Constructors template -gdata::gdata (const int varindex_, const operator_type transport_operator_) - : varindex(varindex_), transport_operator(transport_operator_), - wtime_isend(0.0), wtime_isendwait(0.0), - wtime_irecv(0.0), wtime_irecvwait(0.0), - _has_storage(false), - comm_active(false), - tag(nexttag()) -{ - DECLARE_CCTK_PARAMETERS; - if (barriers) { - MPI_Barrier (dist::comm); - } -} +generic_data::generic_data () + : _has_storage(false) +{ } // Destructors template -gdata::~gdata () -{ - DECLARE_CCTK_PARAMETERS; - if (barriers) { - MPI_Barrier (dist::comm); - } -} - - - -// Data manipulators -template -void gdata::copy_from (comm_state& state, - const gdata* src, const ibbox& box) -{ - switch (state.thestate) { - case state_recv: - copy_from_recv (state, src, box); - break; - case state_send: - copy_from_send (state, src, box); - break; - case state_wait: - copy_from_wait (state, src, box); - break; - default: - assert(0); - } -} - - - -template -void gdata::copy_from_nocomm (const gdata* src, const ibbox& box) -{ - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - if (box.empty()) return; - - assert (proc() == src->proc()); - - // copy on same processor - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - copy_from_innerloop (src, box); - } -} - - - -template -void gdata::copy_from_recv (comm_state& state, - const gdata* src, const ibbox& box) -{ - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - if (box.empty()) return; - - if (proc() == src->proc()) { - // copy on same processor - - } else { - - // copy to different processor - gdata* const tmp = make_typed(varindex, transport_operator); - // TODO: is this efficient? - state.tmps.push_back (tmp); - ++state.current; - tmp->allocate (box, src->proc()); - tmp->change_processor_recv (proc()); - - } -} - - - -template -void gdata::copy_from_send (comm_state& state, - const gdata* src, const ibbox& box) -{ - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - if (box.empty()) return; - - if (proc() == src->proc()) { - // copy on same processor - - copy_from_nocomm (src, box); - - } else { - - // copy to different processor - gdata* const tmp = state.tmps.at(state.current++); - assert (tmp); - tmp->copy_from_nocomm (src, box); - tmp->change_processor_send (proc()); - - } -} - - - -template -void gdata::copy_from_wait (comm_state& state, - const gdata* src, const ibbox& box) -{ - assert (has_storage() && src->has_storage()); - assert (all(box.lower()>=extent().lower() - && box.lower()>=src->extent().lower())); - assert (all(box.upper()<=extent().upper() - && box.upper()<=src->extent().upper())); - assert (all(box.stride()==extent().stride() - && box.stride()==src->extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0 - && (box.lower()-src->extent().lower())%box.stride() == 0)); - - if (box.empty()) return; - - if (proc() == src->proc()) { - // copy on same processor - - } else { - - // copy to different processor - gdata* const tmp = state.tmps.at(state.current++); - assert (tmp); - tmp->change_processor_wait (proc()); - copy_from_nocomm (tmp, box); - delete tmp; - - } -} - - - -template -void gdata -::interpolate_from (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (transport_operator != op_error); - if (transport_operator == op_none) return; - switch (state.thestate) { - case state_recv: - interpolate_from_recv (state, srcs, times, box, time, order_space, order_time); - break; - case state_send: - interpolate_from_send (state, srcs, times, box, time, order_space, order_time); - break; - case state_wait: - interpolate_from_wait (state, srcs, times, box, time, order_space, order_time); - break; - default: - assert(0); - } -} - - - -template -void gdata -::interpolate_from_nocomm (const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (has_storage()); - assert (all(box.lower()>=extent().lower())); - assert (all(box.upper()<=extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs.at(t)->has_storage()); - assert (all(box.lower()>=srcs.at(t)->extent().lower())); - assert (all(box.upper()<=srcs.at(t)->extent().upper())); - } - - assert (! box.empty()); - if (box.empty()) return; - - assert (proc() == srcs.at(0)->proc()); - - assert (transport_operator != op_error); - assert (transport_operator != op_none); - - // interpolate on same processor - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - interpolate_from_innerloop - (srcs, times, box, time, order_space, order_time); - } -} - - - -template -void gdata -::interpolate_from_recv (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (has_storage()); - assert (all(box.lower()>=extent().lower())); - assert (all(box.upper()<=extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs.at(t)->has_storage()); - assert (all(box.lower()>=srcs.at(t)->extent().lower())); - assert (all(box.upper()<=srcs.at(t)->extent().upper())); - } - - assert (! box.empty()); - if (box.empty()) return; - - if (proc() == srcs.at(0)->proc()) { - // interpolate on same processor - - } else { - // interpolate from other processor - - gdata* const tmp = make_typed(varindex, transport_operator); - // TODO: is this efficient? - state.tmps.push_back (tmp); - ++state.current; - tmp->allocate (box, srcs.at(0)->proc()); - tmp->change_processor_recv (proc()); - - } -} +generic_data::~generic_data () { } +// Output template -void gdata -::interpolate_from_send (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (has_storage()); - assert (all(box.lower()>=extent().lower())); - assert (all(box.upper()<=extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs.at(t)->has_storage()); - assert (all(box.lower()>=srcs.at(t)->extent().lower())); - assert (all(box.upper()<=srcs.at(t)->extent().upper())); - } - - assert (! box.empty()); - if (box.empty()) return; - - if (proc() == srcs.at(0)->proc()) { - // interpolate on same processor - - interpolate_from_nocomm (srcs, times, box, time, order_space, order_time); - - } else { - // interpolate from other processor - - gdata* const tmp = state.tmps.at(state.current++); - assert (tmp); - tmp->interpolate_from_nocomm (srcs, times, box, time, order_space, order_time); - tmp->change_processor_send (proc()); - - } +ostream& operator<< (ostream& os, const generic_data& f) { + return f.out(os); } -template -void gdata -::interpolate_from_wait (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) -{ - assert (has_storage()); - assert (all(box.lower()>=extent().lower())); - assert (all(box.upper()<=extent().upper())); - assert (all(box.stride()==extent().stride())); - assert (all((box.lower()-extent().lower())%box.stride() == 0)); - assert (srcs.size() == times.size() && srcs.size()>0); - for (int t=0; t<(int)srcs.size(); ++t) { - assert (srcs.at(t)->has_storage()); - assert (all(box.lower()>=srcs.at(t)->extent().lower())); - assert (all(box.upper()<=srcs.at(t)->extent().upper())); - } - - assert (! box.empty()); - if (box.empty()) return; - - if (proc() == srcs.at(0)->proc()) { - // interpolate on same processor - - } else { - // interpolate from other processor - - gdata* const tmp = state.tmps.at(state.current++); - assert (tmp); - tmp->change_processor_wait (proc()); - copy_from_nocomm (tmp, box); - delete tmp; - - } -} - +#if defined(TMPL_EXPLICIT) +template class generic_data<1>; +template ostream& operator<< (ostream& os, const generic_data<1>& d); +template class generic_data<2>; +template ostream& operator<< (ostream& os, const generic_data<2>& d); -template class comm_state<3>; -template class gdata<3>; +template class generic_data<3>; +template ostream& operator<< (ostream& os, const generic_data<3>& d); +#endif diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh index 9355238dd..423916e3d 100644 --- a/Carpet/CarpetLib/src/gdata.hh +++ b/Carpet/CarpetLib/src/gdata.hh @@ -1,149 +1,118 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.24 2004/04/19 07:56:35 schnetter Exp $ +/*************************************************************************** + gdata.hh - description + ------------------- + begin : Wed Jul 19 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef GDATA_HH #define GDATA_HH -#include -#include - +#include +#include #include #include -#include - -#include "cctk.h" #include "defs.hh" #include "dist.hh" #include "bbox.hh" -#include "operators.hh" +#include "bboxset.hh" #include "vect.hh" -using namespace std; - +// Forward declaration +template class generic_data; +// Output template -class gdata; - - - -// State information for communications -enum astate { state_recv, state_send, state_wait, state_done }; - -template -struct comm_state { - astate thestate; - comm_state (); - void step (); - bool done (); - ~comm_state (); - - vector*> tmps; - size_t current; -}; +ostream& operator<< (ostream& os, const generic_data* f); // A generic data storage without type information template -class gdata { +class generic_data { // Types typedef vect ivect; typedef bbox ibbox; + typedef bboxset ibset; protected: // should be readonly // Fields - int varindex; // Cactus variable index, or -1 - operator_type transport_operator; - - double wtime_isend, wtime_isendwait; - double wtime_irecv, wtime_irecvwait; - bool _has_storage; // has storage associated (on some processor) - bool _owns_storage; // owns the storage - // (only valid if there is storage on this processor; it means that - // the memory is allocated and freed by this class) + ivect _shape, _stride; // shape and index order int _size; // size int _proc; // stored on processor - - ivect _shape, _stride; // shape and index order - + ibbox _extent; // bbox for all data - - bool comm_active; - MPI_Request request; - - int tag; // MPI tag for this object - + public: // Constructors - gdata (const int varindex, - const operator_type transport_operator = op_error); + generic_data (); // Destructors - virtual ~gdata (); + virtual ~generic_data (); // Pseudo constructors - virtual gdata* - make_typed (const int varindex, - const operator_type transport_operator = op_error) const = 0; - - // Processor management - virtual void change_processor (comm_state& state, - const int newproc, void* const mem=0) = 0; - protected: - virtual void change_processor_recv (const int newproc, void* const mem=0) = 0; - virtual void change_processor_send (const int newproc, void* const mem=0) = 0; - virtual void change_processor_wait (const int newproc, void* const mem=0) = 0; - public: - + virtual generic_data* make_typed (const ibbox& extent, const int proc) const + = 0; + // Storage management - virtual void transfer_from (gdata* src) = 0; - - virtual void allocate (const ibbox& extent, const int proc, - void* const mem=0) = 0; + virtual void allocate (const ibbox& extent, const int proc) = 0; virtual void free () = 0; - - // Accessors + virtual void transfer_from (generic_data* src) = 0; + + // Processor management + virtual void change_processor (const int newproc) = 0; + // Accessors bool has_storage () const { return _has_storage; } - bool owns_storage () const { - assert (_has_storage); - return _owns_storage; - } virtual const void* storage () const = 0; virtual void* storage () = 0; - int size () const { + ivect shape () const { assert (_has_storage); - return _size; + return _shape; } - int proc () const { + ivect stride () const { assert (_has_storage); - return _proc; + return _stride; } - - const ivect& shape () const { + + int size () const { assert (_has_storage); - return _shape; + return _size; } - const ivect& stride () const { + int proc () const { assert (_has_storage); - return _stride; + return _proc; } - - const ibbox& extent () const { + + ibbox extent () const { assert (_has_storage); return _extent; } @@ -156,64 +125,75 @@ public: assert (all(ind>=0 && ind<=shape())); return dot(ind, stride()); } - + // Data manipulators - public: - void copy_from (comm_state& state, - const gdata* src, const ibbox& box); - private: - void copy_from_nocomm (const gdata* src, const ibbox& box); - void copy_from_recv (comm_state& state, - const gdata* src, const ibbox& box); - void copy_from_send (comm_state& state, - const gdata* src, const ibbox& box); - void copy_from_wait (comm_state& state, - const gdata* src, const ibbox& box); - public: - void interpolate_from (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); - private: - void interpolate_from_nocomm (const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); - void interpolate_from_recv (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); - void interpolate_from_send (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); - void interpolate_from_wait (comm_state& state, - const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); - public: - + virtual void copy_from (const generic_data* src, + const ibbox& b) = 0; + virtual void interpolate_from (const generic_data* src, + const ibbox& box) = 0; + virtual void interpolate_from (const generic_data* src, const double sfact, + const generic_data* trc, const double tfact, + const ibbox& box) = 0; + + // Output + template + void write_ascii (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const + { + switch (DD) { + case 1: + write_ascii_1 (name, time, *(const vect*)&dirs, tl, rl, c, ml); + break; + case 2: + write_ascii_2 (name, time, *(const vect*)&dirs, tl, rl, c, ml); + break; + case 3: + write_ascii_3 (name, time, *(const vect*)&dirs, tl, rl, c, ml); + break; + default: + abort(); + } + } protected: - virtual void - copy_from_innerloop (const gdata* src, const ibbox& box) = 0; - virtual void - interpolate_from_innerloop (const vector srcs, - const vector times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) = 0; - + virtual void write_ascii_1 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const = 0; + virtual void write_ascii_2 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const = 0; + virtual void write_ascii_3 (const string name, const double time, + const vect& dirs, + const int tl, const int rl, + const int c, const int ml) + const = 0; +// void write_ieee (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; +// void write_hdf (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; +// void write_h5 (const string name, const double time, +// const int tl, const int rl, const int c, const int ml) +// const; +public: + + // Output + friend ostream& operator<< <>(ostream& os, const generic_data* d); + + virtual ostream& out (ostream& os) const = 0; }; +#if defined(TMPL_IMPLICIT) +# include "gdata.cc" +#endif + #endif // GDATA_HH diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc index 69537d537..dcac9fccd 100644 --- a/Carpet/CarpetLib/src/gf.cc +++ b/Carpet/CarpetLib/src/gf.cc @@ -1,47 +1,41 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.21 2004/08/07 19:47:11 schnetter Exp $ +/*************************************************************************** + gf.cc - Grid Function + data for every element of a data hierarchy + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu -#include + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ -#include "cctk.h" + ***************************************************************************/ -#include "defs.hh" +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ -#include "gf.hh" +#include -using namespace std; +#include "defs.hh" + +#if !defined(TMPL_IMPLICIT) || !defined(GF_HH) +# include "gf.hh" +#endif // Constructors template -gf::gf (const int varindex, const operator_type transport_operator, - th& t, dh& d, - const int tmin, const int tmax, const int prolongation_order_time, - const int vectorlength, const int vectorindex, - gf* const vectorleader) - : ggf(varindex, transport_operator, - t, d, tmin, tmax, prolongation_order_time, - vectorlength, vectorindex, vectorleader) +gf::gf (const string name, th& t, dh& d, + const int tmin, const int tmax) + : generic_gf(name, t, d, tmin, tmax) { - // this->recompose (); - this->recompose_crop (); - for (int rl=0; rlh.reflevels(); ++rl) { - this->recompose_allocate (rl); -#if 0 - for (comm_state state; !state.done(); state.step()) { - this->recompose_fill (state, rl, false); - } -#endif - this->recompose_free (rl); -#if 0 - for (comm_state state; !state.done(); state.step()) { - this->recompose_bnd_prolongate (state, rl, false); - } - for (comm_state state; !state.done(); state.step()) { - this->recompose_sync (state, rl, false); - } -#endif - } // for rl + recompose(); } // Destructors @@ -53,39 +47,42 @@ gf::~gf () { } // Access to the data template const data* gf::operator() (int tl, int rl, int c, int ml) const { - assert (tl>=this->tmin && tl<=this->tmax); - assert (rl>=0 && rlh.reflevels()); - assert (c>=0 && ch.components(rl)); - assert (ml>=0 && mlh.mglevels(rl,c)); - return (const data*)this->storage.at(tl-this->tmin).at(rl).at(c).at(ml); + assert (tl>=tmin && tl<=tmax); + assert (rl>=0 && rl=0 && c=0 && ml*)storage[tl-tmin][rl][c][ml]; } template data* gf::operator() (int tl, int rl, int c, int ml) { - assert (tl>=this->tmin && tl<=this->tmax); - assert (rl>=0 && rlh.reflevels()); - assert (c>=0 && ch.components(rl)); - assert (ml>=0 && mlh.mglevels(rl,c)); - return (data*)this->storage.at(tl-this->tmin).at(rl).at(c).at(ml); + assert (tl>=tmin && tl<=tmax); + assert (rl>=0 && rl=0 && c=0 && ml*)storage[tl-tmin][rl][c][ml]; } // Output template -ostream& gf::output (ostream& os) const { - T Tdummy; - os << "gf<" << typestring(Tdummy) << "," << D << ">:" - << this->varindex << "[" << CCTK_VarName(this->varindex) << "]," - << "dt=[" << this->tmin << ":" << this->tmax<< "]"; +ostream& gf::out (ostream& os) const { + os << "gf<" STR(T) "," << D << ">:\"" << name << "\"," + << "dt=[" << tmin << ":" << tmax<< "]"; return os; } +#if defined(TMPL_EXPLICIT) + #define INSTANTIATE(T) \ +template class gf; \ +template class gf; \ template class gf; #include "instantiate" - #undef INSTANTIATE + +#endif diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh index 13fc8107e..298c26589 100644 --- a/Carpet/CarpetLib/src/gf.hh +++ b/Carpet/CarpetLib/src/gf.hh @@ -1,11 +1,29 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.12 2004/03/23 12:40:27 schnetter Exp $ +/*************************************************************************** + gf.hh - Grid Function + data for every element of a data hierarchy + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef GF_HH #define GF_HH -#include -#include - +#include +#include #include #include @@ -18,13 +36,11 @@ #include "th.hh" #include "vect.hh" -using namespace std; - // A real grid function template -class gf: public ggf { +class gf: public generic_gf { // Types typedef vect ivect; @@ -42,11 +58,7 @@ class gf: public ggf { public: // Constructors - gf (const int varindex, const operator_type transport_operator, - th& t, dh& d, - const int tmin, const int tmax, const int prolongation_order_time, - const int vectorlength, const int vectorindex, - gf* const vectorleader); + gf (const string name, th& t, dh& d, const int tmin, const int tmax); // Destructors virtual ~gf (); @@ -57,14 +69,7 @@ public: protected: - virtual gdata* typed_data (int tl, int rl, int c, int ml) - { - return new data(this->varindex, this->transport_operator, - this->vectorlength, this->vectorindex, - this->vectorleader - ? (data*)(*this->vectorleader)(tl,rl,c,ml) - : NULL); - } + virtual generic_data* typed_data() { return new data; } @@ -79,9 +84,13 @@ public: // Output - virtual ostream& output (ostream& os) const; + virtual ostream& out (ostream& os) const; }; +#if defined(TMPL_IMPLICIT) +# include "gf.cc" +#endif + #endif // GF_HH diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index c26b22dbc..665017c99 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -1,59 +1,62 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.46 2004/09/17 16:37:26 schnetter Exp $ - -#include -#include -#include - +/*************************************************************************** + ggf.cc - Generic Grid Function + grid function without type information + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include #include -#include "cctk.h" - #include "defs.hh" #include "dh.hh" #include "th.hh" -#include "ggf.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) || !defined(GGF_HH) +# include "ggf.hh" +#endif // Constructors template -ggf::ggf (const int varindex, const operator_type transport_operator, - th& t, dh& d, - const int tmin, const int tmax, - const int prolongation_order_time, - const int vectorlength, const int vectorindex, - ggf* const vectorleader) - : varindex(varindex), transport_operator(transport_operator), t(t), - tmin(tmin), tmax(tmax), - prolongation_order_time(prolongation_order_time), - h(d.h), d(d), - storage(tmax-tmin+1), - vectorlength(vectorlength), vectorindex(vectorindex), - vectorleader(vectorleader) +generic_gf::generic_gf (const string name, th& t, dh& d, + const int tmin, const int tmax) + : name(name), h(d.h), t(t), d(d), tmin(tmin), tmax(tmax), + storage(tmax-tmin+1) { assert (&t.h == &d.h); - - assert (vectorlength >= 1); - assert (vectorindex >= 0 && vectorindex < vectorlength); - assert ((vectorindex==0 && !vectorleader) - || (vectorindex!=0 && vectorleader)); - + assert (tmin<=tmax+1); + d.add(this); + +// recompose(); } // Destructors template -ggf::~ggf () { +generic_gf::~generic_gf () { d.remove(this); } // Comparison template -bool ggf::operator== (const ggf& f) const { +bool generic_gf::operator== (const generic_gf& f) const { return this == &f; } @@ -61,255 +64,127 @@ bool ggf::operator== (const ggf& f) const { // Modifiers template -void ggf::recompose_crop () -{ - // Free storage that will not be needed - storage.resize(tmax-tmin+1); +void generic_gf::recompose () { + // Retain storage that might be needed + fdata oldstorage(tmax-tmin+1); for (int tl=tmin; tl<=tmax; ++tl) { - for (int rl=h.reflevels(); rl<(int)storage.at(tl-tmin).size(); ++rl) { - for (int c=0; c<(int)storage.at(tl-tmin).at(rl).size(); ++c) { - for (int ml=0; ml<(int)storage.at(tl-tmin).at(rl).at(c).size(); ++ml) { - delete storage.at(tl-tmin).at(rl).at(c).at(ml); + oldstorage[tl-tmin].resize + (min(h.reflevels(), (int)storage[tl-tmin].size())); + for (int rl=0; rl<(int)storage[tl-tmin].size(); ++rl) { + oldstorage[tl-tmin][rl].resize + (min(h.components(rl), (int)storage[tl-tmin][rl].size())); + for (int c=0; c<(int)storage[tl-tmin][rl].size(); ++c) { + oldstorage[tl-tmin][rl][c].resize + (min(h.mglevels(rl,c), (int)storage[tl-tmin][rl][c].size())); + for (int ml=0; ml<(int)storage[tl-tmin][rl][ml].size(); ++ml) { + oldstorage[tl-tmin][rl][c][ml]->transfer_from + (storage[tl-tmin][tl][c][ml]); } // for ml } // for c } // for rl - storage.at(tl-tmin).resize(h.reflevels()); } // for tl -} - -template -void ggf::recompose_allocate (const int rl) -{ - // TODO: restructure storage only when needed - // Retain storage that might be needed - oldstorage.resize(tmax-tmin+1); - for (int tl=tmin; tl<=tmax; ++tl) { - oldstorage.at(tl-tmin).resize(h.reflevels()); - assert (oldstorage.at(tl-tmin).at(rl).size() == 0); - oldstorage.at(tl-tmin).at(rl) = storage.at(tl-tmin).at(rl); - storage.at(tl-tmin).at(rl).resize(0); - } + // Delete storage + storage.clear(); - // Resize structure and allocate storage storage.resize(tmax-tmin+1); for (int tl=tmin; tl<=tmax; ++tl) { - storage.at(tl-tmin).resize(h.reflevels()); - storage.at(tl-tmin).at(rl).resize(h.components(rl)); - for (int c=0; callocate - (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c)); - } // for ml - } // for c + storage[tl-tmin].resize(h.reflevels()); + for (int rl=0; rlallocate + (d.boxes[rl][c][ml].exterior, h.proc(rl,c)); + + // Initialise from coarser level, if possible + // TODO: init only un-copied regions + if (rl>0) { + ref_prolongate (tl,rl,c,ml); + } // if rl + + // Copy from old storage, if possible + if (rl<(int)oldstorage[tl-tmin].size()) { + for (int cc=0; cc<(int)oldstorage[tl-tmin][rl].size(); ++cc) { + if (ml<(int)oldstorage[tl-tmin][rl][cc].size()) { + const ibbox ovlp = + (d.boxes[rl][c][ml].exterior + & oldstorage[tl-tmin][rl][cc][ml]->extent()); + storage[tl-tmin][rl][c][ml]->copy_from + (oldstorage[tl-tmin][rl][cc][ml], ovlp); + } // if ml + } // for cc + } // if rl + + } // for ml + } // for c + + // Delete old storage + if (rl<(int)oldstorage[tl-tmin].size()) { + oldstorage[tl-tmin][rl].clear(); + } + + } // for rl } // for tl -} - -template -void ggf::recompose_fill (comm_state& state, const int rl, - const bool do_prolongate) -{ - // Initialise the new storage - for (int c=0; cextent(); - ovlp.normalize(); - work -= ovlp; - for (typename ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) { - storage.at(tl-tmin).at(rl).at(c).at(ml)->copy_from - (state, oldstorage.at(tl-tmin).at(rl).at(cc).at(ml), *r); - } - } // if ml - } // for cc - } // if rl - - if (do_prolongate) { - // Initialise from coarser level, if possible - if (rl>0) { - if (transport_operator != op_none) { - const int numtl = prolongation_order_time+1; - assert (tmax-tmin+1 >= numtl); - vector tls(numtl); - vector times(numtl); - for (int i=0; i*> gsrcs(numtl); - for (int i=0; iextent() == gsrcs.at(0)->extent()); - } - const CCTK_REAL time = t.time(tl,rl,ml); - - // TODO: choose larger regions first - // TODO: prefer regions from the same processor - const iblist& list - = d.boxes.at(rl).at(c).at(ml).recv_ref_coarse.at(cc); - for (typename iblist::const_iterator iter=list.begin(); iter!=list.end(); ++iter) { - ibset ovlp = work & *iter; - ovlp.normalize(); - work -= ovlp; - for (typename ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) { - storage.at(tl-tmin).at(rl).at(c).at(ml)->interpolate_from - (state, gsrcs, times, *r, time, - d.prolongation_order_space, prolongation_order_time); - } // for r - } // for iter - } // for cc - } // if transport_operator - } // if rl - } // if do_prolongate - - // Note that work need not be empty here; in this case, not - // everything could be initialised. This is okay on outer - // boundaries. - // TODO: check this. - - } // for tl - } // for ml - } // for c -} - -template -void ggf::recompose_free (const int rl) -{ - // Delete old storage + for (int tl=tmin; tl<=tmax; ++tl) { - for (int c=0; c<(int)oldstorage.at(tl-tmin).at(rl).size(); ++c) { - for (int ml=0; ml<(int)oldstorage.at(tl-tmin).at(rl).at(c).size(); ++ml) { - delete oldstorage.at(tl-tmin).at(rl).at(c).at(ml); - } // for ml - } // for c - oldstorage.at(tl-tmin).at(rl).resize(0); - } // for tl -} - -template -void ggf::recompose_bnd_prolongate (comm_state& state, const int rl, - const bool do_prolongate) -{ - if (do_prolongate) { - // Set boundaries - if (rl>0) { + for (int rl=0; rl0) { + ref_bnd_prolongate (tl,rl,c,ml); + } // if rl + } // for ml } // for c - } // if rl - } // if do_prolongate -} - -template -void ggf::recompose_sync (comm_state& state, const int rl, - const bool do_prolongate) -{ - if (do_prolongate) { - // Set boundaries - for (int c=0; c -void ggf::cycle (int rl, int c, int ml) { - assert (rl>=0 && rl=0 && c=0 && ml* tmpdata = storage.at(tmin-tmin).at(rl).at(c).at(ml); - for (int tl=tmin; tl<=tmax-1; ++tl) { - storage.at(tl-tmin).at(rl).at(c).at(ml) = storage.at(tl+1-tmin).at(rl).at(c).at(ml); - } - storage.at(tmax-tmin).at(rl).at(c).at(ml) = tmpdata; -} - -// Flip the time levels by exchanging the data sets -template -void ggf::flip (int rl, int c, int ml) { - assert (rl>=0 && rl=0 && c=0 && ml* tmpdata = storage.at(tl1-tmin).at(rl).at(c).at(ml); - storage.at(tl1-tmin).at(rl).at(c).at(ml) = storage.at(tl2-tmin).at(rl).at(c).at(ml); - storage.at(tl2-tmin).at(rl).at(c).at(ml) = tmpdata; - } + + } // for rl + } // for tl } // Operations -// Copy a region +// Copy region for a component (between time levels) template -void ggf::copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const ibbox dh::dboxes::* recv_box, - int tl2, int rl2, int ml2, - const ibbox dh::dboxes::* send_box) +void generic_gf::copycat (int tl1, int rl1, int c1, int ml1, + const ibbox dh::dboxes::* recv_list, + int tl2, int rl2, int ml2, + const ibbox dh::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1=0 && c1=0 && ml1=tmin && tl2<=tmax); assert (rl2>=0 && rl2copy_from - (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), recv); + storage[tl1-tmin][rl1][c1][ml1]->copy_from + (storage[tl2-tmin][rl2][c2][ml2], recv); } -// Copy regions +// Copy regions for a component (between multigrid levels) template -void ggf::copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblist dh::dboxes::* recv_list, - int tl2, int rl2, int ml2, - const iblist dh::dboxes::* send_list) +void generic_gf::copycat (int tl1, int rl1, int c1, int ml1, + const iblist dh::dboxes::* recv_list, + int tl2, int rl2, int ml2, + const iblist dh::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1::copycat (comm_state& state, assert (rl2>=0 && rl2copy_from - (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r); + storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + (storage[tl2-tmin][rl2][c2][ml2], *r); } } -// Copy regions +// Copy regions for a level (between refinement levels) template -void ggf::copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblistvect dh::dboxes::* recv_listvect, - int tl2, int rl2, int ml2, - const iblistvect dh::dboxes::* send_listvect) +void generic_gf::copycat (int tl1, int rl1, int c1, int ml1, + const iblistvect dh::dboxes::* recv_listvect, + int tl2, int rl2, int ml2, + const iblistvect dh::dboxes::* send_listvect) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1::copycat (comm_state& state, // walk all components for (int c2=0; c2copy_from - (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r); + storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + (storage[tl2-tmin][rl2][c2][ml2], *r); } } } -// Interpolate a region +// Interpolate a component (between time levels) template -void ggf::intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const ibbox dh::dboxes::* recv_list, - const vector tl2s, int rl2, int ml2, - const ibbox dh::dboxes::* send_list, - CCTK_REAL time) +void generic_gf::intercat (int tl1, int rl1, int c1, int ml1, + const ibbox dh::dboxes::* recv_list, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const ibbox dh::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1=0 && c1=0 && ml1=tmin && tl2s.at(i)<=tmax); - } + assert (tl2>=tmin && tl2<=tmax); + assert (tl3>=tmin && tl3<=tmax); assert (rl2>=0 && rl2=0 && ml2*> gsrcs(tl2s.size()); - vector times(tl2s.size()); - for (int i=0; i<(int)gsrcs.size(); ++i) { - assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size()); - assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size()); - assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size()); - gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2); - times.at(i) = t.time(tl2s.at(i),rl2,ml2); - } - - const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list; - const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_list; - assert (all(recv.shape()==send.shape())); + assert (ml2interpolate_from - (state, gsrcs, times, recv, time, - d.prolongation_order_space, prolongation_order_time); + storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + (storage[tl2-tmin][rl2][c2][ml2], fact2, + storage[tl3-tmin][rl2][c2][ml2], fact3, recv); } -// Interpolate regions +// Interpolate a component (between multigrid levels) template -void ggf::intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblist dh::dboxes::* recv_list, - const vector tl2s, int rl2, int ml2, - const iblist dh::dboxes::* send_list, - const CCTK_REAL time) +void generic_gf::intercat (int tl1, int rl1, int c1, int ml1, + const iblist dh::dboxes::* recv_list, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const iblist dh::dboxes::* send_list) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1=0 && c1=0 && ml1=tmin && tl2s.at(i)<=tmax); - } + assert (tl2>=tmin && tl2<=tmax); + assert (tl3>=tmin && tl3<=tmax); assert (rl2>=0 && rl2=0 && ml2*> gsrcs(tl2s.size()); - vector times(tl2s.size()); - for (int i=0; i<(int)gsrcs.size(); ++i) { - assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size()); - assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size()); - assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size()); - gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2); - times.at(i) = t.time(tl2s.at(i),rl2,ml2); - } - - const iblist recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list; - const iblist send = d.boxes.at(rl2).at(c2).at(ml2).*send_list; + assert (ml2interpolate_from - (state, gsrcs, times, *r, time, - d.prolongation_order_space, prolongation_order_time); + storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + (storage[tl2-tmin][rl2][c2][ml2], fact2, + storage[tl3-tmin][rl2][c2][ml2], fact3, *r); } } -// Interpolate regions +// Interpolate a level (between refinement levels) template -void ggf::intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblistvect dh::dboxes::* recv_listvect, - const vector tl2s, int rl2, int ml2, - const iblistvect dh::dboxes::* send_listvect, - const CCTK_REAL time) +void generic_gf::intercat (int tl1, int rl1, int c1, int ml1, + const iblistvect dh::dboxes::* recv_listvect, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const iblistvect dh::dboxes::* send_listvect) { assert (tl1>=tmin && tl1<=tmax); assert (rl1>=0 && rl1=0 && c1=0 && ml1=tmin && tl2s.at(i)<=tmax); - } + assert (tl2>=tmin && tl2<=tmax); assert (rl2>=0 && rl2=tmin && tl3<=tmax); // walk all components for (int c2=0; c2=0 && ml2*> gsrcs(tl2s.size()); - vector times(tl2s.size()); - for (int i=0; i<(int)gsrcs.size(); ++i) { - assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size()); - assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size()); - assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size()); - gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2); - times.at(i) = t.time(tl2s.at(i),rl2,ml2); - } - - const iblist recv = (d.boxes.at(rl1).at(c1).at(ml1).*recv_listvect).at(c2); - const iblist send = (d.boxes.at(rl2).at(c2).at(ml2).*send_listvect).at(c1); + assert (ml2interpolate_from - (state, gsrcs, times, *r, time, - d.prolongation_order_space, prolongation_order_time); + storage[tl1-tmin][rl1][c1][ml1]->interpolate_from + (storage[tl2-tmin][rl2][c2][ml2], fact2, + storage[tl3-tmin][rl2][c2][ml2], fact3, *r); } } } @@ -497,112 +335,99 @@ void ggf::intercat (comm_state& state, // Copy a component from the next time level template -void ggf::copy (comm_state& state, int tl, int rl, int c, int ml) -{ +void generic_gf::copy (int tl, int rl, int c, int ml) { // Copy - copycat (state, - tl ,rl,c,ml, &dh::dboxes::exterior, + copycat (tl ,rl,c,ml, &dh::dboxes::exterior, tl+1,rl, ml, &dh::dboxes::exterior); } // Synchronise the boundaries a component template -void ggf::sync (comm_state& state, int tl, int rl, int c, int ml) -{ +void generic_gf::sync (int tl, int rl, int c, int ml) { // Copy - copycat (state, - tl,rl,c,ml, &dh::dboxes::recv_sync, + copycat (tl,rl,c,ml, &dh::dboxes::recv_sync, tl,rl, ml, &dh::dboxes::send_sync); } // Prolongate the boundaries of a component template -void ggf::ref_bnd_prolongate (comm_state& state, - int tl, int rl, int c, int ml, - CCTK_REAL time) -{ +void generic_gf::ref_bnd_prolongate (int tl, int rl, int c, int ml) { // Interpolate assert (rl>=1); - if (transport_operator == op_none) return; - vector tl2s; - // Interpolation in time - assert (tmax-tmin+1 >= prolongation_order_time+1); - tl2s.resize(prolongation_order_time+1); - for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i; - intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_bnd_coarse, - tl2s,rl-1, ml, &dh::dboxes::send_ref_bnd_fine, - time); + double time = + (t.time(tl,rl,ml) - t.get_time(rl-1,ml)) / (double)t.get_delta(rl-1, ml); + const int tl2 = (int)floor(time); + clog << "### ref_bnd_prolongate tl=" << tl << " rl=" << rl << " c=" << c << " ml=" << ml << " time=" << time << " tl2=" << tl2 << endl; + assert (tl2>=tmin && tl2<=tmax); + if (time==tl2) { + clog << "### (copycat)" << endl; + copycat (tl ,rl ,c,ml, &dh::dboxes::recv_ref_bnd_coarse, + tl2,rl-1, ml, &dh::dboxes::send_ref_bnd_fine); + } else { + int tl3=tl2+1; + assert (tl3>=tmin && tl3<=tmax); + const double fact2 = 1 - (time - tl2); + const double fact3 = 1 - fact2; + clog << "### (intercat) tl3=" << tl3 << " fact2=" << fact2 << " fact3=" << fact3 << endl; + intercat (tl,rl,c,ml, &dh::dboxes::recv_ref_bnd_coarse, + tl2,fact2, tl3,fact3, + rl-1,ml, &dh::dboxes::send_ref_bnd_fine); + } } // Restrict a multigrid level template -void ggf::mg_restrict (comm_state& state, - int tl, int rl, int c, int ml, - CCTK_REAL time) -{ +void generic_gf::mg_restrict (int tl, int rl, int c, int ml) { // Require same times - assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml-1)) - <= 1.0e-8 * abs(t.get_time(rl,ml))); - const vector tl2s(1,tl); - intercat (state, - tl ,rl,c,ml, &dh::dboxes::recv_mg_coarse, - tl2s,rl, ml-1, &dh::dboxes::send_mg_fine, - time); + assert (t.get_time(rl,ml) == t.get_time(rl,ml-1)); + copycat (tl,rl,c,ml, &dh::dboxes::recv_mg_coarse, + tl,rl, ml-1, &dh::dboxes::send_mg_fine); } // Prolongate a multigrid level template -void ggf::mg_prolongate (comm_state& state, - int tl, int rl, int c, int ml, - CCTK_REAL time) -{ +void generic_gf::mg_prolongate (int tl, int rl, int c, int ml) { // Require same times - assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml+1)) - <= 1.0e-8 * abs(t.get_time(rl,ml))); - const vector tl2s(1,tl); - intercat (state, - tl ,rl,c,ml, &dh::dboxes::recv_mg_coarse, - tl2s,rl, ml+1, &dh::dboxes::send_mg_fine, - time); + assert (t.get_time(rl,ml) == t.get_time(rl,ml+1)); + copycat (tl,rl,c,ml, &dh::dboxes::recv_mg_coarse, + tl,rl, ml+1, &dh::dboxes::send_mg_fine); } // Restrict a refinement level template -void ggf::ref_restrict (comm_state& state, - int tl, int rl, int c, int ml, - CCTK_REAL time) -{ +void generic_gf::ref_restrict (int tl, int rl, int c, int ml) { // Require same times - assert (abs(t.get_time(rl,ml) - t.get_time(rl+1,ml)) - <= 1.0e-8 * abs(t.get_time(rl,ml))); - if (transport_operator == op_none) return; - const vector tl2s(1,tl); - intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_fine, - tl2s,rl+1, ml, &dh::dboxes::send_ref_coarse, - time); + assert (t.get_time(rl,ml) == t.get_time(rl+1,ml)); + copycat (tl,rl ,c,ml, &dh::dboxes::recv_ref_fine, + tl,rl+1, ml, &dh::dboxes::send_ref_coarse); } // Prolongate a refinement level template -void ggf::ref_prolongate (comm_state& state, - int tl, int rl, int c, int ml, - CCTK_REAL time) -{ - assert (rl>=1); - if (transport_operator == op_none) return; - vector tl2s; - // Interpolation in time - assert (tmax-tmin+1 >= prolongation_order_time+1); - tl2s.resize(prolongation_order_time+1); - for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i; - intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_coarse, - tl2s,rl-1, ml, &dh::dboxes::send_ref_fine, - time); +void generic_gf::ref_prolongate (int tl, int rl, int c, int ml) { + // Require same times + assert (t.get_time(rl,ml) == t.get_time(rl-1,ml)); + copycat (tl,rl ,c,ml, &dh::dboxes::recv_ref_coarse, + tl,rl-1, ml, &dh::dboxes::send_ref_fine); } -template class ggf<3>; +// Output +template +ostream& operator<< (ostream& os, const generic_gf& f) { + return f.out(os); +} + + + +#if defined(TMPL_EXPLICIT) +template class generic_gf<1>; +template ostream& operator<< (ostream& os, const generic_gf<1>& f); + +template class generic_gf<2>; +template ostream& operator<< (ostream& os, const generic_gf<2>& f); + +template class generic_gf<3>; +template ostream& operator<< (ostream& os, const generic_gf<3>& f); +#endif diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index b6a551cf3..1501aa5dd 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -1,15 +1,30 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.25 2004/08/07 19:47:11 schnetter Exp $ +/*************************************************************************** + ggf.hh - Generic Grid Function + grid function without type information + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef GGF_HH #define GGF_HH -#include - +#include #include #include -#include - -#include "cctk.h" #include "defs.hh" #include "dh.hh" @@ -17,93 +32,63 @@ #include "gh.hh" #include "th.hh" -using namespace std; - // Forward declaration -template class ggf; +template class generic_gf; // Output template -ostream& operator<< (ostream& os, const ggf& f); +ostream& operator<< (ostream& os, const generic_gf& f); // A generic grid function without type information template -class ggf { +class generic_gf { // Types - typedef vect ivect; typedef bbox ibbox; typedef bboxset ibset; typedef list iblist; typedef vector iblistvect; - typedef gdata* tdata; // data ... - typedef vector mdata; // ... for each multigrid level - typedef vector cdata; // ... for each component - typedef vector rdata; // ... for each refinement level - typedef vector fdata; // ... for each time level - + typedef generic_data* tdata; // data ... + typedef vector mdata; // ... for each multigrid level + typedef vector cdata; // ... for each component + typedef vector rdata; // ... for each refinement level + typedef vector fdata; // ... for each time level + public: // should be readonly - + // Fields - int varindex; // Cactus variable index - operator_type transport_operator; - - th &t; // time hierarchy - int tmin, tmax; // timelevels - int prolongation_order_time; // order of temporal prolongation operator - + string name; + gh &h; // grid hierarchy + th &t; // time hierarchy dh &d; // data hierarchy + int tmin, tmax; // timelevels protected: fdata storage; // storage - -public: - int vectorlength; // vector length - int vectorindex; // index of *this - ggf* vectorleader; // first vector element - -private: - fdata oldstorage; - + public: // Constructors - ggf (const int varindex, const operator_type transport_operator, - th& t, dh& d, - const int tmin, const int tmax, - const int prolongation_order_time, - const int vectorlength, const int vectorindex, - ggf* const vectorleader); + generic_gf (const string name, th& t, dh& d, + const int tmin, const int tmax); // Destructors - virtual ~ggf (); + virtual ~generic_gf (); // Comparison - bool operator== (const ggf& f) const; + virtual bool operator== (const generic_gf& f) const; // Modifiers - // void recompose (); - void recompose_crop (); - void recompose_allocate (int rl); - void recompose_fill (comm_state& state, int rl, bool do_prolongate); - void recompose_free (int rl); - void recompose_bnd_prolongate (comm_state& state, int rl, bool do_prolongate); - void recompose_sync (comm_state& state, int rl, bool do_prolongate); - - // Cycle the time levels by rotating the data sets - void cycle (int rl, int c, int ml); - - // Flip the time levels by exchanging the data sets - void flip (int rl, int c, int ml); + virtual void recompose (); @@ -111,7 +96,7 @@ public: protected: - virtual gdata* typed_data (int tl, int rl, int c, int ml) = 0; + virtual generic_data* typed_data() = 0; @@ -119,50 +104,47 @@ protected: protected: - // Copy a region - void copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const ibbox dh::dboxes::* recv_list, - int tl2, int rl2, int ml2, - const ibbox dh::dboxes::* send_list); - - // Copy regions - void copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblist dh::dboxes::* recv_list, - int tl2, int rl2, int ml2, - const iblist dh::dboxes::* send_list); - - // Copy regions - void copycat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblistvect dh::dboxes::* recv_listvect, - int tl2, int rl2, int ml2, - const iblistvect dh::dboxes::* send_listvect); - - // Interpolate a region - void intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const ibbox dh::dboxes::* recv_list, - const vector tl2s, int rl2, int ml2, - const ibbox dh::dboxes::* send_list, - CCTK_REAL time); - - // Interpolate regions - void intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblist dh::dboxes::* recv_list, - const vector tl2s, int rl2, int ml2, - const iblist dh::dboxes::* send_list, - CCTK_REAL time); - - // Interpolate regions - void intercat (comm_state& state, - int tl1, int rl1, int c1, int ml1, - const iblistvect dh::dboxes::* recv_listvect, - const vector tl2s, int rl2, int ml2, - const iblistvect dh::dboxes::* send_listvect, - CCTK_REAL time); + // Copy region for a component (between time levels) + virtual void copycat (int tl1, int rl1, int c1, int ml1, + const ibbox dh::dboxes::* recv_list, + int tl2, int rl2, int ml2, + const ibbox dh::dboxes::* send_list); + + // Copy regions for a component (between multigrid levels) + virtual void copycat (int tl1, int rl1, int c1, int ml1, + const iblist dh::dboxes::* recv_list, + int tl2, int rl2, int ml2, + const iblist dh::dboxes::* send_list); + + // Copy regions for a level (between refinement levels) + virtual void copycat (int tl1, int rl1, int c1, int ml1, + const iblistvect dh::dboxes::* recv_listvect, + int tl2, int rl2, int ml2, + const iblistvect dh::dboxes::* send_listvect); + + // Interpolate a component (between time levels) + virtual void intercat (int tl1, int rl1, int c1, int ml1, + const ibbox dh::dboxes::* recv_list, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const ibbox dh::dboxes::* send_list); + + // Interpolate a component (between multigrid levels) + virtual void intercat (int tl1, int rl1, int c1, int ml1, + const iblist dh::dboxes::* recv_list, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const iblist dh::dboxes::* send_list); + + // Interpolate a level (between refinement levels) + virtual void intercat (int tl1, int rl1, int c1, int ml1, + const iblistvect dh::dboxes::* recv_listvect, + int tl2, const double fact2, + int tl3, const double fact3, + int rl2, int ml2, + const iblistvect dh::dboxes::* send_listvect); @@ -171,55 +153,47 @@ public: // The grid boundaries have to be updated after calling mg_restrict, // mg_prolongate, ref_restrict, or ref_prolongate. - // "Updating" means here that the boundaries have to be - // synchronised. They don't need to be prolongated. - // Copy a component from the next time level - void copy (comm_state& state, int tl, int rl, int c, int ml); + virtual void copy (int tl, int rl, int c, int ml); // Synchronise the boundaries of a component - void sync (comm_state& state, int tl, int rl, int c, int ml); + virtual void sync (int tl, int rl, int c, int ml); // Prolongate the boundaries of a component - void ref_bnd_prolongate (comm_state& state, - int tl, int rl, int c, int ml, CCTK_REAL time); + virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml); // Restrict a multigrid level - void mg_restrict (comm_state& state, - int tl, int rl, int c, int ml, CCTK_REAL time); + virtual void mg_restrict (int tl, int rl, int c, int ml); // Prolongate a multigrid level - void mg_prolongate (comm_state& state, - int tl, int rl, int c, int ml, CCTK_REAL time); + virtual void mg_prolongate (int tl, int rl, int c, int ml); // Restrict a refinement level - void ref_restrict (comm_state& state, - int tl, int rl, int c, int ml, CCTK_REAL time); + virtual void ref_restrict (int tl, int rl, int c, int ml); // Prolongate a refinement level - void ref_prolongate (comm_state& state, - int tl, int rl, int c, int ml, CCTK_REAL time); + virtual void ref_prolongate (int tl, int rl, int c, int ml); // Access to the data - virtual const gdata* operator() (int tl, int rl, int c, int ml) const = 0; + virtual const generic_data* operator() (int tl, int rl, int c, int ml) + const = 0; - virtual gdata* operator() (int tl, int rl, int c, int ml) = 0; + virtual generic_data* operator() (int tl, int rl, int c, int ml) = 0; // Output - virtual ostream& output (ostream& os) const = 0; + friend ostream& operator<< <> (ostream& os, const generic_gf& f); + + virtual ostream& out (ostream& os) const = 0; }; -template -inline ostream& operator<< (ostream& os, const ggf& f) { - return f.output(os); -} - - +#if defined(TMPL_IMPLICIT) +# include "ggf.cc" +#endif #endif // GGF_HH diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index 32df39bf0..5002f5199 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -1,19 +1,36 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.29 2004/08/07 19:47:11 schnetter Exp $ - -#include -#include +/*************************************************************************** + gh.cc - Grid Hierarchy + bounding boxes for each multigrid level of each + component of each refinement level + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include +#include #include -#include "cctk.h" -#include "cctk_Parameters.h" - #include "defs.hh" #include "dh.hh" #include "th.hh" -#include "gh.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) || !defined(GH_HH) +# include "gh.hh" +#endif @@ -21,11 +38,15 @@ using namespace std; template gh::gh (const int reffact, const centering refcent, const int mgfact, const centering mgcent, - const ibbox baseextent) + const ibbox& baseextent) : reffact(reffact), refcent(refcent), mgfact(mgfact), mgcent(mgcent), baseextent(baseextent) { + assert (reffact>=1); + assert (mgfact>=1); + assert (refcent==vertex_centered || refcent==cell_centered); + assert (mgcent==vertex_centered || mgcent==cell_centered); } // Destructors @@ -34,15 +55,8 @@ gh::~gh () { } // Modifiers template -void gh::recompose (const rexts& exts, - const rbnds& outer_bounds, - const rprocs& procs, - const bool do_prolongate) -{ - DECLARE_CCTK_PARAMETERS; - +void gh::recompose (const rexts& exts, const rprocs& procs) { extents = exts; - outer_boundaries = outer_bounds; processors = procs; // Consistency checks @@ -52,10 +66,8 @@ void gh::recompose (const rexts& exts, // Check processor number consistency for (int rl=0; rl::recompose (const rexts& exts, for (int c=0; c0); for (int ml=1; ml::recompose (const rexts& exts, assert (components(rl)>0); for (int c=0; c::recompose (const rexts& exts, // Check base grid extent if (reflevels()>0) { for (int c=0; c all; for (int c=0; c:" - << "reffactor=" << reffact << ",refcentering=" << refcent << "," - << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << "," - << "extents=" << extents << "," - << "outer_boundaries=" << outer_boundaries << "," - << "processors=" << processors << "," + << "reffactor=" << h.reffact << ",refcentering=" << h.refcent << "," + << "mgfactor=" << h.mgfact << ",mgcentering=" << h.mgcent << "," + << "baseextent=" << h.baseextent << "," + << "extents=" << h.extents << "," << "dhs={"; int cnt=0; - for (typename list*>::const_iterator d = dhs.begin(); - d != dhs.end(); ++d) { + for (list*>::const_iterator d = h.dhs.begin(); + d != h.dhs.end(); ++d) { if (cnt++) os << ","; - (*d)->output(os); + os << **d; } os << "}"; return os; @@ -237,4 +229,13 @@ ostream& gh::output (ostream& os) const { +#if defined(TMPL_EXPLICIT) +template class gh<1>; +template ostream& operator<< (ostream& os, const gh<1>& h); + +template class gh<2>; +template ostream& operator<< (ostream& os, const gh<2>& h); + template class gh<3>; +template ostream& operator<< (ostream& os, const gh<3>& h); +#endif diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh index edf149000..1d185bee5 100644 --- a/Carpet/CarpetLib/src/gh.hh +++ b/Carpet/CarpetLib/src/gh.hh @@ -1,10 +1,29 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.18 2004/08/07 19:47:11 schnetter Exp $ +/*************************************************************************** + gh.hh - Grid Hierarchy + bounding boxes for each multigrid level of each + component of each refinement level + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef GH_HH #define GH_HH -#include - +#include #include #include #include @@ -14,14 +33,12 @@ #include "dist.hh" #include "vect.hh" -using namespace std; - // Forward declaration -template class dh; -template class th; template class gh; +template class th; +template class dh; // Output template @@ -34,23 +51,16 @@ ostream& operator<< (ostream& os, const gh& h); template class gh { -public: - // Types typedef vect ivect; typedef bbox ibbox; - typedef vect,D> bvect; - - typedef vector mexts; // ... for each multigrid level - typedef vector cexts; // ... for each component - typedef vector rexts; // ... for each refinement level + typedef vector mexts; // ... for each multigrid level + typedef vector cexts; // ... for each component + typedef vector rexts; // ... for each refinement level - typedef vector cbnds; // ... for each component - typedef vector rbnds; // ... for each refinement level - - typedef vector cprocs; // ... for each component - typedef vector rprocs; // ... for each refinement level + typedef vector cprocs; // ... for each component + typedef vector rprocs; // ... for each refinement level public: // should be readonly @@ -61,16 +71,12 @@ public: // should be readonly int mgfact; // default multigrid factor centering mgcent; // default (vertex or cell centered) - list*> ths; // list of all time hierarchies - - ibbox baseextent; - vector > bases; // [rl][ml] + ibbox baseextent; // bounds (inclusive) of base level - // TODO: invent structure for this - rexts extents; // extents of all grids - rbnds outer_boundaries; // boundary descriptions of all grids + rexts extents; // bounds of all grids rprocs processors; // processor numbers of all grids + list*> ths; // list of all time hierarchies list*> dhs; // list of all data hierarchies public: @@ -78,15 +84,18 @@ public: // Constructors gh (const int reffact, const centering refcent, const int mgfact, const centering mgcent, - const ibbox baseextent); + const ibbox& baseextent); // Destructors - virtual ~gh (); + ~gh (); // Modifiers - void recompose (const rexts& exts, const rbnds& outer_bounds, - const rprocs& procs, - const bool do_prolongate); + void recompose (const rexts& exts, const rprocs& procs); + + // Helpers + rexts make_multigrid_boxes (const vector >& exts, + const int mglevels) + const; // Accessors int reflevels () const { @@ -94,19 +103,20 @@ public: } int components (const int rl) const { - return (int)extents.at(rl).size(); + assert (rl>=0 && rl=0 && rl=0 && c=0 && rl=0 && c* t); void remove (th* t); @@ -126,17 +134,13 @@ public: void remove (dh* d); // Output - virtual ostream& output (ostream& os) const; + friend ostream& operator<< <> (ostream& os, const gh& h); }; -template -inline ostream& operator<< (ostream& os, const gh& h) { - h.output(os); - return os; -} - - +#if defined(TMPL_IMPLICIT) +# include "gh.cc" +#endif #endif // GH_HH diff --git a/Carpet/CarpetLib/src/instantiate b/Carpet/CarpetLib/src/instantiate index 1f4e2d6d7..c87c16eeb 100644 --- a/Carpet/CarpetLib/src/instantiate +++ b/Carpet/CarpetLib/src/instantiate @@ -1,178 +1,33 @@ // Instantiate templates for all available types -*-C++-*- // (C) 2001 Erik Schnetter -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.10 2004/03/01 21:35:13 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.1 2001/03/01 13:40:10 eschnett Exp $ // Usage: // Define the macro INSTANTIATE(T) to instantiate for the type T, // then include this file, // then undefine the macro INSTANTIATE. +#include +#include "cctk.h" -// Decide which types to instantiate - -#ifdef CARPET_ALL -# undef CARPET_BYTE -# undef CARPET_INT -# undef CARPET_REAL -# undef CARPET_COMPLEX -# define CARPET_BYTE -# define CARPET_INT -# define CARPET_REAL -# define CARPET_COMPLEX -#endif - -#ifdef CARPET_ALL_INT -# undef CARPET_INT1 -# undef CARPET_INT2 -# undef CARPET_INT4 -# undef CARPET_INT8 -# define CARPET_INT1 -# define CARPET_INT2 -# define CARPET_INT4 -# define CARPET_INT8 -#endif - -#ifdef CARPET_ALL_REAL -# undef CARPET_REAL4 -# undef CARPET_REAL8 -# undef CARPET_REAL16 -# define CARPET_REAL4 -# define CARPET_REAL8 -# define CARPET_REAL16 -#endif - -#ifdef CARPET_ALL_COMPLEX -# undef CARPET_COMPLEX8 -# undef CARPET_COMPLEX16 -# undef CARPET_COMPLEX32 -# define CARPET_COMPLEX8 -# define CARPET_COMPLEX16 -# define CARPET_COMPLEX32 -#endif - -#if !defined(CARPET_BYTE) && !defined(CARPET_INT) && !defined(CARPET_INT1) && !defined(CARPET_INT2) && !defined(CARPET_INT4) && !defined(CARPET_INT8) && !defined(CARPET_REAL) && !defined(CARPET_REAL4) && !defined(CARPET_REAL8) && !defined(CARPET_REAL16) && !defined(CARPET_COMPLEX) && !defined(CARPET_COMPLEX8) && !defined(CARPET_COMPLEX16) && !defined(CARPET_COMPLEX32) -// Assume the user just wants INT, REAL, and COMPLEX -# undef CARPET_INT -# define CARPET_INT -# undef CARPET_REAL -# define CARPET_REAL -# undef CARPET_COMPLEX -# define CARPET_COMPLEX -#endif - -#ifdef CARPET_INT -# ifdef CCTK_INTEGER_PRECISION_1 -# undef CARPET_INT1 -# define CARPET_INT1 -# endif -# ifdef CCTK_INTEGER_PRECISION_2 -# undef CARPET_INT2 -# define CARPET_INT2 -# endif -# ifdef CCTK_INTEGER_PRECISION_4 -# undef CARPET_INT4 -# define CARPET_INT4 -# endif -# ifdef CCTK_INTEGER_PRECISION_8 -# undef CARPET_INT8 -# define CARPET_INT8 -# endif -#endif -#ifdef CARPET_REAL -# ifdef CCTK_REAL_PRECISION_4 -# undef CARPET_REAL4 -# define CARPET_REAL4 -# endif -# ifdef CCTK_REAL_PRECISION_8 -# undef CARPET_REAL8 -# define CARPET_REAL8 -# endif -# ifdef CCTK_REAL_PRECISION_16 -# undef CARPET_REAL16 -# define CARPET_REAL16 -# endif -#endif -#ifdef CARPET_COMPLEX -# ifdef CCTK_REAL_PRECISION_4 -# undef CARPET_COMPLEX8 -# define CARPET_COMPLEX8 -# endif -# ifdef CCTK_REAL_PRECISION_8 -# undef CARPET_COMPLEX16 -# define CARPET_COMPLEX16 -# endif -# ifdef CCTK_REAL_PRECISION_16 -# undef CARPET_COMPLEX32 -# define CARPET_COMPLEX32 -# endif -#endif - - - -// // Check -// #if !defined(CARPET_BYTE) && !defined(CARPET_INT1) && !defined(CARPET_INT2) && !defined(CARPET_INT4) && !defined(CARPET_INT8) && !defined(CARPET_REAL4) && !defined(CARPET_REAL8) && !defined(CARPET_REAL16) && !defined(CARPET_COMPLEX8) && !defined(CARPET_COMPLEX16) && !defined(CARPET_COMPLEX32) -// # error "You have not defined which grid function types to instantiate." -// #endif - - - -// Instantiate the desired types - -#ifdef CARPET_BYTE INSTANTIATE(CCTK_BYTE) -#endif -#ifdef CARPET_INT1 -# ifdef CCTK_INT1 -INSTANTIATE(CCTK_INT1) -# endif -#endif -#ifdef CARPET_INT2 -# ifdef CCTK_INT2 INSTANTIATE(CCTK_INT2) -# endif -#endif -#ifdef CARPET_INT4 -# ifdef CCTK_INT4 INSTANTIATE(CCTK_INT4) -# endif -#endif -#ifdef CARPET_INT8 -# ifdef CCTK_INT8 +#ifdef CCTK_INT8 INSTANTIATE(CCTK_INT8) -# endif #endif -#ifdef CARPET_REAL4 -# ifdef CCTK_REAL4 INSTANTIATE(CCTK_REAL4) -# endif -#endif -#ifdef CARPET_REAL8 -# ifdef CCTK_REAL8 INSTANTIATE(CCTK_REAL8) -# endif -#endif -#ifdef CARPET_REAL16 -# ifdef CCTK_REAL16 +#ifdef CCTK_REAL16 INSTANTIATE(CCTK_REAL16) -# endif #endif -#ifdef CARPET_COMPLEX8 -# ifdef CCTK_REAL4 -INSTANTIATE(CCTK_COMPLEX8) -# endif -#endif -#ifdef CARPET_COMPLEX16 -# ifdef CCTK_REAL8 -INSTANTIATE(CCTK_COMPLEX16) -# endif -#endif -#ifdef CARPET_COMPLEX32 -# ifdef CCTK_REAL16 -INSTANTIATE(CCTK_COMPLEX32) -# endif +INSTANTIATE(complex) +INSTANTIATE(complex) +#ifdef CCTK_REAL16 +INSTANTIATE(complex) #endif diff --git a/Carpet/CarpetLib/src/io.cc b/Carpet/CarpetLib/src/io.cc new file mode 100644 index 000000000..099fd0890 --- /dev/null +++ b/Carpet/CarpetLib/src/io.cc @@ -0,0 +1,44 @@ +/*************************************************************************** + io.cc - I/O routines + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/io.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include "io.hh" + +#if !defined(TMPL_IMPLICIT) || !defined(GF_HH) +# include "io.hh" +#endif + +template +IObase::DataType datatype(TT dummy) +{ assert(false); return IObase::DataType(); } + +template <> +IObase::DataType datatype(char* dummy) { return IObase::uInt8; } +template <> +IObase::DataType datatype(short* dummy) { return IObase::Int16; } +template <> +IObase::DataType datatype(int* dummy) { return IObase::Int32; } +template <> +IObase::DataType datatype(long long* dummy) { return IObase::Int64; } +template <> +IObase::DataType datatype(float* dummy) { return IObase::Float32; } +template <> +IObase::DataType datatype(double* dummy) { return IObase::Float64; } +template <> +IObase::DataType datatype(long double* dummy) { return IObase::Float64; } diff --git a/Carpet/CarpetLib/src/io.hh b/Carpet/CarpetLib/src/io.hh new file mode 100644 index 000000000..624f3fefe --- /dev/null +++ b/Carpet/CarpetLib/src/io.hh @@ -0,0 +1,61 @@ +/*************************************************************************** + io.hh - I/O routines + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/io.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#ifndef IO_HH +#define IO_HH + +#include +#include + +#include +#include +#include +#include +#include + +#include "data.hh" +#include "defs.hh" + +template +class io { + + enum iotype { ieee, hdf, h5 }; + +// template +// void write_ascii (const data& d, const string name, +// const vect& dirs); + +// static void write_ASCII_1D (const data& d, const string name, +// const int dir); +// static void write_ASCII_2D (const data& d, const string name, +// const int dir1, const int dir2); +// static void write_ASCII_3D (const data& d, const string name); + +// static void write_FlexIO (const data& d, const string name, +// const iotype iot); +}; + + + +#if defined(TMPL_IMPLICIT) +# include "io.cc" +#endif + +#endif // IO_HH diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index 4bfcaec4f..1ff057287 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -1,46 +1,8 @@ # Main make.code.defn file for thorn CarpetLib -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.13 2004/03/03 15:30:40 hawke Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.1 2001/03/01 13:40:10 eschnett Exp $ # Source files in this directory -SRCS = bbox.cc \ - bboxset.cc \ - data.cc \ - defs.cc \ - dh.cc \ - dist.cc \ - gdata.cc \ - gf.cc \ - ggf.cc \ - gh.cc \ - th.cc \ - vect.cc \ - checkindex.F77 \ - copy_3d_complex16.F77 \ - copy_3d_int4.F77 \ - copy_3d_real8.F77 \ - prolongate_3d_real8.F77 \ - prolongate_3d_real8_rf2.F77 \ - prolongate_3d_real8_o3.F77 \ - prolongate_3d_real8_o3_rf2.F77 \ - prolongate_3d_real8_o5.F77 \ - prolongate_3d_real8_2tl.F77 \ - prolongate_3d_real8_2tl_rf2.F77 \ - prolongate_3d_real8_2tl_o3.F77 \ - prolongate_3d_real8_2tl_o3_rf2.F77 \ - prolongate_3d_real8_2tl_o5.F77 \ - prolongate_3d_real8_3tl.F77 \ - prolongate_3d_real8_3tl_rf2.F77 \ - prolongate_3d_real8_3tl_o3.F77 \ - prolongate_3d_real8_3tl_o3_rf2.F77 \ - prolongate_3d_real8_3tl_o5.F77 \ - prolongate_3d_real8_minmod.F77 \ - prolongate_3d_real8_2tl_minmod.F77 \ - prolongate_3d_real8_3tl_minmod.F77 \ - prolongate_3d_real8_eno.F90 \ - prolongate_3d_real8_2tl_eno.F90 \ - prolongate_3d_real8_3tl_eno.F90 \ - restrict_3d_real8.F77 \ - restrict_3d_real8_rf2.F77 +SRCS = bbox.cc bboxset.cc data.cc defs.cc dh.cc dist.cc gdata.cc gf.cc ggf.cc gh.cc io.cc th.cc vect.cc # Subdirectories containing source files SUBDIRS = diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc index 0ad9beca1..1cf81495a 100644 --- a/Carpet/CarpetLib/src/th.cc +++ b/Carpet/CarpetLib/src/th.cc @@ -1,25 +1,39 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.14 2004/03/23 19:30:14 schnetter Exp $ - -#include -#include - +/*************************************************************************** + th.cc - Time Hierarchy + information about time levels + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include -#include "cctk.h" - #include "defs.hh" #include "gh.hh" -#include "th.hh" - -using namespace std; +#if !defined(TMPL_IMPLICIT) || !defined(TH_HH) +# include "th.hh" +#endif // Constructors template -th::th (gh& h, const CCTK_REAL basedelta) - : h(h), delta(basedelta) { +th::th (gh& h, const int basedelta) : h(h), delta(basedelta) { h.add(this); } @@ -35,25 +49,26 @@ void th::recompose () { times.resize(h.reflevels()); deltas.resize(h.reflevels()); for (int rl=0; rl::recompose () { // Output template -void th::output (ostream& os) const { +ostream& operator<< (ostream& os, const th& t) { os << "th<" << D << ">:" << "times={"; - for (int rl=0; rl; +template ostream& operator<< (ostream& os, const th<1>& t); + +template class th<2>; +template ostream& operator<< (ostream& os, const th<2>& t); + template class th<3>; +template ostream& operator<< (ostream& os, const th<3>& t); +#endif diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh index 12c77d782..6089e37cd 100644 --- a/Carpet/CarpetLib/src/th.hh +++ b/Carpet/CarpetLib/src/th.hh @@ -1,20 +1,34 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.11 2004/03/23 19:30:14 schnetter Exp $ +/*************************************************************************** + th.hh - Time Hierarchy + information about time levels + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef TH_HH #define TH_HH -#include - +#include #include #include -#include "cctk.h" - #include "defs.hh" #include "gh.hh" -using namespace std; - // Forward declaration @@ -33,18 +47,18 @@ class th { public: // should be readonly // Fields - gh& h; // hierarchy + gh &h; // hierarchy private: - CCTK_REAL delta; // time step - vector > times; // current times - vector > deltas; // time steps + int delta; // time step + vector > times; // current times + vector > deltas; // time steps public: // Constructors - th (gh& h, const CCTK_REAL basedelta); + th (gh& h, const int basedelta); // Destructors ~th (); @@ -53,52 +67,42 @@ public: void recompose (); // Time management - CCTK_REAL get_time (const int rl, const int ml) const { + int get_time (const int rl, const int ml) const { assert (rl>=0 && rl=0 && ml=0 && rl=0 && ml=0 && rl=0 && ml=0 && rl=0 && ml=0 && rl=0 && ml (ostream& os, const th& d); }; -template -inline ostream& operator<< (ostream& os, const th& t) { - t.output(os); - return os; -} - - +#if defined(TMPL_IMPLICIT) +# include "th.cc" +#endif #endif // TH_HH diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc index 9af8bca1a..d9fd90594 100644 --- a/Carpet/CarpetLib/src/vect.cc +++ b/Carpet/CarpetLib/src/vect.cc @@ -1,59 +1,61 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.14 2004/02/18 15:10:17 schnetter Exp $ - -#include - +/*************************************************************************** + vect.cc - Small inline vectors + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include #include #include "defs.hh" -#include "vect.hh" - -using namespace std; - - - -// Input -template -void vect::input (istream& is) { - skipws (is); - consume (is, '['); - for (int d=0; d> (*this)[d]; - if (d -void vect::output (ostream& os) const { +ostream& operator<< (ostream& os, const vect& a) { os << "["; for (int d=0; d0) os << ","; + os << a[d]; } os << "]"; + return os; } -// Note: We need all dimensions all the time. -template class vect; +#if defined(TMPL_EXPLICIT) template class vect; -template class vect; -template class vect; +// template class vect; +template ostream& operator<< (ostream& os, const vect& a); +// template ostream& operator<< (ostream& os, const vect& a); -template void vect::input (istream& is); -template void vect,3>::input (istream& is); +template class vect; +// template class vect; +template ostream& operator<< (ostream& os, const vect& a); +// template ostream& operator<< (ostream& os, const vect& a); -template void vect::output (ostream& os) const; -template void vect::output (ostream& os) const; -template void vect::output (ostream& os) const; -template void vect,3>::output (ostream& os) const; -template void vect,3>::output (ostream& os) const; +template class vect; +// template class vect; +template ostream& operator<< (ostream& os, const vect& a); +// template ostream& operator<< (ostream& os, const vect& a); +#endif diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index 9160c76aa..7c99d2a39 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.hh @@ -1,310 +1,222 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.29 2004/08/14 07:41:25 schnetter Exp $ +/*************************************************************************** + vect.hh - Small inline vectors + ------------------- + begin : Sun Jun 11 2000 + copyright : (C) 2000 by Erik Schnetter + email : schnetter@astro.psu.edu + + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ #ifndef VECT_HH #define VECT_HH -#include #include #include #include -using namespace std; - // Forward definition template class vect; -// Input/Output -template -istream& operator>> (istream& is, vect& a); +// Output template ostream& operator<< (ostream& os, const vect& a); -/** - * A short vector with a size that is specified at compile time. - */ +// Vect class template class vect { // Fields - - /** Vector elements. */ - T elt[D==0 ? 1 : D]; + T elt[D]; public: // Constructors - - /** Explicit empty constructor. */ explicit vect () { } - /** Copy constructor. */ vect (const vect& a) { for (int d=0; d - /*explicit*/ vect (const vect& a) { + explicit vect (const vect& a) { for (int d=0; d=0 && d=0 && d vect operator[] (const vect& a) const { vect r; - // (*this)[] performs index checking for (int d=0; d=0 && d operator! () const { - vect r; - for (int d=0; d operator&& (const T x) const { - vect r; - for (int d=0; d operator|| (const T x) const { - vect r; - for (int d=0; d operator&& (const vect& a) const { vect r; - for (int d=0; d operator|| (const vect& a) const { vect r; - for (int d=0; d operator== (const vect& a) const { vect r; - for (int d=0; d operator!= (const vect& a) const { vect r; - for (int d=0; d operator< (const vect& a) const { vect r; - for (int d=0; d operator<= (const vect& a) const { vect r; - for (int d=0; d operator> (const vect& a) const { vect r; - for (int d=0; da[d]; + for (int d=0; da[d]; return r; } vect operator>= (const vect& a) const { vect r; - for (int d=0; d=a[d]; + for (int d=0; d=a[d]; return r; } - /** This corresponds to the ?: operator. Return a vector with the - elements set to either a[i] or b[i], depending on whether - (*this)[i] is true or not. */ - template - vect ifthen (const vect& a, const vect& b) const { - vect r; - for (int d=0; d(ostream& os, const vect& a); }; // Operators - -/** This corresponds to the ?: operator. Return a vector with the - elements set to either b[i] or c[i], depending on whether a[i] is - true or not. */ -template -inline vect either (const vect& a, - const vect& b, const vect& c) { - vect r; - for (int d=0; d -inline vect,DD> xpose (vect,D> const & a) { - vect,DD> r; - for (int dd=0; dd inline vect abs (const vect& a) { vect r; @@ -524,23 +412,6 @@ inline vect abs (const vect& a) { return r; } -/** Return the element-wise ceiling. */ -template -inline vect ceil (const vect& a) { - vect r; - for (int d=0; d -inline vect floor (const vect& a) { - vect r; - for (int d=0; d inline vect max (const vect& a, const vect& b) { vect r; @@ -548,7 +419,6 @@ inline vect max (const vect& a, const vect& b) { return r; } -/** Return the element-wise minimum of two vectors. */ template inline vect min (const vect& a, const vect& b) { vect r; @@ -556,19 +426,9 @@ inline vect min (const vect& a, const vect& b) { return r; } -/** Return the element-wise power of two vectors. */ -template -inline vect pow (const vect& a, const vect& b) { - vect r; - for (int d=0; d inline bool any (const vect& a) { bool r(false); @@ -576,7 +436,6 @@ inline bool any (const vect& a) { return r; } -/** Return true iff all of the elements are true (boolean product). */ template inline bool all (const vect& a) { bool r(true); @@ -584,13 +443,11 @@ inline bool all (const vect& a) { return r; } -/** Count the number of elements in the vector. */ template inline int count (const vect& a) { return D; } -/** Return the dot product of two vectors. */ template inline T dot (const vect& a, const vect& b) { T r(0); @@ -598,49 +455,27 @@ inline T dot (const vect& a, const vect& b) { return r; } -/** Return the Euklidean length. */ template inline T hypot (const vect& a) { return sqrt(dot(a,a)); } -/** Return the maximum element. */ template inline T maxval (const vect& a) { assert (D>0); T r(a[0]); - for (int d=1; d inline T minval (const vect& a) { assert (D>0); T r(a[0]); - for (int d=1; d -inline int maxloc (const vect& a) { - assert (D>0); - int r(0); - for (int d=1; da[r]) r=d; + for (int d=0; d -inline int minloc (const vect& a) { - assert (D>0); - int r(0); - for (int d=1; d inline T prod (const vect& a) { T r(1); @@ -648,13 +483,11 @@ inline T prod (const vect& a) { return r; } -/** Return the size (number of elements) of the vector. */ template inline int size (const vect& a) { return D; } -/** Return the sum of the elements. */ template inline T sum (const vect& a) { T r(0); @@ -662,136 +495,10 @@ inline T sum (const vect& a) { return r; } -// Higher order functions - -/** Return a new vector where the function func() has been applied to - all elements. */ -template -inline vect map (U (* const func)(T x), const vect& a) { - vect r; - for (int d=0; d -inline vect zip (U (* const func)(S x, T y), - const vect& a, const vect& b) -{ - vect r; - for (int d=0; d -inline U fold (U (* const func)(U val, T x), U val, const vect& a) -{ - for (int d=0; d -inline U fold1 (U (* const func)(U val, T x), const vect& a) -{ - assert (D>=1); - U val = a[0]; - for (int d=1; d -inline vect scan0 (U (* const func)(U val, T x), U val, - const vect& a) -{ - vect r; - for (int d=0; d -inline vect scan1 (U (* const func)(U val, T x), U val, - const vect& a) -{ - vect r; - for (int d=0; d -inline istream& operator>> (istream& is, vect& a) { - a.input(is); - return is; -} - - - -// Output - -/** Write a vector formatted to a stream. */ -template -inline ostream& operator<< (ostream& os, const vect& a) { - a.output(os); - return os; -} - - - -#if 0 -// Specialise explicit constructors - -/** Constructor for 2-element vectors from 2 elements. */ -template -inline vect::vect (const T x, const T y) { - elt[0]=x; elt[1]=y; -} - -/** Constructor for 3-element vectors from 3 elements. */ -vect (const T x, const T y, const T z) { - assert (D==3); - elt[0]=x; elt[1]=y; elt[2]=z; -} - -/** Constructor for 4-element vectors from 4 elements. */ -vect (const T x, const T y, const T z, const T t) { - assert (D==4); - elt[0]=x; elt[1]=y; elt[2]=z; elt[3]=t; -} +#if defined(TMPL_IMPLICIT) +# include "vect.cc" #endif - - -// Specialise for double - -template<> -inline vect& vect::operator%=(const vect& a) { - for (int d=0; d<3; ++d) { - elt[d]=fmod(elt[d],a[d]); - if (elt[d]>a[d]*double(1.0-1.0e-10)) elt[d]=double(0); - if (elt[d] +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "vect.hh" + +#include "wave.hh" + + + +typedef vect ivect; + + + +wave::wave (const int ref_fact, const int ref_levs, + const int mg_fact, const int mg_levs, + const ivect& lb, const ivect& ub, const ivect& str, const int tstr) + : ref_fact(ref_fact), ref_cent(vertex_centered), ref_levs(ref_levs), + mg_fact(mg_fact), mg_cent(vertex_centered), mg_levs(mg_levs), + global_lb(lb), global_ub(ub), global_str(str), + global_ext(global_lb,global_ub,global_str), + global_lg(1,1,1), global_ug(1,1,1), + time_lb(-2), time_ub(0), time_str(tstr), + hh(ref_fact,ref_cent,mg_fact,mg_cent,global_ext), + tt(hh,time_str), + dd(hh,global_lg,global_ug), + phi("phi",tt,dd,time_lb,time_ub), + anaphi("anaphi",tt,dd,time_lb,time_ub), + output_dir("output") +{ + clog << "WAVE" << endl; + clog << "hh: " << hh << endl; + clog << "tt: " << tt << endl; + clog << "dd: " << dd << endl; + clog << "phi: " << phi << endl; + clog << "anaphi: " << anaphi << endl; +} + + + +void wave::init () { + int size; + MPI_Comm_size (dist::comm, &size); + + vector > bbss(ref_levs); + ivect str = global_str; + ivect lb = global_lb; + ivect ub = global_ub; + for (int rl=0; rl0) { + // refined boxes have smaller stride + assert (all(str%ref_fact==0)); + str /= ref_fact; + // refine (arbitrarily) the center only + lb /= ref_fact; + ub /= ref_fact; + } + vector bbs(size); + for (int c=0; c > > bbsss + = hh.make_multigrid_boxes(bbss, mg_levs); + vector > pss(bbss.size()); + for (int rl=0; rl<(int)bbss.size(); ++rl) { + pss[rl] = vector(bbss[rl].size()); + for (int c=0; c<(int)bbss[rl].size(); ++c) { + pss[rl][c] = c % size; // distribute among processors + } + } + hh.recompose(bbsss, pss); + + clog << "WAVE-INIT" << endl; + clog << "hh: " << hh << endl; + clog << "tt: " << tt << endl; + clog << "dd: " << dd << endl; + clog << "phi: " << phi << endl; + clog << "anaphi: " << anaphi << endl; + + world_xmin = -5; + world_xmax = 5; + world_dx = (world_xmax - world_xmin) / dvect(global_ub - global_lb); + world_tmin = 0; + world_tmax = 2; + cfl_fact = 0.4; + world_dt = cfl_fact * maxval(world_dx); + + amr_init (0); +} + +void wave::amr_init (const int rl) { + const int tl=0; + const int ml=0; + + for (int c=0; cwrite_ascii(name.str(), dirs, t, tl, rl, c, ml); + } + + name.freeze(0); +} + + + +void wave::physics_analytic (const int tl, const int rl, const int c, + const int ml) +{ + if (! hh.is_local(rl,c)) return; + + clog << "analytic tl=" << tl << ", rl=" << rl << ", c=" << c << ", " + << "ml=" << ml << endl; + + const double A = 1; // amplitude + const double radius = 0; // radius + const double sigma = 1; // width + + const double t = tt.time(tl,rl,ml) * world_dt + world_tmin; + const double dt = tt.get_delta(rl,ml) * world_dt; + clog << "t=" << t << " dt=" << dt << endl; + + data& ap = *anaphi(tl,rl,c,ml); + + const ibbox ext = ap.extent(); + clog << " bbox=" << ext << endl; + for (ibbox::iterator it=ext.begin(); it!=ext.end(); ++it) { + const ivect index = *it; + + // equation in cartesian coordinates (t and vector x): + // (d2/dt2 - d2/dx2) phi = 0 + // one solution: + // phi_k(x,t) = f(t + k x) + g(t - k x) + // constraint: + // abs(k) = 1 + // test: + // (d2/dt2 + d2/dx2) phi_k = (1 - k^2) f'' + (1 - k^2) g'' = 0 + + // equation in D-dim spherical coordinates (only t and r, no Omega): + // (d2/dt2 - (D-1)/r d/dr - d2/dr2) phi = 0 + // solution: + // phi(r,t) = A/r f(t+r) + B/r f(t-r) + + const dvect x = dvect(index - global_lb) * world_dx + world_xmin; + const double r = hypot(x); + ap[index] = A * exp(- square((r - t - radius) / sigma)); + } +} + + + +void wave::physics_update (const int tl, const int rl, const int c, + const int ml) +{ + if (! hh.is_local(rl,c)) return; + + clog << "updating tl=" << tl << ", rl=" << rl << ", c=" << c << ", " + << "ml=" << ml << endl; + + const data& pp = *phi(tl-2,rl,c,ml); + const data& p = *phi(tl-1,rl,c,ml); + data& np = *phi(tl ,rl,c,ml); + + const ivect str = p.extent().stride(); + + const double dt = tt.get_delta(rl,ml) * world_dt; + const dvect dx = dvect(str) * world_dx; + const dvect dtdx2 = square(dvect(dt) / dx); + clog << "dt=" << dt << " dx=" << dx << endl; + + const ivect lb = p.extent().lower() + str; + const ivect ub = p.extent().upper() - str; + const ibbox ext(lb,ub,str); + clog << " bbox=" << ext << endl; + for (ibbox::iterator it=ext.begin(); it!=ext.end(); ++it) { + const ivect index = *it; + + double z = 2 * p[index] - pp[index]; + for (int d=0; d ivect; + typedef bbox ibbox; + + // Refinement + int ref_fact; + centering ref_cent; + int ref_levs; + + // Multigrid + int mg_fact; + centering mg_cent; + int mg_levs; + + // Extent + ivect global_lb, global_ub, global_str; + ibbox global_ext; + + // Ghosts + ivect global_lg, global_ug; + + // Time + int time_lb, time_ub, time_str; + + // Data + gh hh; + th tt; + dh dd; + gf phi, anaphi; + + // World + typedef vect dvect; + dvect world_xmin, world_xmax, world_dx; + double cfl_fact; + double world_tmin, world_tmax, world_dt; + + const char* output_dir; + +public: + wave (const int ref_fact, const int ref_levs, + const int mg_fact, const int mg_levs, + const ivect& lb, const ivect& ub, const ivect& str, const int tstr); + void init (); + void update (); + + void output (const int tl, const int rl, const int ml, const char* const dir) + const; + void output_var (const gf& var, + const int tl, const int rl, const int ml, + const char* const dir) + const; + template + void output_var_all_dirs (const gf& var, + const int tl, const int rl, const int ml, + const char* const dir) + const; + template + void output_var_one_dir (const gf& var, + const vect& dirs, + const int tl, const int rl, const int ml, + const char* const dir) + const; + +protected: + void amr_init (const int rl); + void amr_update (const int rl); + + void sync (const int tl, const int rl, const int ml); + + void physics_analytic (const int tl, const int rl, const int c, + const int ml); + void physics_update (const int tl, const int rl, const int c, + const int ml); + void physics_boundary (const int tl, const int rl, const int c, + const int ml); + void physics_setphi (const int tl, const int rl, const int c, + const int ml); + +}; + +#endif // WAVE_HH -- cgit v1.2.3