UWTN: update the time-derivatives note for the composing solvers

articles/symbolic-time-derivatives-in-underworld3/examples/timestepping.py

@@ -0,0 +1,135 @@
+#!/usr/bin/env python3
+"""One field, three time derivatives, the same answer.
+
+A Gaussian blob of temperature is carried once across a periodic-in-effect
+box by a uniform wind while it diffuses. The exact answer is known -- a
+Gaussian advects rigidly and spreads as sqrt(4 k t) -- so each scheme can be
+scored against it rather than against the others.
+
+The point is that the SOLVER is the same in all three runs. Only the history
+manager changes, and with it whether transport is assembled implicitly in the
+residual (Eulerian SUPG), traced along characteristics (semi-Lagrangian), or
+corrected explicitly on the mesh (Eulerian).
+
+What this shows and what it does not: uniform translation of a smooth blob is
+the easiest case there is for a semi-Lagrangian scheme -- the characteristic is
+a straight line and one interpolation per step is nearly exact -- so it wins on
+accuracy here and the cost column is the interesting one. The case that decided
+the default is a convection benchmark, where the flow turns and the
+interpolation error accumulates over a full circuit; see the design note
+`eulerian-supg-transport.md` for that comparison.
+
+Dimensional throughout: metres, seconds, kelvin. The Courant number is set
+deliberately above 1 for the last run, which is where the schemes part company.
+
+Usage:
+    python3 timestepping.py                    # all three, default resolution
+    python3 timestepping.py -uw_cells 48       # finer
+    python3 timestepping.py -uw_courant 2.0    # past the CFL limit
+"""
+
+import time
+
+import numpy as np
+import sympy
+
+import underworld3 as uw
+
+params = uw.Params(
+    cells=uw.Param(32, "cells across the box"),
+    courant=uw.Param(0.5, "Courant number: |v| dt / h"),
+    travel=uw.Param(0.5, "how far the blob is carried, in box widths"),
+)
+
+# --- the physical problem, in units -----------------------------------------
+L = uw.quantity(1.0e6, "m")            # box side, 1000 km
+V = uw.quantity(1.0e-9, "m/s")         # wind, ~3 cm/yr
+KAPPA = uw.quantity(1.0e-6, "m**2/s")  # thermal diffusivity of rock
+T0 = uw.quantity(100.0, "K")           # blob amplitude
+WIDTH = 0.08                           # blob sigma, as a fraction of L
+
+CELLS = int(params.cells)
+h = L / CELLS
+dt = float(params.courant) * h / V     # a time, in seconds
+
+
+def gaussian(x, y, cx, cy, sigma):
+    """The exact blob: amplitude T0, centre (cx, cy), width sigma."""
+    r2 = (x - cx) ** 2 + (y - cy) ** 2
+    return float(T0.magnitude) * np.exp(-r2 / (2.0 * sigma ** 2))
+
+
+def run(flavour):
+    """Advect and diffuse the blob for `steps`, return the L2 error."""
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
+        cellSize=1.0 / CELLS, qdegree=3)
+    T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2)
+    x, y = mesh.X
+
+    # non-dimensional velocity on a unit box: the scaling is L for length and
+    # L/V for time, so the wind is 1 and the diffusivity is 1/Peclet
+    peclet = float((V * L / KAPPA).magnitude)
+    v_fn = sympy.Matrix([[1.0, 0.0]])          # unit wind, +x
+
+    sigma = WIDTH
+    x0, y0 = 0.25, 0.5
+    T.array[:, 0, 0] = gaussian(np.asarray(T.coords)[:, 0],
+                                np.asarray(T.coords)[:, 1], x0, y0, sigma)
+
+    if flavour == "supg":
+        solver = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v_fn)
+    elif flavour == "slcn":
+        solver = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=v_fn)
+    elif flavour == "eulerian":
+        DTdt = uw.systems.Eulerian_DDt(
+            mesh, T, vtype=uw.VarType.SCALAR, degree=T.degree,
+            continuous=True, order=1)
+        solver = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v_fn,
+                                         DuDt=DTdt, order=1)
+    else:
+        raise ValueError(flavour)
+
+    solver.constitutive_model = uw.constitutive_models.DiffusionModel
+    solver.constitutive_model.Parameters.diffusivity = 1.0 / peclet
+    for wall in ("Bottom", "Top", "Left", "Right"):
+        solver.add_dirichlet_bc(0.0, wall)
+
+    # Compare at equal PHYSICAL time, not equal step count: a bigger Courant
+    # number buys fewer steps, which is the whole point of an unconditionally
+    # stable scheme. Comparing at fixed step count flatters whichever scheme
+    # is given the smaller timestep.
+    dt_nd = float(params.courant) / CELLS      # |v|=1 on the unit box
+    n = max(1, int(round(float(params.travel) / dt_nd)))
+    t0 = time.perf_counter()
+    for _step in range(n):
+        solver.solve(timestep=dt_nd)
+    per_step = (time.perf_counter() - t0) / n
+
+    # the exact blob: advected by v*t, spread by the diffusion it has seen
+    t_end = n * dt_nd
+    spread = np.sqrt(sigma ** 2 + 2.0 * t_end / peclet)
+    amp = float(T0.magnitude) * sigma ** 2 / spread ** 2
+    coords = np.asarray(T.coords)
+    exact = amp * np.exp(
+        -((coords[:, 0] - (x0 + t_end)) ** 2 + (coords[:, 1] - y0) ** 2)
+        / (2.0 * spread ** 2))
+    got = T.array[:, 0, 0]
+    err = float(np.sqrt(np.mean((got - exact) ** 2)) / float(T0.magnitude))
+    return err, n, per_step
+
+
+if __name__ == "__main__":
+    uw.pprint(f"box {L}, wind {V}, diffusivity {KAPPA}")
+    uw.pprint(f"{CELLS} cells, Courant {float(params.courant):g}, "
+              f"dt {dt.to('year')}, carried {float(params.travel):g} box widths\n")
+    have_supg = hasattr(uw.systems.ddt, "EulerianSUPG")
+    uw.pprint(f"  {'manager':<24} {'L2 error':>10} {'steps':>7} {'s/step':>9}")
+    for flavour, name in (("supg", "Eulerian SUPG (default)"),
+                          ("slcn", "Semi-Lagrangian"),
+                          ("eulerian", "Eulerian")):
+        if flavour == "supg" and not have_supg:
+            uw.pprint(f"  {name:<24} {'needs EulerianSUPG':>10}")
+            continue
+        err, n, per_step = run(flavour)
+        uw.pprint(f"  {name:<24} {err:>10.4f} {n:>7d} {per_step:>9.3f}")

