Navier-Stokes with Eulerian SUPG momentum transport, and a partition-independent cell size (#687)

docs/advanced/eulerian-advection-diffusion.md

@@ -0,0 +1,153 @@
+# Advection-diffusion composed from a transport manager (Eulerian SUPG by default)
+
+`uw.systems.AdvDiffusion` is the general scalar transport solver. It assembles the
+diffusive flux and the source itself and takes the transport from the history manager
+it holds (`DuDt`); with the default manager, `uw.systems.ddt.EulerianSUPG`, it is the
+Eulerian SUPG scheme this page describes, a drop-in for the semi-Lagrangian solver
+`uw.systems.AdvDiffusionSLCN`. Both solve
+
+$$
+\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.AdvDiffusion(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.
+
+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 term is
+  also weighted by the cell Péclet number, $Pe^2/(Pe^2 + Pe_c^2)$ with
+  $Pe = |\mathbf{u}| h / 2\kappa$ and $Pe_c$ the `peclet_weight` argument (default 4),
+  so the stabilisation is off where a cell is diffusion-dominated and full where advection
+  dominates (pure advection, $\kappa = 0$, is unaffected); `peclet_weight=0` gives the
+  uniform weight.
+- 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`.
+
+## The history manager is the transport plugin
+
+The solver does not assemble its transport itself. Its history manager (`solver.DuDt`)
+contributes three symbolic terms, and the solver composes its residual from them:
+the time derivative of the scheme, the advection, and the stabilisation flux of the
+strong residual. The default manager is `uw.systems.ddt.EulerianSUPG`, which owns the
+advecting velocity (`V_fn` is data on it), the time scheme (`order`, `theta`), and the
+stabilisation knobs (`supg_weight`, `tau_weights`, `tau_shape`, `peclet_weight`); the
+solver's properties of the same names pass through to it.
+
+Any history manager that follows the contract can be supplied instead. A semi-Lagrangian
+manager answers zero for the advection and the stabilisation, because its history is
+already traced back along the characteristics, so the same solver becomes a
+semi-Lagrangian scheme on the field history:
+
+```python
+history = uw.systems.ddt.SemiLagrangian(mesh, T.sym, v.sym, vtype=uw.VarType.SCALAR,
+                                        degree=T.degree, continuous=True, order=1)
+adv = uw.systems.AdvDiffusion(mesh, T, v.sym, DuDt=history)   # no assembled advection
+```
+
+On pure advection this reproduces `AdvDiffusionSLCN` to the solver tolerance; with
+diffusion the two differ in where the diffusive flux history comes from (the traced-back
+field here, the traced-back flux there). The manager works for a vector or tensor unknown
+as well (`vtype`), applying the advection component by component, which is how the
+Navier-Stokes solver and a transported stress use it.
+
+## 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/eulerian-navier-stokes.md

@@ -0,0 +1,81 @@
+# Navier-Stokes composed from a transport manager (Eulerian SUPG by default)
+
+`uw.systems.NavierStokes` solves the incompressible Navier-Stokes equations with
+the momentum transport taken from the history manager it holds. With the default
+manager, `uw.systems.ddt.EulerianSUPG`, the momentum advection is assembled implicitly
+in the Stokes saddle-point residual and stabilised by the streamline-upwind
+Petrov-Galerkin term, which is the scheme this page describes; the semi-Lagrangian
+solver with a stress history is `uw.systems.NavierStokesSLCN`. It is the vector counterpart of {doc}`eulerian-advection-diffusion` and
+takes the same constructor as `uw.systems.Stokes` plus the density and the time
+scheme:
+
+```python
+ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=1)   # Crank-Nicolson
+ns.constitutive_model = uw.constitutive_models.ViscousFlowModel
+ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 / Re
+ns.add_dirichlet_bc((0.0, 0.0), "Bottom")
+...
+for step in range(n):
+    ns.solve(timestep=dt)
+```
+
+`order=1` is the theta rule (Crank-Nicolson at the default `theta=0.5`), `order=2`
+is BDF2. The velocity history lives on the mesh; there is no stress history, the
+viscous stress at an earlier level is rebuilt from the stored velocity. Pressure
+has no history.
+
+## The advecting velocity
+
+The nonlinear term is $(\mathbf{a}\cdot\nabla)\mathbf{u}^{n+1}$ with $\mathbf{a}$
+chosen by `advection=`:
+
+- `"extrapolated"` (default): $\mathbf{a} = 2\mathbf{u}^n - \mathbf{u}^{n-1}$, a
+  second-order lag. Each step is one linear solve through the Stokes fieldsplit.
+- `picard_iterations=n` re-solves up to `n` more times with the latest iterate as
+  $\mathbf{a}$, stopping when the velocity stops changing (`picard_tolerance`).
+  The fixed point is the fully implicit scheme.
+- `"implicit"`: $\mathbf{a} = \mathbf{u}^{n+1}$ and the SNES takes Newton steps on
+  the quadratic term.
+
+The stabilisation parameter is
+$\tau = [(C_t/\Delta t)^2 + (C_u |\mathbf{a}|/h)^2 + (C_\nu \nu/h^2)^2]^{-1/2}$
+with $h$ the local cell size and the three weights in `ns.tau_weights`;
+`ns.supg_weight = 0` gives the plain Galerkin scheme. The term is also weighted by the
+cell Péclet number, $Pe^2/(Pe^2 + Pe_c^2)$ with $Pe = |\mathbf{a}| h / 2\nu$ and
+$Pe_c$ the `peclet_weight` argument (default 4), so the stabilisation is off where a cell
+is diffusion-dominated, where it is not needed and costs a fixed multiple of the Galerkin
+error, and full where advection dominates; `peclet_weight=0` gives the uniform weight.
+`tau_shape` selects the Brooks-Hughes or doubly asymptotic form of $\tau$ in place of the
+inverse sum. The strong residual the
+term acts on carries the time derivative, the advection, the pressure gradient
+and the body force, but not the viscous term (the kernels see first derivatives
+only), so on a smooth, well-resolved flow the Galerkin form is the more accurate
+one and the stabilisation earns its place where the element Reynolds number
+$\rho|\mathbf{a}|h/\eta$ exceeds one.
+
+## Timestep
+
+`ns.estimate_dt()` returns the step at which the velocity changes by a fraction
+(default 0.02) of its range, from the realised rate of the last step; before the
+first solve, and with `basis="resolution"`, it returns the Stokes solver's
+cell-crossing time.
+
+## The history manager is the transport plugin
+
+The momentum transport is not written into the solver. Its history manager
+(`ns.DuDt`, an `EulerianSUPG` on the velocity) contributes the time derivative, the
+implicit advection $\sum_k w_k (\mathbf{a}_k\cdot\nabla)\mathbf{u}^{(k)}$ and the
+stabilisation flux $\tau\,\mathbf{R}\otimes\mathbf{a}$; the solver adds the density,
+the pressure gradient and the body force to form $\mathbf{R}$, the viscous flux of the
+scheme, and the pressure. The advecting velocity is data on the manager: the solver sets
+it to the extrapolated field, the latest Picard iterate, or the unknown itself according
+to `advection`, and names the stored velocity as the carrier of the stored levels
+(`DuDt.V_fn_history`). The stabilisation knobs (`supg_weight`, `tau_weights`,
+`tau_shape`, `peclet_weight`) and `delta_t` live on the manager and the solver's
+properties pass through. The scalar solver `AdvDiffusion` composes the same three
+terms from the same class, and a semi-Lagrangian manager can be supplied to either.
+
+## Further reading
+
+- Design note and measurements: `docs/developer/design/eulerian-supg-transport.md`
+- The semi-Lagrangian Navier-Stokes solver: `uw.systems.NavierStokesSLCN`

docs/advanced/index.md

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

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.AdvDiffusion` 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

docs/api/solvers.md

@@ -46,6 +46,25 @@ Viscoelastic extension of the Stokes solver.
    :show-inheritance:
 ```
 
+### SNES_AdvectionDiffusion_Composed (`uw.systems.AdvDiffusion`)
+
+The scalar transport solver composed from a DDt transport manager; with the
+default `EulerianSUPG` manager it is the implicit Eulerian SUPG scheme.
+
+```{eval-rst}
+.. autoclass:: underworld3.systems.advection_diffusion_eulerian.SNES_AdvectionDiffusion_Composed
+   :members:
+   :show-inheritance:
+```
+
+### EulerianSUPG (the transport manager)
+
+```{eval-rst}
+.. autoclass:: underworld3.systems.ddt.EulerianSUPG
+   :members:
+   :show-inheritance:
+```
+
 ### SNES_Diffusion
 
 ```{eval-rst}
@@ -81,3 +100,14 @@ Viscoelastic extension of the Stokes solver.
    :members:
    :show-inheritance:
 ```
+
+### SNES_NavierStokes_Composed (`uw.systems.NavierStokes`)
+
+Navier-Stokes composed from a DDt transport manager (Eulerian SUPG momentum
+transport by default); the semi-Lagrangian class above is `uw.systems.NavierStokesSLCN`.
+
+```{eval-rst}
+.. autoclass:: underworld3.systems.navier_stokes_eulerian.SNES_NavierStokes_Composed
+   :members:
+   :show-inheritance:
+```

docs/beginner/tutorials/14-Timestepping-with-physical-units.ipynb

@@ -209,7 +209,7 @@
    "outputs": [],
    "source": [
     "# Create advection-diffusion solver\n",
-    "adv_diff = uw.systems.AdvDiffusion(\n",
+    "adv_diff = uw.systems.AdvDiffusionSLCN(\n",
     "    mesh,\n",
     "    u_Field=T,\n",
     "    V_fn=v,\n",

docs/beginner/tutorials/15-Thermal-convection-with-units.ipynb

@@ -247,7 +247,7 @@
    "source": [
     "# Create solver for the energy equation (Advection-Diffusion of temperature)\n",
     "\n",
-    "adv_diff = uw.systems.AdvDiffusion(\n",
+    "adv_diff = uw.systems.AdvDiffusionSLCN(\n",
     "    meshball,\n",
     "    u_Field=t_soln,\n",
     "    V_fn=v_soln,\n",

docs/beginner/tutorials/7-Timestepping-simple.ipynb

@@ -157,7 +157,7 @@
     "T_initial_field = uw.discretisation.MeshVariable(\"T0\", mesh, 1, degree=3)\n",
     "\n",
     "# Create advection-diffusion solver\n",
-    "adv_diff = uw.systems.AdvDiffusion(\n",
+    "adv_diff = uw.systems.AdvDiffusionSLCN(\n",
     "    mesh,\n",
     "    u_Field=T,\n",
     "    V_fn=v,\n",
@@ -265,7 +265,7 @@
     "## Time Stepping\n",
     "\n",
     "In many time-stepping schemes, the time step is constrained by the CFL (Courant-Friedricks-Levy) condition for stability. In the case of the \n",
-    "Semi-Lagrange advection scheme which is used by default by `uw.systems.AdvDiffusion`, the method is implicit and\n",
+    "Semi-Lagrange advection scheme which is used by default by `uw.systems.AdvDiffusionSLCN`, the method is implicit and\n",
     "should work for large timesteps. However, there remains the concept of an element-crossing time that is fundamental\n",
     "in understanding how numerical timestepping operates.\n",
     "\n",

docs/beginner/tutorials/8-Timestepping-coupled.ipynb

@@ -167,7 +167,7 @@
    "source": [
     "# Create solver for the energy equation (Advection-Diffusion of temperature)\n",
     "\n",
-    "adv_diff = uw.systems.AdvDiffusion(\n",
+    "adv_diff = uw.systems.AdvDiffusionSLCN(\n",
     "    meshball,\n",
     "    u_Field=t_soln,\n",
     "    V_fn=v_soln,\n",
@@ -196,7 +196,7 @@
    },
    "outputs": [],
    "source": [
-    "uw.systems.AdvDiffusion.view()"
+    "uw.systems.AdvDiffusionSLCN.view()"
    ]
   },
   {