diff options
-rw-r--r-- | Carpet/CarpetLib/param.ccl | 5 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 68 |
2 files changed, 70 insertions, 3 deletions
diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index 8314a47c6..fa39bfed2 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -36,6 +36,11 @@ BOOLEAN omit_prolongation_points_when_restricting "Do not restrict to points whi { } "no" +INT proper_nesting_distance "Minimum distance (in grid points) between two level interfaces" STEERABLE=always +{ + 0:* :: "any non-negative value is fine; the default value is just a guess" +} 4 + INT print_memstats_every "Report periodically how much memory is used per process" STEERABLE=always diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 937f27410..f2b981fca 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -412,9 +412,6 @@ void dh::setup_refinement_restriction_boxes (dh::dboxes & box, if (omit_prolongation_points_when_restricting) { // remove what is sent during boundary prolongation const int pss = prolongation_stencil_size(); - // TODO: this needs to remove what is sent from other - // processors as well, not only what is sent from processor - // c for (int ccc=0; ccc<h.components(rl); ++ccc) { const dh::dboxes& box2 = boxes.at(ml).at(rl).at(ccc); for (iblistvect::const_iterator slvi = @@ -634,6 +631,71 @@ void dh::check_bboxes (dh::dboxes & box, } } } + + // Check proper nesting + // "Proper nesting" means that prolongation from level L to level + // L+1 does not use any points that are prolongated from level L-1 + // to level L. + // We extend that notion to require a certain distance D in between. + if (c == 0) { + if (rl > 0 and rl < h.reflevels()) { + // Points that are filled by prolongation from level rl-1 + ibset recvs; + for (int cc=0; cc<h.components(rl); ++cc) { + dboxes & box1 = boxes.at(ml).at(rl).at(cc); + for (iblistvect::const_iterator rlvi = box1.recv_ref_bnd_coarse.begin(); + rlvi != box1.recv_ref_bnd_coarse.end(); ++ rlvi) + { + const iblist & recvlist = * rlvi; + for (iblist::const_iterator rli = recvlist.begin(); + rli != recvlist.end(); ++ rli) + { + const ibbox & recv = * rli; + recvs |= recv; + } + } + } + recvs.normalize(); + // Extend the received points by the desired distance + const ivect dist = proper_nesting_distance; + ibset taboo; + for (ibset::const_iterator bi = recvs.begin(); bi != recvs.end(); ++ bi) { + const ibbox & b = * bi; + const ibbox t = b.expand (dist, dist); + taboo |= t; + } + taboo.normalize(); + // Points that are used for prolongation to level rl+1 + ibset sends; + for (int cc=0; cc<h.components(rl); ++cc) { + dboxes & box1 = boxes.at(ml).at(rl).at(cc); + for (iblistvect::const_iterator slvi = box1.send_ref_bnd_fine.begin(); + slvi != box1.send_ref_bnd_fine.end(); ++ slvi) + { + const iblist & sendlist = * slvi; + for (iblist::const_iterator sli = sendlist.begin(); + sli != sendlist.end(); ++ sli) + { + const ibbox & send = * sli; + sends |= send; + } + } + } + sends.normalize(); + // Calculate the overlap + ibset overlap = sends & taboo; + if (not overlap.empty()) { + overlap.normalize(); + cout << "Not properly nested: rl=" << rl << ", " + << "required distance=" << proper_nesting_distance << endl; + cout << "Received from level " << rl-1 << ": " << recvs << endl; + cout << "Taboo region on level " << rl << ": " << taboo << endl; + cout << "Sent to level " << rl+1 << ": " << sends << endl; + cout << "Overlap on level " << rl << ": " << overlap << endl; + CCTK_WARN (1, "Not properly nested"); + } + } + } } void dh::calculate_bases () |