Expose theta on SemiLagrangian DDt for BE/CN choice

src/underworld3/systems/ddt.py

@@ -1139,6 +1139,20 @@ class SemiLagrangian(uw_object):
         :meth:`update_forcing_history` (direct nodal evaluation of
         ``forcing_fn`` — typically the constitutive model's strain-rate
         symbol).
+    theta : float, default=0.5
+        Adams-Moulton θ for the implicit flux integrator at order 1.
+        The order-1 AM coefficients are ``[θ, 1-θ]``:
+
+        - ``θ = 0.5`` → Crank-Nicolson (trapezoidal, second-order
+          accurate, A-stable). Default, matches legacy SLCN behaviour.
+        - ``θ = 1.0`` → Backward Euler (L-stable, monotone for
+          diffusion, first-order accurate). Use for stiff parabolic
+          terms (under-resolved sharp gradients on deformed cells)
+          where CN's lack of L-stability causes sign-flip ringing
+          on stiff modes.
+
+        Settable after construction as a property:
+        ``adv_diff.DuDt.theta = 1.0``.
 
     Notes
     -----
@@ -1173,6 +1187,7 @@ def __init__(
         preserve_moments=False,
         with_forcing_history: bool = False,
         monotone_mode: Optional[str] = None,
+        theta: float = 0.5,
     ):
         super().__init__()
 
@@ -1200,6 +1215,18 @@ def __init__(
         # Settable after construction:
         #   ``adv_diff.DuDt.monotone_mode = "clamp"``
         self.monotone_mode = monotone_mode
+        # Adams-Moulton θ for the implicit flux at order 1.
+        # The order-1 AM coefficients are ``[θ, 1-θ]``:
+        #   θ=0.5  → Crank-Nicolson (A-stable, 2nd order accuracy on
+        #            flux; NOT L-stable — stiff modes get amplification
+        #            factor → -1, can ring on under-resolved sharp
+        #            gradients in deformed cells)
+        #   θ=1.0  → Backward Euler (L-stable, monotone for diffusion;
+        #            1st order accuracy on flux)
+        # Default 0.5 preserves the legacy SLCN behaviour.
+        # Settable after construction:
+        #   ``adv_diff.DuDt.theta = 1.0``
+        self.theta = float(theta)
 
         # Forcing-history storage. Allocated only if requested. Populated
         # each step via update_forcing_history(forcing_fn) — used by ETD-2
@@ -1310,7 +1337,7 @@ def __init__(
         self._exp_coeffs = _create_exp_coefficients(self.instance_number)
         # Initialise to order-1 / viscous values
         _update_bdf_values(self._bdf_coeffs, 1, None, [])
-        _update_am_values(self._am_coeffs, 1, 0.5)
+        _update_am_values(self._am_coeffs, 1, self.theta)
         _update_exp_values(self._exp_coeffs, None, None)
 
         # Working variable that has a potentially different discretisation from psi_star
@@ -1696,7 +1723,7 @@ def update_pre_solve(
 
         # Update coefficient values for current effective_order and dt
         _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history)
-        _update_am_values(self._am_coeffs, self.effective_order, 0.5)
+        _update_am_values(self._am_coeffs, self.effective_order, self.theta)
 
         ## Progress from the oldest part of the history
         # 1. Copy the stored values down the chain in preparation for the next timestep

tests/test_1053_ddt_theta.py

@@ -0,0 +1,93 @@
+"""Regression tests for SemiLagrangian DDt ``theta`` parameter.
+
+``theta`` controls the Adams-Moulton flux integrator's coefficient at
+order 1: AM coefficients are ``[θ, 1-θ]``. Default 0.5 is Crank-
+Nicolson (legacy SLCN); 1.0 is Backward Euler.
+
+These tests verify the constructor parameter is plumbed through to the
+AM coefficient values and that mutation after construction is picked
+up on the next update.
+"""
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+from underworld3.systems import ddt as ddt_module
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+def _make_sl(theta=0.5, order=1):
+    """Build a tiny SemiLagrangian DDt for a scalar field."""
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(4, 4),
+        minCoords=(0.0, 0.0),
+        maxCoords=(1.0, 1.0),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=1)
+    psi = sympy.zeros(1, 1)
+    return ddt_module.SemiLagrangian(
+        mesh,
+        psi_fn=psi,
+        V_fn=v.sym,
+        vtype=uw.VarType.SCALAR,
+        degree=2,
+        continuous=False,
+        order=order,
+        theta=theta,
+    )
+
+
+class TestTheta:
+
+    def test_default_theta_is_cn(self):
+        d = _make_sl()
+        assert d.theta == 0.5
+        # __init__ ran _update_am_values(self._am_coeffs, 1, theta=0.5)
+        assert float(d._am_coeffs[0].sym) == pytest.approx(0.5)
+        assert float(d._am_coeffs[1].sym) == pytest.approx(0.5)
+
+    def test_theta_one_is_backward_euler(self):
+        d = _make_sl(theta=1.0)
+        assert d.theta == 1.0
+        # Order-1 AM coefficients with θ=1.0: [1.0, 0.0]
+        assert float(d._am_coeffs[0].sym) == pytest.approx(1.0)
+        assert float(d._am_coeffs[1].sym) == pytest.approx(0.0)
+
+    def test_theta_zero_is_forward_euler(self):
+        d = _make_sl(theta=0.0)
+        assert d.theta == 0.0
+        # AM coefficients with θ=0.0: [0.0, 1.0]
+        assert float(d._am_coeffs[0].sym) == pytest.approx(0.0)
+        assert float(d._am_coeffs[1].sym) == pytest.approx(1.0)
+
+    def test_theta_attribute_settable_after_construction(self):
+        d = _make_sl(theta=0.5)
+        # Mutate the instance attribute
+        d.theta = 1.0
+        # Re-run the per-step AM-coefficient update path manually
+        # (replicates what update_pre_solve does on each step).
+        ddt_module._update_am_values(
+            d._am_coeffs, d.effective_order, d.theta)
+        assert float(d._am_coeffs[0].sym) == pytest.approx(1.0)
+        assert float(d._am_coeffs[1].sym) == pytest.approx(0.0)
+
+    def test_theta_default_unchanged_when_omitted(self):
+        """Backwards compatibility: omitting theta keeps legacy CN."""
+        mesh = uw.meshing.StructuredQuadBox(
+            elementRes=(4, 4),
+            minCoords=(0.0, 0.0),
+            maxCoords=(1.0, 1.0),
+        )
+        v = uw.discretisation.MeshVariable(
+            "U_legacy", mesh, mesh.dim, degree=1)
+        psi = sympy.zeros(1, 1)
+        d = ddt_module.SemiLagrangian(
+            mesh, psi_fn=psi, V_fn=v.sym,
+            vtype=uw.VarType.SCALAR, degree=2, continuous=False,
+            order=1)
+        assert d.theta == 0.5
+        assert float(d._am_coeffs[0].sym) == pytest.approx(0.5)
+        assert float(d._am_coeffs[1].sym) == pytest.approx(0.5)