From 351d38c7a8b35e6b4250c098e007341b0f8e3001 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 14 Sep 2010 15:56:32 -0500 Subject: CarpetLib: Correct refluxing errors --- Carpet/CarpetLib/src/data.cc | 60 +++-- Carpet/CarpetLib/src/dh.cc | 615 +++++++++++++++++++++--------------------- Carpet/CarpetLib/src/gdata.cc | 18 +- 3 files changed, 361 insertions(+), 332 deletions(-) diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 2c76c6a3b..49229008a 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -56,7 +56,7 @@ call_operator (void #ifndef _OPENMP (* the_operator) (src, srcext, dst, dstext, srcbbox, dstbbox, regbbox); #else -# if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE) +# ifdef CARPET_DEBUG ibset allregbboxes; # endif #pragma omp parallel @@ -86,13 +86,13 @@ call_operator (void regbbox.stride()); if (not myregbbox.empty()) { (* the_operator) (src, srcext, dst, dstext, srcbbox, dstbbox, myregbbox); -# if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE) +# ifdef CARPET_DEBUG #pragma omp critical allregbboxes += myregbbox; # endif } } -# if ! defined (NDEBUG) && ! defined (CARPET_OPTIMISE) +# ifdef CARPET_DEBUG if (not (allregbboxes == ibset (regbbox))) { cout << "allregbboxes=" << allregbboxes << endl << "regbbox=" << regbbox << endl; @@ -408,14 +408,35 @@ copy_from_innerloop (gdata const * const gsrc, assert (proc() == src->proc()); assert (dist::rank() == proc()); + ibbox srcbox, dstbox; + switch (cent) { + case vertex_centered: + srcbox = src->extent(); + dstbox = this->extent(); + break; + case cell_centered: { + ivect const ioff = box.lower() - this->extent().lower(); + ivect const is_centered = ioff % this->extent().stride() == 0; + + // Shift bboxes to be face centred if necessary, since all grid + // functions are stored as if they were cell-centered + // srcbox = src->extent().shift(is_centered-1,2); + srcbox = src->extent(); + dstbox = this->extent().shift(is_centered-1,2); + break; + } + default: + assert (0); + } + #if CARPET_DIM == 3 call_operator (& copy_3d, static_cast (src->storage()), src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); #elif CARPET_DIM == 4 call_operator (& copy_4d, @@ -423,8 +444,8 @@ copy_from_innerloop (gdata const * const gsrc, src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); #else # error "Value for CARPET_DIM not supported" @@ -977,17 +998,22 @@ transfer_restrict (data const * const src, break; case cell_centered: { assert (all (box.stride() == this->extent().stride())); - ivect const izero (0); ivect const ioff = box.lower() - this->extent().lower(); - ivect const is_centered = ioff % this->extent().stride() == izero; + ivect const is_centered = ioff % this->extent().stride() == 0; + + // Shift bboxes to be face centred if necessary, since all grid + // functions are stored as if they were cell-centered + ibbox const srcbox = src->extent().shift(is_centered-1,2); + ibbox const dstbox = this->extent().shift(is_centered-1,2); + if (all(is_centered == ivect(1,1,1))) { call_operator (& restrict_3d_cc_rf2, static_cast (src->storage()), src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); } else if (all(is_centered == ivect(0,1,1))) { call_operator (& restrict_3d_vc_rf2, @@ -995,8 +1021,8 @@ transfer_restrict (data const * const src, src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); } else if (all(is_centered == ivect(1,0,1))) { call_operator (& restrict_3d_vc_rf2, @@ -1004,8 +1030,8 @@ transfer_restrict (data const * const src, src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); } else if (all(is_centered == ivect(1,1,0))) { call_operator (& restrict_3d_vc_rf2, @@ -1013,8 +1039,8 @@ transfer_restrict (data const * const src, src->shape(), static_cast (this->storage()), this->shape(), - src->extent(), - this->extent(), + srcbox, + dstbox, box); } else { assert (0); diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 7ecba6785..4e728717f 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -936,59 +936,317 @@ regrid (bool const do_init) if (rl > 0) { int const orl = rl - 1; - fast_dboxes & fast_olevel = fast_boxes.AT(ml).AT(orl); + full_cboxes const& full_olevel = full_boxes.AT(ml).AT(orl); + fast_dboxes& fast_olevel = fast_boxes.AT(ml).AT(orl); + ibbox const& odomext = h.baseextent(ml,orl); + + // Refinement restriction may fill all coarse interior points, + // and must use all fine active points + + ibset allrestricted; + switch (h.refcent) { + case vertex_centered: + allrestricted = allactive.contracted_for(odomext); + break; + case cell_centered: { + ibset const& source = allactive; + ibbox const& target = odomext; + ibbox const all = allactive.container().expand(10,10); + ibbox const all_target = all.contracted_for(target); + ibset const tmp0 = source; + ibset const tmp1 = tmp0.expanded_for(target); + ibset const tmp2 = all_target - tmp1; + ibset const tmp3 = tmp2.expand(1,1); + ibset const tmp4 = all_target - tmp3; + allrestricted = tmp4; + cout << "source=" << source << "\n" + << "target=" << target << "\n" + << "all=" << all << "\n" + << "all_target=" << all_target << "\n" + << "tmp1=" << tmp1 << "\n" + << "tmp2=" << tmp2 << "\n" + << "tmp3=" << tmp3 << "\n" + << "allrestricted=" << allrestricted << "\n"; + break; + } + default: + assert(0); + } - if (h.components(orl) > 0) { - for (int lc = 0; lc < h.local_components(rl); ++ lc) { - int const c = h.get_component (rl, lc); - - full_dboxes const & box = full_level.AT(c); - full_dboxes const & obox0 = full_boxes.AT(ml).AT(orl).AT(0); - - // Refinement restriction may fill all coarse active - // points, and must use all fine active points + for (int olc = 0; olc < h.local_components(orl); ++ olc) { + int const oc = h.get_component(orl, olc); + full_dboxes const& obox = full_olevel.AT(oc); + + ibset needrecv = allrestricted & obox.owned; + // Cannot restrict into buffer zones + assert ((allrestricted & obox.buffers).empty()); + + for (int c = 0; c < h.components(rl); ++ c) { + full_dboxes const& box = full_level.AT(c); - ibset needrecv (box.active.contracted_for (obox0.interior)); + ibbox const contracted_exterior = + box.exterior.contracted_for(odomext); + ibset const ovlp = needrecv & contracted_exterior; - for (int cc = 0; cc < h.components(orl); ++ cc) { - full_dboxes & obox = full_boxes.AT(ml).AT(orl).AT(cc); - - ibset const contracted_active - (box.active.contracted_for (obox.interior)); - ibset const ovlp = obox.active & contracted_active; + for (ibset::const_iterator + ri = ovlp.begin(); ri != ovlp.end(); ++ ri) + { + ibbox const& recv = *ri; + ibbox const send = recv.expanded_for(box.exterior); + ASSERT_c (send <= box.exterior, + "Refinement restriction: Send region must be contained in exterior"); - for (ibset::const_iterator - ri = ovlp.begin(); ri != ovlp.end(); ++ ri) - { - ibbox const & recv = * ri; - ibbox const send = recv.expanded_for (box.interior); - ASSERT_c (send <= box.active, - "Refinement restriction: Send region must be contained in active part"); - fast_olevel.fast_ref_rest_sendrecv.push_back - (sendrecv_pseudoregion_t (send, c, recv, cc)); - if (not on_this_proc (orl, cc)) { - fast_dboxes & fast_level_otherproc = - fast_level_otherprocs.AT(this_proc(orl, cc)); - fast_level_otherproc.fast_ref_rest_sendrecv.push_back - (sendrecv_pseudoregion_t (send, c, recv, cc)); - } + sendrecv_pseudoregion_t const preg (send, c, recv, oc); + fast_olevel.fast_ref_rest_sendrecv.push_back(preg); + if (not on_this_proc(rl, c)) { + fast_dboxes& fast_level_otherproc = + fast_level_otherprocs.AT(this_proc(rl, c)); + fast_level_otherproc.fast_ref_rest_sendrecv.push_back(preg); } - - needrecv -= ovlp; - - } // for cc + } - // All points must have been received - ASSERT_rl (needrecv.empty(), - "Refinement restriction: All points must have been received"); + needrecv -= ovlp; - } // for lc - } // if orl not empty + } // for c + + // All points must have been received + ASSERT_rl (needrecv.empty(), + "Refinement restriction: All points must have been received"); + + } // for olc } // if rl > 0 timer_comm_refrest.stop(); + + + + // Refinement refluxing: + + static Carpet::Timer timer_comm_reflux + ("CarpetLib::dh::regrid::comm::reflux"); + timer_comm_reflux.start(); + + // If there is no coarser level, do nothing + if (rl > 0) { + int const orl = rl - 1; + local_cboxes & local_olevel = local_boxes.AT(ml).AT(orl); + full_cboxes const& full_olevel = full_boxes.AT(ml).AT(orl); + fast_dboxes & fast_olevel = fast_boxes.AT(ml).AT(orl); + + // This works only with cell centred grids + if (h.refcent == cell_centered) { + + + + // Fine grids + ibset const& fine_level = allactive; + + ivect const izero = ivect (0); + ivect const ione = ivect (1); + + assert (all (h.reffacts.AT(rl) % h.reffacts.AT(orl) == 0)); + ivect const reffact = h.reffacts.AT(rl) / h.reffacts.AT(orl); + + vect,dim> all_fine_boundary; + + for (int dir=0; dir,dim> all_coarse_boundary; + for (int dir=0; dir<3; ++dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + ibbox const coarse_faces = coarse_ext.shift(-idir,2); + ibbox const fine_faces = fine_ext.shift(-idir,2); + for (int face=0; face<2; ++face) { + cout << "REFREF rl=" << rl << " dir=" << dir << " face=" << face << "\n" + << " all_fine_boundary=" << all_fine_boundary[dir][face] << "\n" + << " coarse_ext=" << coarse_ext << "\n" + << " coarse_faces=" << coarse_faces << "\n"; + all_coarse_boundary[dir][face] = + all_fine_boundary[dir][face].contracted_for (coarse_faces); + cout << " all_coarse_boundary=" << all_coarse_boundary[dir][face] << "\n"; + ASSERT_rl (all_coarse_boundary[dir][face].expanded_for(fine_faces) == all_fine_boundary[dir][face], + "Refluxing: Coarse grid boundaries must be consistent with fine grid boundaries"); + } // for face + } // for dir + + + + // Split onto components + for (int lc = 0; lc < h.local_components(rl); ++ lc) { + int const c = h.get_component(rl, lc); + full_dboxes & box = full_level.AT(c); + local_dboxes & local_box = local_level.AT(lc); + + for (int dir = 0; dir < dim; ++ dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + for (int face = 0; face < 2; ++ face) { + + // This is not really used (only for debugging) + local_box.fine_boundary[dir][face] = + box.exterior.shift(-idir,2) & all_fine_boundary[dir][face]; + cout << "REFREF rl=" << rl << " lc=" << lc << " dir=" << dir << " face=" << face << "\n" + << " local.fine_boundary=" << local_box.fine_boundary[dir][face] << "\n"; + + } // for face + } // for dir + + } // for lc + +#ifdef CARPET_DEBUG + for (int dir = 0; dir < dim; ++ dir) { + ivect const idir = ivect::dir(dir); // Unit vector + for (int face = 0; face < 2; ++ face) { + ibset all_fine_boundary_combined; + for (int c = 0; c < h.components(rl); ++ c) { + full_dboxes const & box = full_level.AT(c); + all_fine_boundary_combined += + box.exterior.shift(-idir,2) & all_fine_boundary[dir][face]; + } // for c + ASSERT_rl (all_fine_boundary_combined == all_fine_boundary[dir][face], + "Refluxing: Union of all fine grid boundaries must be the total fine grid boundary "); + } + } +#endif + + for (int lc = 0; lc < h.local_components(orl); ++ lc) { + int const c = h.get_component(orl, lc); + full_dboxes const & obox = full_olevel.AT(c); + local_dboxes & local_obox = local_olevel.AT(lc); + + for (int dir = 0; dir < dim; ++ dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + for (int face = 0; face < 2; ++ face) { + + // This is used for post-processing the fluxes + // (calculating the difference between coarse and fine + // fluxes, adjusting the state) + local_obox.coarse_boundary[dir][face] = + obox.owned.shift(-idir,2) & all_coarse_boundary[dir][face]; + // Cannot reflux into buffer zones + cout << "REFREF orl=" << orl << " lc=" << lc << " dir=" << dir << " face=" << face << "\n" + << " local.coarse_boundary=" << local_obox.coarse_boundary[dir][face] << "\n"; + assert ((obox.buffers.shift(-idir,2) & all_coarse_boundary[dir][face]).empty()); + + } // for face + } // for dir + + } // for lc + +#ifdef CARPET_DEBUG + for (int dir = 0; dir < dim; ++ dir) { + ivect const idir = ivect::dir(dir); // Unit vector + for (int face = 0; face < 2; ++ face) { + ibset all_coarse_boundary_combined; + for (int c = 0; c < h.components(orl); ++ c) { + full_dboxes const & obox = full_olevel.AT(c); + all_coarse_boundary_combined += + obox.owned.shift(face ? +idir : -idir,2) & all_coarse_boundary[dir][face]; + } // for c + ASSERT_rl (all_coarse_boundary_combined == all_coarse_boundary[dir][face], + "Refluxing: Union of all coarse grid boundaries must be the total coarse grid boundary "); + } + } +#endif + + + + // Set up communication schedule + for (int olc = 0; olc < h.local_components(orl); ++ olc) { + int const oc = h.get_component(orl, olc); + local_dboxes & local_obox = local_olevel.AT(olc); + cout << "REF ref_refl ml=" << ml << " rl=" << rl << " olc=" << olc << " oc=" << oc << "\n"; + + for (int dir = 0; dir < dim; ++ dir) { + // Unit vector + ivect const idir = ivect::dir(dir); + ibbox const coarse_faces = coarse_ext.shift(-idir,2); + ibbox const fine_faces = h.baseextent(ml, rl).shift(-idir,2); + for (int face = 0; face < 2; ++ face) { + + srpvect fast_dboxes::* const fast_ref_refl_sendrecv = + fast_dboxes::fast_ref_refl_sendrecv[dir][face]; + + // Refluxing must fill all coarse refluxing boundary + // points, and may use all fine points + + ibset needrecv (local_obox.coarse_boundary[dir][face]); + + for (int c = 0; c < h.components(rl); ++ c) { + full_dboxes const & box = full_level.AT(c); + cout << "REF ref_refl ml=" << ml << " rl=" << rl << " olc=" << olc << " oc=" << oc << " c=" << c << " dir=" << dir << " face=" << face << "\n"; + + ibbox const contracted_exterior = + box.exterior.shift(-idir, 2).contracted_for(coarse_faces); + cout << " exterior=" << box.exterior.shift(-idir, 2) << "\n" + << " contracted=" << contracted_exterior << "\n"; + ibset const ovlp = needrecv & contracted_exterior; + cout << " ovlp=" << ovlp << "\n"; + + for (ibset::const_iterator + ri = ovlp.begin(); ri != ovlp.end(); ++ ri) + { + ibbox const & recv = * ri; + ibbox const send = + recv.expanded_for (box.exterior.shift(-idir, 2)); + ASSERT_c (send <= box.exterior.shift(-idir, 2), + "Refinement restriction: Send region must be contained in exterior"); + + sendrecv_pseudoregion_t const preg (send, c, recv, oc); + cout << "REF ref_refl ml=" << ml << " rl=" << rl << " olc=" << olc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << "\n" + << " preg=" << preg << "\n"; + (fast_olevel.*fast_ref_refl_sendrecv).push_back (preg); + if (not on_this_proc (rl, c)) { + fast_dboxes & fast_level_otherproc = + fast_level_otherprocs.AT(this_proc(rl, c)); + (fast_level_otherproc.*fast_ref_refl_sendrecv). + push_back (preg); + } + } + + needrecv -= ovlp; + + } // for c + + // All points must have been received + if (not needrecv.empty()) { + cout << "needrecv=" << needrecv << endl; + } + ASSERT_rl (needrecv.empty(), + "Refinement refluxing: All points must have been received"); + + } // for face + } // for dir + + } // for olc + + } // if cell centered + } // if rl > 0 + + timer_comm_reflux.stop(); + timer_comm.stop(); @@ -1110,7 +1368,7 @@ regrid (bool const do_init) // Mask for unused points on coarser level (which do not influence the future // evolution provided regridding is done at the right times): static Carpet::Timer timer_overwrittenmask ("CarpetLib::dh::regrid::unusedpoints_mask"); - timer_mask.start(); + timer_overwrittenmask.start(); if (rl > 0) { int const orl = rl - 1; @@ -1148,278 +1406,7 @@ regrid (bool const do_init) } // if reffact != 2 } // if not coarsest level - timer_mask.stop(); - - - - // Refluxing: - static Carpet::Timer timer_reflux ("CarpetLib::dh::regrid::reflux"); - timer_reflux.start(); - - // If there is no coarser level, do nothing - if (rl > 0) { - int const orl = rl - 1; - light_cboxes & light_olevel = light_boxes.AT(ml).AT(orl); - local_cboxes & local_olevel = local_boxes.AT(ml).AT(orl); - full_cboxes const& full_olevel = full_boxes.AT(ml).AT(orl); - fast_dboxes & fast_olevel = fast_boxes.AT(ml).AT(orl); - - // This works only with cell centred grids - if (h.refcent == cell_centered) { - - // Fine grids - ibset const& fine_level = allactive; - - ivect const izero = ivect (0); - ivect const ione = ivect (1); - - assert (all (h.reffacts.AT(rl) % h.reffacts.AT(orl) == 0)); - ivect const reffact = h.reffacts.AT(rl) / h.reffacts.AT(orl); - - vect,dim> all_fine_boundary; - - for (int dir=0; dir,dim> all_coarse_boundary; - for (int dir=0; dir<3; ++dir) { - // Unit vector - ivect const idir = ivect::dir(dir); - ibbox const coarse_faces = coarse_ext.shift(idir,2); - for (int face=0; face<2; ++face) { - cout << "REFREF rl=" << rl << " dir=" << dir << " face=" << face << "\n" - << " all_fine_boundary=" << all_fine_boundary[dir][face] << "\n" - << " coarse_ext=" << coarse_ext << "\n" - << " coarse_faces=" << coarse_faces << "\n"; - all_coarse_boundary[dir][face] = - all_fine_boundary[dir][face].contracted_for (coarse_faces); - cout << " all_coarse_boundary=" << all_coarse_boundary[dir][face] << "\n"; - } // for face - } // for dir - - - - for (int lc = 0; lc < h.local_components(rl); ++ lc) { - int const c = h.get_component(rl, lc); - full_dboxes & box = full_level.AT(c); - local_dboxes & local_box = local_level.AT(lc); - - for (int dir = 0; dir < dim; ++ dir) { - for (int face = 0; face < 2; ++ face) { - - // This is not really used (only for debugging) - local_box.fine_boundary[dir][face] = - box.exterior & all_fine_boundary[dir][face]; - cout << "REFREF rl=" << rl << " lc=" << lc << " dir=" << dir << " face=" << face << "\n" - << " local.fine_boundary=" << local_box.fine_boundary[dir][face] << "\n"; - - } // for face - } // for dir - - } // for lc - - for (int lc = 0; lc < h.local_components(orl); ++ lc) { - int const c = h.get_component(orl, lc); - full_dboxes const & obox = full_olevel.AT(c); - local_dboxes & local_obox = local_olevel.AT(lc); - - for (int dir = 0; dir < dim; ++ dir) { - for (int face = 0; face < 2; ++ face) { - - // This is used for post-processing the fluxes - // (calculating the difference between coarse and fine - // fluxes, adjusting the state) - local_obox.coarse_boundary[dir][face] = - obox.exterior & all_coarse_boundary[dir][face]; - cout << "REFREF orl=" << orl << " lc=" << lc << " dir=" << dir << " face=" << face << "\n" - << " local.coarse_boundary=" << local_obox.coarse_boundary[dir][face] << "\n"; - - } // for face - } // for dir - - } // for lc - - for (int lc = 0; lc < h.local_components(orl); ++ lc) { - int const oc = h.get_component(orl, lc); - light_dboxes & light_obox = light_olevel.AT(oc); - local_dboxes & local_obox = local_olevel.AT(lc); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " oc=" << oc << "\n"; - - for (int dir = 0; dir < dim; ++ dir) { - for (int face = 0; face < 2; ++ face) { - // Unit vector - ivect const idir = ivect::dir(dir); - - srpvect fast_dboxes::* const fast_ref_refl_sendrecv = - fast_dboxes::fast_ref_refl_sendrecv[dir][face]; - - // Refluxing must fill all coarse refluxing boundary - // points, and may use all fine points - - ibset needrecv (local_obox.coarse_boundary[dir][face]); - - for (int c = 0; c < h.components(rl); ++ c) { - full_dboxes const & box = full_level.AT(c); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " oc=" << oc << " c=" << c << " dir=" << dir << " face=" << face << "\n"; - - ibbox const contracted_exterior = - box.exterior.contracted_for (light_obox.exterior); - cout << " exterior=" << box.exterior << "\n" - << " contracted=" << contracted_exterior << "\n"; - ibset const ovlp = needrecv & contracted_exterior; - cout << " ovlp=" << ovlp << "\n"; - - for (ibset::const_iterator - ri = ovlp.begin(); ri != ovlp.end(); ++ ri) - { - ibbox const & recv = * ri; - ibbox const send = recv.expanded_for (box.exterior); - ASSERT_c (send <= box.exterior, - "Refinement restriction: Send region must be contained in exterior"); - - sendrecv_pseudoregion_t const preg (send, c, recv, oc); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << "\n" - << " preg=" << preg << "\n"; - (fast_olevel.*fast_ref_refl_sendrecv).push_back (preg); - if (not on_this_proc (orl, oc)) { - fast_dboxes & fast_level_otherproc = - fast_level_otherprocs.AT(this_proc(orl, oc)); - (fast_level_otherproc.*fast_ref_refl_sendrecv). - push_back (preg); - } - } - - needrecv -= ovlp; - - } // for c - - // All points must have been received - if (not needrecv.empty()) { - cout << "needrecv=" << needrecv << endl; - } - ASSERT_rl (needrecv.empty(), - "Refinement refluxing: All points must have been received"); - - } // for face - } // for dir - - } // for lc - -#if 0 - for (int lc = 0; lc < h.local_components(rl); ++ lc) { - int const c = h.get_component(rl, lc); - full_dboxes const & box = full_level.AT(c); - local_dboxes const & local_box = local_level.AT(lc); - full_dboxes const & obox0 = full_boxes.AT(ml).AT(orl).AT(0); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << "\n"; - - for (int dir = 0; dir < dim; ++ dir) { - for (int face = 0; face < 2; ++ face) { - // Unit vector - ivect const idir = ivect::dir(dir); - - srpvect fast_dboxes::* const fast_ref_refl_sendrecv = - fast_dboxes::fast_ref_refl_sendrecv[dir][face]; - - // Refluxing may fill all coarse active points, and - // must use all fine refluxing boundary points - - // ibset needrecv (box.active.contracted_for (obox0.interior)); - // ibset needrecv - // (local_box.coarse_boundary[dir][face]. - // contracted_for (obox0.interior)); -#error "should this be local_obox instead of local_box?" - ibset needrecv - (local_box.coarse_boundary[dir][face]. - shift (idir*(1-reffact), 2). - contracted_for (obox0.interior)); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << " dir=" << dir << " face=" << face << "\n" - << " coarse_boundary=" << local_box.coarse_boundary[dir][face] << "\n" - << " shifted=" << local_box.coarse_boundary[dir][face].shift (idir*(1-reffact), 2) << "\n" - << " needrecv=" << needrecv << "\n"; - - for (int oc = 0; oc < h.components(orl); ++ oc) { - full_dboxes & obox = full_boxes.AT(ml).AT(orl).AT(oc); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << "\n"; - - ibset const contracted_boundary - (local_box.fine_boundary[dir][face]. - contracted_for (obox.interior)); - cout << " fine_boundary=" << local_box.fine_boundary[dir][face] << "\n"; - cout << " contracted_boundary=" << contracted_boundary << "\n"; - ibset const ovlp = obox.active & contracted_boundary; - cout << " ovlp=" << ovlp << "\n"; - - for (ibset::const_iterator - ri = ovlp.begin(); ri != ovlp.end(); ++ ri) - { - ibbox const & recv = * ri; -#warning "TODO: need to use correct coarse/fine transition" - ibbox const send = recv.expanded_for (box.interior); - ASSERT_c (send <= box.active, - "Refinement restriction: Send region must be contained in active part"); - - // Modify the extents since the flux grid - // functions are staggered - assert (all (send.stride() % 2 == 0)); - ibbox const msend (send.lower() + idir * send.stride()/2, - send.upper() + idir * send.stride()/2, - send.stride()); - assert (all (recv.stride() % 2 == 0)); - ibbox const mrecv (recv.lower() + idir * recv.stride()/2, - recv.upper() + idir * recv.stride()/2, - recv.stride()); - - sendrecv_pseudoregion_t const preg (msend, c, mrecv, oc); - (fast_olevel.*fast_ref_refl_sendrecv).push_back (preg); - cout << "REF ref_refl ml=" << ml << " rl=" << rl << " lc=" << lc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << " preg=" << preg << "\n"; - if (not on_this_proc (orl, oc)) { - fast_dboxes & fast_level_otherproc = - fast_level_otherprocs.AT(this_proc(orl, oc)); - (fast_level_otherproc.*fast_ref_refl_sendrecv). - push_back (preg); - } - } - - needrecv -= ovlp; - - } // for oc - - // All points must have been received - if (not needrecv.empty()) { - cout << "needrecv=" << needrecv << endl; - } - ASSERT_rl (needrecv.empty(), - "Refinement refluxing: All points must have been received"); - - } // for face - } // for dir - - } // for lc -#endif - - } // if cell centered - - } // if rl > 0 - - timer_reflux.stop(); + timer_overwrittenmask.stop(); diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index 0d8980a39..c9f9f8ad2 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -191,12 +191,21 @@ transfer_from (comm_state & state, if (is_src) { // copy the data into the send buffer if (interp_on_src) { + ivect ioffset (0); + if (src->cent == cell_centered) { + assert (all (srcbox.stride() == src->extent().stride())); + ivect const ioff = srcbox.lower() - src->extent().lower(); + ivect const is_centered = ioff % src->extent().stride() == 0; + ioffset = 1 - is_centered; + } + ibbox const bufbox = dstbox.shift(ioffset, 2); + assert (bufbox.size() == dstbox.size()); size_t const sendbufsize = src->c_datatype_size() * dstbox.size(); void * const sendbuf = state.send_buffer (src->c_datatype(), dstproc, dstbox.size()); gdata * const buf = src->make_typed (src->varindex, src->cent, src->transport_operator); - buf->allocate (dstbox, srcproc, sendbuf, sendbufsize); + buf->allocate (bufbox, srcproc, sendbuf, sendbufsize); buf->transfer_from_innerloop (srcs, times, dstbox, time, order_space, order_time); delete buf; @@ -241,6 +250,13 @@ transfer_from (comm_state & state, copy_from_innerloop (buf, dstbox); delete buf; } else { + if (cent == cell_centered) { + assert (all (dstbox.stride() == this->extent().stride())); + ivect const ioff = dstbox.lower() - this->extent().lower(); + ivect const is_centered = ioff % this->extent().stride() == 0; + ivect const ioffset = not is_centered; + assert (all (ioffset == 0)); + } gdata const * const null = NULL; vector bufs (ntimelevels, null); vector timebuf (ntimelevels); -- cgit v1.2.3