fv_ops: WENO3 cell edge reconstruction

Author: bendudson

examples/finite-volume/fluid/README.md

@@ -54,3 +54,7 @@ for the second order MinMod method.
 Or, to use a smooth symmetric limiter:
 
     FV::Div_par<FV::VanAlbada>
+
+Or a higher-order smooth WENO reconstruction:
+
+    FV::Div_par<FV::WENO3>

include/bout/fv_ops.hxx

@@ -208,6 +208,61 @@ struct VanAlbada {
   }
 };
 
+/*!
+   * WENO3-JS (Jiang-Shu) reconstruction to cell faces
+   *
+   * This is a third-order essentially non-oscillatory reconstruction using two
+   * candidate second-order polynomials and smoothness-weighted blending.
+   *
+   * Unlike TVD slope limiters (e.g. ``MC``), WENO reconstruction is generally
+   * smooth (differentiable) for all inputs, but it does not enforce strict
+   * monotonicity.
+   *
+   * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in
+   * replacement anywhere `Stencil1D` is populated with those values.
+   */
+struct WENO3 {
+  void operator()(Stencil1D& n) {
+    // Right face (between c and p): value from cell c (left state at i+1/2)
+    const BoutReal p0_r = 0.5 * (-n.m + 3.0 * n.c);
+    const BoutReal p1_r = 0.5 * (n.c + n.p);
+
+    const BoutReal beta0_r = SQ(n.c - n.m);
+    const BoutReal beta1_r = SQ(n.p - n.c);
+
+    // Left face (between m and c): value from cell c (right state at i-1/2)
+    const BoutReal p0_l = 0.5 * (-n.p + 3.0 * n.c);
+    const BoutReal p1_l = 0.5 * (n.m + n.c);
+
+    const BoutReal beta0_l = beta1_r;
+    const BoutReal beta1_l = beta0_r;
+
+    // Smoothness parameter (scaled to local variation)
+    const BoutReal eps = 1e-12 * (beta0_r + beta1_r) + 1e-30;
+
+    // Linear weights for WENO3-JS
+    constexpr BoutReal d0 = 1.0 / 3.0;
+    constexpr BoutReal d1 = 2.0 / 3.0;
+
+    // Right face weights
+    const BoutReal a0_r = d0 / SQ(eps + beta0_r);
+    const BoutReal a1_r = d1 / SQ(eps + beta1_r);
+    const BoutReal wsum_r = a0_r + a1_r;
+    const BoutReal w0_r = a0_r / wsum_r;
+    const BoutReal w1_r = a1_r / wsum_r;
+
+    // Left face weights (mirrored)
+    const BoutReal a0_l = d0 / SQ(eps + beta0_l);
+    const BoutReal a1_l = d1 / SQ(eps + beta1_l);
+    const BoutReal wsum_l = a0_l + a1_l;
+    const BoutReal w0_l = a0_l / wsum_l;
+    const BoutReal w1_l = a1_l / wsum_l;
+
+    n.R = w0_r * p0_r + w1_r * p1_r;
+    n.L = w0_l * p0_l + w1_l * p1_l;
+  }
+};
+
 /*!
    * Communicate fluxes between processors
    * Takes values in guard cells, and adds them to cells

manual/sphinx/user_docs/differential_operators.rst

@@ -671,6 +671,10 @@ 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.
 
+* ``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.
+
 Useful resources on slope limiters include:
 
 * `Wikipedia's Flux Limiter page <https://en.wikipedia.org/wiki/Flux_limiter>`_

tests/unit/include/bout/test_fv_ops.cxx

@@ -54,3 +54,29 @@ TEST(FVOpsLimiterTest, VanAlbadaLimitsToSmallerGradient) {
   EXPECT_NEAR(s.L, 0.7, 1e-12);
   EXPECT_NEAR(s.R, 1.3, 1e-12);
 }
+
+TEST(FVOpsLimiterTest, WENO3Constant) {
+  FV::Stencil1D s{};
+  s.m = 1.0;
+  s.c = 1.0;
+  s.p = 1.0;
+
+  FV::WENO3 recon;
+  recon(s);
+
+  EXPECT_DOUBLE_EQ(s.L, 1.0);
+  EXPECT_DOUBLE_EQ(s.R, 1.0);
+}
+
+TEST(FVOpsLimiterTest, WENO3Linear) {
+  FV::Stencil1D s{};
+  s.m = 0.0;
+  s.c = 1.0;
+  s.p = 2.0;
+
+  FV::WENO3 recon;
+  recon(s);
+
+  EXPECT_NEAR(s.L, 0.5, 1e-14);
+  EXPECT_NEAR(s.R, 1.5, 1e-14);
+}