Conservative level-set transport on the shared solver interface (NengLu's #657 work)

docs/advanced/eulerian-advection-diffusion.md

@@ -0,0 +1,123 @@
+# 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 = 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 |
+|---|---|---|
+| `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()` | 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 |
+
+`order=3` (BDF3) is available; see below for when it is safe.
+
+## 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.
+
+In parallel the Eulerian step stays seven times cheaper in serial and about five
+times at eight ranks (the departure-point search parallelises perfectly, the
+ILU preconditioner a little less), and its answer is identical to ten digits
+at every rank count, where the semi-Lagrangian answer moves with the partition.
+
+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 (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
+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 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
+
+- 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`

docs/advanced/index.md

@@ -139,6 +139,8 @@ custom-meshes
 curved-boundary-conditions
 mesh-adaptation
 semi-lagrangian-time-integration
+eulerian-advection-diffusion
+level-set-transport
 porous-flow
 snapshot-restore
 troubleshooting

docs/advanced/level-set-transport.md

@@ -0,0 +1,93 @@
+# Conservative level sets
+
+`uw.systems.LevelSetSolver` carries a material interface as the 0.5 contour of a
+smoothed indicator
+
+$$
+\psi = \tfrac12\left(1 + \tanh\frac{\varphi}{2\varepsilon}\right),
+$$
+
+where $\varphi$ is the signed distance to the interface (positive inside) and
+$\varepsilon$ the interface thickness, a fraction of the local cell size. The
+field is transported by an ordinary scalar solver; what makes it a level set is
+what happens after each step:
+
+- **reinitialisation** restores the $\tanh$ profile without moving the 0.5
+  contour (Parameswaran and Mandal 2023, integrated in pseudo-time with SSP-RK3);
+- **mass correction** restores the enclosed volume by a uniform, clipped shift
+  found by bisection (Zhang, Zou and Greaves 2010).
+
+Neither depends on the transport scheme, so the solver takes either the Eulerian
+SUPG solver (the default) or the semi-Lagrangian one.
+
+```python
+from underworld3.systems import level_set
+
+psi = uw.discretisation.MeshVariable("psi", mesh, 1, degree=2)
+eps = level_set.interface_thickness(mesh, psi, scale=0.35)
+level_set.initialise_psi(psi, eps, interface_geometry="polygon",
+                         interface_coordinates=circle_points)   # or signed_distance=...
+
+ls = uw.systems.LevelSetSolver(psi, velocity=v.sym, epsilon=eps)  # advection="slcn" to compare
+for step in range(n_steps):
+    ls.solve(dt)                     # advect, reinitialise when due, restore the volume
+
+viscosity = level_set.material_property_field(psi.sym[0], [eta_outside, eta_inside], "geometric")
+```
+
+## Choices
+
+| argument | meaning |
+|---|---|
+| `advection` | `"supg"` (default) or `"slcn"`; both run pure advection |
+| `order`, `theta` | the transport solver's time scheme; Crank-Nicolson by default, which preserves the profile's amplitude between reinitialisations |
+| `reini_frequency`, `reini_steps`, `reini_dt` | how often, how many pseudo-time steps, and how long each is (half the smallest $\varepsilon$ by default) |
+| `far_field` | the value of $\psi$ imposed on the domain boundary; set it whenever the flow crosses the boundary (an inflow boundary with no value lets mass in) |
+| `conserve_mass` | `"auto"` (default): the global correction is on for `"slcn"`, which loses volume by interpolation, and off for `"supg"`, which conserves it to solver tolerance on its own; the clip to [0, 1] of the ringing at a one-cell band then costs about 0.2% per revolution, which `volume_drift` reports |
+| `adv_solver_bc` | box wall labels on which a zero normal gradient is imposed by copying the neighbouring interior nodes |
+
+**Band thickness.** `interface_thickness(scale=0.35)`, the g-adopt default,
+gives a band well under one cell, which a continuous-Galerkin transport rings
+at. Measured on a rotating circle at 32 cells across, one revolution, SUPG with
+no mass correction:
+
+| `scale` | $\varepsilon / h$ | volume drift |
+|---|---|---|
+| 0.35 | 0.12 | +0.84% (clipped ringing) |
+| 1.0 | 0.36 | +0.28% |
+| 2.0 | 0.71 | -0.18% (ringing gone) |
+| 3.0 | 1.07 | -0.85% (reinitialisation curvature error) |
+
+For the SUPG transport a `scale` of 1.5 to 2, a band of two to three cells, is
+the sensible setting; the thickness trades interface resolution for a clean
+transport.
+
+`initialise_psi` accepts a precomputed signed distance, or a polygon, curve or
+`shapely` geometry (the latter three need the optional `shapely` package).
+`material_property_field` blends a property across one or more level sets with
+a sharp, arithmetic, geometric or harmonic transition.
+
+## Cost
+
+Per step at 64 by 64 (LeVeque flow, Courant 0.5, reinitialisation every fifth
+step): the SUPG advection takes 0.13 s and the SLCN advection 1.24 s; the
+reinitialisation 0.06 to 0.11 s averaged; the mass correction 0.13 to 0.24 s.
+Since the Eulerian transport does not need the correction, its level-set step
+costs about 0.19 s against 1.59 s for the semi-Lagrangian one.
+
+## Which transport solver
+
+On the LeVeque swirling flow at 64 by 64 and Courant 0.5 (period 2), the SUPG
+level set returns with a shape error of 0.028 against 0.051 for the
+semi-Lagrangian one, at half the wall time; the mass correction pins both to
+the same volume. The Eulerian solver's advantage is the same as for any scalar:
+no interpolation loss per step, and cells refined for the Stokes problem cost
+nothing. See {doc}`eulerian-advection-diffusion`. The example
+`docs/examples/convection/advanced/Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py` runs the
+comparison.
+
+## Credit
+
+The level-set pipeline, its SUPG transport and the LeVeque comparison are
+NengLu's contribution (issue #657); this module unifies the two variants of
+that work on the shared solver interface.

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 `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
+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