diff options
Diffstat (limited to 'Carpet/CarpetLib')
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 37 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.hh | 24 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/ggf.cc | 96 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/ggf.hh | 15 |
4 files changed, 161 insertions, 11 deletions
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index fb3a98538..7ea31573e 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -542,6 +542,43 @@ void dh::trim_unsynced_boundaries (dh::dboxes & box, } } +void +dh:: +optimise_field (dboxes & box, + iblistvect const dboxes::* const field, + pvect dboxes::* const field_fast, + int const rl, int const c, int const ml) +{ + size_t num_regions = 0; + for (int cc=0; cc<h.components(rl); ++cc) { + num_regions += (box.*field).AT(cc).size(); + } + assert ((box.*field_fast).empty()); + (box.*field_fast).reserve (num_regions); + for (int cc=0; cc<h.components(rl); ++cc) { + for (iblist::const_iterator + r=(box.*field).AT(cc).begin(); r!=(box.*field).AT(cc).end(); ++r) + { + pseudoregion pr; + pr.extent = * r; + pr.processor = cc; + (box.*field_fast).push_back (pr); + } + } + assert ((box.*field_fast).size() == num_regions); +} + +void +dh:: +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); +} + void dh::check_bboxes (dh::dboxes & box, int const rl, int const c, int const ml) { diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index bddc4360b..86e3c4f1a 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -17,6 +17,17 @@ using namespace std; +// A pseudoregion is almost a region; it is a bbox that lives on a +// certain processor. Pseudoregions are a compact way to store +// information about what processors needs to send data to what other +// processor during synchronisation or regridding. +struct pseudoregion { + ibbox extent; + int processor; +}; + + + // Forward declaration class ggf; class dh; @@ -34,6 +45,9 @@ public: typedef list<ibbox> iblist; typedef vector<iblist> iblistvect; // vector of lists + typedef vector <pseudoregion> pvect; + + // in here, the term "boundary" means both ghost zones and // refinement boundaries, but does not refer to outer (physical) @@ -57,13 +71,17 @@ public: iblistvect send_ref_fine; iblistvect send_ref_coarse; iblistvect recv_ref_fine; + pvect recv_ref_fine_fast; iblistvect recv_ref_coarse; + pvect recv_ref_coarse_fast; iblistvect send_sync; // send while syncing iblistvect send_ref_bnd_fine; // sent to finer grids ibset boundaries; // boundaries iblistvect recv_sync; // received while syncing + pvect recv_sync_fast; iblistvect recv_ref_bnd_coarse; // received from coarser grids + pvect recv_ref_bnd_coarse_fast; ibset sync_not; // not received while syncing (outer // boundary of that level) ibset recv_not; // not received while syncing or @@ -100,6 +118,12 @@ private: void setup_refinement_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 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 check_bboxes (dboxes & b, int rl, int c, int ml); diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 0778d83fe..7009ed0f9 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -375,16 +375,51 @@ void ggf::copycat (comm_state& state, // walk all components static Timer copycat1 ("copycat_listvect_1"); copycat1.start (); - for (int c2=0; c2<h.components(rl2); ++c2) { - const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2); + iblistvect::const_iterator recvi = + (d.boxes.AT(ml1).AT(rl1).AT(c1).*recv_listvect).begin(); + gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1); + cdata::const_iterator srci = storage.AT(ml2).AT(rl2).begin(); + int const nc = h.components(rl2); + for (int c2=0; c2<nc; ++c2) { + iblist const & recv = (*recvi); + gdata const * const src = (*srci).AT(tl2); // walk all boxes for (iblist::const_iterator r=recv.begin(); r!=recv.end(); ++r) { // (use the send boxes for communication) // copy the content - gdata* const dst = storage.at(ml1).at(rl1).at(c1).at(tl1); - gdata* const src = storage.at(ml2).at(rl2).at(c2).at(tl2); dst->copy_from(state, src, *r); } + ++ recvi; + ++ srci; + } + copycat1.stop (0); +} + +// Copy regions +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) +{ + assert (rl1>=0 and rl1<h.reflevels()); + assert (c1>=0 and c1<h.components(rl1)); + assert (ml1>=0 and ml1<h.mglevels()); + assert (tl1>=0 and tl1<timelevels(ml1,rl1)); + assert ( ml2<h.mglevels()); + assert (rl2>=0 and rl2<h.reflevels()); + assert (tl2>=0 and tl2<timelevels(ml2,rl2)); + // 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); + 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; + ibbox const & r = pr.extent; + int const c2 = pr.processor; + gdata const * const src = srcs.AT(c2).AT(tl2); + dst->copy_from(state, src, r); } copycat1.stop (0); } @@ -509,6 +544,45 @@ void ggf::intercat (comm_state& state, intercat1.stop (0); } +// Interpolate regions +void ggf::intercat (comm_state& state, + int tl1, int rl1, int c1, int ml1, + const pvect dh::dboxes::* recv_pvect, + const vector<int> tl2s, int rl2, int ml2, + const CCTK_REAL time) +{ + assert (rl1>=0 and rl1<h.reflevels()); + assert (c1>=0 and c1<h.components(rl1)); + assert (ml1>=0 and ml1<h.mglevels()); + assert (tl1>=0 and tl1<timelevels(ml1,rl1)); + assert (rl2>=0 and rl2<h.reflevels()); + assert (ml2>=0 and ml2<h.mglevels()); + vector<CCTK_REAL> times(tl2s.size()); + for (int i=0; i<(int)tl2s.size(); ++i) { + assert (tl2s.AT(i)>=0 and tl2s.AT(i)<timelevels(ml2,rl2)); + times.AT(i) = t.time(tl2s.AT(i),rl2,ml2); + } + // walk all components + static Timer intercat1 ("intercat_pvect_1"); + intercat1.start (); + int const pos = d.prolongation_order_space; + int const pot = transport_operator == op_copy ? 0 : prolongation_order_time; + gdata * const dst = storage.AT(ml1).AT(rl1).AT(c1).AT(tl1); + cdata const & srcs = storage.AT(ml2).AT(rl2); + vector<const gdata*> gsrcs(tl2s.size()); + 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; + ibbox const & r = pr.extent; + int const c2 = pr.processor; + for (int i=0; i<(int)gsrcs.size(); ++i) { + gsrcs.AT(i) = srcs.AT(c2).AT(tl2s.AT(i)); + } + dst->interpolate_from (state, gsrcs, times, r, time, pos, pot); + } + intercat1.stop (0); +} + // Copy a component from the next time level @@ -530,7 +604,7 @@ void ggf::sync (comm_state& state, int tl, int rl, int c, int ml) static Timer timer ("sync"); timer.start (); copycat (state, - tl,rl,c,ml, &dh::dboxes::recv_sync, + tl,rl,c,ml, &dh::dboxes::recv_sync_fast, tl,rl, ml); timer.stop (0); } @@ -565,7 +639,7 @@ void ggf::ref_bnd_prolongate (comm_state& state, tl2s.AT(0) = 0; } intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_bnd_coarse, + tl ,rl ,c,ml, &dh::dboxes::recv_ref_bnd_coarse_fast, tl2s,rl-1, ml, time); timer.stop (0); } @@ -613,11 +687,11 @@ void ggf::ref_restrict (comm_state& state, 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); static Timer timer ("ref_restrict"); timer.start (); + vector<int> const tl2s(1,tl); intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_fine, + tl ,rl ,c,ml, &dh::dboxes::recv_ref_fine_fast, tl2s,rl+1, ml, time); timer.stop (0); } @@ -629,15 +703,15 @@ void ggf::ref_prolongate (comm_state& state, { assert (rl>=1); if (transport_operator == op_none) return; - vector<int> tl2s; static Timer timer ("ref_prolongate"); timer.start (); + vector<int> tl2s; // Interpolation in time assert (timelevels(ml,rl) >= prolongation_order_time+1); tl2s.resize(prolongation_order_time+1); - for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i; + for (int i=0; i<=prolongation_order_time; ++i) tl2s.AT(i) = i; intercat (state, - tl ,rl ,c,ml, &dh::dboxes::recv_ref_coarse, + tl ,rl ,c,ml, &dh::dboxes::recv_ref_coarse_fast, tl2s,rl-1, ml, time); timer.stop (0); } diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 30752a422..3e113ee01 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -33,6 +33,8 @@ class ggf { typedef list<ibbox> iblist; typedef vector<iblist> iblistvect; + typedef vector <pseudoregion> pvect; + typedef gdata* tdata; // data ... typedef vector<tdata> fdata; // ... for each time level typedef vector<fdata> cdata; // ... for each component @@ -134,6 +136,12 @@ protected: int tl1, int rl1, int c1, int ml1, const iblistvect dh::dboxes::* recv_listvect, int tl2, int rl2, int ml2); + + // Copy regions + 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); // Interpolate a region void intercat (comm_state& state, @@ -156,6 +164,13 @@ protected: const vector<int> tl2s, int rl2, int ml2, CCTK_REAL time); + // Interpolate regions + void intercat (comm_state& state, + int tl1, int rl1, int c1, int ml1, + const pvect dh::dboxes::* recv_pvect, + const vector<int> tl2s, int rl2, int ml2, + CCTK_REAL time); + public: |