diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 2777f4f3..f703d1f0 100644 --- a/src/underworld3/systems/ddt.py +++ b/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 diff --git a/tests/test_1053_ddt_theta.py b/tests/test_1053_ddt_theta.py new file mode 100644 index 00000000..7c5deb70 --- /dev/null +++ b/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)