Expose theta on SemiLagrangian DDt for BE/CN choice#187
Merged
Conversation
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
Contributor
There was a problem hiding this comment.
Pull request overview
This PR exposes the theta parameter on SemiLagrangian DDt so users can choose Crank-Nicolson or Backward Euler behavior for order-1 Adams-Moulton flux integration while preserving the default legacy behavior.
Changes:
- Adds
thetato theSemiLagrangian.__init__signature with default0.5. - Stores
self.thetaand uses it when initializing/updating AM coefficients. - Replaces hardcoded
0.5AM coefficient updates withself.theta.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| smoothing=0.0, | ||
| preserve_moments=False, | ||
| with_forcing_history: bool = False, | ||
| theta: float = 0.5, |
| # 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) |
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
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Exposes the Adams-Moulton
thetaparameter onSemiLagrangian_DDt. Previously hardcoded to0.5(Crank-Nicolson) at two sites in the class. Default is unchanged so existing users see no behaviour difference; settingadv_diff.DuDt.theta = 1.0now gives Backward Euler on the implicit flux integrator.Why
The AM order-1 coefficients are
[θ, 1-θ]:θ = 0.5→ Crank-Nicolson. Trapezoidal. A-stable, second-order accuracy on the flux term. NOT L-stable: for stiff modes the amplification factor(1 − λΔt/2) / (1 + λΔt/2)approaches −1, so the integrator reflects-with-sign-flip rather than damps. On under-resolved sharp gradients in deformed cells this can ignite step-on-step oscillating overshoots in T.θ = 1.0→ Backward Euler. L-stable, monotone for diffusion. First-order accurate on the flux term. The standard fallback for stiff parabolic problems.The other
SemiLagrangian-family classes (Symbolic,Eulerian,Lagrangian,Lagrangian_Swarm) already takethetaas an__init__parameter and useself.thetain their AM update sites. This PR bringsSemiLagrangian(the one used byAdvDiffusionSLCN) into line.Usage
Or via the constructor:
Test plan
Files changed
src/underworld3/systems/ddt.py(+15 lines, -2 lines)Underworld development team with AI support from Claude Code