Add VanAlbadaSmooth, an improved Van Albada slope limiter

include/bout/fv_ops.hxx

@@ -90,6 +90,7 @@ struct MinMod;
 struct MC;
 struct Superbee;
 struct VanAlbada;
+struct VanAlbadaSmooth;
 struct WENO3;
 
 /*!

include/bout/fv_ops_impl.hxx

@@ -21,6 +21,7 @@
 #include "bout/vector2d.hxx"
 #include <bout/yboundary_regions.hxx>
 
+#include <algorithm>
 #include <cmath>
 
 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
    *

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.

src/mesh/fv_ops.cxx

@@ -671,6 +671,24 @@ template Field3D Div_par_fvv_heating<VanAlbada>(const Field3D& f_in, const Field
 template Field3D Div_a_Grad_perp_limit<VanAlbada>(const Field3D& a, const Field3D& g,
                                                   const Field3D& f);
 
+template Field3D Div_par_fvv<VanAlbadaSmooth>(const Field3D& f_in, const Field3D& v_in,
+                                              const Field3D& wave_speed_in, bool fixflux);
+template Field3D Div_par<VanAlbadaSmooth>(const Field3D& f_in, const Field3D& v_in,
+                                          const Field3D& wave_speed_in, bool fixflux);
+template Field3D Div_f_v<VanAlbadaSmooth>(const Field3D& n_in, const Vector3D& v,
+                                          bool bndry_flux);
+template Field3D Div_par_mod<VanAlbadaSmooth>(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<VanAlbadaSmooth>(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<VanAlbadaSmooth>(const Field3D& a,
+                                                        const Field3D& g,
+                                                        const Field3D& f);
+
 template Field3D Div_par_fvv<WENO3>(const Field3D& f_in, const Field3D& v_in,
                                     const Field3D& wave_speed_in, bool fixflux);
 template Field3D Div_par<WENO3>(const Field3D& f_in, const Field3D& v_in,

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

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<FV::Superbee>(input, v, v),
              "_Superbee");
 
+  fv_op_test("FV_Div_par", dump, input, FV::Div_par<FV::VanAlbada>(input, v, v),
+             "_VanAlbada");
+  fv_op_test("FV_Div_par_mod", dump, input,
+             FV::Div_par_mod<FV::VanAlbada>(input, v, v, flow_ylow), "_VanAlbada");
+  fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv<FV::VanAlbada>(input, v, v),
+             "_VanAlbada");
+  fv_op_test("FV_Div_par_mod_stagnation", dump, input,
+             FV::Div_par_mod<FV::VanAlbada>(input, v_stagnation, wave, flow_ylow),
+             "_VanAlbada");
+  fv_op_test("FV_Div_par_fvv_stagnation", dump, input,
+             FV::Div_par_fvv<FV::VanAlbada>(input, v_stagnation, wave), "_VanAlbada");
+
+  fv_op_test("FV_Div_par", dump, input, FV::Div_par<FV::VanAlbadaSmooth>(input, v, v),
+             "_VanAlbadaSmooth");
+  fv_op_test("FV_Div_par_mod", dump, input,
+             FV::Div_par_mod<FV::VanAlbadaSmooth>(input, v, v, flow_ylow),
+             "_VanAlbadaSmooth");
+  fv_op_test("FV_Div_par_fvv", dump, input,
+             FV::Div_par_fvv<FV::VanAlbadaSmooth>(input, v, v), "_VanAlbadaSmooth");
+  fv_op_test("FV_Div_par_mod_stagnation", dump, input,
+             FV::Div_par_mod<FV::VanAlbadaSmooth>(input, v_stagnation, wave, flow_ylow),
+             "_VanAlbadaSmooth");
+  fv_op_test("FV_Div_par_fvv_stagnation", dump, input,
+             FV::Div_par_fvv<FV::VanAlbadaSmooth>(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));

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}")

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

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