diff options
Diffstat (limited to 'Carpet/CarpetLib/src/ggf.cc')
-rw-r--r-- | Carpet/CarpetLib/src/ggf.cc | 608 |
1 files changed, 608 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc new file mode 100644 index 000000000..c26b22dbc --- /dev/null +++ b/Carpet/CarpetLib/src/ggf.cc @@ -0,0 +1,608 @@ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.46 2004/09/17 16:37:26 schnetter Exp $ + +#include <assert.h> +#include <math.h> +#include <stdlib.h> + +#include <iostream> +#include <string> + +#include "cctk.h" + +#include "defs.hh" +#include "dh.hh" +#include "th.hh" + +#include "ggf.hh" + +using namespace std; + + + +// Constructors +template<int D> +ggf<D>::ggf (const int varindex, const operator_type transport_operator, + th<D>& t, dh<D>& 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) +{ + assert (&t.h == &d.h); + + assert (vectorlength >= 1); + assert (vectorindex >= 0 && vectorindex < vectorlength); + assert ((vectorindex==0 && !vectorleader) + || (vectorindex!=0 && vectorleader)); + + d.add(this); +} + +// Destructors +template<int D> +ggf<D>::~ggf () { + d.remove(this); +} + +// Comparison +template<int D> +bool ggf<D>::operator== (const ggf<D>& f) const { + return this == &f; +} + + + +// Modifiers +template<int D> +void ggf<D>::recompose_crop () +{ + // Free storage that will not be needed + storage.resize(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); + } // for ml + } // for c + } // for rl + storage.at(tl-tmin).resize(h.reflevels()); + } // for tl +} + +template<int D> +void ggf<D>::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); + } + + // 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; c<h.components(rl); ++c) { + storage.at(tl-tmin).at(rl).at(c).resize(h.mglevels(rl,c)); + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + storage.at(tl-tmin).at(rl).at(c).at(ml) = typed_data(tl,rl,c,ml); + storage.at(tl-tmin).at(rl).at(c).at(ml)->allocate + (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c)); + } // for ml + } // for c + } // for tl +} + +template<int D> +void ggf<D>::recompose_fill (comm_state<D>& state, const int rl, + const bool do_prolongate) +{ + // Initialise the new storage + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + for (int tl=tmin; tl<=tmax; ++tl) { + + // Find out which regions need to be prolongated + // (Copy the exterior because some variables are not prolongated) + // TODO: do this once in the dh instead of for each variable here + ibset work (d.boxes.at(rl).at(c).at(ml).exterior); + + // Copy from old storage, if possible + // TODO: copy only from interior regions? + if (rl<(int)oldstorage.at(tl-tmin).size()) { + for (int cc=0; cc<(int)oldstorage.at(tl-tmin).at(rl).size(); ++cc) { + if (ml<(int)oldstorage.at(tl-tmin).at(rl).at(cc).size()) { + // TODO: prefer same processor, etc., see dh.cc + ibset ovlp + = work & oldstorage.at(tl-tmin).at(rl).at(cc).at(ml)->extent(); + 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<int> tls(numtl); + vector<CCTK_REAL> times(numtl); + for (int i=0; i<numtl; ++i) { + tls.at(i) = tmax - i; + times.at(i) = t.time(tls.at(i),rl-1,ml); + } + for (int cc=0; cc<(int)storage.at(tl-tmin).at(rl-1).size(); ++cc) { + vector<const gdata<D>*> gsrcs(numtl); + for (int i=0; i<numtl; ++i) { + gsrcs.at(i) + = storage.at(tls.at(i)-tmin).at(rl-1).at(cc).at(ml); + assert (gsrcs.at(i)->extent() == 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<int D> +void ggf<D>::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<int D> +void ggf<D>::recompose_bnd_prolongate (comm_state<D>& state, const int rl, + const bool do_prolongate) +{ + if (do_prolongate) { + // Set boundaries + if (rl>0) { + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + for (int tl=tmin; tl<=tmax; ++tl) { + + // TODO: assert that reflevel 0 boundaries are copied + const CCTK_REAL time = t.time(tl,rl,ml); + ref_bnd_prolongate (state,tl,rl,c,ml,time); + + } // for tl + } // for ml + } // for c + } // if rl + } // if do_prolongate +} + +template<int D> +void ggf<D>::recompose_sync (comm_state<D>& state, const int rl, + const bool do_prolongate) +{ + if (do_prolongate) { + // Set boundaries + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + for (int tl=tmin; tl<=tmax; ++tl) { + + sync (state,tl,rl,c,ml); + + } // for tl + } // for ml + } // for c + } // if do_prolongate +} + + + +// Cycle the time levels by rotating the data sets +template<int D> +void ggf<D>::cycle (int rl, int c, int ml) { + assert (rl>=0 && rl<h.reflevels()); + assert (c>=0 && c<h.components(rl)); + assert (ml>=0 && ml<h.mglevels(rl,c)); + gdata<D>* 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<int D> +void ggf<D>::flip (int rl, int c, int ml) { + assert (rl>=0 && rl<h.reflevels()); + assert (c>=0 && c<h.components(rl)); + assert (ml>=0 && ml<h.mglevels(rl,c)); + for (int t=0; t<(tmax-tmin)/2; ++t) { + const int tl1 = tmin + t; + const int tl2 = tmax - t; + assert (tl1 < tl2); + gdata<D>* 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; + } +} + + + +// Operations + +// Copy a region +template<int D> +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const ibbox dh<D>::dboxes::* recv_box, + int tl2, int rl2, int ml2, + const ibbox dh<D>::dboxes::* send_box) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + assert (tl2>=tmin && tl2<=tmax); + assert (rl2>=0 && rl2<h.reflevels()); + const int c2=c1; + assert (ml2<h.mglevels(rl2,c2)); + const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_box; + const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_box; + assert (all(recv.shape()==send.shape())); + // copy the content + assert (recv==send); + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from + (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), recv); +} + +// Copy regions +template<int D> +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const iblist dh<D>::dboxes::* recv_list, + int tl2, int rl2, int ml2, + const iblist dh<D>::dboxes::* send_list) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + assert (tl2>=tmin && tl2<=tmax); + assert (rl2>=0 && rl2<h.reflevels()); + const int c2=c1; + assert (ml2<h.mglevels(rl2,c2)); + 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 (recv.size()==send.size()); + // walk all boxes + for (typename iblist::const_iterator r=recv.begin(), s=send.begin(); + r!=recv.end(); ++r, ++s) { + // (use the send boxes for communication) + // copy the content + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from + (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r); + } +} + +// Copy regions +template<int D> +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const iblistvect dh<D>::dboxes::* recv_listvect, + int tl2, int rl2, int ml2, + const iblistvect dh<D>::dboxes::* send_listvect) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + assert (tl2>=tmin && tl2<=tmax); + assert (rl2>=0 && rl2<h.reflevels()); + // walk all components + for (int c2=0; c2<h.components(rl2); ++c2) { + assert (ml2<h.mglevels(rl2,c2)); + 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 (recv.size()==send.size()); + // walk all boxes + for (typename iblist::const_iterator r=recv.begin(), s=send.begin(); + r!=recv.end(); ++r, ++s) { + // (use the send boxes for communication) + // copy the content + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from + (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r); + } + } +} + +// Interpolate a region +template<int D> +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const ibbox dh<D>::dboxes::* recv_list, + const vector<int> tl2s, int rl2, int ml2, + const ibbox dh<D>::dboxes::* send_list, + CCTK_REAL time) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax); + } + assert (rl2>=0 && rl2<h.reflevels()); + const int c2=c1; + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); + + vector<const gdata<D>*> gsrcs(tl2s.size()); + vector<CCTK_REAL> 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())); + // interpolate the content + assert (recv==send); + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from + (state, gsrcs, times, recv, time, + d.prolongation_order_space, prolongation_order_time); +} + +// Interpolate regions +template<int D> +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const iblist dh<D>::dboxes::* recv_list, + const vector<int> tl2s, int rl2, int ml2, + const iblist dh<D>::dboxes::* send_list, + const CCTK_REAL time) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax); + } + assert (rl2>=0 && rl2<h.reflevels()); + const int c2=c1; + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); + + vector<const gdata<D>*> gsrcs(tl2s.size()); + vector<CCTK_REAL> 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 (recv.size()==send.size()); + // walk all boxes + for (typename iblist::const_iterator r=recv.begin(), s=send.begin(); + r!=recv.end(); ++r, ++s) { + // (use the send boxes for communication) + // interpolate the content + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from + (state, gsrcs, times, *r, time, + d.prolongation_order_space, prolongation_order_time); + } +} + +// Interpolate regions +template<int D> +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, + const iblistvect dh<D>::dboxes::* recv_listvect, + const vector<int> tl2s, int rl2, int ml2, + const iblistvect dh<D>::dboxes::* send_listvect, + const CCTK_REAL time) +{ + assert (tl1>=tmin && tl1<=tmax); + assert (rl1>=0 && rl1<h.reflevels()); + assert (c1>=0 && c1<h.components(rl1)); + assert (ml1>=0 && ml1<h.mglevels(rl1,c1)); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax); + } + assert (rl2>=0 && rl2<h.reflevels()); + // walk all components + for (int c2=0; c2<h.components(rl2); ++c2) { + assert (ml2>=0 && ml2<h.mglevels(rl2,c2)); + + vector<const gdata<D>*> gsrcs(tl2s.size()); + vector<CCTK_REAL> 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 (recv.size()==send.size()); + // walk all boxes + for (typename iblist::const_iterator r=recv.begin(), s=send.begin(); + r!=recv.end(); ++r, ++s) { + // (use the send boxes for communication) + // interpolate the content + storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from + (state, gsrcs, times, *r, time, + d.prolongation_order_space, prolongation_order_time); + } + } +} + + + +// Copy a component from the next time level +template<int D> +void ggf<D>::copy (comm_state<D>& state, int tl, int rl, int c, int ml) +{ + // Copy + copycat (state, + tl ,rl,c,ml, &dh<D>::dboxes::exterior, + tl+1,rl, ml, &dh<D>::dboxes::exterior); +} + +// Synchronise the boundaries a component +template<int D> +void ggf<D>::sync (comm_state<D>& state, int tl, int rl, int c, int ml) +{ + // Copy + copycat (state, + tl,rl,c,ml, &dh<D>::dboxes::recv_sync, + tl,rl, ml, &dh<D>::dboxes::send_sync); +} + +// Prolongate the boundaries of a component +template<int D> +void ggf<D>::ref_bnd_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, + CCTK_REAL time) +{ + // Interpolate + assert (rl>=1); + if (transport_operator == op_none) return; + vector<int> 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<D>::dboxes::recv_ref_bnd_coarse, + tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine, + time); +} + +// Restrict a multigrid level +template<int D> +void ggf<D>::mg_restrict (comm_state<D>& state, + int tl, int rl, int c, int ml, + CCTK_REAL time) +{ + // 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<int> tl2s(1,tl); + intercat (state, + tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, + tl2s,rl, ml-1, &dh<D>::dboxes::send_mg_fine, + time); +} + +// Prolongate a multigrid level +template<int D> +void ggf<D>::mg_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, + CCTK_REAL time) +{ + // 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<int> tl2s(1,tl); + intercat (state, + tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, + tl2s,rl, ml+1, &dh<D>::dboxes::send_mg_fine, + time); +} + +// Restrict a refinement level +template<int D> +void ggf<D>::ref_restrict (comm_state<D>& state, + int tl, int rl, int c, int ml, + CCTK_REAL time) +{ + // 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<int> tl2s(1,tl); + intercat (state, + tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine, + tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse, + time); +} + +// Prolongate a refinement level +template<int D> +void ggf<D>::ref_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, + CCTK_REAL time) +{ + assert (rl>=1); + if (transport_operator == op_none) return; + vector<int> 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<D>::dboxes::recv_ref_coarse, + tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine, + time); +} + + + +template class ggf<3>; |