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, 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{};