From 7a5f190dbbc5b4fa7d6ef931e7feab87038ca428 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 6 Mar 2024 17:02:30 -0500 Subject: [PATCH 1/5] CarpetX: Test curvature instead of slope for ENO fallback condition --- CarpetX/src/prolongate_3d_rf2_impl.hxx | 121 ++++++++++++++++++------- 1 file changed, 90 insertions(+), 31 deletions(-) diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx index b2296ca99..202e58672 100644 --- a/CarpetX/src/prolongate_3d_rf2_impl.hxx +++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx @@ -1403,62 +1403,121 @@ void prolongate_3d_rf2< // the stencil if constexpr (INTPI == CONS || INTPI == ENO) { - bool need_fallback_i = false; + if constexpr (false) { + // Check whether slopes change sign + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + const CCTK_REAL s0 = + crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) - + crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk); + for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { + const CCTK_REAL s = crse(icrse[0] + di - 0, + icrse[1] + dj, icrse[2] + dk) - + crse(icrse[0] + di - 1, + icrse[1] + dj, icrse[2] + dk); + need_fallback |= s * s0 < 0; + } + } + } + } + // Check whether curvatures change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - const CCTK_REAL s0 = - crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) - - crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk); - for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - - crse(icrse[0] + (di - 1), icrse[1] + dj, - icrse[2] + dk); - need_fallback_i |= s * s0 < 0; + const CCTK_REAL c0 = + crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk) - + 2 * crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) + + crse(icrse[0] + 2, icrse[1] + dj, icrse[2] + dk); + for (int di = sradi[0] + 3; di <= sradi[1]; ++di) { + const CCTK_REAL c = + crse(icrse[0] + di - 2, icrse[1] + dj, + icrse[2] + dk) - + 2 * crse(icrse[0] + di - 1, icrse[1] + dj, + icrse[2] + dk) + + crse(icrse[0] + di - 0, icrse[1] + dj, icrse[2] + dk); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_i; } if constexpr (INTPJ == CONS || INTPJ == ENO) { - bool need_fallback_j = false; + if constexpr (false) { + // Check whether slopes change sign + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const CCTK_REAL s0 = + crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) - + crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk); + for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { + const CCTK_REAL s = + crse(icrse[0] + di, icrse[1] + dj - 0, + icrse[2] + dk) - + crse(icrse[0] + di, icrse[1] + dj - 1, + icrse[2] + dk); + need_fallback |= s * s0 < 0; + } + } + } + } + // Check whether curvatures change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) - - crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk); - for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - - crse(icrse[0] + di, icrse[1] + (dj - 1), - icrse[2] + dk); - need_fallback_j |= s * s0 < 0; + const CCTK_REAL c0 = + crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk) - + 2 * crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) + + crse(icrse[0] + di, icrse[1] + 2, icrse[2] + dk); + for (int dj = sradj[0] + 3; dj <= sradj[1]; ++dj) { + const CCTK_REAL c = + crse(icrse[0] + di, icrse[1] + dj - 2, + icrse[2] + dk) - + 2 * crse(icrse[0] + di, icrse[1] + dj - 1, + icrse[2] + dk) + + crse(icrse[0] + di, icrse[1] + dj - 0, icrse[2] + dk); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_j; } if constexpr (INTPK == CONS || INTPK == ENO) { - bool need_fallback_k = false; + if constexpr (false) { + // Check whether slopes change sign + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const CCTK_REAL s0 = + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) - + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0); + for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { + const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk - 0) - + crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk - 1); + need_fallback |= s * s0 < 0; + } + } + } + } + // Check whether curvatures change sign for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) - - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0); + const CCTK_REAL c0 = + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0) - + 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) + + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2); for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - + const CCTK_REAL c = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + (dk - 1)); - need_fallback_k |= s * s0 < 0; + icrse[2] + dk - 2) - + 2 * crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk - 1) + + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk - 0); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_k; } } + res = need_fallback ? val_lin : val; } // if FB != FB_NONE From 33d76b51b143ad7d74cca619a1881333f0f5e2ff Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 6 Mar 2024 19:28:58 -0500 Subject: [PATCH 2/5] CarpetX: Correct prolongation in per-group operators as well --- CarpetX/src/prolongate_3d_rf2_impl.hxx | 195 ++++++++++++++++--------- 1 file changed, 125 insertions(+), 70 deletions(-) diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx index 202e58672..610bd7b08 100644 --- a/CarpetX/src/prolongate_3d_rf2_impl.hxx +++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx @@ -1398,9 +1398,9 @@ void prolongate_3d_rf2< need_fallback |= val < minval || val > maxval; // Fallback condition 2 (checked for every direction in - // which the operator is conservative): The slope of the - // values at the stencil points changes sign anywhere in - // the stencil + // which the operator is conservative): The slope or + // curvature of the values at the stencil points changes + // sign anywhere in the stencil if constexpr (INTPI == CONS || INTPI == ENO) { if constexpr (false) { @@ -1834,15 +1834,19 @@ void prolongate_3d_rf2< std::array ress; - if constexpr (FB == FB_NONE || ((INTPI != CONS || ORDERI <= 1) && - (INTPJ != CONS || ORDERJ <= 1) && - (INTPK != CONS || ORDERK <= 1))) { + if constexpr (FB == FB_NONE || + (((INTPI != CONS && INTPI != ENO) || ORDERI <= 1) && + ((INTPJ != CONS && INTPJ != ENO) || ORDERJ <= 1) && + ((INTPK != CONS && INTPK != ENO) || ORDERK <= 1))) { + // We are not falling back to linear interpolation for (int comp = 0; comp < ncomps; ++comp) ress[comp] = vals[comp]; } else { + // We might want to fall back to linear interpolation + // Check whether we need to fall back bool need_fallback = false; { const std::array sradi = @@ -1854,6 +1858,8 @@ void prolongate_3d_rf2< const std::array sradk = interp1d().stencil_radius(shift[2], off[2]); + // Fallback condition 1: The interpolated value introduces a new + // extremum CCTK_REAL minval = +1 / CCTK_REAL(0), maxval = -1 / CCTK_REAL(0); for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { @@ -1870,102 +1876,151 @@ void prolongate_3d_rf2< need_fallback |= vals[comp] < minval || vals[comp] > maxval; } - if constexpr (INTPI == CONS) { - bool need_fallback_i = false; + // Fallback condition 2 (checked for every direction in + // which the operator is conservative): The slope or + // curvature of the values at the stencil points changes + // sign anywhere in the stencil + + if constexpr (INTPI == CONS || INTPI == ENO) { + if constexpr (false) { + // Check whether slopes change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + const CCTK_REAL s0 = crse(icrse[0] + 1, icrse[1] + dj, + icrse[2] + dk, comp) - + crse(icrse[0] + 0, icrse[1] + dj, + icrse[2] + dk, comp); + for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { + const CCTK_REAL s = + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, + comp) - + crse(icrse[0] + (di - 1), icrse[1] + dj, + icrse[2] + dk, comp); + need_fallback |= s * s0 < 0; + } + } + } + } + } + // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - // using std::signbit; - // const bool s0 = signbit( - // crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) - - // crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk)); - const CCTK_REAL s0 = - crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk, comp) - - crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk, comp); - for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { - // const bool s = signbit( - // crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - - // crse(icrse[0] + (di - 1), icrse[1] + dj, - // icrse[2] + dk)); - // need_fallback_i |= s != s0; - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, + const CCTK_REAL c0 = + crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk, comp) - + 2 * crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk, + comp) + + crse(icrse[0] + 2, icrse[1] + dj, icrse[2] + dk, comp); + for (int di = sradi[0] + 3; di <= sradi[1]; ++di) { + const CCTK_REAL c = + crse(icrse[0] + di - 2, icrse[1] + dj, icrse[2] + dk, comp) - - crse(icrse[0] + (di - 1), icrse[1] + dj, - icrse[2] + dk, comp); - need_fallback_i |= s * s0 < 0; + 2 * crse(icrse[0] + di - 1, icrse[1] + dj, + icrse[2] + dk, comp) + + crse(icrse[0] + di - 0, icrse[1] + dj, icrse[2] + dk, + comp); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_i; } } - if constexpr (INTPJ == CONS) { - bool need_fallback_j = false; + if constexpr (INTPJ == CONS || INTPJ == ENO) { + if constexpr (false) { + // Check whether slopes change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const CCTK_REAL s0 = crse(icrse[0] + di, icrse[1] + 1, + icrse[2] + dk, comp) - + crse(icrse[0] + di, icrse[1] + 0, + icrse[2] + dk, comp); + for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { + const CCTK_REAL s = + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, + comp) - + crse(icrse[0] + di, icrse[1] + (dj - 1), + icrse[2] + dk, comp); + need_fallback |= s * s0 < 0; + } + } + } + } + } + // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - // using std::signbit; - // const bool s0 = signbit( - // crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) - - // crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk)); - const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk, comp) - - crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk, comp); - for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { - // const bool s = signbit( - // crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - - // crse(icrse[0] + di, icrse[1] + (dj - 1), - // icrse[2] + dk)); - // need_fallback_j |= s != s0; - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, + const CCTK_REAL c0 = + crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk, comp) - + 2 * crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk, + comp) + + crse(icrse[0] + di, icrse[1] + 2, icrse[2] + dk, comp); + for (int dj = sradj[0] + 3; dj <= sradj[1]; ++dj) { + const CCTK_REAL c = + crse(icrse[0] + di, icrse[1] + dj - 2, icrse[2] + dk, comp) - - crse(icrse[0] + di, icrse[1] + (dj - 1), - icrse[2] + dk, comp); - need_fallback_j |= s * s0 < 0; + 2 * crse(icrse[0] + di, icrse[1] + dj - 1, + icrse[2] + dk, comp) + + crse(icrse[0] + di, icrse[1] + dj - 0, icrse[2] + dk, + comp); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_j; } } - if constexpr (INTPK == CONS) { - bool need_fallback_k = false; + if constexpr (INTPK == CONS || INTPK == ENO) { + if constexpr (false) { + // Check whether slopes change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const CCTK_REAL s0 = crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + 1, comp) - + crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + 0, comp); + for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { + const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk, comp) - + crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + (dk - 1), comp); + need_fallback |= s * s0 < 0; + } + } + } + } + } + // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - // using std::signbit; - // const bool s0 = signbit( - // crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) - - // crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0)); - const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1, comp) - - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0, comp); + const CCTK_REAL c0 = + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0, comp) - + 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1, + comp) + + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2, comp); for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { - // const bool s = signbit( - // crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk) - - // crse(icrse[0] + di, icrse[1] + dj, - // icrse[2] + (dk - 1))); - // need_fallback_k |= s != s0; - const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk, comp) - + const CCTK_REAL c = crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk - 2, comp) - + 2 * crse(icrse[0] + di, icrse[1] + dj, + icrse[2] + dk - 1, comp) + crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + (dk - 1), comp); - need_fallback_k |= s * s0 < 0; + icrse[2] + dk - 0, comp); + need_fallback |= c * c0 < 0; } } } - need_fallback |= need_fallback_k; } } } - constexpr int LINORDERI = INTPI == CONS ? 1 : ORDERI; - constexpr int LINORDERJ = INTPJ == CONS ? 1 : ORDERJ; - constexpr int LINORDERK = INTPK == CONS ? 1 : ORDERK; + constexpr int LINORDERI = INTPI == CONS || INTPI == ENO ? 1 : ORDERI; + constexpr int LINORDERJ = INTPJ == CONS || INTPJ == ENO ? 1 : ORDERJ; + constexpr int LINORDERK = INTPK == CONS || INTPK == ENO ? 1 : ORDERK; std::array val_lins; for (int comp = 0; comp < ncomps; ++comp) { val_lins[comp] = call_stencil_3d( From dbbe58de7cadce13917f38d972dcd14dca712de0 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 6 Mar 2024 19:54:39 -0500 Subject: [PATCH 3/5] CarpetX: Correct some loop bounds in fallback prolongation operators --- CarpetX/src/prolongate_3d_rf2_impl.hxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx index 610bd7b08..c404f398c 100644 --- a/CarpetX/src/prolongate_3d_rf2_impl.hxx +++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx @@ -1504,7 +1504,7 @@ void prolongate_3d_rf2< crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0) - 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2); - for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { + for (int dk = sradk[0] + 3; dk <= sradk[1]; ++dk) { const CCTK_REAL c = crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk - 2) - @@ -1893,9 +1893,9 @@ void prolongate_3d_rf2< icrse[2] + dk, comp); for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, - comp) - - crse(icrse[0] + (di - 1), icrse[1] + dj, + crse(icrse[0] + di - 0, icrse[1] + dj, + icrse[2] + dk, comp) - + crse(icrse[0] + di - 1, icrse[1] + dj, icrse[2] + dk, comp); need_fallback |= s * s0 < 0; } @@ -1939,9 +1939,9 @@ void prolongate_3d_rf2< icrse[2] + dk, comp); for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, - comp) - - crse(icrse[0] + di, icrse[1] + (dj - 1), + crse(icrse[0] + di, icrse[1] + dj - 0, + icrse[2] + dk, comp) - + crse(icrse[0] + di, icrse[1] + dj - 1, icrse[2] + dk, comp); need_fallback |= s * s0 < 0; } @@ -1985,9 +1985,9 @@ void prolongate_3d_rf2< icrse[2] + 0, comp); for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk, comp) - + icrse[2] + dk - 0, comp) - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + (dk - 1), comp); + icrse[2] + dk - 1, comp); need_fallback |= s * s0 < 0; } } @@ -2003,7 +2003,7 @@ void prolongate_3d_rf2< 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1, comp) + crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2, comp); - for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { + for (int dk = sradk[0] + 3; dk <= sradk[1]; ++dk) { const CCTK_REAL c = crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk - 2, comp) - 2 * crse(icrse[0] + di, icrse[1] + dj, From 8ef9cf747ade6594896b09e67427a48714e5c5fa Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 6 Mar 2024 20:32:08 -0500 Subject: [PATCH 4/5] CarpetX: Correct indexing in fallback check operators --- CarpetX/src/prolongate_3d_rf2_impl.hxx | 343 +++++++++++-------------- 1 file changed, 154 insertions(+), 189 deletions(-) diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx index c404f398c..4a388857d 100644 --- a/CarpetX/src/prolongate_3d_rf2_impl.hxx +++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx @@ -1217,6 +1217,10 @@ void prolongate_3d_rf2< const vect icrse = ifine >> 1; const vect off = ifine & 0x1; + const auto pcrse = [&](const int di, const int dj, const int dk) { + return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk); + }; + // Choose stencil shift vect shift{0, 0, 0}; constexpr bool any_use_shift = any(use_shift); @@ -1260,9 +1264,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int di) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk); + return pcrse(si + di, sj + dj, sk + dk); }); ddx = fmax(ddx, fabs(dd)); } @@ -1279,9 +1281,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int dj) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk); + return pcrse(si + di, sj + dj, sk + dk); }); ddy = fmax(ddy, fabs(dd)); } @@ -1298,9 +1298,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int dk) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk); + return pcrse(si + di, sj + dj, sk + dk); }); ddz = fmax(ddz, fabs(dd)); } @@ -1322,9 +1320,10 @@ void prolongate_3d_rf2< } // if any_use_shift const CCTK_REAL val = call_stencil_3d( - [&](const int di, const int dj, const int dk) { - return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk); - }, + // [&](const int di, const int dj, const int dk) { + // return pcrse(di, dj, dk); + // }, + pcrse, [&](const auto &crse) { return interp1d()(crse, shift[0], off[0]); }, @@ -1356,9 +1355,10 @@ void prolongate_3d_rf2< constexpr int LINORDERK = INTPK == CONS || INTPK == ENO ? 1 : ORDERK; const CCTK_REAL val_lin = call_stencil_3d( - [&](const int di, const int dj, const int dk) { - return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk); - }, + // [&](const int di, const int dj, const int dk) { + // return pcrse(di, dj, dk); + // }, + pcrse, [&](const auto &crse) { return interp1d()(crse, 0, off[0]); }, @@ -1388,10 +1388,8 @@ void prolongate_3d_rf2< for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { using std::fmax, std::fmin; - minval = fmin(minval, crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk)); - maxval = fmax(maxval, crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk)); + minval = fmin(minval, pcrse(di, dj, dk)); + maxval = fmax(maxval, pcrse(di, dj, dk)); } } } @@ -1403,115 +1401,105 @@ void prolongate_3d_rf2< // sign anywhere in the stencil if constexpr (INTPI == CONS || INTPI == ENO) { - if constexpr (false) { + if (false && sradi[1] - sradi[0] > 2) { // Check whether slopes change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + const int di = sradi[0]; const CCTK_REAL s0 = - crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) - - crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk); - for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { - const CCTK_REAL s = crse(icrse[0] + di - 0, - icrse[1] + dj, icrse[2] + dk) - - crse(icrse[0] + di - 1, - icrse[1] + dj, icrse[2] + dk); + pcrse(di + 1, dj, dk) - pcrse(di + 0, dj, dk); + for (int di = sradi[0] + 1; di <= sradi[1] - 1; ++di) { + const CCTK_REAL s = + pcrse(di + 1, dj, dk) - pcrse(di + 0, dj, dk); need_fallback |= s * s0 < 0; } } } } - // Check whether curvatures change sign - for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { - for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - const CCTK_REAL c0 = - crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk) - - 2 * crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk) + - crse(icrse[0] + 2, icrse[1] + dj, icrse[2] + dk); - for (int di = sradi[0] + 3; di <= sradi[1]; ++di) { - const CCTK_REAL c = - crse(icrse[0] + di - 2, icrse[1] + dj, - icrse[2] + dk) - - 2 * crse(icrse[0] + di - 1, icrse[1] + dj, - icrse[2] + dk) + - crse(icrse[0] + di - 0, icrse[1] + dj, icrse[2] + dk); - need_fallback |= c * c0 < 0; + if (sradi[1] - sradi[0] > 3) { + // Check whether curvatures change sign + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + const int di = sradi[0]; + const CCTK_REAL c0 = pcrse(di + 0, dj, dk) - + 2 * pcrse(di + 1, dj, dk) + + pcrse(di + 2, dj, dk); + for (int di = sradi[0] + 1; di <= sradi[1] - 2; ++di) { + const CCTK_REAL c = pcrse(di + 0, dj, dk) - + 2 * pcrse(di + 1, dj, dk) + + pcrse(di + 2, dj, dk); + need_fallback |= c * c0 < 0; + } } } } } if constexpr (INTPJ == CONS || INTPJ == ENO) { - if constexpr (false) { + if (false && sradj[1] - sradj[0] > 2) { // Check whether slopes change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dj = sradj[0]; const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) - - crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk); - for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { + pcrse(di, dj + 1, dk) - pcrse(di, dj + 0, dk); + for (int dj = sradj[0] + 1; dj <= sradj[1] - 1; ++dj) { const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj - 0, - icrse[2] + dk) - - crse(icrse[0] + di, icrse[1] + dj - 1, - icrse[2] + dk); + pcrse(di, dj + 1, dk) - pcrse(di, dj + 0, dk); need_fallback |= s * s0 < 0; } } } } - // Check whether curvatures change sign - for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { - for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL c0 = - crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk) - - 2 * crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk) + - crse(icrse[0] + di, icrse[1] + 2, icrse[2] + dk); - for (int dj = sradj[0] + 3; dj <= sradj[1]; ++dj) { - const CCTK_REAL c = - crse(icrse[0] + di, icrse[1] + dj - 2, - icrse[2] + dk) - - 2 * crse(icrse[0] + di, icrse[1] + dj - 1, - icrse[2] + dk) + - crse(icrse[0] + di, icrse[1] + dj - 0, icrse[2] + dk); - need_fallback |= c * c0 < 0; + if (sradj[1] - sradj[0] > 3) { + // Check whether curvatures change sign + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dj = sradj[0]; + const CCTK_REAL c0 = pcrse(di, dj + 0, dk) - + 2 * pcrse(di, dj + 1, dk) + + pcrse(di, dj + 2, dk); + for (int dj = sradj[0] + 1; dj <= sradj[1] - 2; ++dj) { + const CCTK_REAL c = pcrse(di, dj + 0, dk) - + 2 * pcrse(di, dj + 1, dk) + + pcrse(di, dj + 2, dk); + need_fallback |= c * c0 < 0; + } } } } } if constexpr (INTPK == CONS || INTPK == ENO) { - if constexpr (false) { + if (false && sradk[1] - sradk[0] > 2) { // Check whether slopes change sign for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dk = sradk[0]; const CCTK_REAL s0 = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) - - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0); - for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { - const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 0) - - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 1); + pcrse(di, dj, dk + 1) - pcrse(di, dj, dk + 0); + for (int dk = sradk[0] + 1; dk <= sradk[1] - 1; ++dk) { + const CCTK_REAL s = + pcrse(di, dj, dk + 1) - pcrse(di, dj, dk + 0); need_fallback |= s * s0 < 0; } } } } - // Check whether curvatures change sign - for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL c0 = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0) - - 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1) + - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2); - for (int dk = sradk[0] + 3; dk <= sradk[1]; ++dk) { - const CCTK_REAL c = - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 2) - - 2 * crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 1) + - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk - 0); - need_fallback |= c * c0 < 0; + if (sradk[1] - sradk[0] > 3) { + // Check whether curvatures change sign + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dk = sradk[0]; + const CCTK_REAL c0 = pcrse(di, dj, dk + 0) - + 2 * pcrse(di, dj, dk + 1) + + pcrse(di, dj, dk + 2); + for (int dk = sradk[0] + 1; dk <= sradk[1] - 2; ++dk) { + const CCTK_REAL c = pcrse(di, dj, dk + 0) - + 2 * pcrse(di, dj, dk + 1) + + pcrse(di, dj, dk + 2); + need_fallback |= c * c0 < 0; + } } } } @@ -1706,6 +1694,11 @@ void prolongate_3d_rf2< const vect icrse = ifine >> 1; const vect off = ifine & 0x1; + const auto pcrse = [&](const int di, const int dj, const int dk, + const int comp) { + return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, comp); + }; + // Choose stencil shift vect shift{0, 0, 0}; constexpr bool any_use_shift = any(use_shift); @@ -1750,9 +1743,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int di) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk, comp); + return pcrse(si + di, sj + dj, sk + dk, comp); }); ddx = fmax(ddx, fabs(dd)); } @@ -1771,9 +1762,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int dj) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk, comp); + return pcrse(si + di, sj + dj, sk + dk, comp); }); ddy = fmax(ddy, fabs(dd)); } @@ -1792,9 +1781,7 @@ void prolongate_3d_rf2< const CCTK_REAL dd = undivided_difference_1d()( [&](const int dk) { - return crse(icrse[0] + si + di, - icrse[1] + sj + dj, - icrse[2] + sk + dk, comp); + return pcrse(si + di, sj + dj, sk + dk, comp); }); ddz = fmax(ddz, fabs(dd)); } @@ -1820,7 +1807,7 @@ void prolongate_3d_rf2< for (int comp = 0; comp < ncomps; ++comp) vals[comp] = call_stencil_3d( [&](const int di, const int dj, const int dk) { - return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, comp); + return pcrse(di, dj, dk, comp); }, [&](const auto &crse) { return interp1d()(crse, shift[0], off[0]); @@ -1866,10 +1853,8 @@ void prolongate_3d_rf2< for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { using std::fmax, std::fmin; - minval = fmin(minval, crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk, comp)); - maxval = fmax(maxval, crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk, comp)); + minval = fmin(minval, pcrse(di, dj, dk, comp)); + maxval = fmax(maxval, pcrse(di, dj, dk, comp)); } } } @@ -1882,45 +1867,38 @@ void prolongate_3d_rf2< // sign anywhere in the stencil if constexpr (INTPI == CONS || INTPI == ENO) { - if constexpr (false) { + if (false && sradi[1] - sradi[0] > 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - const CCTK_REAL s0 = crse(icrse[0] + 1, icrse[1] + dj, - icrse[2] + dk, comp) - - crse(icrse[0] + 0, icrse[1] + dj, - icrse[2] + dk, comp); - for (int di = sradi[0] + 2; di <= sradi[1]; ++di) { - const CCTK_REAL s = - crse(icrse[0] + di - 0, icrse[1] + dj, - icrse[2] + dk, comp) - - crse(icrse[0] + di - 1, icrse[1] + dj, - icrse[2] + dk, comp); + const int di = sradi[0]; + const CCTK_REAL s0 = pcrse(di + 1, dj, dk, comp) - + pcrse(di + 0, dj, dk, comp); + for (int di = sradi[0] + 1; di <= sradi[1] - 1; ++di) { + const CCTK_REAL s = pcrse(di + 1, dj, dk, comp) - + pcrse(di + 0, dj, dk, comp); need_fallback |= s * s0 < 0; } } } } } - // Check whether curvatures change sign - for (int comp = 0; comp < ncomps; ++comp) { - for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { - for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - const CCTK_REAL c0 = - crse(icrse[0] + 0, icrse[1] + dj, icrse[2] + dk, comp) - - 2 * crse(icrse[0] + 1, icrse[1] + dj, icrse[2] + dk, - comp) + - crse(icrse[0] + 2, icrse[1] + dj, icrse[2] + dk, comp); - for (int di = sradi[0] + 3; di <= sradi[1]; ++di) { - const CCTK_REAL c = - crse(icrse[0] + di - 2, icrse[1] + dj, icrse[2] + dk, - comp) - - 2 * crse(icrse[0] + di - 1, icrse[1] + dj, - icrse[2] + dk, comp) + - crse(icrse[0] + di - 0, icrse[1] + dj, icrse[2] + dk, - comp); - need_fallback |= c * c0 < 0; + if (sradi[1] - sradi[0] > 3) { + // Check whether curvatures change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + const int di = sradi[0]; + const CCTK_REAL c0 = pcrse(di + 0, dj, dk, comp) - + 2 * pcrse(di + 1, dj, dk, comp) + + pcrse(di + 2, dj, dk, comp); + for (int di = sradi[0] + 1; di <= sradi[1] - 2; ++di) { + const CCTK_REAL c = pcrse(di + 0, dj, dk, comp) - + 2 * pcrse(di + 1, dj, dk, comp) + + pcrse(di + 2, dj, dk, comp); + need_fallback |= c * c0 < 0; + } } } } @@ -1928,45 +1906,38 @@ void prolongate_3d_rf2< } if constexpr (INTPJ == CONS || INTPJ == ENO) { - if constexpr (false) { + if (false && sradj[1] - sradj[0] > 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL s0 = crse(icrse[0] + di, icrse[1] + 1, - icrse[2] + dk, comp) - - crse(icrse[0] + di, icrse[1] + 0, - icrse[2] + dk, comp); - for (int dj = sradj[0] + 2; dj <= sradj[1]; ++dj) { - const CCTK_REAL s = - crse(icrse[0] + di, icrse[1] + dj - 0, - icrse[2] + dk, comp) - - crse(icrse[0] + di, icrse[1] + dj - 1, - icrse[2] + dk, comp); + const int dj = sradj[0]; + const CCTK_REAL s0 = pcrse(di, dj + 1, dk, comp) - + pcrse(di, dj + 0, dk, comp); + for (int dj = sradj[0] + 1; dj <= sradj[1] - 1; ++dj) { + const CCTK_REAL s = pcrse(di, dj + 1, dk, comp) - + pcrse(di, dj + 0, dk, comp); need_fallback |= s * s0 < 0; } } } } } - // Check whether curvatures change sign - for (int comp = 0; comp < ncomps; ++comp) { - for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { - for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL c0 = - crse(icrse[0] + di, icrse[1] + 0, icrse[2] + dk, comp) - - 2 * crse(icrse[0] + di, icrse[1] + 1, icrse[2] + dk, - comp) + - crse(icrse[0] + di, icrse[1] + 2, icrse[2] + dk, comp); - for (int dj = sradj[0] + 3; dj <= sradj[1]; ++dj) { - const CCTK_REAL c = - crse(icrse[0] + di, icrse[1] + dj - 2, icrse[2] + dk, - comp) - - 2 * crse(icrse[0] + di, icrse[1] + dj - 1, - icrse[2] + dk, comp) + - crse(icrse[0] + di, icrse[1] + dj - 0, icrse[2] + dk, - comp); - need_fallback |= c * c0 < 0; + if (sradj[1] - sradj[0] > 3) { + // Check whether curvatures change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dj = sradj[0]; + const CCTK_REAL c0 = pcrse(di, dj + 0, dk, comp) - + 2 * pcrse(di, dj + 1, dk, comp) + + pcrse(di, dj + 2, dk, comp); + for (int dj = sradj[0] + 1; dj <= sradj[1] - 2; ++dj) { + const CCTK_REAL c = pcrse(di, dj + 0, dk, comp) - + 2 * pcrse(di, dj + 1, dk, comp) + + pcrse(di, dj + 2, dk, comp); + need_fallback |= c * c0 < 0; + } } } } @@ -1974,43 +1945,38 @@ void prolongate_3d_rf2< } if constexpr (INTPK == CONS || INTPK == ENO) { - if constexpr (false) { + if (false && sradk[1] - sradk[0] > 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL s0 = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + 1, comp) - - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + 0, comp); - for (int dk = sradk[0] + 2; dk <= sradk[1]; ++dk) { - const CCTK_REAL s = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 0, comp) - - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 1, comp); + const int dk = sradk[0]; + const CCTK_REAL s0 = pcrse(di, dj, dk + 1, comp) - + pcrse(di, dj, dk + 0, comp); + for (int dk = sradk[0] + 1; dk <= sradk[1] - 1; ++dk) { + const CCTK_REAL s = pcrse(di, dj, dk + 1, comp) - + pcrse(di, dj, dk + 0, comp); need_fallback |= s * s0 < 0; } } } } } - // Check whether curvatures change sign - for (int comp = 0; comp < ncomps; ++comp) { - for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { - for (int di = sradi[0]; di <= sradi[1]; ++di) { - const CCTK_REAL c0 = - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 0, comp) - - 2 * crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 1, - comp) + - crse(icrse[0] + di, icrse[1] + dj, icrse[2] + 2, comp); - for (int dk = sradk[0] + 3; dk <= sradk[1]; ++dk) { - const CCTK_REAL c = crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 2, comp) - - 2 * crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 1, comp) + - crse(icrse[0] + di, icrse[1] + dj, - icrse[2] + dk - 0, comp); - need_fallback |= c * c0 < 0; + if (sradk[1] - sradk[0] > 3) { + // Check whether curvatures change sign + for (int comp = 0; comp < ncomps; ++comp) { + for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { + for (int di = sradi[0]; di <= sradi[1]; ++di) { + const int dk = sradk[0]; + const CCTK_REAL c0 = pcrse(di, dj, dk + 0, comp) - + 2 * pcrse(di, dj, dk + 1, comp) + + pcrse(di, dj, dk + 2, comp); + for (int dk = sradk[0] + 1; dk <= sradk[1] - 2; ++dk) { + const CCTK_REAL c = pcrse(di, dj, dk + 0, comp) - + 2 * pcrse(di, dj, dk + 1, comp) + + pcrse(di, dj, dk + 2, comp); + need_fallback |= c * c0 < 0; + } } } } @@ -2025,8 +1991,7 @@ void prolongate_3d_rf2< for (int comp = 0; comp < ncomps; ++comp) { val_lins[comp] = call_stencil_3d( [&](const int di, const int dj, const int dk) { - return crse(icrse[0] + di, icrse[1] + dj, icrse[2] + dk, - comp); + return pcrse(di, dj, dk, comp); }, [&](const auto &crse) { return interp1d()(crse, 0, off[0]); From a9d937d84af55ab01c1e351a2bea019fb6a29d46 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Thu, 7 Mar 2024 10:48:51 -0500 Subject: [PATCH 5/5] CarpetX: Correct fallback curvature tests --- CarpetX/src/prolongate_3d_rf2_impl.hxx | 32 +++++++++++--------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/CarpetX/src/prolongate_3d_rf2_impl.hxx b/CarpetX/src/prolongate_3d_rf2_impl.hxx index 4a388857d..08d769920 100644 --- a/CarpetX/src/prolongate_3d_rf2_impl.hxx +++ b/CarpetX/src/prolongate_3d_rf2_impl.hxx @@ -1320,9 +1320,6 @@ void prolongate_3d_rf2< } // if any_use_shift const CCTK_REAL val = call_stencil_3d( - // [&](const int di, const int dj, const int dk) { - // return pcrse(di, dj, dk); - // }, pcrse, [&](const auto &crse) { return interp1d()(crse, shift[0], off[0]); @@ -1355,9 +1352,6 @@ void prolongate_3d_rf2< constexpr int LINORDERK = INTPK == CONS || INTPK == ENO ? 1 : ORDERK; const CCTK_REAL val_lin = call_stencil_3d( - // [&](const int di, const int dj, const int dk) { - // return pcrse(di, dj, dk); - // }, pcrse, [&](const auto &crse) { return interp1d()(crse, 0, off[0]); @@ -1401,7 +1395,7 @@ void prolongate_3d_rf2< // sign anywhere in the stencil if constexpr (INTPI == CONS || INTPI == ENO) { - if (false && sradi[1] - sradi[0] > 2) { + if (false && sradi[1] - sradi[0] >= 2) { // Check whether slopes change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { @@ -1416,7 +1410,7 @@ void prolongate_3d_rf2< } } } - if (sradi[1] - sradi[0] > 3) { + if (sradi[1] - sradi[0] >= 3) { // Check whether curvatures change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { @@ -1436,7 +1430,7 @@ void prolongate_3d_rf2< } if constexpr (INTPJ == CONS || INTPJ == ENO) { - if (false && sradj[1] - sradj[0] > 2) { + if (false && sradj[1] - sradj[0] >= 2) { // Check whether slopes change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { @@ -1451,7 +1445,7 @@ void prolongate_3d_rf2< } } } - if (sradj[1] - sradj[0] > 3) { + if (sradj[1] - sradj[0] >= 3) { // Check whether curvatures change sign for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int di = sradi[0]; di <= sradi[1]; ++di) { @@ -1471,7 +1465,7 @@ void prolongate_3d_rf2< } if constexpr (INTPK == CONS || INTPK == ENO) { - if (false && sradk[1] - sradk[0] > 2) { + if (false && sradk[1] - sradk[0] >= 2) { // Check whether slopes change sign for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { @@ -1486,7 +1480,7 @@ void prolongate_3d_rf2< } } } - if (sradk[1] - sradk[0] > 3) { + if (sradk[1] - sradk[0] >= 3) { // Check whether curvatures change sign for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { @@ -1847,8 +1841,8 @@ void prolongate_3d_rf2< off[2]); // Fallback condition 1: The interpolated value introduces a new // extremum - CCTK_REAL minval = +1 / CCTK_REAL(0), maxval = -1 / CCTK_REAL(0); for (int comp = 0; comp < ncomps; ++comp) { + CCTK_REAL minval = +1 / CCTK_REAL(0), maxval = -1 / CCTK_REAL(0); for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { for (int di = sradi[0]; di <= sradi[1]; ++di) { @@ -1867,7 +1861,7 @@ void prolongate_3d_rf2< // sign anywhere in the stencil if constexpr (INTPI == CONS || INTPI == ENO) { - if (false && sradi[1] - sradi[0] > 2) { + if (false && sradi[1] - sradi[0] >= 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { @@ -1884,7 +1878,7 @@ void prolongate_3d_rf2< } } } - if (sradi[1] - sradi[0] > 3) { + if (sradi[1] - sradi[0] >= 3) { // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { @@ -1906,7 +1900,7 @@ void prolongate_3d_rf2< } if constexpr (INTPJ == CONS || INTPJ == ENO) { - if (false && sradj[1] - sradj[0] > 2) { + if (false && sradj[1] - sradj[0] >= 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { @@ -1923,7 +1917,7 @@ void prolongate_3d_rf2< } } } - if (sradj[1] - sradj[0] > 3) { + if (sradj[1] - sradj[0] >= 3) { // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dk = sradk[0]; dk <= sradk[1]; ++dk) { @@ -1945,7 +1939,7 @@ void prolongate_3d_rf2< } if constexpr (INTPK == CONS || INTPK == ENO) { - if (false && sradk[1] - sradk[0] > 2) { + if (false && sradk[1] - sradk[0] >= 2) { // Check whether slopes change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) { @@ -1962,7 +1956,7 @@ void prolongate_3d_rf2< } } } - if (sradk[1] - sradk[0] > 3) { + if (sradk[1] - sradk[0] >= 3) { // Check whether curvatures change sign for (int comp = 0; comp < ncomps; ++comp) { for (int dj = sradj[0]; dj <= sradj[1]; ++dj) {