From 539e5cecb0bc4c41694389160c2790660212f25a Mon Sep 17 00:00:00 2001 From: Christian Reisswig Date: Sat, 11 Jun 2011 11:27:23 -0700 Subject: Select 2nd order cc-eno for 3rd order interpolation since stencil radius is larger. This is consistent with vertex-centered eno. The trouble is that Carpet derives the stencil radius from the interpolation order imposing Lagrange interpolation. I also commented out the code that checks whether we have to switch down to first order. According to Shu, this is not necessary! This way interpolation runs faster. --- Carpet/CarpetLib/src/data.cc | 19 +++++++++++-- Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc | 36 ++++++++++++------------ 2 files changed, 34 insertions(+), 21 deletions(-) (limited to 'Carpet/CarpetLib') diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index cac634dde..58bb2adf1 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -836,12 +836,25 @@ transfer_prolongate (data const * const src, ibbox3 const & restrict regbbox) = { & prolongate_3d_cc_eno_rf2, - & prolongate_3d_cc_eno_rf2 + & prolongate_3d_cc_eno_rf2, // note that we cheat here: order is still 2 even though 3 was requested! + & prolongate_3d_cc_eno_rf2, // note that we cheat here: order is 3 even though 4 was requested! + & prolongate_3d_cc_eno_rf2 // note that we cheat here: order is 3 even though 5 was requested! + // have cheated here for two reasons: first, the ENO prolongation operator stencil radius is larger than Lagrange (and dh.cc assumes that the stencil goes as order_space/2!), + // and second, we want to allow spacetime interpolation to be of higher order while keeping the implemeneted ENO order! }; - if (order_space < 2 or order_space > 3) { + if (order_space < 2 or order_space > 5) { CCTK_WARN (CCTK_WARN_ABORT, - "There is no cell-centred stencil for op=\"ENO\" with order_space not in {2,3}"); + "There is no cell-centred stencil for op=\"ENO\" with order_space not in {2,3,4,5}"); } + /*if (order_space == 3) { + CCTK_WARN (1, + "Currently, in order to keep a consistent stencil radius in the buffer zone, we use 2nd order ENO. For third order ENO set order_space to 4 or 5!"); + } + if (order_space > 3) { + CCTK_WARN (1, + "Currently, in order to keep a consistent stencil radius in the buffer zone, we use 3rd order ENO. There are no higher orders than 3rd order implemented!"); + }*/ + call_operator (the_operators[order_space-2], static_cast (src->storage()), src->shape(), diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc index 0ebb2f03d..aefdeb5d7 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc +++ b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc @@ -289,7 +289,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (lV*rV <= 0) + /*if (lV*rV <= 0) { // switch back to first order! res = 0; @@ -298,7 +298,7 @@ namespace CarpetLib { res += coeffs1::get(i) * f[i-coeffs1d::minimin]; } break; - } + }*/ if (abs(lV) < abs(rV)) { // use left-shifted stencil since it is smoother @@ -338,7 +338,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] <= 0) + /*if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -347,7 +347,7 @@ namespace CarpetLib { res += coeffs1::get(i) * f[i-coeffs1d::minimin]; } break; - } + }*/ int min = 0; for (int i=0; i < 3; ++i) @@ -379,7 +379,7 @@ namespace CarpetLib { } // check that result is reasonable! - if ((res - f[-coeffs1d::minimin-1+di]) * (f[-coeffs1d::minimin+di] - res) < 0) + /*if ((res - f[-coeffs1d::minimin-1+di]) * (f[-coeffs1d::minimin+di] - res) < 0) { res = 0; // switch back to first order @@ -387,7 +387,7 @@ namespace CarpetLib { for (ptrdiff_t i=coeffs1::imin; i::minimin]; } - } + }*/ /*typedef coeffs1d coeffs; for (ptrdiff_t i=coeffs::imin; i::minimin]; } break; - } + }*/ if (abs(lV) < abs(rV)) { // use left-shifted stencil since it is smoother @@ -496,7 +496,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] <= 0) + /*if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -505,7 +505,7 @@ namespace CarpetLib { res += coeffs1::get(i) * f[i-coeffs1d::minimin]; } break; - } + }*/ int min = 0; for (int i=0; i < 3; ++i) @@ -537,7 +537,7 @@ namespace CarpetLib { } // check that result is reasonable! - if ((res - f[-coeffs1d::minimin-1+dj]) * (f[-coeffs1d::minimin+dj] - res) < 0) + /*if ((res - f[-coeffs1d::minimin-1+dj]) * (f[-coeffs1d::minimin+dj] - res) < 0) { res = 0; // switch back to first order @@ -545,7 +545,7 @@ namespace CarpetLib { for (ptrdiff_t i=coeffs1::imin; i::minimin]; } - } + }*/ /*typedef coeffs1d coeffs; for (ptrdiff_t i=coeffs::imin; i::minimin]; } break; - } + }*/ if (abs(lV) < abs(rV)) { // use left-shifted stencil since it is smoother @@ -656,7 +656,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] <= 0) + /*if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -665,7 +665,7 @@ namespace CarpetLib { res += coeffs1::get(i) * f[i-coeffs1d::minimin]; } break; - } + }*/ int min = 0; for (int i=0; i < 3; ++i) @@ -697,7 +697,7 @@ namespace CarpetLib { } // check that result is reasonable! - if ((res - f[-coeffs1d::minimin-1+dk]) * (f[-coeffs1d::minimin+dk] - res) < 0) + /*if ((res - f[-coeffs1d::minimin-1+dk]) * (f[-coeffs1d::minimin+dk] - res) < 0) { res = 0; // switch back to first order @@ -705,7 +705,7 @@ namespace CarpetLib { for (ptrdiff_t i=coeffs1::imin; i::minimin]; } - } + }*/ /*typedef coeffs1d coeffs; for (ptrdiff_t i=coeffs::imin; i