From 379cdde60f0a227af5dc28371edbf0271a1debfb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 14 May 2026 20:55:43 +1000 Subject: [PATCH 1/2] Expose theta on SemiLagrangian DDt for BE/CN choice on AM flux MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Adams-Moulton flux integrator in SemiLagrangian_DDt was hardcoded to theta=0.5 (Crank-Nicolson, 2nd-order accurate, A-stable). For stiff parabolic terms — e.g. advection-diffusion with a thin boundary layer on deformed elements — CN is not L-stable: the amplification factor (1 - λΔt/2)/(1 + λΔt/2) → -1 for stiff modes, so the integrator sign-flip-reflects them rather than damping. Combined with the Picard iteration's coupling, this can produce step-on-step oscillating T overshoots. Adds ``theta`` as an __init__ parameter (default 0.5, preserving the legacy SLCN behaviour) and as an instance attribute settable after construction: adv_diff.DuDt.theta = 1.0 # Backward Euler (L-stable, 1st order) adv_diff.DFDt.theta = 1.0 The two hardcoded ``_update_am_values(..., 0.5)`` sites in SemiLagrangian (initial coefficient setup and per-step refresh in update_pre_solve) now use ``self.theta``. The matching constructs in the Symbolic, Eulerian, Lagrangian and Lagrangian_Swarm classes already used self.theta — this PR brings SemiLagrangian into line. No behaviour change for existing users (default unchanged). Useful as a per-solver knob when CN-driven ringing is observed on stiff problems. Underworld development team with AI support from Claude Code --- src/underworld3/systems/ddt.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index bda85a8a..1d21cc31 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1172,6 +1172,7 @@ def __init__( smoothing=0.0, preserve_moments=False, with_forcing_history: bool = False, + theta: float = 0.5, ): super().__init__() @@ -1185,6 +1186,18 @@ def __init__( self.order = order self.preserve_moments = preserve_moments self.with_forcing_history = with_forcing_history + # 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 @@ -1295,7 +1308,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 @@ -1671,7 +1684,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 From 290951bfc3c4aefc25494496edafb91b9e0ac637 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 14 May 2026 21:35:57 +1000 Subject: [PATCH 2/2] Address Copilot review on PR #187: docstring + regression test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. Add the new ``theta`` constructor parameter to the SemiLagrangian class docstring's Parameters section. Was missing while the other DDt classes (Symbolic, Eulerian, Lagrangian, Lagrangian_Swarm) already documented theirs. 2. Add a regression test (tests/test_1053_ddt_theta.py) covering: - default theta=0.5 yields order-1 AM coefficients [0.5, 0.5] (CN) - theta=1.0 yields [1.0, 0.0] (Backward Euler) - theta=0.0 yields [0.0, 1.0] (Forward Euler — completeness) - mutation after construction (d.theta = 1.0) is picked up on the next AM coefficient refresh - omitting theta preserves legacy CN behaviour Underworld development team with AI support from Claude Code --- src/underworld3/systems/ddt.py | 14 +++++ tests/test_1053_ddt_theta.py | 93 ++++++++++++++++++++++++++++++++++ 2 files changed, 107 insertions(+) create mode 100644 tests/test_1053_ddt_theta.py diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 1d21cc31..dedccf04 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 ----- 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)