diff options
Diffstat (limited to 'Carpet')
| -rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 79 | ||||
| -rw-r--r-- | Carpet/CarpetLib/src/dh.hh | 12 | ||||
| -rw-r--r-- | Carpet/CarpetLib/src/ggf.cc | 132 | ||||
| -rw-r--r-- | Carpet/CarpetLib/src/ggf.hh | 5 |
4 files changed, 142 insertions, 86 deletions
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 2895ac881..08fe03722 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -64,8 +64,9 @@ void dh::regrid () static Timer total ("dh::regrid"); total.start (); - boxes.clear(); - + oldboxes.clear(); + swap (boxes, oldboxes); + { static Timer timer ("dh::regrid::allocate_bboxes"); timer.start (); @@ -125,6 +126,12 @@ void dh::regrid () foreach_reflevel_component_mglevel (&dh::trim_unsynced_boundaries); timer.stop (0); } + { + static Timer timer ("dh::regrid::setup_old2new"); + timer.start (); + foreach_reflevel_component_mglevel (&dh::setup_old2new); + timer.stop (0); + } { static Timer timer ("dh::regrid::optimise_fields"); @@ -193,6 +200,10 @@ void dh::recompose (const int rl, const bool do_prolongate) (*f)->recompose_free_old (rl); timer.stop (0); } + // Omit prolongation and synchronisation. This is supposed to + // happen in the postregrid bin as part of the boundary conditions + // which are applied there. +#if 0 { static Timer timer ("dh::recompose::bnd_prolongate"); timer.start (); @@ -209,6 +220,7 @@ void dh::recompose (const int rl, const bool do_prolongate) } timer.stop (0); } +#endif } // for all grid functions of same vartype total.stop (0); @@ -711,6 +723,62 @@ void dh::trim_unsynced_boundaries (dh::dboxes & box, void dh:: +setup_old2new (dh::dboxes & box, + int const rl, int const c, int const ml) +{ + // Find out which regions need to be prolongated + // (Copy the exterior because some variables are not prolongated) + + if (oldboxes.empty()) return; + + ibset work = box.exterior; + + // Copy from old storage, if possible + // TODO: copy only from interior regions? + if (rl < (int)oldboxes.AT(ml).size()) { + for (int cc = 0; cc < (int)oldboxes.AT(ml).AT(rl).size(); ++ cc) { + // TODO: prefer same processor, etc. + ibset ovlp = work & oldboxes.AT(ml).AT(rl).AT(cc).exterior; + ovlp.normalize(); + work -= ovlp; + for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) { + box.old2new_recv_sync.AT(cc).push_back (* ri); + } + } // for cc + } // if rl + + // Initialise from coarser level, if possible + if (rl>0) { + for (int cc=0; cc<(int)boxes.AT(ml).AT(rl-1).size(); ++cc) { + + // TODO: choose larger regions first + // TODO: prefer regions from the same processor + iblist const & clist = box.recv_ref_coarse.AT(cc); + ibset cset; + for (iblist::const_iterator + iter = clist.begin();iter != clist.end(); ++iter) + { + cset += * iter; + } + ibset ovlp = work & cset; + ovlp.normalize(); + work -= ovlp; + for (ibset::const_iterator ri = ovlp.begin(); ri != ovlp.end(); ++ ri) { + box.old2new_recv_ref_coarse.AT(cc).push_back (* ri); + } + } // for cc + } // if rl + + // 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. +} + + + +void +dh:: optimise_field (dboxes & box, iblistvect const dboxes::* const field, pvect dboxes::* const field_fast, @@ -737,13 +805,14 @@ optimise_field (dboxes & box, void dh:: -optimise_fields (dboxes & box, - int const rl, int const c, int const ml) +optimise_fields (dboxes & box, int const rl, int const c, int const ml) { optimise_field (box, &dboxes::recv_ref_fine, &dboxes::recv_ref_fine_fast, rl, c, ml); optimise_field (box, &dboxes::recv_ref_coarse, &dboxes::recv_ref_coarse_fast, rl, c, ml); optimise_field (box, &dboxes::recv_sync, &dboxes::recv_sync_fast, rl, c, ml); optimise_field (box, &dboxes::recv_ref_bnd_coarse, &dboxes::recv_ref_bnd_coarse_fast, rl, c, ml); + optimise_field (box, &dboxes::old2new_recv_sync, &dboxes::old2new_recv_sync_fast, rl, c, ml); + optimise_field (box, &dboxes::old2new_recv_ref_coarse, &dboxes::old2new_recv_ref_coarse_fast, rl, c, ml); } void dh::do_check_bboxes (dh::dboxes & box, @@ -1157,6 +1226,8 @@ void dh::do_output_bboxes(dh::dboxes & box, cout << "recv_ref_bnd_coarse=" << box.recv_ref_bnd_coarse << endl; cout << "sync_not=" << box.sync_not << endl; cout << "recv_not=" << box.recv_not << endl; + cout << "old2new_recv_sync=" << box.old2new_recv_sync << endl; + cout << "old2new_recv_ref_coarse=" << box.old2new_recv_ref_coarse << endl; } void dh::output_bases () diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index ad3024dc1..e0f0749c5 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -91,6 +91,13 @@ public: ibset recv_not; // not received while syncing or // prolongating (globally outer // boundary) + + // Information for regridding, i.e., for copying data from the old + // to the new hierarchy + iblistvect old2new_recv_sync; + pvect old2new_recv_sync_fast; + iblistvect old2new_recv_ref_coarse; + pvect old2new_recv_ref_coarse_fast; }; private: @@ -123,13 +130,14 @@ private: void prepare_refinement_boundary_prolongation_boxes (dboxes & b, int rl, int c, int ml); void setup_refinement_boundary_prolongation_boxes (dboxes & b, int rl, int c, int ml); void setup_refinement_restriction_boxes (dboxes & b, int rl, int c, int ml); + void setup_old2new (dboxes & b, int rl, int c, int ml); + void trim_unsynced_boundaries (dboxes & b, int rl, int c, int ml); void optimise_field (dboxes & b, iblistvect const dboxes::* field, pvect dboxes::* field_fast, int rl, int c, int ml); void optimise_fields (dboxes & b, int rl, int c, int ml); - void trim_unsynced_boundaries (dboxes & b, int rl, int c, int ml); void do_output_bboxes (dboxes & b, int rl, int c, int ml); void do_check_bboxes (dboxes & b, int rl, int c, int ml); @@ -152,6 +160,8 @@ public: // should be readonly mbases bases; // bounding boxes around the grid // hierarchy + mboxes oldboxes; // old grid hierarchy, used during regridding + list<ggf*> gfs; // list of all grid functions public: diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 2a4e3857f..7c6ed319b 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -138,87 +138,55 @@ void ggf::recompose_allocate (const int rl) } // for ml } -void ggf::recompose_fill (comm_state& state, const int rl, - const bool do_prolongate) +void ggf::recompose_fill (comm_state & state, int const rl, + bool const do_prolongate) { // Initialise the new storage - for (int ml=0; ml<h.mglevels(); ++ml) { - for (int c=0; c<h.components(rl); ++c) { - for (int tl=0; tl<timelevels(ml,rl); ++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(ml).AT(rl).AT(c).exterior); - - // Copy from old storage, if possible - // TODO: copy only from interior regions? - if (rl<(int)oldstorage.AT(ml).size()) { - for (int cc=0; cc<(int)oldstorage.AT(ml).AT(rl).size(); ++cc) { - // TODO: prefer same processor, etc., see dh.cc - ibset ovlp - = work & oldstorage.AT(ml).AT(rl).AT(cc).AT(tl)->extent(); - ovlp.normalize(); - work -= ovlp; - for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) { - storage.AT(ml).AT(rl).AT(c).AT(tl)->copy_from - (state, oldstorage.AT(ml).AT(rl).AT(cc).AT(tl), *r); - } - } // for cc + for (int ml = 0; ml < h.mglevels(); ++ ml) { + + int const numtl = timelevels (ml, rl); + vector <int> tls (numtl); + vector <CCTK_REAL> times (numtl); + for (int tl = 0; tl < numtl; ++ tl) { + tls.AT(tl) = tl; + times.AT(tl) = t.time (tls.AT(tl), rl - 1, ml); + } + + for (int c = 0; c < h.components (rl); ++c) { + + // Initialise from the same level of the old hierarchy, where + // possible + if (rl < (int)oldstorage.AT(ml).size()) { + for (int tl = 0; tl < timelevels (ml, rl); ++tl) { + copycat (state, + tl, rl, c, ml, + & dh::dboxes::old2new_recv_sync_fast, + tl, rl, ml, + & oldstorage); + } // for tl + } // if rl + + if (do_prolongate) { + // Initialise from a coarser level of the new hierarchy, where + // possible + if (rl > 0) { + if (transport_operator != op_none) { + for (int tl = 0; tl < timelevels (ml, rl); ++tl) { + intercat (state, + tl, rl, c, ml, + & dh::dboxes::old2new_recv_ref_coarse_fast, + tls, rl - 1, ml, + times.AT(tl)); + } // for tl + } // if transport_operator } // if rl - - if (do_prolongate) { - // Initialise from coarser level, if possible - if (rl>0) { - if (transport_operator != op_none) { - const int pos = d.prolongation_order_space; - const int pot = (transport_operator != op_copy - ? prolongation_order_time - : 0); - const int numtl = pot+1; - assert (timelevels(ml,rl) >= numtl); - vector<int> tls(numtl); - vector<CCTK_REAL> times(numtl); - for (int i=0; i<numtl; ++i) { - tls.AT(i) = i; - times.AT(i) = t.time(tls.AT(i),rl-1,ml); - } - for (int cc=0; cc<(int)storage.AT(ml).AT(rl-1).size(); ++cc) { - vector<const gdata*> gsrcs(numtl); - for (int i=0; i<numtl; ++i) { - gsrcs.AT(i) = storage.AT(ml).AT(rl-1).AT(cc).AT(tls.AT(i)); - 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(ml).AT(rl).AT(c).recv_ref_coarse.AT(cc); - for (iblist::const_iterator iter=list.begin(); - iter!=list.end(); ++iter) - { - ibset ovlp = work & *iter; - ovlp.normalize(); - work -= ovlp; - for (ibset::const_iterator r=ovlp.begin(); - r!=ovlp.end(); ++r) - { - storage.AT(ml).AT(rl).AT(c).AT(tl)->interpolate_from - (state, gsrcs, times, *r, time, pos, pot); - } // 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 + } // 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 c } // for ml } @@ -236,6 +204,7 @@ void ggf::recompose_free_old (const int rl) } // for ml } +#if 0 void ggf::recompose_bnd_prolongate (comm_state& state, const int rl, const bool do_prolongate) { @@ -273,6 +242,7 @@ void ggf::recompose_sync (comm_state& state, const int rl, } // for ml } // if do_prolongate } +#endif @@ -402,7 +372,8 @@ void ggf::copycat (comm_state& state, void ggf::copycat (comm_state& state, int tl1, int rl1, int c1, int ml1, const pvect dh::dboxes::* recv_pvect, - int tl2, int rl2, int ml2) + int tl2, int rl2, int ml2, + mdata * const srcstorage_) { assert (rl1>=0 and rl1<h.reflevels()); assert (c1>=0 and c1<h.components(rl1)); @@ -411,11 +382,12 @@ void ggf::copycat (comm_state& state, assert ( ml2<h.mglevels()); assert (rl2>=0 and rl2<h.reflevels()); assert (tl2>=0 and tl2<timelevels(ml2,rl2)); + mdata & srcstorage = srcstorage_ ? * srcstorage_ : storage; // walk all components static Timer copycat1 ("copycat_pvect_1"); copycat1.start (); gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1); - cdata const & srcs = storage.AT(ml2).AT(rl2); + cdata const & srcs = srcstorage.AT(ml2).AT(rl2); pvect const & prs = d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_pvect; for (pvect::const_iterator ipr=prs.begin(); ipr!=prs.end(); ++ipr) { pseudoregion const & pr = * ipr; diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 3e113ee01..8908c30f7 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -96,8 +96,10 @@ public: void recompose_allocate (int rl); void recompose_fill (comm_state& state, int rl, bool do_prolongate); void recompose_free_old (int rl); +#if 0 void recompose_bnd_prolongate (comm_state& state, int rl, bool do_prolongate); void recompose_sync (comm_state& state, int rl, bool do_prolongate); +#endif // Cycle the time levels by rotating the data sets void cycle (int rl, int c, int ml); @@ -141,7 +143,8 @@ protected: void copycat (comm_state& state, int tl1, int rl1, int c1, int ml1, const pvect dh::dboxes::* recv_pvect, - int tl2, int rl2, int ml2); + int tl2, int rl2, int ml2, + mdata * srcstorage = 0); // Interpolate a region void intercat (comm_state& state, |