articles/symbolic-time-derivatives-in-underworld3/metadata.yml

@@ -9,7 +9,7 @@ authors:
     orcid: 0000-0003-3685-174X
     affiliation: Australian National University
 publication_date: 2026-04-16
-version: 1.0.0
+version: 1.1.0
 legacy_doi: null
 archive_doi: 10.6084/m9.figshare.33193596
 license: CC-BY-4.0

articles/symbolic-time-derivatives-in-underworld3/symbolic-time-derivatives-in-underworld3.md

@@ -23,7 +23,7 @@ exports:
     template: ../../templates/pdf
     output: symbolic-time-derivatives-in-underworld3.pdf
     article_id: UWTN 2026-007
-    article_version: 1.0.0
+    article_version: 1.1.0
 parts:
   abstract: "In Underworld3, the time derivative is a symbolic object. It appears in the solver's strong form as a SymPy expression, alongside the constitutive stress and the body force."
 doi: 10.6084/m9.figshare.33193596
@@ -44,17 +44,19 @@ In many finite element codes, these choices are baked into the solver implementa
 
 ## The DDt Hierarchy
 
-UW3 provides four implementations of the time derivative, all sharing the same calling interface and working for scalar, vector and tensor quantities.
+UW3 provides five implementations of the time derivative, all sharing the same calling interface and working for scalar, vector and tensor quantities.
 
 **Eulerian** stores history on the mesh. The material derivative is approximated by finite differences in time at fixed grid points, with an advection correction. This is the classical approach: simple, mesh-based, but subject to CFL stability constraints when advection dominates.
 
