From f7d1a216492756281ada84373f9ca2eaefb2dde0 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 01/14] Expose the BDF and Adams-Moulton coefficient symbols on the DDt managers A solver that assembles its own weighted sum of history terms (an Eulerian scheme applying a multistep rule to a spatial operator) needs the constants-routed coefficient expressions, not just their current values. Read-only accessors; no behaviour change. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- src/underworld3/systems/ddt.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e198aad57..a37533cbe 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -666,6 +666,29 @@ def bdf_coefficients(self): """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + @property + def bdf_coefficient_expressions(self): + r"""The BDF coefficient symbols :math:`[c_0, c_1, \dots]` as UWexpressions. + + For a solver that assembles its own weighted sum of history terms + (an Eulerian scheme applying the multistep rule to a spatial + operator, say). The symbols are routed through PETSc's + ``constants[]`` array, so their values follow ``effective_order`` + and the timestep without a recompile; ``bdf_coefficients`` gives + the current values. + """ + return list(self._bdf_coeffs) + + @property + def am_coefficient_expressions(self): + r"""The Adams-Moulton coefficient symbols :math:`[a_0, a_1, \dots]` as UWexpressions. + + :math:`a_0` weights the new state, :math:`a_k` the history slot + ``psi_star[k-1]``. Same constants-routing as + :attr:`bdf_coefficient_expressions`. + """ + return list(self._am_coeffs) + def _history_syms(self): """History terms as sympy expressions for the weighted sums. From e84dea98e0c6d4aaeac3810526ea075c92923ee7 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 02/14] Pack and index auxiliary fields by DM field, not by position in mesh.vars A MeshVariable that is dropped and garbage-collected (the default Model holds the only strong reference; uw.reset_default_model() releases it, and the statistics helpers delete temporaries deliberately) leaves its PETSc field in the DM. Mesh.update_lvec zipped mesh.vars.values() against the field decomposition by position, and the JIT's petsc_a[] offsets were a running count over the live variables, so every later variable was packed into, and read from, the wrong slots. Measured: a P0 cell-size field landing in a P2 slot as garbage, NaN residuals in one run and a subtly wrong answer in the next, depending on when the collector ran. update_lvec now packs by field name and zeroes an orphaned field; the JIT reads component offsets from the DM's own field list and patches each variable from its field_id. Regression test: 2 of its 3 checks fail without the fix. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../discretisation/discretisation_mesh.py | 20 +++- src/underworld3/utilities/_jitextension.py | 34 +++++- ...st_1058_dropped_meshvariable_aux_layout.py | 113 ++++++++++++++++++ 3 files changed, 159 insertions(+), 8 deletions(-) create mode 100644 tests/test_1058_dropped_meshvariable_aux_layout.py diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d3c67611b..c0e65cf82 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3896,13 +3896,21 @@ def update_lvec(self, swarm_sync=True): # The field decomposition seems to fail if coarse DMs are present names, isets, dms = self.dm.createFieldDecomposition() - # traverse subdms, taking user generated data in the subdm - # local vec, pushing it into a global sub vec - for var, subiset, subdm in zip(self.vars.values(), isets, dms): - # var.vec lazily creates the PETSc local vector on first access - lvec = var.vec + # Traverse the DM's fields BY NAME. `self.vars` holds its + # variables weakly, so a dropped-and-collected variable leaves + # a field behind in the DM; a positional zip would then pack + # every later variable into the wrong field (measured: the + # cell-size field landing in a P2 slot as garbage, NaN + # residuals in a solver that reads it). An orphaned field is + # zeroed so nothing stale can reach a kernel. + for name, subiset, subdm in zip(names, isets, dms): + var = self.vars.get(name) subvec = a_global.getSubVector(subiset) - subdm.localToGlobal(lvec, subvec, addv=False) + if var is None: + subvec.set(0.0) + else: + # var.vec lazily creates the PETSc local vector on first access + subdm.localToGlobal(var.vec, subvec, addv=False) a_global.restoreSubVector(subiset, subvec) for iset in isets: diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 01e1a0582..21d190348 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -773,6 +773,24 @@ def getext( @timing.routine_timer_decorator +def _aux_component_offsets(mesh): + """Component offset of every field of the mesh DM, keyed by field id. + + Read from the DM itself, not from ``mesh.vars``: a MeshVariable that + was dropped and collected leaves its PETSc field in the DM (a DMPlex + cannot shed a field), and PETSc lays the auxiliary arrays out over + ALL fields in field order. The offsets therefore have to count the + orphaned fields too. + """ + offsets = {} + total = 0 + for field_id in range(mesh.dm.getNumFields()): + fe, _label = mesh.dm.getField(field_id) + offsets[field_id] = total + total += fe.getNumComponents() + return offsets + + def generate_c_source( name, mesh: underworld3.discretisation.Mesh, @@ -822,7 +840,7 @@ def generate_c_source( count_bd_residual_sig, count_bd_jacobian_sig = callbacks.counts # `_ccode` patching - def ccode_patch_fns(varlist, prefix_str): + def ccode_patch_fns(varlist, prefix_str, component_offsets=None): """ This function patches uw functions with the necessary ccode routines for the code printing. @@ -848,11 +866,22 @@ def ccode_patch_fns(varlist, prefix_str): ordered according to their `field_id`. prefix_str: str The string prefix to write. + component_offsets: dict, optional + Component offset of every field in the DM, by ``field_id`` + (see ``_aux_component_offsets``). When given, each variable + is patched from ITS OWN field's offset instead of a running + count over ``varlist``: a field whose Python variable has + been dropped stays in the DM and still occupies its slots, + so a running count would shift every later variable onto + the wrong data. """ u_i = 0 # variable increment u_x_i = 0 # variable gradient increment lambdafunc = lambda self, printer: self._ccodestr for var in varlist: + if component_offsets is not None: + u_i = component_offsets[var.field_id] + u_x_i = u_i * mesh.cdim if var.vtype == VarType.SCALAR: # monkey patch this guy into the function type(var.fn)._ccodestr = f"{prefix_str}[{u_i}]" @@ -898,7 +927,8 @@ def ccode_patch_fns(varlist, prefix_str): # is important, as the secondary call will overwrite # those patched in the first call. - ccode_patch_fns(_stable_sorted(mesh.vars.values()), "petsc_a") + ccode_patch_fns(_stable_sorted(mesh.vars.values()), "petsc_a", + component_offsets=_aux_component_offsets(mesh)) ccode_patch_fns(primary_field_list, "petsc_u") # Also patch `BaseScalar` types. Nothing fancy - patch the overall type, diff --git a/tests/test_1058_dropped_meshvariable_aux_layout.py b/tests/test_1058_dropped_meshvariable_aux_layout.py new file mode 100644 index 000000000..cdc9d3e6e --- /dev/null +++ b/tests/test_1058_dropped_meshvariable_aux_layout.py @@ -0,0 +1,113 @@ +"""A dropped MeshVariable must not corrupt the auxiliary data of later solves. + +`mesh.vars` holds variables weakly, but a DMPlex cannot shed a field: a +variable that is dropped and garbage-collected leaves its PETSc field in +the DM. Two places used to assume the registry and the DM field list line +up by position: + +- `Mesh.update_lvec` zipped `mesh.vars.values()` against the DM's field + decomposition, so every later variable was packed into the wrong field + (the orphan's slot) and its own slot stayed at whatever it held; +- the JIT's `petsc_a[]` offsets were a running count over the live + variables, skipping the orphan's components. + +Measured before the fix: a cell-size (P0) field landing in a P2 slot as +garbage, NaN residuals (`DIVERGED_FUNCTION_NANORINF`) in one run and a +subtly wrong answer in the next, depending on when the collector ran. The +default Model holds the only strong reference to a variable (the mesh +outlives the model it was created under), so `uw.reset_default_model()`, +which the test suite runs between tests, releases every variable a script +no longer names; the variable-statistics +helpers also delete temporaries from the registry on purpose. The orphan is +an ordinary state, not a misuse. + +Run: pixi run python -m pytest tests/test_1058_dropped_meshvariable_aux_layout.py -v +""" +import gc + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + + +def _poisson_with_field_coefficient(mesh, tag): + """A Poisson solve whose answer depends on an auxiliary field (the + diffusivity is a MeshVariable), so mis-packed aux data changes it.""" + x, y = mesh.X + kappa = uw.discretisation.MeshVariable(f"kappa_{tag}", mesh, 1, degree=1) + kappa.array[:, 0, 0] = uw.function.evaluate(1.0 + 4.0 * x * y, kappa.coords).reshape(-1) + u = uw.discretisation.MeshVariable(f"u_{tag}", mesh, 1, degree=2) + poisson = uw.systems.Poisson(mesh, u) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = kappa.sym[0] + poisson.f = 1.0 + for b in ("Left", "Right", "Top", "Bottom"): + poisson.add_dirichlet_bc(0.0, b) + poisson.solve() + return np.array(u.array), kappa, u + + +def test_dropped_variable_leaves_an_orphaned_field(): + """The premise: dropping a variable does not shrink the DM.""" + mesh = _mesh() + n_fields = mesh.dm.getNumFields() + # The mesh keeps the model it was created under alive; a variable + # registers with the CURRENT default model, so a reset before and + # after creating it is what releases it (the suite's per-test reset). + uw.reset_default_model() + uw.discretisation.MeshVariable("temporary", mesh, 2, degree=2) + uw.reset_default_model() + gc.collect() + assert "temporary" not in mesh.vars + assert mesh.dm.getNumFields() == n_fields + 1 + + +def test_solve_after_a_dropped_variable_matches_a_clean_mesh(): + reference, _k, _u = _poisson_with_field_coefficient(_mesh(), "ref") + + mesh = _mesh() + uw.reset_default_model() + uw.discretisation.MeshVariable("dropped_vector", mesh, 2, degree=2) + uw.discretisation.MeshVariable("dropped_scalar", mesh, 1, degree=1) + uw.reset_default_model() + gc.collect() + assert mesh.dm.getNumFields() > len(mesh.vars) + + answer, _k, _u = _poisson_with_field_coefficient(mesh, "orphan") + assert np.allclose(answer, reference, rtol=0, atol=1e-10) + + +def test_packed_aux_vector_lands_in_the_named_fields(): + mesh = _mesh() + uw.reset_default_model() + uw.discretisation.MeshVariable("dropped", mesh, 2, degree=1) + uw.reset_default_model() + gc.collect() + assert "dropped" not in mesh.vars + x, y = mesh.X + a = uw.discretisation.MeshVariable("a_live", mesh, 1, degree=1) + a.array[:, 0, 0] = uw.function.evaluate(x + 2 * y, a.coords).reshape(-1) + + mesh.update_lvec() + names, isets, _dms = mesh.dm.createFieldDecomposition() + g = mesh.dm.getGlobalVec() + mesh.dm.localToGlobal(mesh.lvec, g) + packed = {} + for name, iset in zip(names, isets): + sub = g.getSubVector(iset) + packed[name] = (sub.min()[1], sub.max()[1]) + g.restoreSubVector(iset, sub) + mesh.dm.restoreGlobalVec(g) + + assert packed["dropped"] == (0.0, 0.0) + lo, hi = packed["a_live"] + assert lo == pytest.approx(0.0) and hi == pytest.approx(3.0) From ccfc23bbb6abdd633e3263c96959ca3c880468a2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 03/14] Add the RotatingGaussian transport oracle; fix the integral-norm error for scalar variables A Gaussian carried round the origin by rigid rotation while diffusing is exact at every time (rotation commutes with the Laplacian), so a transport scheme's error can be measured directly and the round trip after one revolution is an absolute check. AnalyticSolution.error(norm='integral') added a 1x1 Matrix symbol to a scalar expression and had never been exercised on a scalar variable. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- src/underworld3/analytic/__init__.py | 4 +- src/underworld3/analytic/_base.py | 8 ++- src/underworld3/analytic/transport.py | 86 +++++++++++++++++++++++++++ 3 files changed, 96 insertions(+), 2 deletions(-) diff --git a/src/underworld3/analytic/__init__.py b/src/underworld3/analytic/__init__.py index 6ccd3eb7d..f1f4c2102 100644 --- a/src/underworld3/analytic/__init__.py +++ b/src/underworld3/analytic/__init__.py @@ -37,7 +37,7 @@ from .inclusion import EllipticalInclusion from .kramer import CylindricalStokes from .richards import GardnerSteady, GardnerTransient -from .transport import AdvectedFront, ErfcDiffusion, Poisson1D, TwoLayerDarcy +from .transport import AdvectedFront, ErfcDiffusion, Poisson1D, RotatingGaussian, TwoLayerDarcy from .velic import ( SolA, SolB, @@ -67,6 +67,7 @@ "GardnerSteady", "GardnerTransient", "Poisson1D", + "RotatingGaussian", "SolA", "SolB", "SolC", @@ -100,6 +101,7 @@ "GardnerSteady": GardnerSteady, "GardnerTransient": GardnerTransient, "Poisson1D": Poisson1D, + "RotatingGaussian": RotatingGaussian, "SolA": SolA, "SolB": SolB, "SolC": SolC, diff --git a/src/underworld3/analytic/_base.py b/src/underworld3/analytic/_base.py index fbec7ec33..9a54e4b08 100644 --- a/src/underworld3/analytic/_base.py +++ b/src/underworld3/analytic/_base.py @@ -483,7 +483,13 @@ def error(self, field, meshvar, norm="l2"): else sympy.S.Zero ) magnitude = uw.maths.L2_norm(zero, exact, self.mesh) - return float(uw.maths.L2_norm(meshvar.sym, exact, self.mesh) / magnitude) + computed = meshvar.sym + if (isinstance(computed, sympy.MatrixBase) and computed.shape == (1, 1) + and not isinstance(exact, sympy.MatrixBase)): + # A scalar variable's symbol is a 1x1 Matrix; the exact + # scalar is not. Compare like with like. + computed = computed[0] + return float(uw.maths.L2_norm(computed, exact, self.mesh) / magnitude) if norm != "l2": raise ValueError(f"norm must be 'l2' or 'integral'; got {norm!r}") diff --git a/src/underworld3/analytic/transport.py b/src/underworld3/analytic/transport.py index 153044409..dc6872503 100644 --- a/src/underworld3/analytic/transport.py +++ b/src/underworld3/analytic/transport.py @@ -231,6 +231,92 @@ def __init__(self, mesh, kappa=1.0e-3, speed=1.0, x0=0.1, x1=0.3): ) +class RotatingGaussian(_Transport): + r"""A Gaussian carried round the origin by rigid rotation while it diffuses. + + The velocity :math:`\mathbf{u} = \omega(-y, x)` is solenoidal and rigid, so + it commutes with the Laplacian: the exact field is the free-space diffusing + Gaussian with its centre following the rotation, + + .. math:: + \phi(\mathbf{x}, t) = \frac{\sigma^2}{\sigma^2 + 2\kappa t} + \exp\!\left(-\frac{|\mathbf{x} - \mathbf{c}(t)|^2} + {2(\sigma^2 + 2\kappa t)}\right), + \qquad + \mathbf{c}(t) = R\,(\cos(\omega t + \varphi_0),\ \sin(\omega t + \varphi_0)). + + The transport test with a known answer at every time: after one + revolution, :math:`t = 2\pi/\omega`, a pure-advection field must return + to its initial state, so the round-trip error is an absolute measure and + the quarter-turn errors give the growth in between. With + :math:`\kappa = 0` the solution is regular at :math:`t = 0` and a + benchmark may start there. + + The domain is whatever mesh is supplied; the solution is exact on the + plane, so the walls should sit where the field is negligible (a few + :math:`\sigma` from the orbit) and carry :math:`\phi = 0`. + + Parameters + ---------- + mesh : Mesh + A 2D mesh containing the orbit. + sigma : float + Standard deviation of the initial Gaussian. + centre_radius : float + Orbit radius :math:`R`. + omega : float + Angular velocity; the period is :math:`2\pi/\omega`. + diffusivity : float + :math:`\kappa \ge 0`; zero is pure advection. + phase : float + Initial angular position :math:`\varphi_0` of the centre. + """ + + reference = ( + "Rigid rotation of a diffusing Gaussian; classical (e.g. the rotating " + "cone/Gaussian tests of Zalesak 1979 and LeVeque 1996, here in closed form)." + ) + eqn_solution = ( + r"\frac{\sigma^2}{\sigma^2 + 2\kappa t}" + r"\exp\left(-\frac{|\mathbf{x}-\mathbf{c}(t)|^2}{2(\sigma^2+2\kappa t)}\right)" + ) + singular_at_origin = False + + def __init__(self, mesh, sigma=0.12, centre_radius=0.5, omega=1.0, + diffusivity=0.0, phase=0.0): + super().__init__(mesh) + + if float(sigma) <= 0.0: + raise ValueError("sigma must be positive.") + if float(diffusivity) < 0.0: + raise ValueError("diffusivity must not be negative.") + + self.sigma = float(sigma) + self.centre_radius = float(centre_radius) + self.omega = float(omega) + self.diffusivity = float(diffusivity) + self.kappa = float(diffusivity) + self.phase = float(phase) + self.t = sympy.Symbol("t", positive=True) + + x, y = mesh.X + angle = self.omega * self.t + self.phase + cx = self.centre_radius * sympy.cos(angle) + cy = self.centre_radius * sympy.sin(angle) + variance = self.sigma ** 2 + 2 * self.kappa * self.t + profile = (self.sigma ** 2 / variance) * sympy.exp( + -((x - cx) ** 2 + (y - cy) ** 2) / (2 * variance)) + + self.set_scalar_field( + profile, coefficient=self.kappa, source=0, + advection=(-self.omega * y, self.omega * x)) + + @property + def period(self): + r"""Time of one revolution, :math:`2\pi/\omega`.""" + return 2.0 * sympy.pi.evalf() / self.omega + + class TwoLayerDarcy(_Transport): r"""Steady Darcy flow through two layers of different permeability. From 6d8e2b752528aa98d9c3392b641ef5674ecb1f64 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 04/14] Extract the per-element timestep estimate shared by the advection-diffusion solvers The cell-crossing / diffusion-time reduction (isotropic or direction-aware, minimum or percentile) becomes a module-level helper so the Eulerian solver can call it rather than carrying a copy. SLCN behaviour unchanged. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- src/underworld3/systems/solvers.py | 204 +++++++++++++++-------------- 1 file changed, 105 insertions(+), 99 deletions(-) diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index e8b82d845..2044c8571 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -321,6 +321,104 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True): return vel +def _advective_diffusive_dt(constitutive_K, V_fn, mesh, direction_aware=False, + percentile=0.0): + r"""Per-element resolution timestep, reduced to one global value. + + The minimum over cells of the advective crossing time :math:`h/|v|` and + the diffusive time :math:`h^2/\kappa`, nondimensional. Shared by the + semi-Lagrangian and the Eulerian advection-diffusion solvers: for both + it is a *resolution* estimate, not a stability limit. The semi-Lagrangian + scheme is unconditionally stable and the implicit Eulerian scheme is + stable at any cell Courant number; what bounds either one is accuracy + on the feature being transported, which the mesh cannot know. + + Parameters + ---------- + constitutive_K : sympy expression or number + Diffusivity (the constitutive model's unified ``K``). + V_fn : sympy Matrix + Advecting velocity, evaluated at cell centroids. + mesh : Mesh + direction_aware : bool, default False + Use the per-cell extent along the local velocity instead of the + isotropic radius (triangles only; falls back otherwise). + percentile : float, default 0.0 + ``0`` takes the strict global minimum; ``> 0`` takes that global + percentile of the per-element timesteps, so a few sliver cells + cannot collapse the estimate. + + Returns + ------- + (dt, dt_adv, dt_diff) : floats + The estimate and its two components; ``inf`` where a component does + not apply (zero velocity, zero diffusivity). + """ + from mpi4py import MPI + + comm = uw.mpi.comm + + diffusivity_glob = _global_max_diffusivity(constitutive_K, mesh) + vel = _centroid_velocities_nd(V_fn, mesh) + vel_magnitudes = np.linalg.norm(vel, axis=1) + element_radii = mesh._radii + + def _reduce_dt(per_elem): + fin = per_elem[np.isfinite(per_elem)] if len(per_elem) else per_elem + if percentile and percentile > 0: + gathered = comm.allgather(np.ascontiguousarray(fin, dtype=float)) + allv = (np.concatenate([a for a in gathered if a.size]) + if any(a.size for a in gathered) else np.empty(0)) + return float(np.percentile(allv, percentile)) if allv.size else np.inf + loc = float(np.min(fin)) if len(fin) else np.inf + return comm.allreduce(loc, op=MPI.MIN) + + if diffusivity_glob > 0: + dt_diff_per_element = (element_radii ** 2) / diffusivity_glob + else: + dt_diff_per_element = np.array([np.inf]) + + if direction_aware: + from underworld3.meshing.smoothing import _tri_cells + tris = _tri_cells(mesh.dm) + if tris is None: + h_per_element = element_radii + else: + coords = np.asarray(mesh.X.coords) + centroids = coords[tris].mean(axis=1) + vhat = np.where( + vel_magnitudes[:, None] > 0, + vel / np.maximum(vel_magnitudes[:, None], 1.0e-30), + 0.0) + D = coords[tris] - centroids[:, None, :] + # Signed projections of the cell vertices along v-hat: the + # extent material actually traverses through the cell. + s = np.einsum('cvd,cd->cv', D, vhat) + h_per_element = np.maximum(s.max(axis=1) - s.min(axis=1), 0.0) + else: + h_per_element = element_radii + + with np.errstate(divide='ignore', invalid='ignore'): + dt_adv_per_element = np.where( + vel_magnitudes > 0, h_per_element / vel_magnitudes, np.inf) + + dt_diff = _reduce_dt(dt_diff_per_element) + dt_adv = _reduce_dt(dt_adv_per_element) + return min(dt_diff, dt_adv), dt_adv, dt_diff + + +def _dimensionalise_dt(dt_estimate): + """Return a timestep estimate with physical time units when a model with + reference scales is active, otherwise as a plain nondimensional scalar.""" + try: + return uw.dimensionalise(np.squeeze(dt_estimate), {'[time]': 1}) + except Exception: + # Sanctioned fallback: no active scaling model. _as_scalar because + # np.squeeze promotes a Python float to a 0-d array, which is not a + # number any caller expects (see _apply_unit_aware_scaling). + return _as_scalar(np.squeeze(dt_estimate)) + + def _invalidate_solution_cache(u): """Drop the cached data view of a just-solved variable. @@ -4302,111 +4400,19 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): with reference scales is available, otherwise nondimensional. """ - ### required modules - from mpi4py import MPI - - comm = uw.mpi.comm - - ## global max diffusivity (unified .K property: diffusivity for - ## diffusion models) - diffusivity_glob = _global_max_diffusivity( - self.constitutive_model.K, self.mesh) - - ### velocity values at element centroids (nondimensional) - vel = _centroid_velocities_nd(self.V_fn, self.mesh) - - # Get per-element velocity magnitudes - vel_magnitudes = np.linalg.norm(vel, axis=1) - - # Get per-element radii (characteristic element size) - element_radii = self.mesh._radii - - ## estimate dt of adv and diff components using per-element approach - ## dt_adv_i = h_i / |v_i| for advection - ## dt_diff_i = h_i^2 / κ for diffusion (using global κ for now) - - # Reduce per-element dt to one global value. Default (percentile=0) = - # strict global MINIMUM — one cell sets the limit. percentile>0 takes the - # Nth global percentile (50 = median) of the per-element dt instead, so a - # few anisotropic SLIVER cells (velocity ACROSS a thin cell) don't collapse - # dt. SLCN is unconditionally stable, and ``direction_aware`` already - # credits cells stretched ALONG the flow — together they give the - # orientation-aware + sliver-robust timestep. - def _reduce_dt(per_elem): - fin = per_elem[np.isfinite(per_elem)] if len(per_elem) else per_elem - if percentile and percentile > 0: - gathered = comm.allgather(np.ascontiguousarray(fin, dtype=float)) - allv = (np.concatenate([a for a in gathered if a.size]) - if any(a.size for a in gathered) else np.empty(0)) - return float(np.percentile(allv, percentile)) if allv.size else np.inf - loc = float(np.min(fin)) if len(fin) else np.inf - return comm.allreduce(loc, op=MPI.MIN) - - # Per-element diffusive timestep (all elements use same diffusivity) - if diffusivity_glob > 0: - dt_diff_per_element = (element_radii ** 2) / diffusivity_glob - else: - dt_diff_per_element = np.array([np.inf]) - - # Per-element advective timestep — either isotropic - # (mesh._radii / |v|) or direction-aware (v-aligned cell - # extent / |v|). - if direction_aware: - # Per-cell vertex indices (triangle / tet). - from underworld3.meshing.smoothing import _tri_cells - tris = _tri_cells(self.mesh.dm) - if tris is None: - # Fall back to isotropic for non-triangle meshes. - h_per_element = element_radii - else: - coords = np.asarray(self.mesh.X.coords) - centroids = coords[tris].mean(axis=1) - # v-hat per cell (use centroid v we already have) - vhat = np.where( - vel_magnitudes[:, None] > 0, - vel / np.maximum(vel_magnitudes[:, None], - 1.0e-30), - 0.0) - D = coords[tris] - centroids[:, None, :] - # Signed projections along v̂ per cell vertex - s = np.einsum('cvd,cd->cv', D, vhat) - h_per_element = s.max(axis=1) - s.min(axis=1) - # Sanity-floor — for zero-velocity cells s=0 - # ⇒ h_eff=0 ⇒ dt_adv=inf via the where below - h_per_element = np.maximum( - h_per_element, 0.0) - else: - h_per_element = element_radii - - with np.errstate(divide='ignore', invalid='ignore'): - dt_adv_per_element = np.where( - vel_magnitudes > 0, - h_per_element / vel_magnitudes, - np.inf - ) - # Global reduction — strict min (percentile=0) or Nth percentile (median). - min_dt_diff_glob = _reduce_dt(dt_diff_per_element) - min_dt_adv_glob = _reduce_dt(dt_adv_per_element) + dt_estimate, dt_adv, dt_diff = _advective_diffusive_dt( + self.constitutive_model.K, self.V_fn, self.mesh, + direction_aware=direction_aware, percentile=percentile) # Store for user inspection - self.dt_adv = min_dt_adv_glob if not np.isinf(min_dt_adv_glob) else 0.0 - self.dt_diff = min_dt_diff_glob if not np.isinf(min_dt_diff_glob) else 0.0 + self.dt_adv = dt_adv if not np.isinf(dt_adv) else 0.0 + self.dt_diff = dt_diff if not np.isinf(dt_diff) else 0.0 - # Take overall minimum (respecting infinity for zero velocity/diffusivity cases) - dt_estimate = min(min_dt_diff_glob, min_dt_adv_glob) - - # If both are infinite (no velocity and no diffusivity), return infinity + # Both infinite (no velocity and no diffusivity): nothing to bound if np.isinf(dt_estimate): return np.inf - # Dimensionalise the result to physical time - try: - return uw.dimensionalise(np.squeeze(dt_estimate), {'[time]': 1}) - except Exception: - # Fallback: return plain nondimensional number. _as_scalar because - # np.squeeze promotes a Python float to a 0-d array, which is not - # a number any caller expects (see _apply_unit_aware_scaling). - return _as_scalar(np.squeeze(dt_estimate)) + return _dimensionalise_dt(dt_estimate) @timing.routine_timer_decorator def solve( From 9c8a125b5e24a55a3c299ada0b2fef765c716c27 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 05/14] Skip the mesh-owned multigrid pickup for a solver that owns its preconditioner A solver with no managed option block (_pc_option_prefix is None) sets its own PC; installing the adapt child's PCMG hierarchy on it segfaulted inside PETSc (additive-Schwarz PC, PCMG calls). The gate now treats that state as the explicit choice it is, alongside preconditioner='gamg' and the user override latch. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- src/underworld3/utilities/custom_mg.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index e463e31c8..4c57b1ec0 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -1783,8 +1783,13 @@ def build_transfers(solver, field_id=None): # `return` here is what turned the gate into a TypeError at the call site # when this hunk migrated from auto_inject_custom_mg (which returns nothing) # during the #488 x #471 merge. + # A solver with no managed option block (`_pc_option_prefix is None`) + # owns its PC outright, so the pickup would install a PCMG hierarchy + # on a PC of another type (measured: SEGV in _configure_pcmg with an + # additive-Schwarz PC on an adapt child). if (getattr(solver, "_preconditioner", "auto") == "gamg" - or getattr(solver, "_pc_user_override", False)): + or getattr(solver, "_pc_user_override", False) + or getattr(solver, "_pc_option_prefix", "") is None): return None, None level_tail = list(coarse) builder = getattr(solver.mesh, "_custom_mg_builder", "barycentric") From b68653069fda5acbd909df3be4676878febd155a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:29:58 -0700 Subject: [PATCH 06/14] Eulerian advection-diffusion with SUPG: BDF and Adams-Moulton orders from the symbolic history uw.systems.AdvDiffusionSUPG(mesh, T, V_fn, order=N, integrator='bdf'|'am') assembles the implicit weak form from the Eulerian DDt history: the BDF stencil or the Adams-Moulton weights on the advective and diffusive terms at every stored time level, plus the SUPG flux tau R u with the strong residual of the same scheme. Timestep, multistep coefficients and the tau weights are runtime constants of the compiled kernels, so a change of dt costs nothing (the issue #657 prototype recompiled on every change). Diffusivity comes from the constitutive model like every scalar solver. Measured on the rotating Gaussian: stable at any cell Courant number, error set by u dt against the feature width (dt^2 for the second-order schemes), unchanged to three digits by a band refined to h/9 at local Courant 13; Crank-Nicolson reproduces the prototype's numbers to four digits. Tests: API and no-recompile contract, temporal convergence (slopes 0.8/0.9 for BDF1, 1.9 for BDF2), band invariance, round trip, np=2 = serial. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../Ex_AdvectionDiffusionSUPG_RotationTest.py | 171 ++++++ src/underworld3/systems/__init__.py | 3 + .../systems/advection_diffusion_eulerian.py | 489 ++++++++++++++++++ .../test_1077_advdiff_supg_parallel.py | 49 ++ tests/test_1055_advdiff_supg_api.py | 165 ++++++ ...est_1100_advdiff_supg_rotating_gaussian.py | 121 +++++ 6 files changed, 998 insertions(+) create mode 100644 docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py create mode 100644 src/underworld3/systems/advection_diffusion_eulerian.py create mode 100644 tests/parallel/test_1077_advdiff_supg_parallel.py create mode 100644 tests/test_1055_advdiff_supg_api.py create mode 100644 tests/test_1100_advdiff_supg_rotating_gaussian.py diff --git a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py new file mode 100644 index 000000000..d5a0db03e --- /dev/null +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py @@ -0,0 +1,171 @@ +# --- +# jupyter: +# jupytext: +# formats: py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% [markdown] +""" +# Eulerian SUPG Advection-Diffusion Rotation Test + +**PHYSICS:** convection +**DIFFICULTY:** advanced + +## Description + +A Gaussian anomaly carried round the origin by rigid rotation, solved with +the fully implicit Eulerian solver `uw.systems.AdvDiffusionSUPG`. The exact +solution is known at every time (`uw.analytic.RotatingGaussian`), so the +error is measured directly rather than inferred from a picture. + +The scheme is stable at any cell Courant number; what limits the timestep +is how far the anomaly moves per step relative to its own width. Try +`-uw_courant 4` to see the accuracy fall off as `dt**2` while the solve +stays perfectly stable, and `-uw_order 2` to see the second-order scheme. + +## Key Concepts + +- **Implicit Eulerian transport**: no trace-back, no departure points; the + timestep is a runtime constant of the compiled kernels. +- **SUPG stabilisation**: the streamline-upwind test-function perturbation + written as a flux, so PETSc needs no modified test space. +- **Multistep order**: `order=1, 2, 3` with `integrator="bdf"` or `"am"`. + +## Parameters + +- `uw_res`: cells across the box +- `uw_courant`: timestep as a multiple of the cell-crossing time +- `uw_order`, `uw_integrator`, `uw_theta`: the time scheme +- `uw_diffusivity`: thermal diffusivity (0 is pure advection) +""" + +# %% +import numpy as np +import sympy +import underworld3 as uw + +# %% [markdown] +""" +## Configurable Parameters + +Override from the command line: + +```bash +python Ex_AdvectionDiffusionSUPG_RotationTest.py -uw_courant 4 -uw_order 2 +``` +""" + +# %% +params = uw.Params( + uw_res=32, + uw_courant=1.0, + uw_order=2, + uw_integrator="bdf", + uw_theta=1.0, + uw_diffusivity=0.0, + uw_sigma=0.12, +) + +# %% [markdown] +""" +## Mesh, exact solution and the transported field +""" + +# %% +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=2.0 / params.uw_res, qdegree=3) +x, y = mesh.X + +exact = uw.analytic.RotatingGaussian( + mesh, sigma=params.uw_sigma, centre_radius=0.5, omega=1.0, + diffusivity=params.uw_diffusivity) + +T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) +T.array[:, 0, 0] = uw.function.evaluate(exact.at(0.0), T.coords).reshape(-1) + +# Rigid rotation about the origin, one revolution in 2 pi +velocity = sympy.Matrix([[-y, x]]) + +# %% [markdown] +""" +## The solver + +Diffusivity is set on the constitutive model, as for every scalar solver. The +walls carry T = 0, which is exact to rounding a few sigma from the orbit. +""" + +# %% +adv_diff = uw.systems.AdvDiffusionSUPG( + mesh, T, velocity, order=params.uw_order, + integrator=params.uw_integrator, theta=params.uw_theta) +adv_diff.constitutive_model.Parameters.diffusivity = params.uw_diffusivity +for boundary in ("Left", "Right", "Top", "Bottom"): + adv_diff.add_dirichlet_bc(0.0, boundary) + +# %% [markdown] +""" +## Time loop + +`estimate_dt` returns the cell-crossing time. It is a resolution guide, not a +stability limit, so the timestep is a chosen multiple of it. For a multistep +scheme the exact history is planted so the first step already runs at full +order. +""" + +# %% +period = float(exact.period) +dt_cell = float(adv_diff.estimate_dt()) +n_steps = int(np.ceil(period / (params.uw_courant * dt_cell))) +dt = period / n_steps + +if params.uw_order > 1: + history = [uw.function.evaluate(exact.at(-k * dt), T.coords).reshape(-1, 1, 1) + for k in range(params.uw_order)] + adv_diff.DuDt.set_initial_history(history, dt=dt) + +t = 0.0 +for step in range(n_steps): + adv_diff.solve(timestep=dt) + t += dt + if step % max(1, n_steps // 4) == 0 or step == n_steps - 1: + err = exact.error(exact.at(t), T, norm="integral") + uw.pprint(f"step {step:4d} t = {t:6.3f} relative L2 error = {err:.3e}") + +# %% [markdown] +""" +## Result + +After one revolution the field should match its initial state. At a Courant +number of one half the round-trip error is below one per cent on this mesh; +it grows as `dt**2` from there. +""" + +# %% +round_trip = exact.error(exact.at(t), T, norm="integral") +uw.pprint(f"round-trip relative L2 error: {round_trip:.3e} " + f"(min {float(T.array.min()):.3f}, max {float(T.array.max()):.3f})") + +# %% +if uw.mpi.size == 1: + import pyvista as pv + import underworld3.visualisation as vis + + pvmesh = vis.mesh_to_pv_mesh(mesh) + pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, T.sym) + pvmesh.point_data["T_exact"] = vis.scalar_fn_to_pv_points(pvmesh, exact.at(t)) + pvmesh.point_data["error"] = pvmesh.point_data["T"] - pvmesh.point_data["T_exact"] + + pl = pv.Plotter(window_size=(900, 450), shape=(1, 2)) + pl.subplot(0, 0) + pl.add_mesh(pvmesh, scalars="T", cmap="RdBu_r", clim=(0, 1), show_edges=False) + pl.subplot(0, 1) + pl.add_mesh(pvmesh, scalars="error", cmap="RdBu_r", show_edges=False) + pl.show(cpos="xy") diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index ea6823ab5..ed7ee58fc 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -19,6 +19,8 @@ L2 projection of fields onto mesh variables. AdvDiffusion : class Advection-diffusion with semi-Lagrangian transport. +AdvDiffusionSUPG : class + Advection-diffusion, implicit Eulerian with SUPG stabilisation. NavierStokes : class Navier-Stokes equations with inertia. Diffusion : class @@ -64,6 +66,7 @@ # These are now implemented the same way using the ddt module from .solvers import SNES_AdvectionDiffusion as AdvDiffusionSLCN from .solvers import SNES_AdvectionDiffusion as AdvDiffusion +from .advection_diffusion_eulerian import SNES_AdvectionDiffusion_SUPG as AdvDiffusionSUPG # import diffusion-only solver from .solvers import SNES_Diffusion as Diffusion diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py new file mode 100644 index 000000000..2dc77ab2f --- /dev/null +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -0,0 +1,489 @@ +r"""Fully implicit Eulerian advection-diffusion with SUPG stabilisation. + +The scalar transport equation + +.. math:: + \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f + +discretised on the mesh with a linear multistep rule in time and a +streamline-upwind Petrov-Galerkin (SUPG) term in space. Every time level +is a mesh variable held by an :class:`~underworld3.systems.ddt.Eulerian` +history manager, so the scheme's order is a construction argument and the +timestep and multistep coefficients are runtime constants of the compiled +kernels: neither changes the generated code. + +The companion of :class:`~underworld3.systems.solvers.SNES_AdvectionDiffusion` +(semi-Lagrangian). The Eulerian scheme is stable at any cell Courant number +and its accuracy is set by how far the transported feature moves in one +step; the semi-Lagrangian scheme's accuracy is set by how far a characteristic +turns in one step. See ``docs/developer/design/eulerian-supg-transport.md``. +""" + +import numpy as np +import sympy +from typing import Optional + +import underworld3 as uw +import underworld3.timing as timing +from underworld3.systems import SNES_Scalar +from underworld3.utilities._api_tools import Template +from underworld3.function import expression as public_expression +from underworld3.systems.ddt import Eulerian as Eulerian_DDt +from underworld3.systems.solvers import ( + _advective_diffusive_dt, + _dimensionalise_dt, + _invalidate_solution_cache, + _nondimensionalise_timestep, +) + + +def _as_row_vector(V_fn, dim): + """Coerce a velocity expression to a ``(1, dim)`` sympy row Matrix.""" + if isinstance(V_fn, uw.discretisation.MeshVariable): + V_fn = V_fn.sym + if isinstance(V_fn, sympy.MatrixBase): + if V_fn.shape == (1, dim): + return V_fn + if V_fn.shape == (dim, 1): + return V_fn.T + raise ValueError( + f"V_fn has shape {V_fn.shape} but the mesh is {dim}-D; expected a " + f"(1, {dim}) row vector such as `v.sym` of a vector MeshVariable." + ) + raise ValueError( + f"V_fn must be a (1, {dim}) sympy Matrix or a vector MeshVariable, " + f"not {type(V_fn).__name__}." + ) + + +class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): + r"""Eulerian advection-diffusion solver, implicit in time, SUPG in space. + + .. math:: + \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f + + Two families of time integration are built from the same stored history + :math:`\phi^{n}, \phi^{n-1}, \dots` (real mesh variables, so their + gradients are available inside the kernels): + + ``integrator="bdf"`` (backward differentiation, order 1-3) + + .. math:: + \frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k} + + \mathbf{u}\cdot\nabla\phi^{n+1} + - \nabla\cdot(\kappa\nabla\phi^{n+1}) = f + + ``integrator="am"`` (Adams-Moulton, order 1-3; ``theta`` at order 1) + + .. math:: + \frac{\phi^{n+1}-\phi^{n}}{\Delta t} + + \sum_{k=0}^{N} a_k\left[\mathbf{u}\cdot\nabla\phi^{n+1-k} + - \nabla\cdot(\kappa\nabla\phi^{n+1-k})\right] = f + + Order 1 with ``theta=1`` is backward Euler in both families; + ``integrator="am", order=1, theta=0.5`` is Crank-Nicolson. The BDF + coefficients :math:`c_k` and Adams-Moulton weights :math:`a_k` are the + ones the :class:`~underworld3.systems.ddt.Eulerian` manager maintains; + both ramp from first order over the opening steps unless a history is + planted with ``solver.DuDt.set_initial_history``. A BDF3 request falls + back to variable-step BDF2 whenever consecutive timesteps differ by more + than 5%. + + **Weak form.** With the strong residual of the chosen scheme + :math:`R(\phi)` (time derivative and advection; see below) the residual + assembled through PETSc's pointwise interface is + + .. math:: + f_0 = R(\phi), \qquad + \mathbf{f}_1 = \sum_k w_k\,\kappa\nabla\phi^{n+1-k} + + \tau\,R(\phi)\,\mathbf{u}, + + where :math:`w_k` are the weights of the spatial operator (:math:`w_0 = 1` + for BDF, :math:`w_k = a_k` for Adams-Moulton). The SUPG contribution is + the Petrov-Galerkin test-function perturbation + :math:`\tau\,\mathbf{u}\cdot\nabla w` written as a flux against + :math:`\nabla w`, so PETSc needs no modified test space. The strong + residual carries no diffusion term because the pointwise kernels see + first derivatives only; for linear elements that term vanishes + identically, for higher orders it is the usual inconsistency of SUPG + without a Laplacian reconstruction. + + **Stabilisation parameter.** + + .. math:: + \tau = \left[\left(\frac{2 c_0}{\Delta t}\right)^2 + + \left(\frac{2|\mathbf{u}|}{h}\right)^2 + + \left(\frac{4\kappa}{h^2}\right)^2\right]^{-1/2} + + with :math:`h` the local cell size (``mesh.cell_size()``, the + :math:`\mathrm{volume}^{1/d}` equivalent radius) and :math:`c_0` the + leading multistep coefficient. The three weights are runtime constants + (``tau_weights``) and ``supg_weight`` scales the whole term, so a + Galerkin baseline needs no rebuild. + + **What limits the timestep.** Nothing, for stability. The implicit + scheme is stable at any cell Courant number, including on cells refined + for a Stokes problem that the scalar does not need. Accuracy is set by + how far the transported feature moves per step relative to its own + width: the error grows as :math:`(\mathbf{u}\Delta t)^2` for the + second-order schemes, and the transient term of :math:`\tau` cannot + hide that. :meth:`estimate_dt` returns the cell-crossing time as a + resolution guide only. + + Parameters + ---------- + mesh : Mesh + u_Field : MeshVariable + Continuous scalar field :math:`\phi`. + V_fn : MeshVariable or sympy Matrix + Advecting velocity, ``(1, dim)``. + order : int, default 1 + Order of the time integration, 1 to 3. + integrator : {"bdf", "am"}, default "bdf" + theta : float, default 1.0 + Adams-Moulton blend at order 1 only (0.5 is Crank-Nicolson). Must be + 1.0 for BDF and for Adams-Moulton above order 1. + verbose : bool, default False + DuDt : Eulerian, optional + A pre-built history manager (order at least ``order``, no ``V_fn``). + + Notes + ----- + The diffusivity is set through the constitutive model, as for the other + scalar solvers; the solver starts with a + :class:`~underworld3.constitutive_models.DiffusionModel` at + :math:`\kappa = 0` (pure advection):: + + adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=2) + adv.constitutive_model.Parameters.diffusivity = 1.0e-3 + adv.add_dirichlet_bc(0.0, "Left") + adv.solve(timestep=dt) + + The linear system is nonsymmetric, so the preconditioner defaults are + GMRES with an additive-Schwarz ILU preconditioner rather than the + algebraic multigrid the symmetric scalar solvers use. + """ + + _INTEGRATORS = ("bdf", "am") + + @timing.routine_timer_decorator + def __init__( + self, + mesh: uw.discretisation.Mesh, + u_Field: uw.discretisation.MeshVariable, + V_fn, + order: int = 1, + integrator: str = "bdf", + theta: float = 1.0, + verbose: bool = False, + DuDt: Optional[Eulerian_DDt] = None, + ): + if not u_Field.continuous: + raise ValueError( + "u_Field must be a continuous MeshVariable: the SUPG weak form " + "is continuous Galerkin." + ) + if integrator not in self._INTEGRATORS: + raise ValueError( + f"integrator must be one of {self._INTEGRATORS}, not {integrator!r}." + ) + order = int(order) + if order not in (1, 2, 3): + raise ValueError(f"order must be 1, 2 or 3, not {order}.") + theta = float(theta) + if theta != 1.0 and not (integrator == "am" and order == 1): + raise ValueError( + "theta applies to integrator='am' at order 1 only " + "(0.5 is Crank-Nicolson); higher orders and BDF take theta=1." + ) + + super().__init__(mesh, u_Field, u_Field.degree, verbose, DuDt=DuDt, DFDt=None) + + self.f = sympy.Matrix.zeros(1, 1) + self._integrator = integrator + self._time_order = order + self._theta = theta + self._V_fn = _as_row_vector(V_fn, mesh.dim) + + tag = self.instance_number + self._delta_t = public_expression( + rf"\Delta t_{{{tag}}}", 1.0, "Eulerian advection-diffusion timestep") + self._last_timestep = None + + # SUPG on/off and the three tau weights are runtime constants: the + # compiled kernels read them from PETSc's constants[] array. + self._supg_weight = public_expression( + rf"w^{{\mathrm{{SUPG}}}}_{{{tag}}}", 1.0, "SUPG term weight (0 = Galerkin)") + self._tau_weights = [ + public_expression(rf"C^{{\tau}}_{{t,{tag}}}", 2.0, "tau transient weight"), + public_expression(rf"C^{{\tau}}_{{u,{tag}}}", 2.0, "tau advective weight"), + public_expression(rf"C^{{\tau}}_{{\kappa,{tag}}}", 4.0, "tau diffusive weight"), + ] + + if DuDt is None: + self.Unknowns.DuDt = Eulerian_DDt( + self.mesh, + u_Field, + vtype=uw.VarType.SCALAR, + degree=u_Field.degree, + continuous=u_Field.continuous, + V_fn=None, + theta=theta, + varsymbol=u_Field.symbol, + verbose=verbose, + bcs=self.essential_bcs, + order=order, + smoothing=0.0, + ) + else: + if DuDt.order < order: + raise ValueError( + f"DuDt supplied is order {DuDt.order} but order {order} was requested." + ) + if getattr(DuDt, "V_fn", None) is not None: + raise ValueError( + "DuDt must be built with V_fn=None: advection is assembled " + "implicitly by this solver, not as an explicit history correction." + ) + self.Unknowns.DuDt = DuDt + + # Diffusivity lives on the constitutive model, as for every scalar + # solver; kappa = 0 until the user sets it. + self.constitutive_model = uw.constitutive_models.DiffusionModel + self.constitutive_model.Parameters.diffusivity = 0.0 + + # Nonsymmetric operator: opt out of the managed GAMG/FMG block and + # use GMRES with an additive-Schwarz ILU preconditioner. RCM + # ordering improves the ILU fill on convection-dominated operators. + self._pc_option_prefix = None + self.petsc_options["ksp_type"] = "gmres" + self.petsc_options["ksp_gmres_restart"] = 200 + self.petsc_options["pc_type"] = "asm" + self.petsc_options["sub_pc_type"] = "ilu" + self.petsc_options["sub_pc_factor_mat_ordering_type"] = "rcm" + self.petsc_options["snes_rtol"] = 1.0e-8 + self.petsc_options["snes_max_it"] = 20 + + # ------------------------------------------------------------------ + # Scheme description + # ------------------------------------------------------------------ + + @property + def integrator(self) -> str: + """``"bdf"`` or ``"am"``.""" + return self._integrator + + @property + def order(self) -> int: + """Requested order of the time integration.""" + return self._time_order + + @property + def theta(self) -> float: + """Adams-Moulton blend at order 1 (1.0 backward Euler, 0.5 Crank-Nicolson).""" + return self._theta + + @property + def delta_t(self): + r"""The timestep :math:`\Delta t` as a UW expression (set by :meth:`solve`).""" + return self._delta_t + + @property + def V_fn(self): + """Advecting velocity, ``(1, dim)``.""" + return self._V_fn + + @V_fn.setter + def V_fn(self, value): + self._V_fn = _as_row_vector(value, self.mesh.dim) + self.is_setup = False + + @property + def f(self): + """Volumetric source term.""" + return self._f + + @f.setter + def f(self, value): + self._f = sympy.Matrix((value,)) + self._needs_function_rewire = True + + @property + def supg_weight(self) -> float: + """Scale of the SUPG term: 1 (default) or 0 for plain Galerkin. No rebuild.""" + return float(self._supg_weight.sym) + + @supg_weight.setter + def supg_weight(self, value): + self._supg_weight.sym = float(value) + + @property + def tau_weights(self): + r"""The weights :math:`(C_t, C_u, C_\kappa)` of the three terms in :math:`\tau`.""" + return tuple(float(w.sym) for w in self._tau_weights) + + @tau_weights.setter + def tau_weights(self, values): + ct, cu, ck = (float(v) for v in values) + self._tau_weights[0].sym = ct + self._tau_weights[1].sym = cu + self._tau_weights[2].sym = ck + + # ------------------------------------------------------------------ + # Residual pieces (raw field symbols only, so the Jacobian sees them) + # ------------------------------------------------------------------ + + def _states(self): + r"""``[phi^{n+1}, phi^{n}, phi^{n-1}, ...]`` as scalar field symbols.""" + return [self.u.sym[0]] + [ps.sym[0] for ps in self.DuDt.psi_star] + + def _spatial_weights(self): + """Weight of the spatial operator at each time level of ``_states``.""" + n = len(self.DuDt.psi_star) + if self._integrator == "bdf": + return [sympy.Integer(1)] + [sympy.Integer(0)] * n + return self.DuDt.am_coefficient_expressions[: n + 1] + + def _time_derivative(self): + if self._integrator == "bdf": + return self.DuDt.bdf()[0] / self._delta_t + phi_new, phi_old = self._states()[:2] + return (phi_new - phi_old) / self._delta_t + + def _advection(self): + dim = self.mesh.dim + u = self._V_fn + total = sympy.Integer(0) + for w, phi in zip(self._spatial_weights(), self._states()): + if w == 0: + continue + grad = self.mesh.vector.gradient(phi) + total = total + w * sum(u[0, i] * grad[0, i] for i in range(dim)) + return total + + def _diffusive_flux(self): + r"""``(1, dim)`` flux :math:`\sum_k w_k\,\nabla\phi^{(k)}\cdot\kappa` from the constitutive tensor.""" + dim = self.mesh.dim + c = self.constitutive_model.c + total = sympy.zeros(1, dim) + for w, phi in zip(self._spatial_weights(), self._states()): + if w == 0: + continue + grad = self.mesh.vector.gradient(phi) + total = total + w * (grad * c) + return total + + def _strong_residual(self): + return self._time_derivative() + self._advection() - self._f[0] + + def _scalar_diffusivity(self): + kappa = self.constitutive_model.Parameters.diffusivity + if isinstance(kappa, sympy.MatrixBase): + raise ValueError( + "The SUPG parameter needs a scalar diffusivity; anisotropic " + "diffusion is not supported by this solver." + ) + return kappa + + def _tau(self): + dim = self.mesh.dim + u = self._V_fn + u_mag2 = sum(u[0, i] ** 2 for i in range(dim)) + h = self.mesh.cell_size() + kappa = self._scalar_diffusivity() + if self._integrator == "bdf": + c0 = self.DuDt.bdf_coefficient_expressions[0] + else: + c0 = sympy.Integer(1) + ct, cu, ck = self._tau_weights + transient = (ct * c0 / self._delta_t) ** 2 + advective = (cu * sympy.sqrt(u_mag2) / h) ** 2 + diffusive = (ck * kappa / h ** 2) ** 2 + return self._supg_weight / sympy.sqrt(transient + advective + diffusive + 1.0e-30) + + F0 = Template( + r"f_0(\phi)", + lambda self: sympy.Matrix([[self._strong_residual()]]), + "Strong residual of the time scheme: time derivative, advection and source.", + ) + F1 = Template( + r"\mathbf{F}_1(\phi)", + lambda self: self._diffusive_flux() + self._tau() * self._strong_residual() * self._V_fn, + "Diffusive flux of the time scheme plus the SUPG flux tau R u.", + ) + + # ------------------------------------------------------------------ + # Timestep and solve + # ------------------------------------------------------------------ + + @timing.routine_timer_decorator + def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): + r"""Cell-crossing timestep, as a resolution guide. + + The minimum over cells of :math:`h/|\mathbf{u}|` and + :math:`h^2/\kappa`, exactly as for the semi-Lagrangian solver. It is + not a stability limit for this scheme, and on a mesh refined for + another problem it is far smaller than the timestep the transported + field needs. Choose the timestep from the feature being transported: + :math:`|\mathbf{u}|\Delta t` should be a fraction of its width. + + Parameters + ---------- + direction_aware : bool, default False + Use the per-cell extent along the local velocity. + percentile : float, default 0.0 + Global percentile of the per-cell timesteps instead of the minimum. + """ + dt_estimate, dt_adv, dt_diff = _advective_diffusive_dt( + self.constitutive_model.K, self._V_fn, self.mesh, + direction_aware=direction_aware, percentile=percentile) + self.dt_adv = dt_adv if not np.isinf(dt_adv) else 0.0 + self.dt_diff = dt_diff if not np.isinf(dt_diff) else 0.0 + if np.isinf(dt_estimate): + return np.inf + return _dimensionalise_dt(dt_estimate) + + def solve( + self, + *, + timestep=None, + zero_init_guess: Optional[bool] = None, + verbose: bool = False, + _force_setup: bool = False, + divergence_retries: int = 0, + ): + r"""Advance :math:`\phi` by one step of size ``timestep``. + + ``timestep`` is required and keyword-only. Changing it between calls + updates a runtime constant of the compiled kernels; nothing is + recompiled. + """ + if timestep is None: + raise ValueError( + "solve() requires timestep=
; there is no default timestep." + ) + dt = float(_nondimensionalise_timestep(timestep)) + if dt <= 0.0: + raise ValueError(f"timestep must be positive, not {dt}.") + if dt != self._last_timestep: + self._delta_t.sym = dt + self._last_timestep = dt + + if _force_setup: + self._needs_function_rewire = True + if not self.constitutive_model._solver_is_setup: + self._needs_function_rewire = True + if not self.is_setup: + self._setup_pointwise_functions(verbose) + self._setup_discretisation(verbose) + self._setup_solver(verbose) + + self.DuDt.update_pre_solve(dt, verbose=verbose) + super().solve(zero_init_guess, _force_setup, divergence_retries=divergence_retries) + _invalidate_solution_cache(self.u) + self.DuDt.update_post_solve(dt, verbose=verbose) + + self.is_setup = True + self.constitutive_model._solver_is_setup = True diff --git a/tests/parallel/test_1077_advdiff_supg_parallel.py b/tests/parallel/test_1077_advdiff_supg_parallel.py new file mode 100644 index 000000000..ae9938fd9 --- /dev/null +++ b/tests/parallel/test_1077_advdiff_supg_parallel.py @@ -0,0 +1,49 @@ +"""The Eulerian SUPG solver gives the serial answer on any number of ranks. + +The scheme has no rank-local step: history is a mesh variable, the residual +is assembled by PETSc, the timestep is a runtime constant. So the integral +error against the rotating-Gaussian oracle after a few steps must match a +serial reference to solver tolerance, whatever the partition. + +Run: mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1077_advdiff_supg_parallel.py +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a, pytest.mark.mpi] + +# Serial reference, res 16, BDF2, dt 0.05, 8 steps (recorded with this file; +# np=2 reproduced it to 1.4e-12). +SERIAL_ERROR = 0.0301522514 + + +def _run(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, + qdegree=3, regular=False) + x, y = mesh.X + sol = uw.analytic.RotatingGaussian(mesh, sigma=0.12, centre_radius=0.5, omega=1.0) + T = uw.discretisation.MeshVariable("T1077", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), order=2) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + dt = 0.05 + adv.DuDt.set_initial_history( + [uw.function.evaluate(sol.at(-k * dt), T.coords).reshape(-1, 1, 1) for k in range(2)], + dt=dt) + for _ in range(8): + adv.solve(timestep=dt) + return sol.error(sol.at(8 * dt), T, norm="integral") + + +def test_error_is_partition_independent(): + err = _run() + assert np.isfinite(err) and err < 0.05, err + gathered = uw.mpi.comm.allgather(err) + assert max(gathered) - min(gathered) < 1e-12, gathered + if SERIAL_ERROR is not None: + assert abs(err - SERIAL_ERROR) < 1e-8, (err, SERIAL_ERROR) diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py new file mode 100644 index 000000000..2173fec2e --- /dev/null +++ b/tests/test_1055_advdiff_supg_api.py @@ -0,0 +1,165 @@ +"""API contract of the Eulerian SUPG advection-diffusion solver. + +Structural checks that run in seconds: the export, argument validation, the +scheme assembled from the history manager, and the rule that a change of +timestep is a change of a runtime constant, never a recompile. + +Run: pixi run python -m pytest tests/test_1055_advdiff_supg_api.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +@pytest.fixture(scope="module") +def mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + + +def _solver(mesh, tag, **kwargs): + x, y = mesh.X + T = uw.discretisation.MeshVariable(f"T_{tag}", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate( + sympy.exp(-((x - 0.5) ** 2 + y ** 2) / 0.03), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), **kwargs) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + return adv, T + + +def test_exported_and_constructs(mesh): + adv, _T = _solver(mesh, "a") + assert type(adv).__name__ == "SNES_AdvectionDiffusion_SUPG" + assert adv.integrator == "bdf" and adv.order == 1 + assert isinstance(adv.DuDt, uw.systems.ddt.Eulerian) + assert adv.DuDt.V_fn is None, "advection is implicit, not a history correction" + + +@pytest.mark.parametrize("tag, kwargs, message", [ + ("v0", dict(order=4), "order must be"), + ("v1", dict(integrator="rk4"), "integrator must be"), + ("v2", dict(integrator="bdf", theta=0.5), "theta applies"), + ("v3", dict(integrator="am", order=2, theta=0.5), "theta applies"), +]) +def test_scheme_arguments_are_validated(mesh, tag, kwargs, message): + with pytest.raises(ValueError, match=message): + _solver(mesh, tag, **kwargs) + + +def test_timestep_is_required(mesh): + adv, _T = _solver(mesh, "b") + with pytest.raises(ValueError, match="requires timestep"): + adv.solve() + + +def test_bdf1_diffusive_flux_is_the_constitutive_flux(mesh): + """At order 1 the assembled diffusive flux is exactly the constitutive + model's own flux of the new state; the history weights are inert.""" + adv, _T = _solver(mesh, "c") + adv.constitutive_model.Parameters.diffusivity = 0.7 + difference = adv._diffusive_flux() - adv.constitutive_model.flux.T + assert all(sympy.simplify(e) == 0 for e in difference) + + +def test_am_order2_uses_all_three_time_levels(mesh): + adv, _T = _solver(mesh, "d", integrator="am", order=2) + weights = adv._spatial_weights() + assert len(weights) == 3 + states = adv._states() + assert len(states) == 3 + # every history state appears (through its derivatives) in the advection operator + names = {str(atom.func) for atom in adv._advection().atoms(sympy.Function)} + for s in states[1:]: + assert any(str(s.func) in n for n in names), (s, names) + + +def test_timestep_change_is_a_constant_update_not_a_recompile(mesh): + adv, _T = _solver(mesh, "e", order=2) + adv.solve(timestep=0.01) + key = adv._current_jit_cache_key + names = [getattr(c, "name", str(c)) for c in adv.constants_manifest] + assert any(r"\Delta t" in n for n in names), names + assert any("BDF" in n for n in names), names + adv.solve(timestep=0.013) + assert adv._current_jit_cache_key == key + assert float(adv.delta_t.sym) == 0.013 + + +def test_timestep_change_reaches_the_kernels(mesh): + """A solver stepped 0.01 then 0.02 gives the same field as a fresh solver + stepped 0.02 from the same state: the constant is really updated.""" + adv1, T1 = _solver(mesh, "f1") + adv1.solve(timestep=0.01) + state = np.array(T1.array) + adv1.solve(timestep=0.02) + + adv2, T2 = _solver(mesh, "f2") + T2.array[...] = state + adv2.DuDt.initialise_history() + adv2.solve(timestep=0.02) + # to the linear-solver tolerance (measured 2e-11 against a 2e-2 control) + assert np.allclose(np.asarray(T1.array), np.asarray(T2.array), rtol=0, atol=1e-8) + + # negative control: a different timestep gives a visibly different field + adv3, T3 = _solver(mesh, "f3") + T3.array[...] = state + adv3.DuDt.initialise_history() + adv3.solve(timestep=0.01) + assert np.abs(np.asarray(T2.array) - np.asarray(T3.array)).max() > 1e-3 + + +def test_order_ramps_from_one_unless_history_is_planted(mesh): + adv, T = _solver(mesh, "g", order=2) + adv.solve(timestep=0.01) + assert adv.DuDt.effective_order == 1 + adv.solve(timestep=0.01) + assert adv.DuDt.effective_order == 2 + + adv2, T2 = _solver(mesh, "h", order=2) + adv2.DuDt.set_initial_history([np.array(T2.array), np.array(T2.array)], dt=0.01) + adv2.solve(timestep=0.01) + assert adv2.DuDt.effective_order == 2 + + +def test_galerkin_baseline_needs_no_rebuild(mesh): + adv, _T = _solver(mesh, "i") + adv.solve(timestep=0.01) + key = adv._current_jit_cache_key + adv.supg_weight = 0.0 + adv.solve(timestep=0.01) + assert adv._current_jit_cache_key == key + assert adv.supg_weight == 0.0 + + +def test_solves_on_an_adapt_child_with_its_own_preconditioner(): + """An adapt child carries a mesh-owned multigrid hierarchy that the + solver base installs opportunistically. This solver owns its (additive + Schwarz) preconditioner, so the pickup must be skipped: installing a + PCMG hierarchy on a non-MG preconditioner segfaulted inside PETSc.""" + base = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3, + refinement=1) + x, y = base.X + + def metric(pts): + h = np.where(np.abs(pts[:, 0]) < 0.1, 0.03, 0.125) + return 1.0 / h ** 2 + + child = base.adapt(metric, max_levels=2) + xc, yc = child.X + T = uw.discretisation.MeshVariable("T_child", child, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate( + sympy.exp(-((xc - 0.5) ** 2 + yc ** 2) / 0.03), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(child, T, sympy.Matrix([[-yc, xc]])) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + adv.solve(timestep=0.02) + assert adv.snes.getKSP().getPC().getType() == "asm" + assert adv._custom_mg is None + data = np.asarray(T.array[:, 0, 0]) + assert np.isfinite(data).all() and 0.9 < data.max() < 1.01 diff --git a/tests/test_1100_advdiff_supg_rotating_gaussian.py b/tests/test_1100_advdiff_supg_rotating_gaussian.py new file mode 100644 index 000000000..5607d22a6 --- /dev/null +++ b/tests/test_1100_advdiff_supg_rotating_gaussian.py @@ -0,0 +1,121 @@ +"""The Eulerian SUPG solver against the rotating Gaussian. + +Three properties measured on ``uw.analytic.RotatingGaussian`` (rigid rotation, +exact at every time): + +1. temporal order: the error at a quarter turn falls as dt (BDF1) and dt^2 + (BDF2) when the timestep is halved, with the exact history planted so + the multistep scheme runs at full order from the first step; +2. mesh refinement the scalar does not need leaves the answer alone: a band + refined to h/8 across the orbit, at the same timestep, gives the same + error to three digits even though its cells sit at a local Courant + number of several; +3. the round trip: after one revolution the field returns to its initial + state to a few per cent at a Courant number of one half. + +Run: pixi run python -m pytest tests/test_1100_advdiff_supg_rotating_gaussian.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +SIGMA = 0.12 + + +def _box(res, refinement=0): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=2.0 / res, + qdegree=3, regular=False, refinement=refinement) + + +def _problem(mesh, tag, order, integrator="bdf", theta=1.0, kappa=0.0): + x, y = mesh.X + sol = uw.analytic.RotatingGaussian(mesh, sigma=SIGMA, centre_radius=0.5, + omega=1.0, diffusivity=kappa) + T = uw.discretisation.MeshVariable(f"T_{tag}", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), + order=order, integrator=integrator, theta=theta) + adv.constitutive_model.Parameters.diffusivity = kappa + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + return sol, T, adv + + +def _run(sol, T, adv, dt, t_end, plant=True): + nsteps = int(round(t_end / dt)) + dt = t_end / nsteps + if plant and adv.order > 1: + values = [uw.function.evaluate(sol.at(-k * dt), T.coords).reshape(-1, 1, 1) + for k in range(adv.order)] + adv.DuDt.set_initial_history(values, dt=dt) + for _ in range(nsteps): + adv.solve(timestep=dt) + return sol.error(sol.at(t_end), T, norm="integral") + + +@pytest.mark.parametrize("order, timesteps, expected_slope", [ + (1, (0.02, 0.01, 0.005), 1.0), + (2, (0.04, 0.02, 0.01), 2.0), +]) +def test_temporal_convergence_order(order, timesteps, expected_slope): + """Halving dt divides the quarter-turn error by 2 (BDF1) or 4 (BDF2). + + The timesteps sit where the temporal error dominates the fixed spatial + error but is still in its asymptotic range (backward Euler at + u dt > sigma/2 is already saturated), which is why the slope is checked + with a tolerance. + """ + mesh = _box(32) + t_end = float(sympy.pi) / 2 + errors = [] + for i, dt in enumerate(timesteps): + sol, T, adv = _problem(mesh, f"c{order}{i}", order) + errors.append(_run(sol, T, adv, dt, t_end)) + slopes = np.log2(np.array(errors[:-1]) / np.array(errors[1:])) + print(f"order {order}: errors {errors} slopes {slopes}") + assert slopes.min() > expected_slope - 0.35, (order, errors, slopes) + + +def test_refinement_the_scalar_does_not_need_leaves_the_error_alone(): + """A band at h/8 across the orbit, same dt as the uniform mesh.""" + dt = 0.0433 + t_end = float(sympy.pi) / 2 + + uniform = _box(32) + sol, T, adv = _problem(uniform, "u", 2) + err_uniform = _run(sol, T, adv, dt, t_end) + + base = _box(16, refinement=1) + fault = uw.meshing.Surface("band", base, + np.array([[0.0, -1.0, 0.0], [0.0, 1.0, 0.0]]), symbol="F") + fault.discretize() + h = 1.0 / 16 + + def metric(pts, _f=fault, _hn=h / 8, _hf=h, _core=0.03, _ramp=0.06): + d = _f.unsigned_distance(pts) + hh = np.where(d < _core, _hn, np.minimum(_hn + (_hf - _hn) * (d - _core) / _ramp, _hf)) + return 1.0 / hh ** 2 + + child = base.adapt(metric, max_levels=3) + assert float(np.min(child._radii)) < 0.3 * float(np.min(uniform._radii)) + + sol_c, T_c, adv_c = _problem(child, "b", 2) + err_band = _run(sol_c, T_c, adv_c, dt, t_end) + + # the band cells are at a local Courant number well above one + assert dt / float(adv_c.estimate_dt()) > 4.0 + assert abs(err_band - err_uniform) < 0.15 * err_uniform, (err_uniform, err_band) + + +def test_round_trip_at_moderate_courant(): + mesh = _box(32) + sol, T, adv = _problem(mesh, "r", 2) + err = _run(sol, T, adv, 0.5 * float(adv.estimate_dt()), float(sol.period)) + assert err < 0.03, err + data = np.asarray(T.array[:, 0, 0]) + assert data.min() > -0.02 and data.max() < 1.02, (data.min(), data.max()) From 451efe3fac16bb17caa3270c90f5e83d91d6cbc3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 17:34:29 -0700 Subject: [PATCH 07/14] Design note for the Eulerian SUPG solver; BDF2 becomes the default from the integrator study Rotating-Gaussian study at res 32, Courant 0.25 to 8, pure advection and kappa 1e-3: Adams-Moulton above order 1 blows up from Courant 1 (bounded stability region), BDF3 fails from Courant 4, Crank-Nicolson is three to four times more accurate than BDF2 at the same timestep but rings once the feature is under-resolved in time, backward Euler carries 20-40% error at any practical timestep. Cost per step is the same for every scheme. BDF2 is the robust default; the note records the alternatives and when to pick them. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../semi-lagrangian-time-integration.md | 13 ++ .../design/eulerian-supg-transport.md | 199 ++++++++++++++++++ docs/developer/index.md | 1 + .../systems/advection_diffusion_eulerian.py | 16 +- tests/test_1055_advdiff_supg_api.py | 2 +- 5 files changed, 227 insertions(+), 4 deletions(-) create mode 100644 docs/developer/design/eulerian-supg-transport.md diff --git a/docs/advanced/semi-lagrangian-time-integration.md b/docs/advanced/semi-lagrangian-time-integration.md index 23a935db2..7812b7c78 100644 --- a/docs/advanced/semi-lagrangian-time-integration.md +++ b/docs/advanced/semi-lagrangian-time-integration.md @@ -112,6 +112,19 @@ $[\theta,\,1-\theta]$: `theta` is settable after construction: `adv_diff.DFDt.theta = 1.0`. +## The Eulerian alternative + +`uw.systems.AdvDiffusionSUPG` solves the same equation without a trace-back: +all terms are assembled on the mesh, implicit in time, with SUPG +stabilisation. Its `order=` and `integrator="bdf"|"am"` arguments select the +multistep scheme, built from the same stored history as above; `order=1, +integrator="am", theta=0.5` is Crank-Nicolson. The scheme is stable at any +cell Courant number, so cells refined for a Stokes problem never limit the +transport timestep; its accuracy is set by how far the transported feature +moves per step. The semi-Lagrangian scheme's accuracy is instead set by how +far a characteristic turns per step. The measurements behind that split are +in `docs/developer/design/eulerian-supg-transport.md`. + ## Related options - **`monotone_mode`** (`"clamp"` / `"pick"`) bounds the semi-Lagrangian diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md new file mode 100644 index 000000000..268c108a6 --- /dev/null +++ b/docs/developer/design/eulerian-supg-transport.md @@ -0,0 +1,199 @@ +# Eulerian SUPG transport: design and measurements + +**Status**: implemented on `feature/eulerian-supg-transport` (2026-09-02), static mesh. +Supersedes the Crank-Nicolson prototype of issue #657 as the implementation route +while keeping its weak-form idea. + +## Why an Eulerian scheme + +Underworld3 meshes are usually refined for the momentum problem: faults, viscosity +jumps, boundary layers. A transported scalar rarely needs that resolution, so a +scheme whose timestep is bounded by the smallest cell pays for cells it does not +use. The semi-Lagrangian solver (`AdvDiffusionSLCN`) escapes that bound but pays +for departure points, which are expensive per step and irregular in parallel, and +its moving-mesh staging needs a lagged copy of the previous geometry. + +An implicit Eulerian scheme has no stability bound at all. Its cost is a +nonsymmetric solve per step, and its accuracy is bounded by how far the transported +feature moves in one step. The measurements below say when each is the better tool. + +## The scheme + +The equation is + +$$ +\frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f . +$$ + +Every past time level $\phi^{n}, \phi^{n-1}, \dots$ is a mesh variable held by an +`Eulerian` history manager, so first derivatives of past states are available in +the kernels and two multistep families share one code path: + +| `integrator` | time derivative | spatial operator | +|---|---|---| +| `bdf`, order $N$ | $\frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k}$ | at $n+1$ only | +| `am`, order $N$ | $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$ | $\sum_{k=0}^{N} a_k\,S(\phi^{n+1-k})$ | + +with $S(\phi) = \mathbf{u}\cdot\nabla\phi - \nabla\cdot(\kappa\nabla\phi)$ and the +coefficients those the history manager already maintains (`theta` is the +Adams-Moulton weight at order 1; 0.5 is Crank-Nicolson). Both families ramp from +first order over the opening steps unless `solver.DuDt.set_initial_history` plants +the history. The pointwise residual is + +$$ +f_0 = R(\phi), \qquad +\mathbf{f}_1 = \sum_k w_k\,\kappa\nabla\phi^{n+1-k} + \tau\,R(\phi)\,\mathbf{u}, +$$ + +where $R$ is the strong residual of the chosen scheme (time derivative, advection +and source) and $w_k$ the spatial weights of the family. The SUPG term is the +Petrov-Galerkin test-function perturbation $\tau\,\mathbf{u}\cdot\nabla w$ written +as a flux against $\nabla w$, so PETSc needs no modified test space. + +$$ +\tau = \left[\left(\frac{2 c_0}{\Delta t}\right)^2 + + \left(\frac{2|\mathbf{u}|}{h}\right)^2 + + \left(\frac{4\kappa}{h^2}\right)^2\right]^{-1/2}, +\qquad h = \texttt{mesh.cell\_size()} . +$$ + +### Decisions and their reasons + +- **No diffusion in the strong residual.** PETSc's pointwise kernels see first + derivatives only, so $-\nabla\cdot(\kappa\nabla\phi)$ cannot appear in $R$. For + linear elements it vanishes identically; for higher orders this is the usual + inconsistency of SUPG without a Laplacian reconstruction. Diffusion enters as the + Galerkin flux only. +- **Every knob is a runtime constant.** The timestep, the multistep coefficients, + the three weights in $\tau$ and the overall SUPG weight are UW expressions routed + through PETSc's `constants[]` array. A change of timestep costs nothing; the + prototype recompiled its kernels on every change (1.2 s against 0.03 s for a step). +- **Diffusivity on the constitutive model**, as for every scalar solver, starting at + $\kappa = 0$. The prototype carried a float attribute with a warning bridge. +- **Own preconditioner.** The operator is nonsymmetric, so GMRES with an + additive-Schwarz ILU preconditioner replaces the managed GAMG block. The solver + sets `_pc_option_prefix = None`, and the mesh-owned multigrid pickup on adapt + children now respects that (it segfaulted otherwise). +- **Moving meshes, phase 1.** The unknown and its history stay on the default + `REMAP` transfer policy with the material velocity. The remap re-interpolates old + states onto the new nodes, so the Eulerian form is already correct to + interpolation accuracy. The `CARRY` + $\mathbf{u} - \mathbf{u}_\text{mesh}$ form + is phase 2 and must not be mixed with `REMAP`. +- **Not yet:** discontinuity capturing (the prototype's residual omitted the time + derivative and added first-order diffusion everywhere; a correct lagged residual + needs $\phi^{n-1}$), a streamline element length from a mesh-owned metric tensor, + the ALE hook. + +## Measurements + +Rotating Gaussian (`uw.analytic.RotatingGaussian`, $\sigma = 0.12$, orbit radius +0.5), P2 field, unstructured simplex box, one revolution; relative $L_2$ error at +the end. "Courant" is on the cell size. Study scripts and CSVs are in +`~/+Simulations/supg_vs_slcn_657/`. + +### Eulerian against semi-Lagrangian (the #657 prototype, Crank-Nicolson) + +| mesh | Courant | SUPG CN | SLCN | cost per step SUPG : SLCN | +|---|---|---|---|---| +| uniform 32 | 0.5 | 0.6% | 21% | 1 : 6.3 | +| uniform 32 | 2 | 9.8% | 7.7% | 1 : 6.4 | +| uniform 32 | 8 | 66%, min $-0.35$ | 8.8% | 1 : 5.6 | +| uniform 32 | 32 | 113% | 93%, mass $-32$% | 1 : 5.7 | +| uniform 64 | 2 | 2.5% | 2.2% | 1 : 3.6 | +| uniform 64 | 8 | 31% | 2.2% | 1 : 3.6 | +| band $h/9$ at $x = 0$ | 0.5 / 2 | 0.6% / 9.8% | 18% / 6.5% | 1 : 5.6 | + +Three facts follow. + +1. The implicit scheme is stable at any cell Courant number, and cells the scalar + does not need are free: the band refined to $h/9$ sits at local Courant 13 and + changes the error in the third digit only. +2. Its accuracy is set by $\mathbf{u}\Delta t$ against the feature width. The error + scales as $\Delta t^2$ for Crank-Nicolson, which is A-stable but not L-stable + and rings once the feature is under-resolved in time. +3. SLCN's error is flat in $\Delta t$ but accumulates at small Courant (one + interpolation per step), so it is the worse scheme exactly where it is not meant + to run; its limit is the arc a characteristic turns per step, about 10 degrees + for the RK2 trace-back, a property of the flow rather than the mesh. + +The new class reproduces the prototype's Crank-Nicolson numbers to four digits +(0.5993% and 9.777% at Courant 0.5 and 2 on the uniform mesh). + +### BDF against Adams-Moulton + +`time_integrator_study.py`: the same rotating Gaussian, res 32, every scheme +the class offers, at Courant 0.25 to 8; relative $L_2$ error after one +revolution, "X" where the run blew up (with the step). Pure advection first, +then $\kappa = 10^{-3}$ (cell Peclet about 40). + +| scheme | C 0.25 | 0.5 | 1 | 2 | 4 | 8 | +|---|---|---|---|---|---|---| +| BDF1 = backward Euler | 19% | 30% | 44% | 57% | 68% | 77% | +| BDF2 | 0.6% | 2.4% | 9.3% | 28% | 53% | 73% | +| BDF3 | 0.32% | 0.28% | 2.7% | 18% | X | X | +| Crank-Nicolson (`am`, 1, theta 0.5) | 0.27% | 0.6% | 2.5% | 9.8% | 31% | 66% | +| Adams-Moulton 2 (third order) | 0.28% | 0.24% | 0.24% | X@68 | X@41 | X@32 | +| Adams-Moulton 3 (fourth order) | 0.28% | 0.25% | X@155 | X@32 | X@22 | X@19 | + +| scheme, $\kappa = 10^{-3}$ | C 0.25 | 0.5 | 1 | 2 | 4 | 8 | +|---|---|---|---|---|---|---| +| BDF1 = backward Euler | 12% | 20% | 31% | 45% | 58% | 69% | +| BDF2 | 0.27% | 0.71% | 3.3% | 13% | 35% | 59% | +| BDF3 | 0.31% | 0.45% | 0.87% | 4.4% | 51% | X | +| Crank-Nicolson | 0.38% | 0.51% | 0.63% | 2.5% | 13% | 42% | +| Adams-Moulton 2 | 0.42% | 0.71% | 1.3% | X | X | X | +| Adams-Moulton 3 | 0.42% | 0.71% | X | X | X | X | + +Cost per step is the same for every scheme (0.058 to 0.068 s at res 32): the +history terms are extra kernel inputs, not extra solves. BDF1 and backward Euler +agree to every digit, which checks that the two families are assembled +consistently. + +What the table says: + +- **Adams-Moulton above order 1 is unusable for advection.** Its stability region + is bounded and covers only a short segment of the imaginary axis, so on a pure + advection operator it blows up once the Courant number reaches about 1, and + diffusion at this Peclet number does not rescue it. It is kept in the class for + the record and for diffusion-dominated use, with that warning in the docstring. +- **BDF3 is the most accurate scheme below Courant 1** (0.3%, on the spatial + floor) but it is not A-stable either, and it fails from Courant 4. +- **Crank-Nicolson is three to four times more accurate than BDF2 at the same + timestep** across the usable range, because it does not damp; the price is + ringing once the feature is under-resolved in time (minimum $-0.35$ at Courant 8 + against $-0.20$ for BDF2), and no damping of stiff modes at all. +- **BDF2 is the robust choice**: stable at every Courant number, damped, second + order, and the error is still set by $\mathbf{u}\Delta t$ against the feature + width. + +**Default: `order=2, integrator="bdf"`.** Defaults err toward robustness; a user +with a smooth field at Courant 2 or below gets the better answer from +`integrator="am", order=1, theta=0.5` (Crank-Nicolson), and below Courant 1 from +`order=3`. Backward Euler (order 1) is a first-order scheme with 20 to 40% error +at any practical timestep and is not a sensible default for transport. + +### Temporal convergence (tests/test_1100) + +Quarter-turn error on the uniform res-32 mesh with the exact history planted: +BDF1 slopes 0.80 and 0.88 between $\Delta t$ = 0.02, 0.01, 0.005; BDF2 slopes above +1.65 between 0.04, 0.02, 0.01. + +## What the timestep estimate means + +`estimate_dt` returns the cell-crossing time, the same resolution estimate the +semi-Lagrangian solver reports, because that is the only quantity the mesh knows. +It is not a stability limit for either scheme. Choose the Eulerian timestep from +the transported feature: $|\mathbf{u}|\Delta t$ a fraction of its width. For SLCN +the honest limit is the trace-back arc, $\Delta t \lesssim 0.25 / \max|\nabla\mathbf{u}|$, +which is a separate change to that solver. + +## A defect found on the way + +The API test was flaky only after a test that dropped mesh variables. The cause +is general and predates this work: `mesh.vars` holds variables weakly, a +garbage-collected variable leaves its PETSc field in the DM, and both +`Mesh.update_lvec` and the JIT's auxiliary-field offsets assumed the registry and +the DM fields line up by position. Every later variable was then packed into, and +read from, the wrong slots. Fixed in the same branch (pack by field name, offsets +from the DM's field list) with `tests/test_1058_dropped_meshvariable_aux_layout.py`. diff --git a/docs/developer/index.md b/docs/developer/index.md index c32976e1f..26a102965 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -165,6 +165,7 @@ design/TURBULENCE_MODEL_DESIGN design/declined-coord-units-proposal design/nonlinear-solver-homotopy-warmstart design/fault-zone-hybrid-architecture +design/eulerian-supg-transport ``` ```{toctree} diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 2dc77ab2f..25bdc0209 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -139,9 +139,19 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): Continuous scalar field :math:`\phi`. V_fn : MeshVariable or sympy Matrix Advecting velocity, ``(1, dim)``. - order : int, default 1 - Order of the time integration, 1 to 3. + order : int, default 2 + Order of the time integration, 1 to 3. BDF2 is the default because it + is stable at every Courant number and damped. Measured on a rotating + Gaussian (``docs/developer/design/eulerian-supg-transport.md``): + Crank-Nicolson is three to four times more accurate than BDF2 at the + same timestep below Courant 2 but rings once the feature is + under-resolved in time; BDF3 is the most accurate scheme below + Courant 1 and fails from Courant 4; backward Euler (order 1) carries + 20 to 40% error at any practical timestep. integrator : {"bdf", "am"}, default "bdf" + Adams-Moulton above order 1 has a bounded stability region and blows + up on an advection operator from about Courant 1; it is provided for + diffusion-dominated problems and for comparison. theta : float, default 1.0 Adams-Moulton blend at order 1 only (0.5 is Crank-Nicolson). Must be 1.0 for BDF and for Adams-Moulton above order 1. @@ -174,7 +184,7 @@ def __init__( mesh: uw.discretisation.Mesh, u_Field: uw.discretisation.MeshVariable, V_fn, - order: int = 1, + order: int = 2, integrator: str = "bdf", theta: float = 1.0, verbose: bool = False, diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index 2173fec2e..bb95e3fd4 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -35,7 +35,7 @@ def _solver(mesh, tag, **kwargs): def test_exported_and_constructs(mesh): adv, _T = _solver(mesh, "a") assert type(adv).__name__ == "SNES_AdvectionDiffusion_SUPG" - assert adv.integrator == "bdf" and adv.order == 1 + assert adv.integrator == "bdf" and adv.order == 2 assert isinstance(adv.DuDt, uw.systems.ddt.Eulerian) assert adv.DuDt.V_fn is None, "advection is implicit, not a history correction" From f076781dfc6ca83599254ceab769aa0ddd18f46c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 18:13:03 -0700 Subject: [PATCH 08/14] Integrator study, res 64: BDF3 grows slowly on pure advection at any Courant number BDF2 and Crank-Nicolson track their res-32 errors at the same u dt. BDF3's stability region misses the imaginary axis near the origin, so the low-frequency modes of a finer mesh grow: 31x the exact field after 590 steps at Courant 1. Safe only with diffusion, below Courant 2. Note and docstring updated; the BDF2 default stands. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../design/eulerian-supg-transport.md | 31 ++++++++++++++++--- .../systems/advection_diffusion_eulerian.py | 3 +- 2 files changed, 28 insertions(+), 6 deletions(-) diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 268c108a6..7410e80ab 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -145,7 +145,27 @@ then $\kappa = 10^{-3}$ (cell Peclet about 40). | Adams-Moulton 2 | 0.42% | 0.71% | 1.3% | X | X | X | | Adams-Moulton 3 | 0.42% | 0.71% | X | X | X | X | -Cost per step is the same for every scheme (0.058 to 0.068 s at res 32): the +At res 64 (pure advection, Courant 1 to 8, 590 to 74 steps per revolution): + +| scheme, res 64 | C 1 | 2 | 4 | 8 | +|---|---|---|---|---| +| BDF1 = backward Euler | 30% | 43% | 57% | 68% | +| BDF2 | 2.5% | 9.1% | 27% | 53% | +| BDF3 | 3100% (slow growth) | 1.9% | 17% | 130% | +| Crank-Nicolson | 0.62% | 2.5% | 9.5% | 31% | +| Adams-Moulton 2 | 310% (slow growth) | X@76 | X@49 | X@38 | +| Adams-Moulton 3 | X@116 | X@37 | X@25 | X@22 | + +BDF2 and Crank-Nicolson track their res-32 values at the same $\mathbf{u}\Delta t$ +(the error is set by the timestep, not the mesh). BDF3 is not safe for pure +advection at any Courant number: its stability region misses the imaginary axis +near the origin, so the low-frequency modes a finer mesh carries grow slowly (31 +times the exact field after 590 steps at Courant 1, where the coarser mesh with +half the steps still looked fine); with $\kappa = 10^{-3}$ it behaved. Use it +only with diffusion and below Courant 2. + +Cost per step is the same for every scheme (0.058 to 0.068 s at res 32, 0.32 to +0.36 s at res 64): the history terms are extra kernel inputs, not extra solves. BDF1 and backward Euler agree to every digit, which checks that the two families are assembled consistently. @@ -157,8 +177,9 @@ What the table says: advection operator it blows up once the Courant number reaches about 1, and diffusion at this Peclet number does not rescue it. It is kept in the class for the record and for diffusion-dominated use, with that warning in the docstring. -- **BDF3 is the most accurate scheme below Courant 1** (0.3%, on the spatial - floor) but it is not A-stable either, and it fails from Courant 4. +- **BDF3 is the most accurate scheme below Courant 1 with diffusion present** + (0.3%, on the spatial floor) but it is not A-stable, fails from Courant 4, and + on pure advection grows slowly at any Courant number (the res-64 rows). - **Crank-Nicolson is three to four times more accurate than BDF2 at the same timestep** across the usable range, because it does not damp; the price is ringing once the feature is under-resolved in time (minimum $-0.35$ at Courant 8 @@ -169,8 +190,8 @@ What the table says: **Default: `order=2, integrator="bdf"`.** Defaults err toward robustness; a user with a smooth field at Courant 2 or below gets the better answer from -`integrator="am", order=1, theta=0.5` (Crank-Nicolson), and below Courant 1 from -`order=3`. Backward Euler (order 1) is a first-order scheme with 20 to 40% error +`integrator="am", order=1, theta=0.5` (Crank-Nicolson), and with diffusion below +Courant 1 from `order=3`. Backward Euler (order 1) is a first-order scheme with 20 to 40% error at any practical timestep and is not a sensible default for transport. ### Temporal convergence (tests/test_1100) diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 25bdc0209..6c8e23f41 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -146,7 +146,8 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): Crank-Nicolson is three to four times more accurate than BDF2 at the same timestep below Courant 2 but rings once the feature is under-resolved in time; BDF3 is the most accurate scheme below - Courant 1 and fails from Courant 4; backward Euler (order 1) carries + Courant 1 when diffusion is present, fails from Courant 4, and on + pure advection grows slowly at any Courant number; backward Euler (order 1) carries 20 to 40% error at any practical timestep. integrator : {"bdf", "am"}, default "bdf" Adams-Moulton above order 1 has a bounded stability region and blows From 2b9941c1f58d7bb1917e0ad09fa729cdf1de6727 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 21:16:48 -0700 Subject: [PATCH 09/14] AdvDiffusionSUPG takes the semi-Lagrangian solver's interface: a drop-in replacement The constructor, order, theta, f, V_fn, constitutive_model, delta_t, estimate_dt and solve keep the meaning they have for AdvDiffusionSLCN, so a script changes the class name and nothing else. order=1 with theta=0.5 is Crank-Nicolson and the default, as for SLCN; order=2 takes theta=1 (BDF2) unless 0.5 is asked for explicitly, which is refused for the reason the SLCN documentation gives. The trace-back-only arguments (restore_points_func, monotone_mode, old_frame_traceback, DFDt) are accepted and ignored with a warning. integrator is inferred and only needs setting to reach the higher Adams-Moulton rules. delta_t is settable and solve() reuses it; the notebook viewer reports the scheme. User page docs/advanced/eulerian-advection-diffusion.md with the swap table and the when-to-use-which guidance. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- docs/advanced/eulerian-advection-diffusion.md | 105 ++++++++ docs/advanced/index.md | 1 + .../design/eulerian-supg-transport.md | 16 +- .../Ex_AdvectionDiffusionSUPG_RotationTest.py | 13 +- .../systems/advection_diffusion_eulerian.py | 240 ++++++++++++------ tests/test_1055_advdiff_supg_api.py | 32 ++- 6 files changed, 312 insertions(+), 95 deletions(-) create mode 100644 docs/advanced/eulerian-advection-diffusion.md diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md new file mode 100644 index 000000000..9e0951526 --- /dev/null +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -0,0 +1,105 @@ +# Eulerian advection-diffusion (SUPG): a drop-in for SLCN + +`uw.systems.AdvDiffusionSUPG` solves the same scalar transport equation as the +semi-Lagrangian solver `uw.systems.AdvDiffusionSLCN`, + +$$ +\frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f , +$$ + +but assembles every term on the mesh, implicit in time, with streamline-upwind +(SUPG) stabilisation. There is no trace-back and no departure point. The two +classes share their interface, so switching is one line: + +```python +adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=1) # was AdvDiffusionSLCN +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = 1.0e-3 +adv.add_dirichlet_bc(1.0, "Bottom") +adv.add_dirichlet_bc(0.0, "Top") + +dt = 0.5 * adv.estimate_dt() +adv.solve(timestep=dt) +``` + +## What carries over + +| SLCN | SUPG | note | +|---|---|---| +| `order=1, theta=0.5` | same | Crank-Nicolson, the default for both | +| `order=1, theta=1.0` | same | backward Euler | +| `order=2, theta=1.0` | same | SL-BDF2 becomes BDF2 | +| `order=2, theta=0.5` | refused | refused for the same reason: a BDF stencil does not pair with a centred flux | +| `f`, `V_fn`, `constitutive_model`, `delta_t` | same | | +| `estimate_dt(direction_aware, percentile)` | same | the cell-crossing time, a resolution guide for both | +| `solve(zero_init_guess, timestep, ...)` | same | | +| `DuDt.set_initial_history(values, dt)` | same | plant an exact history to start at full order | +| `restore_points_func`, `monotone_mode`, `old_frame_traceback`, `DFDt` | ignored, with a warning | they configure the trace-back | + +`order=3` (BDF3) is available; see below for when it is safe. `integrator="am"` +above order 1 reaches the higher Adams-Moulton rules, which are for +diffusion-dominated problems only. + +## When to use which + +Both solvers are free of any stability limit on the timestep, so cells refined +for the Stokes problem never dictate the transport step. They differ in what +bounds their accuracy and in what a step costs. + +**Eulerian SUPG.** The error is set by how far the transported feature moves per +step relative to its own width, as $(\mathbf{u}\Delta t)^2$ for the second-order +schemes. It does not depend on the cell size at all: on a rotating Gaussian a band +refined to $h/9$, with its cells at a local Courant number of 13, changes the error +in the third digit only. A step costs one nonsymmetric solve, four to six times +less than a semi-Lagrangian step in serial, and it needs no departure points in +parallel. On a moving mesh the field and its history are re-interpolated by the +ordinary remesh transfer, so no special staging is needed. + +**Semi-Lagrangian.** The error is nearly independent of the timestep but +accumulates one interpolation per step, so at small Courant numbers it is the +worse scheme (21% against 0.6% after one revolution at Courant 0.5 on the same +mesh). Its limit is the arc a characteristic turns per step, about 10 degrees for +the RK2 trace-back, a property of the flow rather than the mesh. Above roughly +Courant 2 on the feature's own scale it keeps its accuracy where the Eulerian +scheme loses it. + +A practical rule: if the timestep is chosen so that the temperature field itself +is resolved in time (a fraction of a feature width per step), the Eulerian solver +is cheaper and more accurate; if the step is deliberately long relative to the +transported features, the semi-Lagrangian solver is the one that survives it. + +## Choosing the time scheme + +Measured on a rotating Gaussian, one revolution, relative $L_2$ error; the full +tables are in the design note. + +| scheme | behaviour | +|---|---| +| Crank-Nicolson (`order=1`) | three to four times more accurate than BDF2 at the same timestep below Courant 2; rings once the feature is under-resolved in time | +| BDF2 (`order=2`) | damped and stable at every Courant number; the choice for sharp or under-resolved fields | +| BDF3 (`order=3`) | the most accurate scheme below Courant 1 when diffusion is present; on pure advection it grows slowly at any Courant number, so use it only with diffusion | +| backward Euler (`order=1, theta=1.0`) | 20 to 40% error at any practical timestep; not for transport | +| Adams-Moulton 2, 3 (`integrator="am"`) | third and fourth order below Courant 1; blow up on advection from about Courant 1 | + +All schemes cost the same per step: the history terms are extra kernel inputs, +not extra solves. Changing the timestep between steps changes a runtime constant +of the compiled kernels; nothing is recompiled. + +## Details that differ from SLCN + +- The strong residual used in the SUPG term carries the time derivative and the + advection but no diffusion term, because PETSc's pointwise kernels see first + derivatives only. For linear elements the missing term is identically zero. +- The stabilisation parameter uses the local cell size (`mesh.cell_size()`) and + three weights that are runtime constants (`solver.tau_weights`); + `solver.supg_weight = 0` gives the plain Galerkin scheme for comparison. +- The linear system is nonsymmetric, so the solver defaults to GMRES with an + additive-Schwarz ILU preconditioner instead of algebraic multigrid. Every + option can be overridden through `solver.petsc_options`. + +## Further reading + +- Design note and measurements: `docs/developer/design/eulerian-supg-transport.md` +- The semi-Lagrangian schemes: {doc}`semi-lagrangian-time-integration` +- Example: `docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py` diff --git a/docs/advanced/index.md b/docs/advanced/index.md index cb47f8dea..b569bf0a7 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -139,6 +139,7 @@ custom-meshes curved-boundary-conditions mesh-adaptation semi-lagrangian-time-integration +eulerian-advection-diffusion porous-flow snapshot-restore troubleshooting diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 7410e80ab..580870ebf 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -188,11 +188,17 @@ What the table says: order, and the error is still set by $\mathbf{u}\Delta t$ against the feature width. -**Default: `order=2, integrator="bdf"`.** Defaults err toward robustness; a user -with a smooth field at Courant 2 or below gets the better answer from -`integrator="am", order=1, theta=0.5` (Crank-Nicolson), and with diffusion below -Courant 1 from `order=3`. Backward Euler (order 1) is a first-order scheme with 20 to 40% error -at any practical timestep and is not a sensible default for transport. +**Interface and default.** The class is a drop-in replacement for the +semi-Lagrangian solver: the same constructor, and `order` and `theta` with the same +meaning (`order=1, theta=0.5` is Crank-Nicolson and the default, as for SLCN; +`order=2, theta=1.0` is BDF2, the counterpart of SL-BDF2; `order=2, theta=0.5` is +refused for the reason the SLCN documentation gives). `integrator` is inferred and +only needs setting to reach Adams-Moulton above order 1. The choice of +Crank-Nicolson as the default follows the drop-in contract and the table: it is +the more accurate scheme wherever the answer is good, and where it rings the +answer is already wrong for every scheme. A user who wants damping asks for +`order=2`; below Courant 1 with diffusion, `order=3`. Backward Euler is not a +sensible choice for transport. ### Temporal convergence (tests/test_1100) diff --git a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py index d5a0db03e..5987ac3bc 100644 --- a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py @@ -37,13 +37,14 @@ timestep is a runtime constant of the compiled kernels. - **SUPG stabilisation**: the streamline-upwind test-function perturbation written as a flux, so PETSc needs no modified test space. -- **Multistep order**: `order=1, 2, 3` with `integrator="bdf"` or `"am"`. +- **Drop-in for SLCN**: the same constructor, `order`, `theta`, `estimate_dt` + and `solve`; change the class name and nothing else. ## Parameters - `uw_res`: cells across the box - `uw_courant`: timestep as a multiple of the cell-crossing time -- `uw_order`, `uw_integrator`, `uw_theta`: the time scheme +- `uw_order`, `uw_theta`: the time scheme, with the semi-Lagrangian solver's meaning - `uw_diffusivity`: thermal diffusivity (0 is pure advection) """ @@ -67,9 +68,8 @@ params = uw.Params( uw_res=32, uw_courant=1.0, - uw_order=2, - uw_integrator="bdf", - uw_theta=1.0, + uw_order=1, # 1 with theta 0.5 is Crank-Nicolson; 2 with theta 1.0 is BDF2 + uw_theta=0.5, uw_diffusivity=0.0, uw_sigma=0.12, ) @@ -104,8 +104,7 @@ # %% adv_diff = uw.systems.AdvDiffusionSUPG( - mesh, T, velocity, order=params.uw_order, - integrator=params.uw_integrator, theta=params.uw_theta) + mesh, T, velocity, order=params.uw_order, theta=params.uw_theta) adv_diff.constitutive_model.Parameters.diffusivity = params.uw_diffusivity for boundary in ("Left", "Right", "Top", "Bottom"): adv_diff.add_dirichlet_bc(0.0, boundary) diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 6c8e23f41..924aae9aa 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -20,9 +20,11 @@ turns in one step. See ``docs/developer/design/eulerian-supg-transport.md``. """ +import warnings + import numpy as np import sympy -from typing import Optional +from typing import Callable, Optional, Union import underworld3 as uw import underworld3.timing as timing @@ -64,35 +66,75 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi - \nabla\cdot(\kappa\nabla\phi) = f - Two families of time integration are built from the same stored history - :math:`\phi^{n}, \phi^{n-1}, \dots` (real mesh variables, so their - gradients are available inside the kernels): + A drop-in replacement for :class:`~underworld3.systems.solvers.SNES_AdvectionDiffusion` + (``uw.systems.AdvDiffusionSLCN``): the constructor, ``order``, ``theta``, + ``f``, ``V_fn``, ``constitutive_model``, ``delta_t``, ``estimate_dt`` and + ``solve`` all keep the semi-Lagrangian solver's meaning, so a script changes + the class name and nothing else:: + + adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=1) # was AdvDiffusionSLCN + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = 1.0e-3 + adv.add_dirichlet_bc(0.0, "Left") + adv.solve(timestep=dt) + + The arguments that only make sense for a trace-back + (``restore_points_func``, ``monotone_mode``, ``old_frame_traceback``, + ``DFDt``) are accepted and ignored with a warning. + + **Time schemes.** ``order`` and ``theta`` select the same schemes as for + the semi-Lagrangian solver: + + ========== ======= ===================================================== + ``order`` ``theta`` scheme + ========== ======= ===================================================== + 1 0.5 Crank-Nicolson (default; the SLCN convention) + 1 1.0 backward Euler + 2 1.0 BDF2, all spatial terms at :math:`n+1` (the SL-BDF2 convention) + 3 1.0 BDF3 + ========== ======= ===================================================== - ``integrator="bdf"`` (backward differentiation, order 1-3) + ``order=2`` with ``theta=0.5`` is refused, as the semi-Lagrangian + documentation says: a BDF stencil pairs with terms at :math:`n+1`, not + with a centred flux. Every past time level is a mesh variable held by an + :class:`~underworld3.systems.ddt.Eulerian` history manager, so gradients + of past states are available in the kernels and both families come from + one code path: + + ``integrator="bdf"`` (order :math:`N`) .. math:: \frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k} + \mathbf{u}\cdot\nabla\phi^{n+1} - \nabla\cdot(\kappa\nabla\phi^{n+1}) = f - ``integrator="am"`` (Adams-Moulton, order 1-3; ``theta`` at order 1) + ``integrator="am"`` (order :math:`N`; ``theta`` at order 1) .. math:: \frac{\phi^{n+1}-\phi^{n}}{\Delta t} + \sum_{k=0}^{N} a_k\left[\mathbf{u}\cdot\nabla\phi^{n+1-k} - \nabla\cdot(\kappa\nabla\phi^{n+1-k})\right] = f - Order 1 with ``theta=1`` is backward Euler in both families; - ``integrator="am", order=1, theta=0.5`` is Crank-Nicolson. The BDF - coefficients :math:`c_k` and Adams-Moulton weights :math:`a_k` are the - ones the :class:`~underworld3.systems.ddt.Eulerian` manager maintains; - both ramp from first order over the opening steps unless a history is - planted with ``solver.DuDt.set_initial_history``. A BDF3 request falls - back to variable-step BDF2 whenever consecutive timesteps differ by more - than 5%. + ``integrator`` is inferred from ``order`` and ``theta`` (Adams-Moulton at + order 1, BDF above) and only needs setting to reach Adams-Moulton above + order 1, which is provided for diffusion-dominated problems: its bounded + stability region makes it blow up on an advection operator from about + Courant 1. Both families ramp from first order over the opening steps + unless a history is planted with ``solver.DuDt.set_initial_history``. A + BDF3 request falls back to variable-step BDF2 whenever consecutive + timesteps differ by more than 5%. + + **Which scheme.** Measured on a rotating Gaussian + (``docs/developer/design/eulerian-supg-transport.md``): Crank-Nicolson is + three to four times more accurate than BDF2 at the same timestep below + Courant 2 on the feature scale, and rings once the feature is + under-resolved in time; BDF2 is damped and stable at every Courant + number; BDF3 is the most accurate scheme below Courant 1 when diffusion + is present but grows slowly on pure advection; backward Euler carries 20 + to 40% error at any practical timestep. **Weak form.** With the strong residual of the chosen scheme - :math:`R(\phi)` (time derivative and advection; see below) the residual + :math:`R(\phi)` (time derivative, advection, source) the residual assembled through PETSc's pointwise interface is .. math:: @@ -117,20 +159,22 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): + \left(\frac{2|\mathbf{u}|}{h}\right)^2 + \left(\frac{4\kappa}{h^2}\right)^2\right]^{-1/2} - with :math:`h` the local cell size (``mesh.cell_size()``, the - :math:`\mathrm{volume}^{1/d}` equivalent radius) and :math:`c_0` the - leading multistep coefficient. The three weights are runtime constants - (``tau_weights``) and ``supg_weight`` scales the whole term, so a - Galerkin baseline needs no rebuild. + with :math:`h` the local cell size (``mesh.cell_size()``) and + :math:`c_0` the leading multistep coefficient. The three weights are + runtime constants (``tau_weights``) and ``supg_weight`` scales the whole + term, so a Galerkin baseline needs no rebuild. - **What limits the timestep.** Nothing, for stability. The implicit + **What limits the timestep.** Nothing, for stability: the implicit scheme is stable at any cell Courant number, including on cells refined for a Stokes problem that the scalar does not need. Accuracy is set by how far the transported feature moves per step relative to its own - width: the error grows as :math:`(\mathbf{u}\Delta t)^2` for the - second-order schemes, and the transient term of :math:`\tau` cannot - hide that. :meth:`estimate_dt` returns the cell-crossing time as a - resolution guide only. + width, as :math:`(\mathbf{u}\Delta t)^2` for the second-order schemes. + :meth:`estimate_dt` returns the cell-crossing time as a resolution guide, + exactly as the semi-Lagrangian solver does. Against that solver: the + semi-Lagrangian error is flat in the timestep but accumulates one + interpolation per step, and its limit is the arc a characteristic turns + per step; the Eulerian solve costs four to six times less per step in + serial and needs no departure points in parallel. Parameters ---------- @@ -139,42 +183,31 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): Continuous scalar field :math:`\phi`. V_fn : MeshVariable or sympy Matrix Advecting velocity, ``(1, dim)``. - order : int, default 2 - Order of the time integration, 1 to 3. BDF2 is the default because it - is stable at every Courant number and damped. Measured on a rotating - Gaussian (``docs/developer/design/eulerian-supg-transport.md``): - Crank-Nicolson is three to four times more accurate than BDF2 at the - same timestep below Courant 2 but rings once the feature is - under-resolved in time; BDF3 is the most accurate scheme below - Courant 1 when diffusion is present, fails from Courant 4, and on - pure advection grows slowly at any Courant number; backward Euler (order 1) carries - 20 to 40% error at any practical timestep. - integrator : {"bdf", "am"}, default "bdf" - Adams-Moulton above order 1 has a bounded stability region and blows - up on an advection operator from about Courant 1; it is provided for - diffusion-dominated problems and for comparison. - theta : float, default 1.0 - Adams-Moulton blend at order 1 only (0.5 is Crank-Nicolson). Must be - 1.0 for BDF and for Adams-Moulton above order 1. + order : int, default 1 + Time-integration order, 1 to 3 (see the table above). + theta : float, optional + Crank-Nicolson blend at order 1: 0.5 (the default there) is + Crank-Nicolson, 1.0 is backward Euler. Above order 1 the only + consistent value is 1.0, which is taken when ``theta`` is not given + and refused when 0.5 is asked for explicitly. + integrator : {"bdf", "am"}, optional + Inferred from ``order`` and ``theta`` when omitted. verbose : bool, default False DuDt : Eulerian, optional A pre-built history manager (order at least ``order``, no ``V_fn``). + restore_points_func, monotone_mode, old_frame_traceback, DFDt + Semi-Lagrangian arguments, accepted for drop-in compatibility and + ignored with a warning: there is no trace-back here. Notes ----- - The diffusivity is set through the constitutive model, as for the other - scalar solvers; the solver starts with a + The diffusivity is set through the constitutive model, as for every + scalar solver; the solver starts with a :class:`~underworld3.constitutive_models.DiffusionModel` at - :math:`\kappa = 0` (pure advection):: - - adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=2) - adv.constitutive_model.Parameters.diffusivity = 1.0e-3 - adv.add_dirichlet_bc(0.0, "Left") - adv.solve(timestep=dt) - - The linear system is nonsymmetric, so the preconditioner defaults are - GMRES with an additive-Schwarz ILU preconditioner rather than the - algebraic multigrid the symmetric scalar solvers use. + :math:`\kappa = 0` (pure advection). The linear system is nonsymmetric, + so the solver defaults to GMRES with an additive-Schwarz ILU + preconditioner rather than the algebraic multigrid the symmetric scalar + solvers use; every option is overridable through ``petsc_options``. """ _INTEGRATORS = ("bdf", "am") @@ -185,29 +218,52 @@ def __init__( mesh: uw.discretisation.Mesh, u_Field: uw.discretisation.MeshVariable, V_fn, - order: int = 2, - integrator: str = "bdf", - theta: float = 1.0, + order: int = 1, + theta: Optional[float] = None, + integrator: Optional[str] = None, verbose: bool = False, DuDt: Optional[Eulerian_DDt] = None, + DFDt=None, + restore_points_func: Optional[Callable] = None, + monotone_mode: Optional[str] = None, + old_frame_traceback: bool = False, ): if not u_Field.continuous: raise ValueError( "u_Field must be a continuous MeshVariable: the SUPG weak form " "is continuous Galerkin." ) - if integrator not in self._INTEGRATORS: - raise ValueError( - f"integrator must be one of {self._INTEGRATORS}, not {integrator!r}." + ignored = [name for name, value in ( + ("restore_points_func", restore_points_func), + ("monotone_mode", monotone_mode), + ("old_frame_traceback", old_frame_traceback), + ("DFDt", DFDt), + ) if value] + if ignored: + warnings.warn( + f"AdvDiffusionSUPG ignores {', '.join(ignored)}: these configure " + "the semi-Lagrangian trace-back and the Eulerian scheme has none.", + stacklevel=2, ) order = int(order) if order not in (1, 2, 3): raise ValueError(f"order must be 1, 2 or 3, not {order}.") - theta = float(theta) + # theta means what it means for the semi-Lagrangian solver: the + # Crank-Nicolson blend at order 1. Left unset, order 2 and 3 take the + # only consistent value; set explicitly to 0.5 there, it is refused. + theta = float(theta) if theta is not None else (0.5 if order == 1 else 1.0) + if integrator is None: + integrator = "am" if order == 1 else "bdf" + if integrator not in self._INTEGRATORS: + raise ValueError( + f"integrator must be one of {self._INTEGRATORS}, not {integrator!r}." + ) if theta != 1.0 and not (integrator == "am" and order == 1): raise ValueError( - "theta applies to integrator='am' at order 1 only " - "(0.5 is Crank-Nicolson); higher orders and BDF take theta=1." + "theta applies at order 1 only (0.5 is Crank-Nicolson, 1.0 is " + "backward Euler); order 2 and 3 take theta=1.0, the same rule as " + "the semi-Lagrangian solver (a BDF stencil pairs with terms at n+1, " + "not with a centred flux)." ) super().__init__(mesh, u_Field, u_Field.degree, verbose, DuDt=DuDt, DFDt=None) @@ -277,6 +333,19 @@ def __init__( self.petsc_options["snes_rtol"] = 1.0e-8 self.petsc_options["snes_max_it"] = 20 + def _object_viewer(self): + from IPython.display import Latex, display + + super()._object_viewer() + scheme = {("am", 1): f"Adams-Moulton order 1, theta = {self._theta}", + ("bdf", 1): "backward Euler"}.get( + (self._integrator, self._time_order), + f"{self._integrator.upper()} order {self._time_order}") + display(Latex(r"$\quad\mathrm{u} = $ " + self.u.sym._repr_latex_())) + display(Latex(r"$\quad\mathbf{v} = $ " + self._V_fn._repr_latex_())) + display(Latex(r"$\quad\Delta t = $ " + self._delta_t._repr_latex_())) + display(Latex(rf"$\quad$ time scheme: {scheme}")) + # ------------------------------------------------------------------ # Scheme description # ------------------------------------------------------------------ @@ -298,9 +367,24 @@ def theta(self) -> float: @property def delta_t(self): - r"""The timestep :math:`\Delta t` as a UW expression (set by :meth:`solve`).""" + r"""The timestep :math:`\Delta t` as a UW expression. + + Set by :meth:`solve`, or assign it directly (a number or a quantity + with time units) and call ``solve()`` without ``timestep``, as with + the semi-Lagrangian solver. A new value updates a runtime constant of + the compiled kernels; nothing is recompiled. + """ return self._delta_t + @delta_t.setter + def delta_t(self, value): + dt = float(_nondimensionalise_timestep(value)) + if dt <= 0.0: + raise ValueError(f"timestep must be positive, not {dt}.") + if dt != self._last_timestep: + self._delta_t.sym = dt + self._last_timestep = dt + @property def V_fn(self): """Advecting velocity, ``(1, dim)``.""" @@ -458,29 +542,27 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): def solve( self, - *, - timestep=None, zero_init_guess: Optional[bool] = None, - verbose: bool = False, + timestep=None, _force_setup: bool = False, + _evalf: bool = False, + verbose: bool = False, divergence_retries: int = 0, ): - r"""Advance :math:`\phi` by one step of size ``timestep``. + r"""Advance :math:`\phi` by one step. - ``timestep`` is required and keyword-only. Changing it between calls - updates a runtime constant of the compiled kernels; nothing is - recompiled. + Same signature as the semi-Lagrangian solver. ``timestep`` sets + :attr:`delta_t`; omit it to reuse the value already set. Changing it + between calls updates a runtime constant of the compiled kernels; + nothing is recompiled. """ - if timestep is None: + if timestep is not None: + self.delta_t = timestep + elif self._last_timestep is None: raise ValueError( - "solve() requires timestep=
; there is no default timestep." + "solve() needs a timestep: pass timestep=
or set solver.delta_t first." ) - dt = float(_nondimensionalise_timestep(timestep)) - if dt <= 0.0: - raise ValueError(f"timestep must be positive, not {dt}.") - if dt != self._last_timestep: - self._delta_t.sym = dt - self._last_timestep = dt + dt = self._last_timestep if _force_setup: self._needs_function_rewire = True diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index bb95e3fd4..d1e045f10 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -32,14 +32,38 @@ def _solver(mesh, tag, **kwargs): return adv, T -def test_exported_and_constructs(mesh): +def test_exported_and_constructs_with_the_slcn_defaults(mesh): adv, _T = _solver(mesh, "a") assert type(adv).__name__ == "SNES_AdvectionDiffusion_SUPG" - assert adv.integrator == "bdf" and adv.order == 2 + # order 1, theta 0.5: Crank-Nicolson, the semi-Lagrangian solver's default + assert adv.integrator == "am" and adv.order == 1 and adv.theta == 0.5 assert isinstance(adv.DuDt, uw.systems.ddt.Eulerian) assert adv.DuDt.V_fn is None, "advection is implicit, not a history correction" +def test_slcn_order_theta_pairs_select_the_documented_schemes(mesh): + assert _solver(mesh, "p1", order=1, theta=1.0)[0].integrator == "am" # backward Euler + assert _solver(mesh, "p2", order=2, theta=1.0)[0].integrator == "bdf" # SL-BDF2's counterpart + assert _solver(mesh, "p3", order=2)[0].integrator == "bdf" # theta 0.5 only bites at order 1 + with pytest.raises(ValueError, match="theta applies"): + _solver(mesh, "p4", order=2, theta=0.5) + + +def test_semi_lagrangian_only_arguments_are_ignored_with_a_warning(mesh): + with pytest.warns(UserWarning, match="monotone_mode, old_frame_traceback"): + adv, _T = _solver(mesh, "q", monotone_mode="clamp", old_frame_traceback=True) + adv.solve(timestep=0.01) + + +def test_solve_takes_the_slcn_signature_and_delta_t(mesh): + adv, T = _solver(mesh, "s") + adv.solve(False, 0.01) # positional, as SLCN allows + adv.delta_t = 0.02 # set once ... + adv.solve() # ... and reuse + assert float(adv.delta_t.sym) == 0.02 + assert np.isfinite(np.asarray(T.array)).all() + + @pytest.mark.parametrize("tag, kwargs, message", [ ("v0", dict(order=4), "order must be"), ("v1", dict(integrator="rk4"), "integrator must be"), @@ -53,14 +77,14 @@ def test_scheme_arguments_are_validated(mesh, tag, kwargs, message): def test_timestep_is_required(mesh): adv, _T = _solver(mesh, "b") - with pytest.raises(ValueError, match="requires timestep"): + with pytest.raises(ValueError, match="needs a timestep"): adv.solve() def test_bdf1_diffusive_flux_is_the_constitutive_flux(mesh): """At order 1 the assembled diffusive flux is exactly the constitutive model's own flux of the new state; the history weights are inert.""" - adv, _T = _solver(mesh, "c") + adv, _T = _solver(mesh, "c", order=1, theta=1.0, integrator="bdf") adv.constitutive_model.Parameters.diffusivity = 0.7 difference = adv._diffusive_flux() - adv.constitutive_model.flux.T assert all(sympy.simplify(e) == 0 for e in difference) From 252165e791a4456fda158121a4beb7c6a4e33d86 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 21:48:05 -0700 Subject: [PATCH 10/14] Drop the integrator argument: order and theta already reach every safe scheme The only schemes the argument added were Adams-Moulton at orders 2 and 3, which the integrator study shows blowing up on advection from Courant 1. The multistep family now follows the order (the theta rule at order 1, BDF above); the higher Adams-Moulton assembly stays in the code, reachable only by switching the family on the instance, which is how the study measured it. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- docs/advanced/eulerian-advection-diffusion.md | 6 +-- .../semi-lagrangian-time-integration.md | 6 +-- .../design/eulerian-supg-transport.md | 16 ++++---- .../systems/advection_diffusion_eulerian.py | 39 ++++++++----------- tests/test_1055_advdiff_supg_api.py | 22 ++++++----- ...est_1100_advdiff_supg_rotating_gaussian.py | 6 +-- 6 files changed, 46 insertions(+), 49 deletions(-) diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md index 9e0951526..a18f18925 100644 --- a/docs/advanced/eulerian-advection-diffusion.md +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -37,9 +37,7 @@ adv.solve(timestep=dt) | `DuDt.set_initial_history(values, dt)` | same | plant an exact history to start at full order | | `restore_points_func`, `monotone_mode`, `old_frame_traceback`, `DFDt` | ignored, with a warning | they configure the trace-back | -`order=3` (BDF3) is available; see below for when it is safe. `integrator="am"` -above order 1 reaches the higher Adams-Moulton rules, which are for -diffusion-dominated problems only. +`order=3` (BDF3) is available; see below for when it is safe. ## When to use which @@ -80,7 +78,7 @@ tables are in the design note. | BDF2 (`order=2`) | damped and stable at every Courant number; the choice for sharp or under-resolved fields | | BDF3 (`order=3`) | the most accurate scheme below Courant 1 when diffusion is present; on pure advection it grows slowly at any Courant number, so use it only with diffusion | | backward Euler (`order=1, theta=1.0`) | 20 to 40% error at any practical timestep; not for transport | -| Adams-Moulton 2, 3 (`integrator="am"`) | third and fourth order below Courant 1; blow up on advection from about Courant 1 | +| Adams-Moulton 2, 3 (not offered) | third and fourth order below Courant 1 but blow up on advection from about Courant 1, which is why there is no knob for them | All schemes cost the same per step: the history terms are extra kernel inputs, not extra solves. Changing the timestep between steps changes a runtime constant diff --git a/docs/advanced/semi-lagrangian-time-integration.md b/docs/advanced/semi-lagrangian-time-integration.md index 7812b7c78..f48c36ed7 100644 --- a/docs/advanced/semi-lagrangian-time-integration.md +++ b/docs/advanced/semi-lagrangian-time-integration.md @@ -116,9 +116,9 @@ $[\theta,\,1-\theta]$: `uw.systems.AdvDiffusionSUPG` solves the same equation without a trace-back: all terms are assembled on the mesh, implicit in time, with SUPG -stabilisation. Its `order=` and `integrator="bdf"|"am"` arguments select the -multistep scheme, built from the same stored history as above; `order=1, -integrator="am", theta=0.5` is Crank-Nicolson. The scheme is stable at any +stabilisation. Its `order=` and `theta=` arguments mean what they mean here: +`order=1, theta=0.5` is Crank-Nicolson, `order=2` is BDF2, built from the same +stored history as above. The scheme is stable at any cell Courant number, so cells refined for a Stokes problem never limit the transport timestep; its accuracy is set by how far the transported feature moves per step. The semi-Lagrangian scheme's accuracy is instead set by how diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 580870ebf..8b98803c5 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -30,10 +30,10 @@ Every past time level $\phi^{n}, \phi^{n-1}, \dots$ is a mesh variable held by a `Eulerian` history manager, so first derivatives of past states are available in the kernels and two multistep families share one code path: -| `integrator` | time derivative | spatial operator | +| family | time derivative | spatial operator | |---|---|---| -| `bdf`, order $N$ | $\frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k}$ | at $n+1$ only | -| `am`, order $N$ | $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$ | $\sum_{k=0}^{N} a_k\,S(\phi^{n+1-k})$ | +| BDF, order $N$ (`order=2, 3`) | $\frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k}$ | at $n+1$ only | +| theta rule (`order=1`; Adams-Moulton of $N$ steps internally) | $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$ | $\sum_{k=0}^{N} a_k\,S(\phi^{n+1-k})$ | with $S(\phi) = \mathbf{u}\cdot\nabla\phi - \nabla\cdot(\kappa\nabla\phi)$ and the coefficients those the history manager already maintains (`theta` is the @@ -175,8 +175,8 @@ What the table says: - **Adams-Moulton above order 1 is unusable for advection.** Its stability region is bounded and covers only a short segment of the imaginary axis, so on a pure advection operator it blows up once the Courant number reaches about 1, and - diffusion at this Peclet number does not rescue it. It is kept in the class for - the record and for diffusion-dominated use, with that warning in the docstring. + diffusion at this Peclet number does not rescue it. The assembly code handles + it, but no public argument reaches it. - **BDF3 is the most accurate scheme below Courant 1 with diffusion present** (0.3%, on the spatial floor) but it is not A-stable, fails from Courant 4, and on pure advection grows slowly at any Courant number (the res-64 rows). @@ -192,8 +192,10 @@ What the table says: semi-Lagrangian solver: the same constructor, and `order` and `theta` with the same meaning (`order=1, theta=0.5` is Crank-Nicolson and the default, as for SLCN; `order=2, theta=1.0` is BDF2, the counterpart of SL-BDF2; `order=2, theta=0.5` is -refused for the reason the SLCN documentation gives). `integrator` is inferred and -only needs setting to reach Adams-Moulton above order 1. The choice of +refused for the reason the SLCN documentation gives). There is no `integrator` +argument: the family follows the order, and the only schemes that argument would +have added, Adams-Moulton at orders 2 and 3, are the ones the table rules out. +The study reached them by switching the family on the instance. The choice of Crank-Nicolson as the default follows the drop-in contract and the table: it is the more accurate scheme wherever the answer is good, and where it rings the answer is already wrong for every scheme. A user who wants damping asks for diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 924aae9aa..476a0b5db 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -101,28 +101,26 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): of past states are available in the kernels and both families come from one code path: - ``integrator="bdf"`` (order :math:`N`) + backward differentiation (order :math:`N \ge 2`) .. math:: \frac{1}{\Delta t}\sum_{k=0}^{N} c_k\,\phi^{n+1-k} + \mathbf{u}\cdot\nabla\phi^{n+1} - \nabla\cdot(\kappa\nabla\phi^{n+1}) = f - ``integrator="am"`` (order :math:`N`; ``theta`` at order 1) + the :math:`\theta` rule (order 1; Adams-Moulton of one step) .. math:: \frac{\phi^{n+1}-\phi^{n}}{\Delta t} + \sum_{k=0}^{N} a_k\left[\mathbf{u}\cdot\nabla\phi^{n+1-k} - \nabla\cdot(\kappa\nabla\phi^{n+1-k})\right] = f - ``integrator`` is inferred from ``order`` and ``theta`` (Adams-Moulton at - order 1, BDF above) and only needs setting to reach Adams-Moulton above - order 1, which is provided for diffusion-dominated problems: its bounded - stability region makes it blow up on an advection operator from about - Courant 1. Both families ramp from first order over the opening steps - unless a history is planted with ``solver.DuDt.set_initial_history``. A - BDF3 request falls back to variable-step BDF2 whenever consecutive - timesteps differ by more than 5%. + The higher Adams-Moulton rules are assembled by the same code but are + not offered: their bounded stability region blows up on an advection + operator from about Courant 1 (see the design note). Both families ramp + from first order over the opening steps unless a history is planted with + ``solver.DuDt.set_initial_history``. A BDF3 request falls back to + variable-step BDF2 whenever consecutive timesteps differ by more than 5%. **Which scheme.** Measured on a rotating Gaussian (``docs/developer/design/eulerian-supg-transport.md``): Crank-Nicolson is @@ -190,8 +188,6 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): Crank-Nicolson, 1.0 is backward Euler. Above order 1 the only consistent value is 1.0, which is taken when ``theta`` is not given and refused when 0.5 is asked for explicitly. - integrator : {"bdf", "am"}, optional - Inferred from ``order`` and ``theta`` when omitted. verbose : bool, default False DuDt : Eulerian, optional A pre-built history manager (order at least ``order``, no ``V_fn``). @@ -210,8 +206,6 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): solvers use; every option is overridable through ``petsc_options``. """ - _INTEGRATORS = ("bdf", "am") - @timing.routine_timer_decorator def __init__( self, @@ -220,7 +214,6 @@ def __init__( V_fn, order: int = 1, theta: Optional[float] = None, - integrator: Optional[str] = None, verbose: bool = False, DuDt: Optional[Eulerian_DDt] = None, DFDt=None, @@ -252,13 +245,13 @@ def __init__( # Crank-Nicolson blend at order 1. Left unset, order 2 and 3 take the # only consistent value; set explicitly to 0.5 there, it is refused. theta = float(theta) if theta is not None else (0.5 if order == 1 else 1.0) - if integrator is None: - integrator = "am" if order == 1 else "bdf" - if integrator not in self._INTEGRATORS: - raise ValueError( - f"integrator must be one of {self._INTEGRATORS}, not {integrator!r}." - ) - if theta != 1.0 and not (integrator == "am" and order == 1): + # The multistep family follows the order: the Adams-Moulton (theta) + # rule at order 1, backward differentiation above. Adams-Moulton at + # orders 2 and 3 is assembled by the same code but is not offered: + # its bounded stability region blows up on an advection operator + # from about Courant 1 (design note, integrator study). + integrator = "am" if order == 1 else "bdf" + if theta != 1.0 and order != 1: raise ValueError( "theta applies at order 1 only (0.5 is Crank-Nicolson, 1.0 is " "backward Euler); order 2 and 3 take theta=1.0, the same rule as " @@ -352,7 +345,7 @@ def _object_viewer(self): @property def integrator(self) -> str: - """``"bdf"`` or ``"am"``.""" + """The multistep family in use: ``"am"`` (the theta rule) at order 1, ``"bdf"`` above.""" return self._integrator @property diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index d1e045f10..590dc2b61 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -66,9 +66,9 @@ def test_solve_takes_the_slcn_signature_and_delta_t(mesh): @pytest.mark.parametrize("tag, kwargs, message", [ ("v0", dict(order=4), "order must be"), - ("v1", dict(integrator="rk4"), "integrator must be"), - ("v2", dict(integrator="bdf", theta=0.5), "theta applies"), - ("v3", dict(integrator="am", order=2, theta=0.5), "theta applies"), + ("v1", dict(order=0), "order must be"), + ("v2", dict(order=2, theta=0.5), "theta applies"), + ("v3", dict(order=3, theta=0.5), "theta applies"), ]) def test_scheme_arguments_are_validated(mesh, tag, kwargs, message): with pytest.raises(ValueError, match=message): @@ -81,17 +81,21 @@ def test_timestep_is_required(mesh): adv.solve() -def test_bdf1_diffusive_flux_is_the_constitutive_flux(mesh): - """At order 1 the assembled diffusive flux is exactly the constitutive - model's own flux of the new state; the history weights are inert.""" - adv, _T = _solver(mesh, "c", order=1, theta=1.0, integrator="bdf") +def test_bdf_diffusive_flux_is_the_constitutive_flux(mesh): + """For the BDF family the assembled diffusive flux is exactly the + constitutive model's own flux of the new state; no history enters it.""" + adv, _T = _solver(mesh, "c", order=2) adv.constitutive_model.Parameters.diffusivity = 0.7 difference = adv._diffusive_flux() - adv.constitutive_model.flux.T assert all(sympy.simplify(e) == 0 for e in difference) -def test_am_order2_uses_all_three_time_levels(mesh): - adv, _T = _solver(mesh, "d", integrator="am", order=2) +def test_multistep_weights_reach_every_stored_time_level(mesh): + # The theta rule at higher order is assembled by the same code; it is not + # offered publicly (unstable for advection), so the family is switched + # on the instance here to cover the weighted-sum path. + adv, _T = _solver(mesh, "d", order=2) + adv._integrator = "am" weights = adv._spatial_weights() assert len(weights) == 3 states = adv._states() diff --git a/tests/test_1100_advdiff_supg_rotating_gaussian.py b/tests/test_1100_advdiff_supg_rotating_gaussian.py index 5607d22a6..336ea74e4 100644 --- a/tests/test_1100_advdiff_supg_rotating_gaussian.py +++ b/tests/test_1100_advdiff_supg_rotating_gaussian.py @@ -32,14 +32,14 @@ def _box(res, refinement=0): qdegree=3, regular=False, refinement=refinement) -def _problem(mesh, tag, order, integrator="bdf", theta=1.0, kappa=0.0): +def _problem(mesh, tag, order, theta=None, kappa=0.0): x, y = mesh.X sol = uw.analytic.RotatingGaussian(mesh, sigma=SIGMA, centre_radius=0.5, omega=1.0, diffusivity=kappa) T = uw.discretisation.MeshVariable(f"T_{tag}", mesh, 1, degree=2) T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), - order=order, integrator=integrator, theta=theta) + order=order, theta=theta) adv.constitutive_model.Parameters.diffusivity = kappa for b in ("Left", "Right", "Top", "Bottom"): adv.add_dirichlet_bc(0.0, b) @@ -74,7 +74,7 @@ def test_temporal_convergence_order(order, timesteps, expected_slope): t_end = float(sympy.pi) / 2 errors = [] for i, dt in enumerate(timesteps): - sol, T, adv = _problem(mesh, f"c{order}{i}", order) + sol, T, adv = _problem(mesh, f"c{order}{i}", order, theta=1.0) errors.append(_run(sol, T, adv, dt, t_end)) slopes = np.log2(np.array(errors[:-1]) / np.array(errors[1:])) print(f"order {order}: errors {errors} slopes {slopes}") From 417f7b88504fb7075e9f670478a7c41dbf8aba68 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 2 Sep 2026 22:05:51 -0700 Subject: [PATCH 11/14] An accuracy-based timestep for the Eulerian solver; credit NengLu in the module and note estimate_dt now returns the step at which the field changes by a fraction (0.02) of its range: from the advective rate |u . grad phi| at the vertices before the first solve, and from the rate the last step actually produced after it. The cell-crossing time the semi-Lagrangian solver reports is not a stability limit for this scheme and says nothing about its accuracy; it stays available as basis='resolution'. The estimate is mesh-independent, which the band test now checks (the resolution estimate collapses 3x on the refined child, the accuracy estimate moves under 25%), and at the default fraction Crank-Nicolson completes the rotating-Gaussian round trip under one per cent. The advective rate uses the vertex Clement gradient rather than a point evaluation of a derivative expression, which fails on a mesh carrying many variables. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- docs/advanced/eulerian-advection-diffusion.md | 14 +- .../design/eulerian-supg-transport.md | 34 +++-- .../Ex_AdvectionDiffusionSUPG_RotationTest.py | 34 ++--- .../systems/advection_diffusion_eulerian.py | 127 +++++++++++++++--- tests/test_1055_advdiff_supg_api.py | 27 ++++ ...est_1100_advdiff_supg_rotating_gaussian.py | 23 +++- 6 files changed, 205 insertions(+), 54 deletions(-) diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md index a18f18925..d62a6d416 100644 --- a/docs/advanced/eulerian-advection-diffusion.md +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -19,10 +19,20 @@ adv.constitutive_model.Parameters.diffusivity = 1.0e-3 adv.add_dirichlet_bc(1.0, "Bottom") adv.add_dirichlet_bc(0.0, "Top") -dt = 0.5 * adv.estimate_dt() +dt = adv.estimate_dt() # accuracy-based: 2% of the field's range per step adv.solve(timestep=dt) ``` +The one deliberate difference is the timestep estimate. The semi-Lagrangian +`estimate_dt` reports the cell-crossing time, which for this solver is neither +a stability limit nor an accuracy one. The Eulerian solver's `estimate_dt` +instead returns the step at which the field changes by a given fraction of its +range (0.02 by default), from the advective rate before the first solve and +from the rate the last step actually produced after it. It does not depend on +the mesh, so cells refined for the Stokes problem do not shrink it. A script +that sizes its step in Courant numbers can still ask for +`estimate_dt(basis="resolution")`. + ## What carries over | SLCN | SUPG | note | @@ -32,7 +42,7 @@ adv.solve(timestep=dt) | `order=2, theta=1.0` | same | SL-BDF2 becomes BDF2 | | `order=2, theta=0.5` | refused | refused for the same reason: a BDF stencil does not pair with a centred flux | | `f`, `V_fn`, `constitutive_model`, `delta_t` | same | | -| `estimate_dt(direction_aware, percentile)` | same | the cell-crossing time, a resolution guide for both | +| `estimate_dt()` | accuracy-based by default | the field may change by `fraction` (0.02) of its range per step; `basis="resolution"` returns the cell-crossing time SLCN reports | | `solve(zero_init_guess, timestep, ...)` | same | | | `DuDt.set_initial_history(values, dt)` | same | plant an exact history to start at full order | | `restore_points_func`, `monotone_mode`, `old_frame_traceback`, `DFDt` | ignored, with a warning | they configure the trace-back | diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 8b98803c5..ff894b4f2 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -1,8 +1,14 @@ # Eulerian SUPG transport: design and measurements **Status**: implemented on `feature/eulerian-supg-transport` (2026-09-02), static mesh. -Supersedes the Crank-Nicolson prototype of issue #657 as the implementation route -while keeping its weak-form idea. + +**Credit.** The SUPG weak form used here (the test-function perturbation written as +a flux, so PETSc needs no modified test space), its first working implementation on +PetscDS with P2 elements, the LeVeque swirling-flow comparison against SLCN and the +conservative level-set pipeline that motivated it are NengLu's, on the `levelset` +branch of issue #657. This note builds on that prototype: same formulation and +stabilisation parameter, time integration moved onto the symbolic history +machinery, and the measurements added. ## Why an Eulerian scheme @@ -210,12 +216,24 @@ BDF1 slopes 0.80 and 0.88 between $\Delta t$ = 0.02, 0.01, 0.005; BDF2 slopes ab ## What the timestep estimate means -`estimate_dt` returns the cell-crossing time, the same resolution estimate the -semi-Lagrangian solver reports, because that is the only quantity the mesh knows. -It is not a stability limit for either scheme. Choose the Eulerian timestep from -the transported feature: $|\mathbf{u}|\Delta t$ a fraction of its width. For SLCN -the honest limit is the trace-back arc, $\Delta t \lesssim 0.25 / \max|\nabla\mathbf{u}|$, -which is a separate change to that solver. +The cell-crossing time is not a stability limit for either scheme and says +nothing about this one's accuracy, so the Eulerian solver's `estimate_dt` measures +the field instead: + +$$ +\Delta t = f\,\frac{\max\phi - \min\phi}{\max|\dot\phi|}, +$$ + +with $\dot\phi$ the advective rate $|\mathbf{u}\cdot\nabla\phi|$ before the first +solve and the realised rate $|\phi^{n+1}-\phi^{n}|/\Delta t$ after it (diffusion +and sources included). On the rotating Gaussian the fraction at Courant 0.5 on +the res-32 mesh is about 0.03 (0.6% Crank-Nicolson error) and at Courant 1 about +0.07 (2.5%); the default $f = 0.02$ therefore sits at a few tenths of a per cent. +The estimate is mesh-independent by construction, which is the property the +transport note's section 1 asks for; `basis="resolution"` still returns the +semi-Lagrangian solver's cell-crossing time. For SLCN the honest limit is the +trace-back arc, $\Delta t \lesssim 0.25 / \max|\nabla\mathbf{u}|$, which is a separate +change to that solver. ## A defect found on the way diff --git a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py index 5987ac3bc..61da93c87 100644 --- a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py @@ -27,9 +27,10 @@ error is measured directly rather than inferred from a picture. The scheme is stable at any cell Courant number; what limits the timestep -is how far the anomaly moves per step relative to its own width. Try -`-uw_courant 4` to see the accuracy fall off as `dt**2` while the solve -stays perfectly stable, and `-uw_order 2` to see the second-order scheme. +is how far the anomaly moves per step relative to its own width, which is +what the solver's own `estimate_dt` measures. Try `-uw_dt_fraction 0.1` to +see the accuracy fall off as `dt**2` while the solve stays perfectly +stable, and `-uw_order 2` for the damped second-order scheme. ## Key Concepts @@ -43,7 +44,7 @@ ## Parameters - `uw_res`: cells across the box -- `uw_courant`: timestep as a multiple of the cell-crossing time +- `uw_dt_fraction`: allowed change of the field per step (the timestep follows) - `uw_order`, `uw_theta`: the time scheme, with the semi-Lagrangian solver's meaning - `uw_diffusivity`: thermal diffusivity (0 is pure advection) """ @@ -60,14 +61,14 @@ Override from the command line: ```bash -python Ex_AdvectionDiffusionSUPG_RotationTest.py -uw_courant 4 -uw_order 2 +python Ex_AdvectionDiffusionSUPG_RotationTest.py -uw_dt_fraction 0.1 -uw_order 2 ``` """ # %% params = uw.Params( uw_res=32, - uw_courant=1.0, + uw_dt_fraction=0.02, # allowed change of T per step, as a fraction of its range uw_order=1, # 1 with theta 0.5 is Crank-Nicolson; 2 with theta 1.0 is BDF2 uw_theta=0.5, uw_diffusivity=0.0, @@ -113,16 +114,19 @@ """ ## Time loop -`estimate_dt` returns the cell-crossing time. It is a resolution guide, not a -stability limit, so the timestep is a chosen multiple of it. For a multistep -scheme the exact history is planted so the first step already runs at full -order. +`estimate_dt` returns an accuracy-based step: the field may change by +`uw_dt_fraction` of its range per step. It does not depend on the mesh; the +cell-crossing time the semi-Lagrangian solver reports is available with +`basis="resolution"` and is printed for comparison. For a multistep scheme the +exact history is planted so the first step already runs at full order. """ # %% period = float(exact.period) -dt_cell = float(adv_diff.estimate_dt()) -n_steps = int(np.ceil(period / (params.uw_courant * dt_cell))) +dt_accuracy = float(adv_diff.estimate_dt(fraction=params.uw_dt_fraction)) +dt_cell = float(adv_diff.estimate_dt(basis="resolution")) +uw.pprint(f"accuracy-based dt {dt_accuracy:.4g}, cell-crossing dt {dt_cell:.4g}") +n_steps = int(np.ceil(period / dt_accuracy)) dt = period / n_steps if params.uw_order > 1: @@ -142,9 +146,9 @@ """ ## Result -After one revolution the field should match its initial state. At a Courant -number of one half the round-trip error is below one per cent on this mesh; -it grows as `dt**2` from there. +After one revolution the field should match its initial state. At the +default fraction the round-trip error is a few tenths of a per cent on this +mesh; it grows as `dt**2` with the fraction. """ # %% diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 476a0b5db..a6c3d36dc 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -18,6 +18,11 @@ and its accuracy is set by how far the transported feature moves in one step; the semi-Lagrangian scheme's accuracy is set by how far a characteristic turns in one step. See ``docs/developer/design/eulerian-supg-transport.md``. + +The SUPG weak form, the Petrov-Galerkin test-function perturbation written +as a flux so that PETSc needs no modified test space, and its first +implementation on PetscDS are NengLu's (issue #657, branch ``levelset``); +this module keeps that formulation and its stabilisation parameter. """ import warnings @@ -167,8 +172,10 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): for a Stokes problem that the scalar does not need. Accuracy is set by how far the transported feature moves per step relative to its own width, as :math:`(\mathbf{u}\Delta t)^2` for the second-order schemes. - :meth:`estimate_dt` returns the cell-crossing time as a resolution guide, - exactly as the semi-Lagrangian solver does. Against that solver: the + :meth:`estimate_dt` therefore returns an accuracy-based step, the + allowed change of the field per step as a fraction of its range, and + only reports the cell-crossing time on request + (``basis="resolution"``). Against the semi-Lagrangian solver: the semi-Lagrangian error is flat in the timestep but accumulates one interpolation per step, and its limit is the arc a characteristic turns per step; the Eulerian solve costs four to six times less per step in @@ -271,6 +278,7 @@ def __init__( self._delta_t = public_expression( rf"\Delta t_{{{tag}}}", 1.0, "Eulerian advection-diffusion timestep") self._last_timestep = None + self._last_change_rate = None # SUPG on/off and the three tau weights are runtime constants: the # compiled kernels read them from PETSc's constants[] array. @@ -507,31 +515,99 @@ def _tau(self): # ------------------------------------------------------------------ @timing.routine_timer_decorator - def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): - r"""Cell-crossing timestep, as a resolution guide. - - The minimum over cells of :math:`h/|\mathbf{u}|` and - :math:`h^2/\kappa`, exactly as for the semi-Lagrangian solver. It is - not a stability limit for this scheme, and on a mesh refined for - another problem it is far smaller than the timestep the transported - field needs. Choose the timestep from the feature being transported: - :math:`|\mathbf{u}|\Delta t` should be a fraction of its width. + def estimate_dt(self, fraction: float = 0.02, basis: str = "accuracy", + direction_aware: bool = False, percentile: float = 0.0): + r"""A timestep for this scheme, chosen for accuracy. + + The implicit scheme has no stability limit, so the cell-crossing time + the semi-Lagrangian solver reports says nothing about how large a step + this solver can take. What bounds the error is how much the field + changes per step, and that is what the default estimate measures: + + .. math:: + \Delta t = f\,\frac{\max\phi - \min\phi} + {\max\left|\dot\phi\right|} + + with :math:`\dot\phi` the rate of change of the field. Before the + first solve that rate is the advective one, :math:`|\mathbf{u}\cdot + \nabla\phi|` at the mesh vertices; after a solve it is the rate the + last step actually produced, :math:`|\phi^{n+1}-\phi^{n}|/\Delta t`, + which includes diffusion and sources. The estimate is independent of + the mesh, so a band of cells refined for another problem does not + shrink it; it does shrink for a feature that is genuinely + under-resolved, which is the honest answer. + + On the rotating Gaussian (``docs/developer/design/eulerian-supg-transport.md``) + ``fraction=0.02`` gives Crank-Nicolson a round-trip error of a few + tenths of a per cent after one revolution and BDF2 about 1.5%; + ``fraction=0.07`` gives 2.5% and 9%. Parameters ---------- - direction_aware : bool, default False - Use the per-cell extent along the local velocity. - percentile : float, default 0.0 - Global percentile of the per-cell timesteps instead of the minimum. + fraction : float, default 0.02 + Allowed change of the field per step as a fraction of its range. + basis : {"accuracy", "resolution"} + ``"resolution"`` returns the cell-crossing / diffusion time the + semi-Lagrangian solver's ``estimate_dt`` returns, for scripts that + size the step in Courant numbers. + direction_aware, percentile + Forwarded to the resolution estimate; ignored otherwise. + + Returns + ------- + pint.Quantity or float + With physical time units if a model with reference scales is + active, otherwise nondimensional. ``inf`` if nothing changes. """ - dt_estimate, dt_adv, dt_diff = _advective_diffusive_dt( - self.constitutive_model.K, self._V_fn, self.mesh, - direction_aware=direction_aware, percentile=percentile) - self.dt_adv = dt_adv if not np.isinf(dt_adv) else 0.0 - self.dt_diff = dt_diff if not np.isinf(dt_diff) else 0.0 - if np.isinf(dt_estimate): + from mpi4py import MPI + + if basis == "resolution": + dt_estimate, dt_adv, dt_diff = _advective_diffusive_dt( + self.constitutive_model.K, self._V_fn, self.mesh, + direction_aware=direction_aware, percentile=percentile) + self.dt_adv = dt_adv if not np.isinf(dt_adv) else 0.0 + self.dt_diff = dt_diff if not np.isinf(dt_diff) else 0.0 + if np.isinf(dt_estimate): + return np.inf + return _dimensionalise_dt(dt_estimate) + if basis != "accuracy": + raise ValueError(f"basis must be 'accuracy' or 'resolution', not {basis!r}.") + + comm = uw.mpi.comm + values = np.asarray(self.u.array).reshape(-1) + lo = comm.allreduce(float(values.min()) if values.size else np.inf, op=MPI.MIN) + hi = comm.allreduce(float(values.max()) if values.size else -np.inf, op=MPI.MAX) + field_range = hi - lo + + if self._last_change_rate is not None: + rate = self._last_change_rate + else: + rate = self._advective_rate() + self.dt_accuracy = fraction * field_range / rate if rate > 0.0 else np.inf + if np.isinf(self.dt_accuracy) or field_range <= 0.0: return np.inf - return _dimensionalise_dt(dt_estimate) + return _dimensionalise_dt(self.dt_accuracy) + + def _advective_rate(self): + r"""Global maximum of :math:`|\mathbf{u}\cdot\nabla\phi|` at the mesh vertices. + + The gradient is the Clement recovery at the vertices (no point + location, so it is safe on a mesh carrying many variables) and the + velocity is evaluated at the same points. + """ + from mpi4py import MPI + from underworld3.function.gradient_evaluation import compute_clement_gradient_at_nodes + + coords = np.asarray(self.mesh.X.coords) + n = coords.shape[0] + if n: + grad = np.asarray(compute_clement_gradient_at_nodes(self.u), dtype=float).reshape(n, -1) + vel = uw.function.evaluate(self._V_fn, coords) + vel = np.asarray(getattr(vel, "magnitude", vel), dtype=float).reshape(n, -1) + local = float(np.abs((vel[:, :grad.shape[1]] * grad).sum(axis=1)).max()) + else: + local = 0.0 + return uw.mpi.comm.allreduce(local, op=MPI.MAX) def solve( self, @@ -569,6 +645,13 @@ def solve( self.DuDt.update_pre_solve(dt, verbose=verbose) super().solve(zero_init_guess, _force_setup, divergence_retries=divergence_retries) _invalidate_solution_cache(self.u) + # The realised rate of change of the field over this step feeds the + # accuracy-based estimate_dt; psi_star[0] still holds phi^n here. + from mpi4py import MPI + change = np.abs(np.asarray(self.u.array).reshape(-1) + - np.asarray(self.DuDt.psi_star[0].array).reshape(-1)) + local = float(change.max()) if change.size else 0.0 + self._last_change_rate = uw.mpi.comm.allreduce(local, op=MPI.MAX) / dt self.DuDt.update_post_solve(dt, verbose=verbose) self.is_setup = True diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index 590dc2b61..69a8f6657 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -191,3 +191,30 @@ def metric(pts): assert adv._custom_mg is None data = np.asarray(T.array[:, 0, 0]) assert np.isfinite(data).all() and 0.9 < data.max() < 1.01 + + +def test_estimate_dt_is_accuracy_based_and_resolution_on_request(mesh): + """The default estimate follows the field, not the mesh; the resolution + basis reproduces the semi-Lagrangian solver's cell-crossing time.""" + adv, T = _solver(mesh, "t") + dt_acc = float(adv.estimate_dt()) + dt_res = float(adv.estimate_dt(basis="resolution")) + assert np.isfinite(dt_acc) and dt_acc > 0 and np.isfinite(dt_res) and dt_res > 0 + # a tighter fraction is a proportionally smaller step + assert float(adv.estimate_dt(fraction=0.01)) == pytest.approx(0.5 * dt_acc) + + x, y = mesh.X + T2 = uw.discretisation.MeshVariable("T_t2", mesh, 1, degree=2) + slcn = uw.systems.AdvDiffusionSLCN(mesh, T2, sympy.Matrix([[-y, x]])) + slcn.constitutive_model = uw.constitutive_models.DiffusionModel + slcn.constitutive_model.Parameters.diffusivity = 0.0 + assert dt_res == pytest.approx(float(slcn.estimate_dt()), rel=1e-12) + + # after a step the estimate uses the realised rate of change + adv.solve(timestep=dt_acc) + assert adv._last_change_rate > 0 + dt_after = float(adv.estimate_dt()) + assert np.isfinite(dt_after) and 0.2 * dt_acc < dt_after < 5 * dt_acc + + with pytest.raises(ValueError, match="basis must be"): + adv.estimate_dt(basis="courant") diff --git a/tests/test_1100_advdiff_supg_rotating_gaussian.py b/tests/test_1100_advdiff_supg_rotating_gaussian.py index 336ea74e4..4415bdc03 100644 --- a/tests/test_1100_advdiff_supg_rotating_gaussian.py +++ b/tests/test_1100_advdiff_supg_rotating_gaussian.py @@ -10,8 +10,8 @@ refined to h/8 across the orbit, at the same timestep, gives the same error to three digits even though its cells sit at a local Courant number of several; -3. the round trip: after one revolution the field returns to its initial - state to a few per cent at a Courant number of one half. +3. the round trip: at the solver's own accuracy-based timestep the field + returns to its initial state after one revolution to under one per cent. Run: pixi run python -m pytest tests/test_1100_advdiff_supg_rotating_gaussian.py -v """ @@ -108,14 +108,23 @@ def metric(pts, _f=fault, _hn=h / 8, _hf=h, _core=0.03, _ramp=0.06): err_band = _run(sol_c, T_c, adv_c, dt, t_end) # the band cells are at a local Courant number well above one - assert dt / float(adv_c.estimate_dt()) > 4.0 + assert dt / float(adv_c.estimate_dt(basis="resolution")) > 4.0 assert abs(err_band - err_uniform) < 0.15 * err_uniform, (err_uniform, err_band) + # the accuracy-based estimate follows the field, so the band does not + # shrink it, while the resolution estimate collapses with the cells + dt_acc_uniform = float(adv.estimate_dt()) + dt_acc_band = float(adv_c.estimate_dt()) + assert abs(dt_acc_band - dt_acc_uniform) < 0.25 * dt_acc_uniform, (dt_acc_uniform, dt_acc_band) + assert float(adv.estimate_dt(basis="resolution")) > 3.0 * float(adv_c.estimate_dt(basis="resolution")) -def test_round_trip_at_moderate_courant(): + +def test_round_trip_at_the_default_timestep(): + """The solver's own defaults: Crank-Nicolson at the accuracy-based step + (2% of the range per step). BDF2 at the same step lands near 1.5%.""" mesh = _box(32) - sol, T, adv = _problem(mesh, "r", 2) - err = _run(sol, T, adv, 0.5 * float(adv.estimate_dt()), float(sol.period)) - assert err < 0.03, err + sol, T, adv = _problem(mesh, "r", 1) + err = _run(sol, T, adv, float(adv.estimate_dt()), float(sol.period)) + assert err < 0.01, err data = np.asarray(T.array[:, 0, 0]) assert data.min() > -0.02 and data.max() < 1.02, (data.min(), data.max()) From bd0b1a5ce8637dc1fd0f8e9392b1faadb348b79c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 22:28:32 -0700 Subject: [PATCH 12/14] Match the Krylov tolerance to the SNES tolerance, and make preconditioner="fmg" a real switch on the SUPG solver The Eulerian SUPG step took two Newton iterations on a linear operator: the Krylov default (rtol 1e-5) does not reach the SNES tolerance (1e-8), and the second Jacobian assembly cost more than every linear solve of the step. The Krylov tolerance is now 1e-9 and a step is one Newton iteration: 1.54 s to 0.91 s per step at 256^2 in serial. Measured against geometric multigrid at matched tolerances (design note, "Preconditioner"), GMRES with additive-Schwarz ILU is the cheaper linear solve at every Courant number from 1/2 to 32 and its iteration count is the same on one and eight ranks; the multigrid's cycle count grows with the Courant number nearly as fast, and a cycle costs about three Schwarz iterations. Schwarz stays the default on every mesh. preconditioner = "fmg" now hands the block to the managed multigrid route (custom-P transfers over the refinement hierarchy or an adapt child's coarse tail, flexible GMRES outside) for the rank count where a one-level method runs out of coarse space. The solver's solve() builds through the base _build, where a preconditioner choice is resolved; the pre-run of the three setup stages marked the solver set up first, so the request was silently inert. The semi-Lagrangian solvers share that pattern and the defect (#683). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- docs/advanced/eulerian-advection-diffusion.md | 11 +- .../design/eulerian-supg-transport.md | 64 +++++++++- .../systems/advection_diffusion_eulerian.py | 112 +++++++++++++++--- tests/test_1055_advdiff_supg_api.py | 30 +++++ 4 files changed, 194 insertions(+), 23 deletions(-) diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md index d62a6d416..950ea84cb 100644 --- a/docs/advanced/eulerian-advection-diffusion.md +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -102,9 +102,14 @@ of the compiled kernels; nothing is recompiled. - The stabilisation parameter uses the local cell size (`mesh.cell_size()`) and three weights that are runtime constants (`solver.tau_weights`); `solver.supg_weight = 0` gives the plain Galerkin scheme for comparison. -- The linear system is nonsymmetric, so the solver defaults to GMRES with an - additive-Schwarz ILU preconditioner instead of algebraic multigrid. Every - option can be overridden through `solver.petsc_options`. +- The linear system is nonsymmetric, so the solver uses GMRES with an + additive-Schwarz ILU preconditioner, with the Krylov tolerance matched to the + SNES tolerance so that a step is one Newton iteration. Measured, this is the + cheaper solve at every Courant number up to eight ranks and its iteration + count does not grow with the rank count. `solver.preconditioner = "fmg"` + switches to geometric multigrid over the mesh's refinement hierarchy + (`refinement >= 1`) for very large rank counts. Every option can be + overridden through `solver.petsc_options`. ## Further reading diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index ff894b4f2..8b95b3220 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -77,10 +77,24 @@ $$ prototype recompiled its kernels on every change (1.2 s against 0.03 s for a step). - **Diffusivity on the constitutive model**, as for every scalar solver, starting at $\kappa = 0$. The prototype carried a float attribute with a warning bridge. -- **Own preconditioner.** The operator is nonsymmetric, so GMRES with an - additive-Schwarz ILU preconditioner replaces the managed GAMG block. The solver - sets `_pc_option_prefix = None`, and the mesh-owned multigrid pickup on adapt - children now respects that (it segfaulted otherwise). +- **Additive-Schwarz ILU, one Newton iteration per step.** The operator is + nonsymmetric, so the smoother and the outer Krylov solver have to be safe for + one. Measured (below), GMRES with an additive-Schwarz ILU preconditioner is the + cheaper linear solve at every Courant number from 1/2 to 32 and its iteration + count does not change between one and eight ranks; geometric multigrid's cycle + count grows with the Courant number nearly as fast, and a cycle costs about + three Schwarz iterations. The linear solve is under a tenth of a step either + way; assembly is the rest. What did matter was the tolerance pair: the Krylov + default (1e-5) does not reach the SNES tolerance (1e-8), so the SNES took a + second Newton step on a linear operator, and that Jacobian assembly cost more + than every linear solve of the step. The Krylov tolerance is now 1e-9. + `preconditioner = "fmg"` hands the block to the managed multigrid route + (custom-P transfers over the refinement hierarchy or an adapt child's coarse + tail, flexible GMRES outside), for the rank count where a one-level method + runs out of coarse space. The solver's `solve()` builds through the base + `_build`, which is where a preconditioner choice is resolved; the + semi-Lagrangian solvers run the three setup stages directly and their + `preconditioner` property is inert as a result (#683). - **Moving meshes, phase 1.** The unknown and its history stay on the default `REMAP` transfer policy with the material velocity. The remap re-interpolates old states onto the new nodes, so the Eulerian form is already correct to @@ -214,6 +228,48 @@ Quarter-turn error on the uniform res-32 mesh with the exact history planted: BDF1 slopes 0.80 and 0.88 between $\Delta t$ = 0.02, 0.01, 0.005; BDF2 slopes above 1.65 between 0.04, 0.02, 0.01. +### Preconditioner + +Level-set advection step (`uw.systems.level_set`, a two-cell band, P2, +Crank-Nicolson) on a structured quad box built with a refinement hierarchy, so +every solver sees the same finest operator; the vortex velocity field of the +level-set study. Wall time per step over ten steps after a warm-up step, on a +sixteen-core workstation. Script and logs: +`~/+Simulations/supg_vs_slcn_657/parallel/fmg_timing.py`, `fmg.log`. + +**Schwarz against geometric multigrid at matched tolerances** (Krylov 1e-9, +SNES 1e-8; one Newton iteration per step for both), 256², three levels: + +| Courant | GMRES + ASM-ILU, its (np 1 / 8) | s/step (np 1 / 8) | fgmres + FMG, cycles (np 1 / 8) | s/step (np 1 / 8) | +|---|---|---|---|---| +| 1/2 | 5 / 5 | 0.913 / 0.121 | 1 / 1 | 0.943 / 0.145 | +| 2 | 8.9 / 8.6 | 0.925 / 0.141 | 3.4 / 3.6 | 1.079 / 0.178 | +| 8 | 16.6 / 16.5 | 0.971 / 0.146 | 12.8 / 12.8 | 1.657 / 0.299 | +| 32 | 37 / 37.8 | 1.128 / 0.172 | 23.8 / 24.1 | 2.347 / 0.457 | + +The multigrid smoother is the managed bundle's gmres/4 + SOR with Galerkin coarse +operators, which inherit the fine-grid $\tau$; four levels instead of three +changes nothing at Courant 1/2 (one cycle, 0.935 s either way), so the coarse +operators are not under-stabilised there. Above Courant 8 the scheme rings (the +range of $\phi$ reaches $-0.29$ to $1.29$ at Courant 8), so the rows where +multigrid's cycle count is closest to the Schwarz count are rows nobody runs. + +**Where the step goes** (`-log_view`, np 1, Courant 1/2, eleven solves): residual +evaluation 4.0 s, Jacobian evaluation 4.4 s, `KSPSolve` 0.36 s under Schwarz and +0.95 s under multigrid. With the Krylov tolerance left at its default of 1e-5 the +Schwarz solver stopped at three iterations, the SNES took a second Newton step +(22 Jacobian assemblies over eleven solves), and the step cost 1.54 s; one +multigrid cycle happens to reduce the residual below the SNES tolerance, so it +took one. That looked like a 1.65x win for multigrid and was a Jacobian +assembly. + +**Controls** (Krylov tolerance at its default, 256², np 1 / 8): algebraic +multigrid (the managed GAMG bundle) 5 iterations, 2.07 / 0.245 s; the "fast" +smoother (richardson/3 + SOR) 0.933 s, the same as gmres/4; gmres/2 needs two +cycles and costs 1.61 s; an ILU smoother 1.62 s. At 512² with four levels the +unmatched rows read 6.12 / 0.84 s (Schwarz, two Newton steps) against 3.71 / +0.58 s (multigrid). + ## What the timestep estimate means The cell-crossing time is not a stability limit for either scheme and says diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index a6c3d36dc..2b7cc2c9e 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -208,9 +208,14 @@ class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): scalar solver; the solver starts with a :class:`~underworld3.constitutive_models.DiffusionModel` at :math:`\kappa = 0` (pure advection). The linear system is nonsymmetric, - so the solver defaults to GMRES with an additive-Schwarz ILU - preconditioner rather than the algebraic multigrid the symmetric scalar - solvers use; every option is overridable through ``petsc_options``. + so the solver uses GMRES with an additive-Schwarz ILU preconditioner, the + Krylov tolerance matched to the SNES tolerance so that a step is one + Newton iteration. ``preconditioner = "fmg"`` hands the linear solve to + geometric multigrid over the mesh's refinement hierarchy (a flexible GMRES + outer solver, Galerkin coarse operators); measured, the Schwarz solve is + cheaper at every Courant number to eight ranks, and multigrid is there for + the rank count where a one-level method runs out of coarse space. Every + option is overridable through ``petsc_options``. """ @timing.routine_timer_decorator @@ -322,18 +327,91 @@ def __init__( self.constitutive_model = uw.constitutive_models.DiffusionModel self.constitutive_model.Parameters.diffusivity = 0.0 - # Nonsymmetric operator: opt out of the managed GAMG/FMG block and - # use GMRES with an additive-Schwarz ILU preconditioner. RCM - # ordering improves the ILU fill on convection-dominated operators. - self._pc_option_prefix = None - self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["ksp_gmres_restart"] = 200 - self.petsc_options["pc_type"] = "asm" - self.petsc_options["sub_pc_type"] = "ilu" - self.petsc_options["sub_pc_factor_mat_ordering_type"] = "rcm" + # Linear solver: additive-Schwarz ILU by default, the managed multigrid + # block on request (see ``preconditioner``). One Newton iteration per + # step: the operator is linear in phi, so the Krylov tolerance must + # reach the SNES tolerance or the SNES takes a second step, and a + # second Jacobian assembly costs more than every linear solve of the + # step (design note, "Preconditioner"). + self._set_linear_solver(multigrid=False) self.petsc_options["snes_rtol"] = 1.0e-8 + self.petsc_options["ksp_rtol"] = 1.0e-9 self.petsc_options["snes_max_it"] = 20 + # ------------------------------------------------------------------ + # Linear solver + # ------------------------------------------------------------------ + + _SCHWARZ_OPTIONS = { + "ksp_type": "gmres", + "ksp_gmres_restart": 200, + "pc_type": "asm", + "sub_pc_type": "ilu", + # RCM ordering improves the ILU fill on a convection-dominated operator. + "sub_pc_factor_mat_ordering_type": "rcm", + } + + def _set_linear_solver(self, multigrid: bool): + """Own the linear solver (GMRES + additive-Schwarz ILU) or hand it to + the managed multigrid block. + + Measured on the level-set advection step at 256^2 and 512^2 (design + note, "Preconditioner"): with the Krylov tolerance matched to the + SNES tolerance, additive Schwarz with ILU is the cheaper linear solve + at every Courant number from 1/2 to 32, its iteration count does not + change between one and eight ranks, and the geometric multigrid's + cycle count grows with the Courant number nearly as fast as the + Schwarz iteration count while each cycle costs about three Schwarz + iterations. The linear solve is under a tenth of the step either way; + assembly is the rest. Multigrid keeps its coarse space for a rank + count where a one-level method runs out of one, which is what + ``preconditioner = "fmg"`` is for. + """ + from underworld3.utilities import multigrid_options + + opts = self.petsc_options + bundle_keys = set() + for bundle in (multigrid_options.gamg_bundle(), + multigrid_options.geometric_mg_bundle()): + bundle_keys |= set(bundle.settings) | set(bundle.stale) + if multigrid: + # The managed block starts from the scalar solver's own keys + # (GMRES + the GAMG bundle) and _apply_preconditioner_options + # resolves the request against the mesh hierarchy at build time. + self._pc_option_prefix = "" + for key in self._SCHWARZ_OPTIONS: + opts.delValue(key) + self._push_managed_option("ksp_type", "gmres") + for key, value in multigrid_options.gamg_bundle().settings.items(): + self._push_managed_option(key, value) + else: + self._pc_option_prefix = None + for key in bundle_keys | {"ksp_type"}: + opts.delValue(key) + self._managed_pc_options.pop(self.petsc_options_prefix + key, None) + for key, value in self._SCHWARZ_OPTIONS.items(): + opts[key] = value + + @property + def preconditioner(self): + """Linear preconditioner: ``"auto"`` (default), ``"fmg"`` or ``"gamg"``. + + ``"auto"`` is GMRES with an additive-Schwarz ILU preconditioner, the + measured choice for this operator (see :meth:`_set_linear_solver`). + ``"fmg"`` hands the block to the managed geometric-multigrid route: + custom-P transfers over the mesh's refinement hierarchy or an adapt + child's coarse tail, installed on the live PC at the first solve, + under a flexible GMRES outer solver; without a hierarchy it warns and + degrades to GAMG. ``"gamg"`` is algebraic multigrid. Setting the + property rebuilds the solver at the next solve. + """ + return self._preconditioner + + @preconditioner.setter + def preconditioner(self, value): + SNES_Scalar.preconditioner.fset(self, value) + self._set_linear_solver(multigrid=self._preconditioner != "auto") + def _object_viewer(self): from IPython.display import Latex, display @@ -637,10 +715,12 @@ def solve( self._needs_function_rewire = True if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - if not self.is_setup: - self._setup_pointwise_functions(verbose) - self._setup_discretisation(verbose) - self._setup_solver(verbose) + # The base ``_build`` resolves the preconditioner choice against the + # mesh hierarchy before the SNES reads its options. Running the three + # setup stages directly here (the semi-Lagrangian solvers' pattern) + # marks the solver set up, so ``_build`` returned early and the + # geometric-multigrid request was silently inert (#683). + self._build(verbose) self.DuDt.update_pre_solve(dt, verbose=verbose) super().solve(zero_init_guess, _force_setup, divergence_retries=divergence_retries) diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index 69a8f6657..4d53ffeb5 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -164,6 +164,36 @@ def test_galerkin_baseline_needs_no_rebuild(mesh): assert adv.supg_weight == 0.0 +def test_multigrid_is_one_switch_away_on_a_refinement_hierarchy(): + """The default linear solver is GMRES + additive-Schwarz ILU on any mesh, + one Newton iteration per step. ``preconditioner = "fmg"`` on a mesh with + a refinement hierarchy hands the block to geometric multigrid: custom-P + transfers over ``mesh.dm_hierarchy`` installed on the live PC at the next + solve, under a flexible outer Krylov solver; the two agree to the solve + tolerance, and switching back rebuilds the Schwarz solver.""" + refined = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.5, qdegree=3, + refinement=2) + schwarz, T_s = _solver(refined, "schwarz") + schwarz.solve(timestep=0.01) + assert schwarz.snes.getKSP().getPC().getType() == "asm" + assert schwarz.snes.getIterationNumber() == 1 + + multigrid, T_m = _solver(refined, "multigrid") + multigrid.preconditioner = "fmg" + multigrid.solve(timestep=0.01) + ksp = multigrid.snes.getKSP() + assert ksp.getType() == "fgmres" + assert ksp.getPC().getType() == "mg" + assert ksp.getPC().getMGLevels() == len(refined.dm_hierarchy) == 3 + a, b = np.array(T_s.array[:, 0, 0]), np.array(T_m.array[:, 0, 0]) + assert np.abs(a - b).max() < 1e-6 * np.abs(a).max() + + multigrid.preconditioner = "auto" + multigrid.solve(timestep=0.01) + assert multigrid.snes.getKSP().getPC().getType() == "asm" + + def test_solves_on_an_adapt_child_with_its_own_preconditioner(): """An adapt child carries a mesh-owned multigrid hierarchy that the solver base installs opportunistically. This solver owns its (additive From ce783b1b1ac50b5f4ba6184d729ab9cda435e756 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 3 Sep 2026 22:33:07 -0700 Subject: [PATCH 13/14] Design note: the 512^2 rows at matched tolerance --- docs/developer/design/eulerian-supg-transport.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 8b95b3220..5ac561e2b 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -268,7 +268,8 @@ multigrid (the managed GAMG bundle) 5 iterations, 2.07 / 0.245 s; the "fast" smoother (richardson/3 + SOR) 0.933 s, the same as gmres/4; gmres/2 needs two cycles and costs 1.61 s; an ILU smoother 1.62 s. At 512² with four levels the unmatched rows read 6.12 / 0.84 s (Schwarz, two Newton steps) against 3.71 / -0.58 s (multigrid). +0.58 s (multigrid); matched, with the shipped defaults, 3.51 / 0.48 s (Schwarz, +5 iterations) against 3.62 / 0.54 s (multigrid, one cycle). ## What the timestep estimate means From 4295af77e28d276e7faf580e4c05f2470d8955d9 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 4 Sep 2026 12:08:48 -0700 Subject: [PATCH 14/14] Let theta be set after construction, as the semi-Lagrangian solver allows The shipped convection examples set adv_diff.theta = 0.5 after building the solver; the Eulerian drop-in refused it. The blend is a runtime constant refreshed from the history manager before every solve, so the setter updates it without a recompile (order 1 only, the constructor's rule). Vector and tensor unknowns join the design note's deferred list: the solver is scalar, where the semi-Lagrangian trace-back carries them. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL --- .../design/eulerian-supg-transport.md | 3 ++- .../systems/advection_diffusion_eulerian.py | 18 +++++++++++++++++- tests/test_1055_advdiff_supg_api.py | 15 +++++++++++++++ 3 files changed, 34 insertions(+), 2 deletions(-) diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md index 5ac561e2b..55999bd3e 100644 --- a/docs/developer/design/eulerian-supg-transport.md +++ b/docs/developer/design/eulerian-supg-transport.md @@ -103,7 +103,8 @@ $$ - **Not yet:** discontinuity capturing (the prototype's residual omitted the time derivative and added first-order diffusion everywhere; a correct lagged residual needs $\phi^{n-1}$), a streamline element length from a mesh-owned metric tensor, - the ALE hook. + the ALE hook, and vector or tensor unknowns: the solver is scalar, where the + semi-Lagrangian trace-back carries vectors and tensors through the same machinery. ## Measurements diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 2b7cc2c9e..d294bdb35 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -441,9 +441,25 @@ def order(self) -> int: @property def theta(self) -> float: - """Adams-Moulton blend at order 1 (1.0 backward Euler, 0.5 Crank-Nicolson).""" + """Adams-Moulton blend at order 1 (1.0 backward Euler, 0.5 Crank-Nicolson). + + Settable after construction, as on the semi-Lagrangian solver: the + blend is a runtime constant of the compiled kernels, refreshed from + the history manager before every solve, so nothing is recompiled. + """ return self._theta + @theta.setter + def theta(self, value): + value = float(value) + if value != 1.0 and self._time_order != 1: + raise ValueError( + "theta applies at order 1 only (0.5 is Crank-Nicolson, 1.0 is " + "backward Euler); order 2 and 3 take theta=1.0." + ) + self._theta = value + self.DuDt.theta = value + @property def delta_t(self): r"""The timestep :math:`\Delta t` as a UW expression. diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py index 4d53ffeb5..9b39dc983 100644 --- a/tests/test_1055_advdiff_supg_api.py +++ b/tests/test_1055_advdiff_supg_api.py @@ -49,6 +49,21 @@ def test_slcn_order_theta_pairs_select_the_documented_schemes(mesh): _solver(mesh, "p4", order=2, theta=0.5) +def test_theta_is_settable_after_construction_as_on_slcn(mesh): + """The convection examples set ``adv_diff.theta = 0.5`` after constructing + the semi-Lagrangian solver; the drop-in accepts the same, refreshing the + Adams-Moulton weights at the next solve without a recompile.""" + adv, _T = _solver(mesh, "th") + adv.solve(timestep=0.01) + key = adv._current_jit_cache_key + adv.theta = 1.0 + adv.solve(timestep=0.01) + assert adv.theta == 1.0 and adv.DuDt.theta == 1.0 + assert adv._current_jit_cache_key == key + with pytest.raises(ValueError, match="theta applies"): + _solver(mesh, "th2", order=2)[0].theta = 0.5 + + def test_semi_lagrangian_only_arguments_are_ignored_with_a_warning(mesh): with pytest.warns(UserWarning, match="monotone_mode, old_frame_traceback"): adv, _T = _solver(mesh, "q", monotone_mode="clamp", old_frame_traceback=True)