diff options
author | Roland Haas <roland.haas@physics.gatech.edu> | 2012-06-04 12:51:35 -0700 |
---|---|---|
committer | Barry Wardell <barry.wardell@gmail.com> | 2012-09-11 18:23:22 +0100 |
commit | 1636f53b32fdfc44fc2289acc842e0347343a1c7 (patch) | |
tree | e25525db07000d2fe040fbf3f0fa14aa7589db39 /Carpet | |
parent | 37f6677d699278935fcb93d9c5d6f7b67b17db70 (diff) |
CarpetLib: clean up higher order restriciton some
* rename controlling parameter to use_higher_order_restriction
* introduce parameter restriction_order_space to control which operator
is used (currently order 1 and 3 are suppoted)
* include some comments on what the operator does
* change the way the restrictable region is computed in dh.cc/regrid to
be based on exterior.shrink(stencil_width) rather that the interior
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/CarpetLib/param.ccl | 8 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/data.cc | 36 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 16 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc | 9 |
4 files changed, 53 insertions, 16 deletions
diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index b0a825171..a259fff2d 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -33,10 +33,16 @@ BOOLEAN use_dgfe "Use DGFE operators instead of Lagrange operators" -BOOLEAN use_cc_o3 "Use third order cell centered restriction operators instead of first order" +BOOLEAN use_higher_order_restriction "Use third order cell centered restriction operators instead of first order" { } "no" +INT restriction_order_space "Order of restriction operator to use with use_higher_order_restriction" +{ + 1 :: "linear interpolation, this is Carpet's original implementation" + 3 :: "third order accurate restriction for grid functions where prolongation is not (W)ENO" +} 3 + BOOLEAN output_bboxes "Output bounding box information to the screen" STEERABLE=always diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 9a26eec8a..ba1a438ce 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -1141,17 +1141,35 @@ transfer_restrict (data const * const src, srcregbox, dstregbox, NULL); break; } - if (use_cc_o3 and + if (use_higher_order_restriction and transport_operator != op_WENO and transport_operator != op_ENO) { // HACK - // Don't use call_operator, because we parallelise ourselves - restrict_3d_cc_o3_rf2(static_cast <T const *> (src->storage()), - src->shape(), - static_cast <T *> (this->storage()), - this->shape(), - srcbox, - dstbox, - srcregbox, dstregbox, NULL); + switch (restriction_order_space) { + case 1: + // Don't use call_operator, because we parallelise ourselves + restrict_3d_cc_rf2(static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + srcregbox, dstregbox, NULL); + break; + case 3: + // Don't use call_operator, because we parallelise ourselves + restrict_3d_cc_o3_rf2(static_cast <T const *> (src->storage()), + src->shape(), + static_cast <T *> (this->storage()), + this->shape(), + srcbox, + dstbox, + srcregbox, dstregbox, NULL); + break; + default: + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "There is no restriction stencil with restriction_order_space==%d", restriction_order_space); + break; + } break; } // Don't use call_operator, because we parallelise ourselves diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 572fe2e75..1380a762a 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1034,7 +1034,7 @@ regrid (bool const do_init) full_dboxes const& obox = full_olevel.AT(oc); ibset needrecv = allrestricted & obox.owned; - if(use_cc_o3) { + if(use_higher_order_restriction) { // NOTE: change in behaviour (affects only outer boundaries I think)!!!! // NOTE: b/c of this we need a low-level sync after the restrict needrecv = allrestricted & obox.interior; @@ -1045,10 +1045,14 @@ regrid (bool const do_init) for (int c = 0; c < h.components(rl); ++ c) { full_dboxes const& box = full_level.AT(c); - // HORRIBLE HACK - ibbox const contracted_exterior = - use_cc_o3 ? box.interior.expand(ivect(int(h.refcent==cell_centered))).contracted_for(odomext) : - box.exterior.contracted_for(odomext); + // If we make a mistake expanding the domain of dependence here, it + // _should_ be caught be by the expand()ed is_contained_in(srcbbox) + // test in the actual operator. + int const shrink_by = + use_higher_order_restriction and (h.refcent == cell_centered) ? + restriction_order_space/2 : 0; + ibbox const contracted_exterior = + box.exterior.expand(ivect(-shrink_by)).contracted_for(odomext); ibset const ovlp = needrecv & contracted_exterior; for (ibset::const_iterator @@ -1058,7 +1062,7 @@ regrid (bool const do_init) ibbox const send = recv.expanded_for(box.exterior); ASSERT_c (send <= box.exterior, "Refinement restriction: Send region must be contained in exterior"); - if(use_cc_o3) { + if(use_higher_order_restriction) { ASSERT_c (send <= box.interior.expand(ivect(int(h.refcent==cell_centered))), "Refinement restriction: Send region must be contained in interior"); } diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc index c8ecc040f..991734c85 100644 --- a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc +++ b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc @@ -25,6 +25,15 @@ namespace CarpetLib { dstiext, dstjext, dstkext) +// This operator offers third-order accurate restriction operators for cell +// centered grids when use_higher_order_restriction is set. This interpolation +// is done for samples at the cell centres (so this is not a ppm scheme or +// anything like that). It really only fits a polynomial of the form +// f(x,y,z) = \sum_{i,j,k=0}^3 a_{i,j,k} x^i y^j z^k +// to the fine cells and evaluates at x=y=z=0. So it is good for the metric, +// but bad for matter (since it will destroy the conservation). Because of +// this it should not be used for grid functions whose transport operator is +// not WENO or ENO which hopefully excludes all matter variables. template <typename T> void |