+**Eulerian SUPG** also stores history on the mesh, but it does not correct for advection after the fact. It hands the solver an advection term to assemble *inside* the residual, stabilised by SUPG, so the transport is solved implicitly along with everything else. There is no CFL limit. The plain Eulerian flavour above and this one are different answers to the same problem: one corrects the history explicitly and pays a stability condition, the other puts the transport in the equation and does not.
+
 **Semi-Lagrangian** traces characteristics backward in time. At each mesh node, it asks: where was this material parcel at the previous timestep? It then interpolates the previous solution at that departure point. This is unconditionally stable because the material derivative is evaluated along the characteristic, not at a fixed grid point. There is no CFL constraint, though the user needs to be conscious of accuracy trade-offs inherent in the scheme.
 
 **Lagrangian** follows particles through the flow. The history is stored on swarm variables and advected with the particles. This is the natural choice when material history matters physically, as in viscoelastic stress transport where the stress tensor must be advected and rotated with the material. Accuracy of this scheme depends upon the quality of the particle-layout (density of particles, presence of gaps).
 
 **Symbolic** provides pure symbolic history without mesh or swarm storage. It is used internally by the constitutive model system for building expressions that involve time derivatives.
 
-All four produce the same thing: a symbolic expression that can be embedded in a solver's weak form. The solver does not need to know which implementation is providing the time derivative. It sees a SymPy expression for $D\phi/Dt$ and includes those terms when it differentiates the equation system to determine the Jacobians.
+All five produce the same thing: a symbolic expression that can be embedded in a solver's weak form. The solver does not need to know which implementation is providing the time derivative. It sees a SymPy expression for $D\phi/Dt$ and includes those terms when it differentiates the equation system to determine the Jacobians.
 
 ## BDF Schemes: Multi-Step Time Integration
 
@@ -90,6 +92,8 @@ At order 0, the flux is evaluated purely at the current time (fully implicit). A
 
 In UW3, the solver's flux time derivative (`DFDt`) provides an `adams _ moulton _ flux()` method that returns the appropriately weighted combination of the current flux and previous flux values as symbolic forms backed by stored evaluations. This expression then appears in the solver's $F _ 1$ template as a symbolic expression. Like the BDF coefficients, the AM weights are UWexpressions that update between timesteps.
 
+The composing solver takes these level weights from `DuDt.spatial _ weights()` rather than from a separate `DFDt`, so one manager decides both how the field is transported and how the flux is weighted in time.
+
 The combination of BDF for the time derivative and AM for the flux evaluation gives a family of time integration schemes. BDF-1 with AM-0 is backward Euler. BDF-2 with AM-1 gives second-order accuracy in both the time derivative and the flux evaluation. The user controls this through the `order` parameter when creating the solver.
 
 ## Order Ramping at Startup
@@ -120,13 +124,22 @@ adv_diff.f = H
 adv_diff.solve(timestep=dt)
 ```
 
-The `order` parameter controls the time discretisation. The solver builds its weak form from two template expressions, $F _ 0$ (force-like, paired with the test function) and $F _ 1$ (flux-like, paired with the test function gradient):
+The `order` parameter controls the time discretisation. The solver builds its weak form from two template expressions, $F _ 0$ (force-like, paired with the test function) and $F _ 1$ (flux-like, paired with the test function gradient). It asks the history manager for each piece:
+
+```python
+F0 = DuDt.time_derivative() + DuDt.advection() - H
+F1 = sum_k w_k * k * grad(T_k) + DuDt.stabilisation_flux(R)
+```
+
+The manager decides how transport is done, and the solver does not know. A semi-Lagrangian manager has already moved the field along its characteristics, so it answers zero for `advection()` and for `stabilisation_flux()`, and the pair reduces to the older form:
 
 ```python
 F0 = DuDt.bdf() / delta_t - H
 F1 = DFDt.adams_moulton_flux()
 ```
 
+which is exactly what `uw.systems.AdvDiffusionSLCN` assembles. The Eulerian SUPG manager instead answers with an advection term and a stabilisation flux, and the same solver becomes an implicit Eulerian one.
+
 At **order 1** (backward Euler / fully implicit), these expand to:
 
 $$
@@ -170,28 +183,35 @@ The solve sequence for each timestep is:
 The choice of DDt type is a one-parameter decision at solver construction:
 
 ```python
-# Default for advection-diffusion: Semi-Lagrangian (unconditionally stable)
+# Default for advection-diffusion: Eulerian SUPG (implicit, no CFL limit)
 adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v)
 
