From a08b4adad53a1b60c52c69e84d1472c2a4c8a916 Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 9 Nov 2026 16:27:38 +0000 Subject: [PATCH 1/2] Add VanAlbadaSmooth slope limiter Smoothly tends to 2nd order when cell differences are very small. Based on Van Albada et al 1982 but with the addition of a clipping function present in our other limiters which flattens peaks and dips. The clipping is also subject to the tolerance and is smoothly disabled at very small cell differences. The tolerance is set to 1e-3 of the RMS of the three cell values used in the limiter calculation. This prevents huge performance issues due to Newton noise causing limiter reconstruction chatter at stagnation points. Original VA limiter is preserved for now. --- include/bout/fv_ops.hxx | 1 + include/bout/fv_ops_impl.hxx | 67 +++++++++++++++++++ .../user_docs/differential_operators.rst | 11 +++ src/mesh/fv_ops.cxx | 18 +++++ 4 files changed, 97 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 75c60abdc2..638fe57b5a 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -90,6 +90,7 @@ struct MinMod; struct MC; struct Superbee; struct VanAlbada; +struct VanAlbadaSmooth; struct WENO3; /*! diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 065e467b57..6d62288463 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -21,6 +21,7 @@ #include "bout/vector2d.hxx" #include +#include #include namespace FV { @@ -215,6 +216,72 @@ struct VanAlbada { } }; +/*! + * Smooth symmetric Van Albada second order slope limiter + * + * Uses the averaging formula (Eq. 37) of van Albada, van Leer and Roberts + * (1982), A&A 108, 76, multiplied by a smooth clip that sets the slope to + * zero at sharp extrema: + * + * slope = clip(dl*dr) * ((dr^2 + e2)*dl + (dl^2 + e2)*dr) / (dl^2 + dr^2 + 2*e2) + * + * where `dl = c - m`, `dr = p - c` and `e2 = (rel_eps * RMS(m, c, p))^2`. + * When both differences are small compared to sqrt(e2) the slope tends to + * the central difference, so the limiter is differentiable everywhere. + * This helps nonlinear solvers. + * + * Unlike the paper, e2 does not scale with the grid spacing, and it depends + * on any constant offset in the field. The limiter is not TVD at small + * extrema, where it can overshoot by about rel_eps * RMS. + */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx +struct VanAlbadaSmooth { + explicit VanAlbadaSmooth(BoutReal rel_eps = 1e-3) : rel_eps(rel_eps) {} + + BoutReal rel_eps; + + void operator()(Stencil1D& n) const { + const BoutReal dl = n.c - n.m; + const BoutReal dr = n.p - n.c; + + const BoutReal e2 = tolerance2(rel_eps, n.m, n.c, n.p); + const BoutReal slope = smoothClip(dl, dr, e2) * average(dl, dr, e2); + + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } + +private: + /// Eq. 37 of van Albada et al. (1982). Tends to (a + b)/2 when a and b + /// are small compared to sqrt(e2) + static BoutReal average(BoutReal a, BoutReal b, BoutReal e2) { + return ((((b * b) + e2) * a) + (((a * a) + e2) * b)) + / ((a * a) + (b * b) + (2.0 * e2) + 1e-30); + } + + /// Smooth clip: 1 if a*b >= 0, 1/2 if a*b = -e2, and tends to 0 for + /// a*b << -e2. The ratio is squared, rather than a*b and e2 separately, + /// to avoid overflow. The 1e-150 avoids division by zero if e2 = 0 + static BoutReal smoothClip(BoutReal a, BoutReal b, BoutReal e2) { + const BoutReal ratio = std::min(a * b, 0.0) / (e2 + 1e-150); + return 1.0 / (1.0 + (ratio * ratio)); + } + + /// (rel_eps * RMS(m, c, p))^2, normalised by the largest absolute value + /// to avoid overflow + static BoutReal tolerance2(BoutReal rel_eps, BoutReal m, BoutReal c, BoutReal p) { + const BoutReal scale = std::max({std::abs(m), std::abs(c), std::abs(p)}); + if (rel_eps == 0.0 or scale == 0.0) { + return 0.0; + } + const BoutReal mean_square = + (((m / scale) * (m / scale)) + ((c / scale) * (c / scale)) + + ((p / scale) * (p / scale))) + / 3.0; + return (rel_eps * scale) * (rel_eps * scale) * mean_square; + } +}; + /*! * WENO3-JS (Jiang-Shu) reconstruction to cell faces * diff --git a/manual/sphinx/user_docs/differential_operators.rst b/manual/sphinx/user_docs/differential_operators.rst index de1d370f81..75f4c8e6e9 100644 --- a/manual/sphinx/user_docs/differential_operators.rst +++ b/manual/sphinx/user_docs/differential_operators.rst @@ -800,6 +800,17 @@ values. Several slope limiters are defined in ``fv_ops.hxx``: avoids piecewise branches at extrema. This can be useful for nonlinear solvers and finite-difference Jacobian calculations. +* ``VanAlbadaSmooth`` - An improved version of the VanAlbada limiter which + smoothly disables itself (i.e. smoothly tends to second order) when cell + differences are very small to prevent non-differentiable chatter due to + Newton noise. Based on Eq. 37 of van Albada, van Leer and + Roberts (1982), A&A 108, 76, but with an addition of a smooth clipping + function which smoothly disables the limiter also at small peaks and dips. + Both the reconstruction and clip tolerance is controlled by the constructor + argument ``rel_eps`` (default ``1e-3``) times the RMS of the three cell + values. The tolerance is not grid dependent and so the limiter + will tend towards second-order at extremely fine grids. + * ``WENO3`` - A third-order smooth WENO (Jiang-Shu) cell-face reconstruction using a three-point stencil. This is typically less dissipative than TVD slope limiters, but is not strictly monotonicity-preserving. diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index ff854b6cfa..d3fcebcb0a 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -671,6 +671,24 @@ template Field3D Div_par_fvv_heating(const Field3D& f_in, const Field template Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f); +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, + bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux, + bool dissipative); +template Field3D Div_par_fvv_heating(const Field3D& f_in, + const Field3D& v_in, + const Field3D& wave_speed_in, + Field3D& flow_ylow, bool fixflux); +template Field3D Div_a_Grad_perp_limit(const Field3D& a, + const Field3D& g, + const Field3D& f); + template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux); template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, From e4439a093242b46c0bc5030aa97d677bc216756f Mon Sep 17 00:00:00 2001 From: Mike Kryjak Date: Mon, 9 Nov 2026 16:30:18 +0000 Subject: [PATCH 2/2] VanAlbadaSmooth MMS tests Add MMS tests for original and new Van Albada limiters. Include tests on a new velocity function which has a stagnation point to emulate a midplane stagnation point in a tokamak grid. do tests on more grids to show that the new limiter transitions to 2nd order eventually, as finer grids reduce differences between the cells. --- tests/MMS/spatial/finite-volume/data/BOUT.inp | 4 ++ tests/MMS/spatial/finite-volume/fv_mms.cxx | 33 ++++++++- tests/MMS/spatial/finite-volume/mms.py | 11 ++- tests/MMS/spatial/finite-volume/runtest | 26 ++++++- tests/unit/include/bout/test_fv_ops.cxx | 71 +++++++++++++++++++ 5 files changed, 141 insertions(+), 4 deletions(-) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp index 029011e437..11d8df463a 100644 --- a/tests/MMS/spatial/finite-volume/data/BOUT.inp +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -5,6 +5,10 @@ FV_Div_par_solution = -0.188495559215388*sin(3.0*y) FV_Div_par_K_Grad_par_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) FV_Div_par_K_Grad_par_mod_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) v = (0.1*cos(3.0*y) + 1)/(0.1*sin(2.0*y) + 1) +wave_speed = 2 +v_stagnation = 0.5*sin(y) +FV_Div_par_mod_stagnation_solution = 0.314159265358979*(0.1*sin(2.0*y) + 1)*cos(y) + 0.0628318530717959*sin(y)*cos(2.0*y) +FV_Div_par_fvv_stagnation_solution = 0.314159265358979*(0.1*sin(2.0*y) + 1)*sin(y)*cos(y) + 0.0314159265358979*sin(y)^2*cos(2.0*y) [mesh] MXG = 0 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 1de6052e89..3f7711edb9 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -39,9 +39,14 @@ int main(int argc, char** argv) { Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; + // Velocity profile with stagnation points at 0 and pi + Field3D v_stagnation{ + FieldFactory::get()->create3D("v_stagnation", Options::getRoot(), mesh)}; + + Field3D wave{FieldFactory::get()->create3D("wave_speed", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform. - mesh->communicate(input, v); + mesh->communicate(input, v, v_stagnation, wave); Options dump; // Add mesh geometry variables @@ -81,6 +86,32 @@ int main(int argc, char** argv) { fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_Superbee"); + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_VanAlbada"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_VanAlbada"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_VanAlbada"); + fv_op_test("FV_Div_par_mod_stagnation", dump, input, + FV::Div_par_mod(input, v_stagnation, wave, flow_ylow), + "_VanAlbada"); + fv_op_test("FV_Div_par_fvv_stagnation", dump, input, + FV::Div_par_fvv(input, v_stagnation, wave), "_VanAlbada"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_VanAlbadaSmooth"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), + "_VanAlbadaSmooth"); + fv_op_test("FV_Div_par_fvv", dump, input, + FV::Div_par_fvv(input, v, v), "_VanAlbadaSmooth"); + fv_op_test("FV_Div_par_mod_stagnation", dump, input, + FV::Div_par_mod(input, v_stagnation, wave, flow_ylow), + "_VanAlbadaSmooth"); + fv_op_test("FV_Div_par_fvv_stagnation", dump, input, + FV::Div_par_fvv(input, v_stagnation, wave), + "_VanAlbadaSmooth"); + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, Div_par_K_Grad_par_mod(v, input, flow_ylow)); diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py index a95ecc1328..d4bcb4116f 100755 --- a/tests/MMS/spatial/finite-volume/mms.py +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -9,7 +9,7 @@ from boututils.boutwarnings import AlwaysWarning from boutdata.data import BoutOptionsFile from boutdata.mms import exprToStr, y, Grad_par, Div_par, Metric -from sympy import sin, cos, Expr +from sympy import sin, cos, Expr, Integer warnings.simplefilter("ignore", AlwaysWarning) @@ -27,9 +27,13 @@ # Turn solution into real x and z coordinates replace = [(y, metric.y * 2 * pi / Ly)] +# Velocity profile with stagnation points at y = 0 and y = pi +v_stagnation = 0.5 * sin(y) + f = f.subs(replace) fv = fv.subs(replace) v = fv / f +v_stagnation = v_stagnation.subs(replace) # Substitute back to get input y coordinates replace = [(metric.y, y * Ly / (2 * pi))] @@ -56,6 +60,11 @@ def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, v)), ("FV_Div_par_mod_solution", Div_par(f * v)), ("FV_Div_par_fvv_solution", Div_par(f * v * v)), + ("v_stagnation", v_stagnation), + # Constant wave speed for the stagnation tests, larger than |v| and |v_stagnation| + ("wave_speed", Integer(2)), + ("FV_Div_par_mod_stagnation_solution", Div_par(f * v_stagnation)), + ("FV_Div_par_fvv_stagnation_solution", Div_par(f * v_stagnation * v_stagnation)), ): expr_str = exprToStr(expr.subs(replace)) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index bcd4672545..15254bdea0 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -34,11 +34,24 @@ OPERATORS = { "FV_Div_par_Superbee": 1.5, "FV_Div_par_mod_Superbee": 1.5, "FV_Div_par_fvv_Superbee": 1.5, + # Van Albada is limited at extrema, like the limiters above + "FV_Div_par_VanAlbada": 1.5, + "FV_Div_par_mod_VanAlbada": 1.5, + "FV_Div_par_fvv_VanAlbada": 1.5, + "FV_Div_par_mod_stagnation_VanAlbada": 1.5, + "FV_Div_par_fvv_stagnation_VanAlbada": 1.5, + # VanAlbadaSmooth is smooth at extrema, so reaches 2nd order + # once the grid resolves them (above about 128 cells) + "FV_Div_par_VanAlbadaSmooth": 2, + "FV_Div_par_mod_VanAlbadaSmooth": 2, + "FV_Div_par_fvv_VanAlbadaSmooth": 2, + "FV_Div_par_mod_stagnation_VanAlbadaSmooth": 2, + "FV_Div_par_fvv_stagnation_VanAlbadaSmooth": 2, "FV_Div_par_K_Grad_par": 2, "FV_Div_par_K_Grad_par_mod": 2, } -# Resolution in y and z -NLIST = [8, 16, 32, 64] +# Resolution in y and z. The order is measured between the last two. +NLIST = [8, 16, 32, 64, 128, 256, 512, 1024] dx = 1.0 / array(NLIST) @@ -179,8 +192,17 @@ if __name__ == "__main__": help="Output file to dump errors as JSON", default="fv_operator_errors.json", ) + parser.add_argument( + "--nlist", + type=int, + nargs="+", + help="Resolutions to run (default: %(default)s)", + default=NLIST, + ) args = parser.parse_args() + NLIST = args.nlist + dx = 1.0 / array(NLIST) error2, failures = test_fv_operators() success = len(failures) == 0 diff --git a/tests/unit/include/bout/test_fv_ops.cxx b/tests/unit/include/bout/test_fv_ops.cxx index 5064b076b9..4b21d3bfb6 100644 --- a/tests/unit/include/bout/test_fv_ops.cxx +++ b/tests/unit/include/bout/test_fv_ops.cxx @@ -54,6 +54,77 @@ TEST(FVOpsLimiterTest, VanAlbadaLimitsToSmallerGradient) { EXPECT_NEAR(s.L, 0.7, 1e-12); EXPECT_NEAR(s.R, 1.3, 1e-12); } +// Flat profile gives zero slope +TEST(FVOpsLimiterTest, VanAlbadaSmoothConstant) { + FV::Stencil1D s{}; + s.m = 1.0; + s.c = 1.0; + s.p = 1.0; + + FV::VanAlbadaSmooth limiter{1e-3}; + limiter(s); + + EXPECT_DOUBLE_EQ(s.L, 1.0); + EXPECT_DOUBLE_EQ(s.R, 1.0); +} + +// Linear profile gives no limitation +TEST(FVOpsLimiterTest, VanAlbadaSmoothLinear) { + FV::Stencil1D s{}; + s.m = 0.0; + s.c = 1.0; + s.p = 2.0; + + FV::VanAlbadaSmooth limiter{1e-3}; + limiter(s); + + EXPECT_NEAR(s.L, 0.5, 1e-12); + EXPECT_NEAR(s.R, 1.5, 1e-12); +} + +// A sharp peak gets clipped to zero slope +TEST(FVOpsLimiterTest, VanAlbadaSmoothSharpExtremumGivesZeroSlope) { + FV::Stencil1D s{}; + s.m = 0.0; + s.c = 1.0; + s.p = 0.5; + + FV::VanAlbadaSmooth limiter{1e-3}; + limiter(s); + + EXPECT_NEAR(s.L, 1.0, 1e-12); + EXPECT_NEAR(s.R, 1.0, 1e-12); +} + +// A very small peak below the tolerance tends towards central +// slope (preserving 2nd order but causing small overshoot). +// Note: this makes it not TVD when differences are very small. +TEST(FVOpsLimiterTest, VanAlbadaSmoothSmallExtremumTendsToCentralSlope) { + FV::Stencil1D s{}; + s.m = 1.0; + s.c = 1.0 + 1e-5; + s.p = 1.0 + 0.5e-5; + + FV::VanAlbadaSmooth limiter{1e-3}; + limiter(s); + + EXPECT_NEAR(s.R - s.L, 2.5e-6, 1e-9); + EXPECT_GT(s.R, s.c); +} + +// Very large values: a sharp peak still gets zero slope, with no overflow +TEST(FVOpsLimiterTest, VanAlbadaSmoothLargeValues) { + FV::Stencil1D s{}; + s.m = 1e100; + s.c = 2e100; + s.p = 1e100; + + FV::VanAlbadaSmooth limiter{1e-3}; + limiter(s); + + EXPECT_DOUBLE_EQ(s.L, 2e100); + EXPECT_DOUBLE_EQ(s.R, 2e100); +} TEST(FVOpsLimiterTest, WENO3Constant) { FV::Stencil1D s{};