+# The semi-Lagrangian solver keeps its own name
+adv_diff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=v)
+
 # Override with Lagrangian (particle-based, requires a swarm)
 DTdt = uw.systems.Lagrangian_Swarm_DDt(
     swarm, psi_fn=T.sym,
     vtype=uw.VarType.SCALAR, degree=T.degree,
     continuous=True, order=2
 )
-adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v, DuDt=DTdt)
+adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v,
+                                   DuDt=DTdt, order=2)
 
-# Or Eulerian (mesh-based, for problems without strong advection)
+# Or Eulerian (mesh-based, explicit advection correction)
 DTdt = uw.systems.Eulerian_DDt(
     mesh, T, vtype=uw.VarType.SCALAR,
     degree=T.degree, continuous=True, order=2
 )
-adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v, DuDt=DTdt)
+adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v,
+                                   DuDt=DTdt, order=2)
 ```
 
-The solver does not need to be made aware of which DDt type you chose. It calls `bdf()` and `adams _ moulton _ flux()` and gets SymPy expressions. The physics of the time discretisation is encapsulated in the DDt object. The numerics of the spatial discretisation are encapsulated in the solver. They communicate through symbolic expressions.
+A supplied manager fixes the order, so the solver has to be given the same one; a mismatch raises rather than quietly using two different schemes.
 
-Each solver type has a sensible default. Advection-diffusion and Stokes default to Semi-Lagrangian. Pure diffusion defaults to Eulerian. Viscoelastic solvers use the DFDt infrastructure for stress history on particles. You only need to override the default when your problem requires it.
+The solver does not need to be made aware of which DDt type you chose. It asks the manager for a time derivative, an advection term and a stabilisation flux, and gets SymPy expressions. The physics of the time discretisation is encapsulated in the DDt object. The numerics of the spatial discretisation are encapsulated in the solver. They communicate through symbolic expressions.
+
+Each solver type has a sensible default. Advection-diffusion and Navier-Stokes default to the Eulerian SUPG manager; the semi-Lagrangian solvers keep their `SLCN` names and their place at large Courant numbers, where tracing a characteristic beats stabilising a residual. What decided it was cost and stability rather than accuracy: assembling the transport is several times cheaper per step than tracing characteristics and interpolating, and it does not care what the Courant number is. On smooth translation the semi-Lagrangian scheme is still the more accurate of the two, which is why it keeps its place; the convection benchmarks, where the flow turns and the error accumulates over a full circuit, are where the two draw level on accuracy and the cost difference decides. There is more on that comparison in [Two Ways to Move a Field](/two-ways-to-move-a-field/). Stokes' viscoelastic stress history is still semi-Lagrangian, and pure diffusion is still Eulerian. You only need to override the default when your problem requires it.
 
 ## Why This Matters
 
@@ -201,6 +221,21 @@ In UW3, the time derivative is an object you can create, configure, inspect, and
 
 This is the same design principle we described in the [constitutive models post](/constitutive-models-in-symbolic-form/): separate the physics from the numerics, connect them through symbolic expressions, and make both sides inspectable. For constitutive models, the boundary is the stress tensor. For time derivatives, it is the BDF/AM expression. In both cases, the solver sees a SymPy expression and does not need to know how it was constructed.
 
+## History
+
+- **1.1.0** — 2026-09-08
+  Updated for the composing solvers. `uw.systems.AdvDiffusion` and
+  `NavierStokes` now take their transport from the history manager and default
+  to `EulerianSUPG`; the semi-Lagrangian classes keep the `SLCN` names. The
+  DDt hierarchy gains a fifth flavour and the solver's template expressions are
+  written in the composing form. A runnable example was added
+  (`examples/timestepping.py`), and the reason the default moved is stated as
+  cost and stability rather than accuracy, which is what the measurements show.
+  The scheme theory — BDF, Adams-Moulton, order ramping — is unchanged, and so
+  is the viscoelastic stress history, which is still semi-Lagrangian.
+- **1.0.0** — 2026-04-16 · [10.6084/m9.figshare.33193596.v1](https://doi.org/10.6084/m9.figshare.33193596.v1)
+  First published.
+
 <!-- uwtn-acknowledgement -->
 
 *The Underworld project is supported by AuScope and the Australian Government through the National Collaborative Research Infrastructure Strategy (NCRIS). Source code:* [*github.com/underworldcode/underworld3*](https://github.com/underworldcode/underworld3)