diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md new file mode 100644 index 000000000..2b7a1bafc --- /dev/null +++ b/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` diff --git a/docs/advanced/eulerian-navier-stokes.md b/docs/advanced/eulerian-navier-stokes.md new file mode 100644 index 000000000..957b04114 --- /dev/null +++ b/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` diff --git a/docs/advanced/index.md b/docs/advanced/index.md index cb47f8dea..3a2605939 100644 --- a/docs/advanced/index.md +++ b/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 diff --git a/docs/advanced/semi-lagrangian-time-integration.md b/docs/advanced/semi-lagrangian-time-integration.md index 23a935db2..68b12b609 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.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 diff --git a/docs/api/solvers.md b/docs/api/solvers.md index 6b5c464f3..55027e56e 100644 --- a/docs/api/solvers.md +++ b/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: +``` diff --git a/docs/beginner/tutorials/14-Timestepping-with-physical-units.ipynb b/docs/beginner/tutorials/14-Timestepping-with-physical-units.ipynb index 21930707d..f7527ec1b 100644 --- a/docs/beginner/tutorials/14-Timestepping-with-physical-units.ipynb +++ b/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", diff --git a/docs/beginner/tutorials/15-Thermal-convection-with-units.ipynb b/docs/beginner/tutorials/15-Thermal-convection-with-units.ipynb index 762c45a68..2efcd6b14 100644 --- a/docs/beginner/tutorials/15-Thermal-convection-with-units.ipynb +++ b/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", diff --git a/docs/beginner/tutorials/7-Timestepping-simple.ipynb b/docs/beginner/tutorials/7-Timestepping-simple.ipynb index 85d5094ef..fd90fe3bb 100644 --- a/docs/beginner/tutorials/7-Timestepping-simple.ipynb +++ b/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", diff --git a/docs/beginner/tutorials/8-Timestepping-coupled.ipynb b/docs/beginner/tutorials/8-Timestepping-coupled.ipynb index 45ed4cdb2..045cb6ad9 100644 --- a/docs/beginner/tutorials/8-Timestepping-coupled.ipynb +++ b/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()" ] }, { diff --git a/docs/beginner/tutorials/9-Unsteady_Flow.ipynb b/docs/beginner/tutorials/9-Unsteady_Flow.ipynb index 9f4aba69c..d5ba45a8b 100644 --- a/docs/beginner/tutorials/9-Unsteady_Flow.ipynb +++ b/docs/beginner/tutorials/9-Unsteady_Flow.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "id": "f7a4cbb2-6265-48bd-a646-0e1df6c569de", "metadata": {}, - "source": "# Notebook 7: Unsteady Flow\n\n[![Run this notebook on Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/underworldcode/uw3-binder-launcher/development?labpath=docs%2Fbeginner%2Ftutorials%2F9-Unsteady_Flow.ipynb)\n\n\n
\n\n![](media/CompositeImage.png)\n\n_Flow in a pipe with inflow at the left boundary\n after 50, 100, 150 timesteps (top to bottom) showing the\n progression of the impulsive initial condition. For details,\n see the notebook code._\n\n
\n\nWe'll look at tracking an unsteady flow using a swarm of particle flow-tracers. In this case, the flow is unsteady because we solve the Navier-Stokes equation (that is, the flow has inertia) and we impose an impulsive, initial boundary velocity. \n\nTo begin with, the set up follows the same path as all previous notebooks:\n - Create a mesh\n - Add some variables\n - Create the solver we need (`NavierStokes` this time)\n - Add boundary conditions and constitutive properties.\n\nWe also add a projection solver to compute the vorticity of the flow as we did in Notebook 4 when we needed to compute a heat-flux (thermal gradient) term.\n\nTo track the time evolution of the flow, we introduce a \"passive\" particle\nswarm. Passive, here, refers to the fact that the flow is not changed by the \npresence of the marker particles. \n\nIn the time-loop we have to update the particle locations and we keep this\nas an explicity operation, in general, because it provide the opportunity for\nyou to make changes or perform analyses. In this case, we are adding new particles near the inflow to track the flow.\n\nTo learn more about flow goverened by the Navier-Stokes equation, it may be he helpful to read an elementary fluid dynamics textbook and reproduce some of the simple \"toy\" examples. For example, Acheson, 1990. \n" + "source": "# Notebook 7: Unsteady Flow\n\n[![Run this notebook on Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/underworldcode/uw3-binder-launcher/development?labpath=docs%2Fbeginner%2Ftutorials%2F9-Unsteady_Flow.ipynb)\n\n\n
\n\n![](media/CompositeImage.png)\n\n_Flow in a pipe with inflow at the left boundary\n after 50, 100, 150 timesteps (top to bottom) showing the\n progression of the impulsive initial condition. For details,\n see the notebook code._\n\n
\n\nWe'll look at tracking an unsteady flow using a swarm of particle flow-tracers. In this case, the flow is unsteady because we solve the Navier-Stokes equation (that is, the flow has inertia) and we impose an impulsive, initial boundary velocity. \n\nTo begin with, the set up follows the same path as all previous notebooks:\n - Create a mesh\n - Add some variables\n - Create the solver we need (`NavierStokesSLCN` this time, the semi-Lagrangian Navier-Stokes solver)\n - Add boundary conditions and constitutive properties.\n\nWe also add a projection solver to compute the vorticity of the flow as we did in Notebook 4 when we needed to compute a heat-flux (thermal gradient) term.\n\nTo track the time evolution of the flow, we introduce a \"passive\" particle\nswarm. Passive, here, refers to the fact that the flow is not changed by the \npresence of the marker particles. \n\nIn the time-loop we have to update the particle locations and we keep this\nas an explicity operation, in general, because it provide the opportunity for\nyou to make changes or perform analyses. In this case, we are adding new particles near the inflow to track the flow.\n\nTo learn more about flow goverened by the Navier-Stokes equation, it may be he helpful to read an elementary fluid dynamics textbook and reproduce some of the simple \"toy\" examples. For example, Acheson, 1990. \n" }, { "cell_type": "code", @@ -83,7 +83,7 @@ "metadata": {}, "outputs": [], "source": [ - "navier_stokes = uw.systems.NavierStokes(\n", + "navier_stokes = uw.systems.NavierStokesSLCN(\n", " mesh,\n", " velocityField=v_soln,\n", " pressureField=p_soln,\n", @@ -126,7 +126,7 @@ "metadata": {}, "outputs": [], "source": [ - "uw.systems.NavierStokes.view()" + "uw.systems.NavierStokesSLCN.view()" ] }, { diff --git a/docs/developer/design/SOLVER_UNIFICATION_DESIGN.md b/docs/developer/design/SOLVER_UNIFICATION_DESIGN.md index 2d6f65dcb..b5e2ab3b5 100644 --- a/docs/developer/design/SOLVER_UNIFICATION_DESIGN.md +++ b/docs/developer/design/SOLVER_UNIFICATION_DESIGN.md @@ -23,7 +23,7 @@ creates it lazily. |--------|-----------------|-------------------|-------------------| | `Stokes` | — | — | Viscous, VP | | `VE_Stokes` | — | SemiLagrangian (stress history) | VEP | -| `NavierStokes` | SemiLagrangian (velocity) | SemiLagrangian (AM flux) | Viscous, VP | +| `NavierStokesSLCN` (was `NavierStokes`; the generic name is now the composing Eulerian solver, 2026-09) | SemiLagrangian (velocity) | SemiLagrangian (AM flux) | Viscous, VP | | `VE_NavierStokes` | does not exist | — | — | Problem: user must choose the correct solver class based on the constitutive model. diff --git a/docs/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md new file mode 100644 index 000000000..391464311 --- /dev/null +++ b/docs/developer/design/eulerian-supg-transport.md @@ -0,0 +1,774 @@ +# Eulerian SUPG transport: design and measurements + +**Status**: implemented on `feature/eulerian-supg-transport` (2026-09-02), static mesh. + +**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 + +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: + +| family | time derivative | spatial operator | +|---|---|---| +| 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 +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. +- **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 + 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, and vector or tensor unknowns: the solver is scalar, where the + semi-Lagrangian trace-back carries vectors and tensors through the same machinery. + +## 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 | + +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. + +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. 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). +- **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. + +**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). 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 +`order=2`; below Courant 1 with diffusion, `order=3`. Backward Euler is not a +sensible choice 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. + +### 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); matched, with the shipped defaults, 3.51 / 0.48 s (Schwarz, +5 iterations) against 3.62 / 0.54 s (multigrid, one cycle). + +## Navier-Stokes with SUPG momentum transport + +`uw.systems.NavierStokes` (`systems/navier_stokes_eulerian.py`) is the vector +form of the scalar solver on the Stokes saddle-point class: the momentum advection +is assembled implicitly and the streamline term stabilises it. The residual is + +$$ +\mathbf{f}_0 = \rho\,\big(\dot{\mathbf{u}} + \textstyle\sum_k w_k (\mathbf{a}_k\cdot\nabla)\mathbf{u}^{(k)}\big) - \mathbf{f}, +\qquad +\mathbf{F}_1 = \textstyle\sum_k w_k\,\boldsymbol{\tau}(\mathbf{u}^{(k)}) - p_\mathrm{mech}\mathbf{I} + \tau_s\,\mathbf{R}\otimes\mathbf{a}, +$$ + +with $\mathbf{R} = \mathbf{f}_0 + \nabla p$ the strong residual the SUPG term sees, +$\mathbf{a}$ the advecting velocity at the new level and $\mathbf{a}_k = \mathbf{u}^{(k)}$ +at the stored ones, $w_k$ the weights of the spatial operator (Adams-Moulton at +order 1, all on n+1 for BDF2), and $\tau_s$ the scalar formula with $\nu = \eta/\rho$. +The pressure equation is the Stokes constraint; Taylor-Hood needs no pressure +stabilisation. Since 2026-09-06 the term carries the cell-Péclet weight +$Pe^2/(Pe^2 + Pe_c^2)$ with $Pe_c = 4$ by default (Louis: "the code is still 100% local, +so we should probably just switch to this strategy right away"), measured in "The weight +by cell Péclet number" below; every table before that subsection was made with the +uniform weight (`peclet_weight=0`), and the Pe_c = 4 column there gives the change. The +scalar transport solver carries the same weight (its convection rows are in that +subsection). Decisions, and what they rest on: + +- **No stress history.** The semi-Lagrangian solver carries a stress history + because its Crank-Nicolson viscous term needs the old flux at the departure + points. On the grid the old flux is needed where it was formed, so + $\boldsymbol{\tau}(\mathbf{u}^n) = 2\eta\,\dot\varepsilon(\mathbf{u}^n)$ is rebuilt + from the stored velocity level with the current effective viscosity (exact for + a constant viscosity; use BDF2 with a strain-rate dependent one). Pressure has + no history. A history-dependent stress is the constitutive model's business. +- **The advecting velocity is pluggable.** `advection="extrapolated"` (default), + $\mathbf{a} = 2\mathbf{u}^n - \mathbf{u}^{n-1}$, makes each step one linear Oseen + solve through the Stokes fieldsplit, with a second-order lag and no explicit + stability limit; `picard_iterations=n` re-solves with the latest iterate for the + fully implicit fixed point without a tangent; `advection="implicit"` puts + $\mathbf{u}^{n+1}$ in the term and the SNES takes Newton steps with the symbolic + tangent. At a steady state all three coincide, which the Kovasznay rows below + confirm to every digit; the cylinder wake is where they differ. +- **The pressure gradient belongs in the SUPG residual.** Without it $\mathbf{R}$ is + O(1) at the exact solution and the stabilisation injects an O($\tau$) error: + Kovasznay at 1/16 read 1.9e-3 against 6.6e-4 with it (three times). The viscous + term needs second derivatives the kernels do not see and is the remaining + inconsistency, O($h^2$) for P2 velocity in diffusion-limited cells; a recovered + Laplacian would close it and is deferred. +- **`mesh.cell_size()` was partition-dependent** (#687): the kd-tree radius picked + the nearest centroid among the rank's own cells, so $\tau$ differed across a + partition seam (two-rank Kovasznay error 5e-4 off serial, 1e-15 with a + constant $h$), and after a deform it read stale vertex coordinates against fresh + centroids. The field now reports each cell's own radius from the DM's + coordinates; the kd-tree radii still feed `get_min_radius`. +- **`solve()` builds through `_build`**, one Newton iteration per step at matched + tolerances, as for the scalar solver. + +### Kovasznay flow (Re 40) + +Exact steady Navier-Stokes on $[-0.5, 1] \times [-0.5, 0.5]$, Dirichlet velocity +from the exact solution, P2-P1, 40 steps at Courant 1 from the exact solution +(or 80 from rest); relative $L_2$ velocity error at the end +(`~/+Simulations/navier_stokes_supg/kovasznay/`). + +| h | SUPG, CN | Galerkin, CN | SUPG, BDF2 | SLCN | s/step SUPG / SLCN | +|---|---|---|---|---|---| +| 1/16 | 6.6e-4 | 1.1e-4 | 6.3e-4 | 5.8e-3 | 0.27 / 1.9 | +| 1/32 | 2.6e-4 | 1.6e-5 | 2.5e-4 | 2.9e-3 | 1.5 / 3.7 | +| 1/64 | 6.7e-5 | | | | 6.8 / | + +Galerkin converges at third order here (the interpolation error), SUPG at 1.4 +rising to 2.0 (the missing viscous term), SLCN at first order. At Re 40 the +element Reynolds number is below three on every mesh and the stabilisation is not +needed; it costs a factor of six to sixteen against Galerkin and is still nine +times more accurate than the semi-Lagrangian scheme at seven times less cost per +step. Newton (`advection="implicit"`), two Picard passes, Courant 4, and the +from-rest starts all reach the same steady state (6.60e-4 at 1/16); BDF2 sits at +an exact fixed point (step change 0) where Crank-Nicolson keeps a 5e-5 flicker. +Two ranks reproduce the serial error to 1e-7 (test_1078). + +### Lid-driven cavity + +Unit square, no-slip walls, unit lid (singular at the corners), P2-P1 on an +unstructured mesh, marched from rest; centreline extrema (u on x = 0.5, v on +y = 0.5) against Ghia, Ghia and Shin (1982). `~/+Simulations/navier_stokes_supg/cavity/`. + +| Re | h | scheme | Courant | Picard | u_min | v_max | v_min | steps | s/step | +|---|---|---|---|---|---|---|---|---|---| +| 100 | Ghia | | | | -0.2109 | 0.1753 | -0.2453 | | | +| 100 | 1/32 | SUPG | 1 | 0 | -0.2025 | 0.1710 | -0.2437 | 543 (fixed point) | 0.73 | +| 100 | 1/32 | SLCN | 1 | | -0.1977 | 0.1695 | -0.2365 | 1000 (still moving) | 2.9 | +| 400 | Ghia | | | | -0.3273 | 0.3020 | -0.4499 | | | +| 400 | 1/48 | SUPG | 2 | 0 | -0.3076 | 0.2832 | -0.4288 | 1500 | 2.3 | +| 400 | 1/48 | SUPG | 2 | 1 | -0.3076 | 0.2831 | -0.4288 | 976 (fixed point) | 2.8 | +| 400 | 1/48 | SUPG | 1 | 0 | -0.3075 | 0.2828 | -0.4284 | 1500 (change 3e-6) | 2.1 | +| 1000 | Ghia | | | | -0.3829 | 0.3709 | -0.5155 | | | +| 1000 | 1/64, 3-level FMG | SUPG | 1 | 0 | -0.3413 | 0.0613 | -0.4687 | 1200 (t = 19, still moving) | 3.3 (np 4) | +| 1000 | 1/64, 3-level FMG | Galerkin | 1 | 0 | -0.1437 | 0.0695 | -0.2031 | 300 (t = 4.7) | 3.7 (np 4) | +| 1000 | 1/64, 3-level FMG | SUPG | 2 | 1 | -0.3620 | 0.3491 | -0.4940 | 2281 (t = 71, steady) | 3.3 (np 4) | + +At Re 100 SUPG is within 4% of Ghia on every extremum on a 1/32 mesh and +reaches an exact fixed point; SLCN on the same mesh sits a little further out and +has not settled after 1000 steps at four times the cost. At Re 400 on 1/48 both +runs give the same extrema, 5 to 6% below Ghia (the mesh, not the scheme: the +extrema are steady to four digits), but the extrapolated step alone never +becomes stationary: the max-norm change per step grows to 0.1 and saturates, an +alternating mode of the lagged coefficient fed by the lid singularity while the +interior sits still. One Picard pass removes it (step change exactly zero) for +20% more per step, and so does Courant 1 without any pass. That is the regime the +Picard option was built for: Courant 2 with an element Reynolds number near eight. + +At Re 1000 (element Reynolds number 16) on a 1/64 mesh built with a two-level +refinement so the velocity block runs geometric multigrid, the extrapolated step +takes one Newton and one Krylov iteration per step at 3.3 s on four ranks, and +the Galerkin form runs just as stably for its 300 steps: neither oscillates on +this mesh. The Courant-2, one-Picard run reaches the steady tolerance at step 2281 +(t = 71) with the three extrema at 94 to 96% of Ghia and their positions within 0.01 +(u_min at y 0.175, v_max at x 0.163, v_min at x 0.907), the same shortfall as Re 400 +on 1/48, and the flow (primary vortex, both bottom-corner eddies) as the reference +shows it. The v_max of 0.05 to 0.07 the earlier rows print is the driver reading rank +0's own `evaluate` on four ranks (the left-wall upflow sits in another partition); the +driver now reduces the extrema across ranks, and the row above is a serial +re-evaluation of the final checkpoint. Two earlier four-rank attempts stalled at +their first logged step, which was the driver calling the collective centreline +evaluation on rank 0 only, and a third was killed by the hang watchdog on a rank +that never prints; none of those said anything about the solver. + +### Cylinder wake (DFG 2D-2, Re 100) + +Channel 2.2 by 0.41, cylinder of radius 0.05 at (0.2, 0.2), parabolic inflow with +mean velocity 1, $\nu = 10^{-3}$; mesh 1/20 in the channel and 1/80 on the +cylinder, P2-P1, Courant 1 on the cylinder cells (dt 0.0083), twelve time units +from the parabolic profile; drag and lift from the traction integral on the +cylinder, the Strouhal number from the lift zero crossings over the last three +units. Reference (Schaefer and Turek 1996): $C_D$ max 3.22 to 3.24, $C_L$ max +0.99 to 1.01, St 0.295 to 0.305, $\Delta p$ 2.46 to 2.50. +`~/+Simulations/navier_stokes_supg/cylinder/`. + +Every drag value in the first table below is the PRESSURE drag only: the boundary +integral of the traction dropped the viscous part, because the viscosity is a runtime +expression and the integral kernels read every expression as zero (#695, found through +this benchmark and fixed on this branch). The lift and the Strouhal number were never +affected (the lift is pressure-dominated and the frequency does not go through an integral). + +| scheme | advecting velocity | St | $C_L$ max | $C_D$ max (pressure part only, #695) | $\Delta p$ | s/step | +|---|---|---|---|---|---|---| +| SUPG | extrapolated | 0.298 | 0.82 | 2.33 | 2.41 | 0.73 | +| SUPG | extrapolated + 1 Picard pass | 0.296 | 0.75 | 2.30 | 2.39 | 1.02 | +| SUPG | implicit (Newton) | 0.295 | 0.76 | 2.30 | 2.39 | 0.99 | +| SLCN | (trace-back) | 0.259 | 0.68 | 2.72 | 2.30 | 2.36 | +| SUPG, mesh 1/40 and 1/160, np 4 | extrapolated | 0.304 | 0.89 | 2.48 | | 1.0 (np 4) | + +The shedding frequency and the pressure difference are on the reference at both meshes. +The two fully implicit forms agree with each other to three digits, and the extrapolated +step differs from them by 1% in frequency and 8% on the lift peak: at Courant 1 the lag is +visible on a time-dependent wake but small, and a single Picard pass, or Newton, removes +it at 40% more per step. The semi-Lagrangian solver on the same mesh and step has the +shedding 13% too slow (St 0.259) at three times the cost; the frequency is the quantity +the time integration owns, and there the Eulerian scheme is the accurate one. + +**The drag deficit was a measurement, not the scheme.** The drag read 28% low on the +1/20 mesh and 23% on 1/40, and did not move with the SUPG weights: with the velocity +block solved by LU (the GAMG fallback on this gmsh mesh spins at weak stabilisation, so +the "Galerkin cannot run" of the first attempt was the preconditioner, not the +discretisation), the SUPG weight from 1 to 0 and the tau weights over a factor of four +moved the peak by 2.6%, with the reaction-form drag (the momentum residual integrated +against a hat function on the cylinder nodes) 6% above the traction integral throughout. +The log then showed the total traction drag equal to its pressure part to four digits. +With the integrals fixed, the channel mesh held at 1/20 and only the cylinder cells +refined through gmsh (Louis's prescription: refine the cylinder, keep the step), all at +dt 0.0083 (Courant 1 on the 1/80 cylinder cells, 8 on the 1/640 ones), velocity block LU, +serial; the last row is the whole mesh at 1/40 on four ranks at its own Courant-1 step: + +| cylinder cells | Courant at the cylinder | $C_D$ max traction / reaction | $C_L$ max | $\Delta p$ | St | steps | s/step | +|---|---|---|---|---|---|---|---| +| 1/80 (SUPG) | 1 | 3.046 / 3.115 | 0.897 | 2.41 | 0.298 | 1440 | 0.38 | +| 1/80 (Galerkin) | 1 | 3.098 / 3.168 | 0.909 | 2.41 | 0.295 | 1440 | 0.38 | +| 1/160 | 2 | 3.134 / 3.155 | 0.866 | 2.43 | 0.300 | 1440 | 0.64 | +| 1/160, dt 0.0042 | 1 | 3.131 / 3.153 | 0.841 | 2.43 | 0.298 | 2880 | 0.47 | +| 1/320 | 4 | 3.198 / 3.204 | 0.969 | 2.48 | 0.300 | 1440 | 1.0 | +| 1/640 | 8 | 3.237 / 3.252 | 1.067 | 2.49 | 0.299 | 1440 | 1.6 | +| 1/640, one Picard pass | 8 | 3.218 / 3.220 | 1.018 | 2.48 | 0.296 | 1440 | 1.6 | +| whole mesh 1/40 (cylinder 1/160), np 4 | 1 | 3.182 / 3.204 | 0.979 | | 0.304 | 2880 | 1.5 (np 4) | +| FMG: base 1/10 refined once (cylinder 1/80) | 1 | 3.108 / 3.156 | 0.976 | 2.49 | 0.303 | 1440 | 1.6 | +| FMG: base 1/10 refined once (cylinder 1/160) | 2 | 3.181 / 3.202 | 1.006 | 2.49 | 0.302 | 1440 | 2.4 | +| reference | | 3.22 to 3.24 | 0.99 to 1.01 | 2.46 to 2.50 | 0.295 to 0.305 | | | + +The drag, the pressure difference and the frequency converge onto the reference bands as +the cylinder cells alone are refined, the two force measurements close on each other (6% +apart with the wall shear in one cell, 0.2% at 1/320), and the time step does not enter: +the 1/160 rows at Courant 1 and 2 give the same drag to three digits. Refining the whole +mesh to 1/40 (four ranks, twice the steps) buys less than the 1/320 cylinder cells do on +one core with the 1/20 channel. The lift peak converges from below and overshoots the band +by 6% at 1/640 (Courant 8 on those cells); one Picard pass on the same mesh brings it to +1.018 with the drag at 3.218 and the two force measurements 0.1% apart, so the overshoot is +the extrapolated advecting velocity's lag at that local Courant number, not the mesh. That +is the regime the Picard option exists for, and at Courant 8 on the cells that set the +forces it is worth its 40% per step. Earlier reads of this benchmark (drag "23 to 28% low, not +closing with the mesh, not moving with tau") were the missing viscous traction (#695): the +SUPG weight from 1 to 0 and the tau weights over a factor of four move the drag peak by +2.6%, and the Galerkin form that "could not run" was the GAMG fallback spinning inside the +Schur complement at weak stabilisation (native stack), not the discretisation. + +FMG on this gmsh mesh: building the base mesh at 1/10 and refining once through the circle +callback (`-uw_refinement 1`, the callback snaps the new vertices to the circle) gives the +velocity block its geometric hierarchy, one Krylov iteration per Newton step and no +fallback (the two FMG rows). The refined mesh also gives a better lift and pressure +difference than the directly meshed 1/20 channel with the same cylinder cell: at 1/160 on +the cylinder every quantity but the drag (1% low) is inside the reference band. The +per-step times were taken with twelve cores busy and are not a like-for-like comparison +with LU. LU on the velocity block is serial-only: on more than one rank PETSc's native +factorisation has no parallel path and the run dies in the first solve, so the multigrid +hierarchy is the parallel route on this mesh. + +The same cylinder-only refinement through FMG (base 1/10 refined once, channel 1/20, +fixed dt 0.0083, no LU), the route that scales: + +| cylinder cells | Courant there | advecting velocity | $C_D$ max traction / reaction | $C_L$ max | $\Delta p$ | St | s/step | +|---|---|---|---|---|---|---|---| +| 1/320 | 4 | extrapolated | 3.208 / 3.215 | 1.004 | 2.48 | 0.300 | 4.9 | +| 1/320 | 4 | one Picard pass | 3.187 / 3.194 | 0.954 | 2.47 | 0.297 | 6.1 | +| 1/640 | 8 | one Picard pass | 3.204 / 3.206 | 0.969 | 2.47 | 0.297 | 10.9 | +| 1/640, np 4 | 8 | one Picard pass | 3.205 / 3.206 | 0.969 | | 0.297 | 6.2 (np 4) | +| 1/640 | 8 | Newton | 3.204 / 3.205 | 0.969 | 2.47 | 0.296 | 7.4 | +| 1/640 | 8 | one Picard pass, BDF2 | 3.196 / 3.198 | 0.939 | 2.47 | 0.295 | 8.6 | +| reference | | | 3.22 to 3.24 | 0.99 to 1.01 | 2.46 to 2.50 | 0.295 to 0.305 | | + +(Times with the machine shared by five runs.) The drag and the pressure difference sit +within 1% of the bands with the two force measurements 0.05% apart at 1/640; the +frequency is in band throughout. The lift is the sensitive quantity: the extrapolated step +reads 1.004 at Courant 4 on the cylinder cells and 1.067 at Courant 8 on the unrefined +mesh, the implicit forms 0.954 to 0.969, and BDF2 0.939, so at these local Courant numbers +the extrapolation's lag and BDF2's damping each move the lift peak by 3 to 5% and the +Crank-Nicolson implicit forms are the ones to compare with the reference. One Picard +pass and Newton agree to three digits at 1/640 and Newton is the cheaper of the two. +Serial and four ranks agree to four digits (3.204 / 3.205, 0.9692 / 0.9689), the +partition independence the assembled operator should give, at 1.8x on four ranks with +the machine loaded. `figures/13_cylinder_Re100_supg_c320_picard_tracers.mp4` is the wake +at the 1/320-cell, one-Picard setup with tracers released in the central band. + +Parallel tracers (#693) work with the empty-rank guard from #680: the two further +failures reported there were the driver's (an advection before the first release, on a +swarm that had never been populated and so carries the DMSwarm local size of −1, which +fails in serial in the same way; and a timing variable shadowed by a rank-local array). + +### Vortex decay (Taylor-Green) + +The exact unsteady solution on $[0,\pi]^2$, $\mathbf{u} = (-\sin x\cos y,\ \cos x\sin y)\,e^{-2\nu t}$, +$p = \tfrac14(\cos 2x + \cos 2y)\,e^{-4\nu t}$, has no normal flow and no tangential +stress on the walls, so free-slip walls (the normal component fixed) are exact and carry +no time dependence. Relative $L_2$ velocity error at $t = 1$ against the exact solution, +$\nu = 0.01$, P2-P1 on a regular simplex mesh, velocity block by LU, from the exact +initial state (`~/+Simulations/navier_stokes_supg/vortex_decay/`, `scripts/taylor_green.py`). +The interpolation error of the exact field is 1.7e-5 on the 1/32 mesh and 2.2e-6 on 1/64. + +| dt (mesh 1/32) | SUPG, CN | Galerkin, CN | SUPG, BDF2 | Galerkin, BDF2 | +|---|---|---|---|---| +| 0.2 | 4.3e-4 | 6.8e-5 | 2.5e-4 | 5.3e-5 | +| 0.1 | 2.1e-4 | 4.8e-5 | 2.0e-4 | 4.9e-5 | +| 0.05 | 1.7e-4 | 4.9e-5 | 1.4e-4 | 4.9e-5 | +| 0.025 | 1.2e-4 | 4.9e-5 | 9.2e-5 | 4.9e-5 | +| 0.0125 | 7.8e-5 | 4.9e-5 | 6.5e-5 | 4.9e-5 | + +| h (dt 0.0125) | SUPG | Galerkin | interpolation | +|---|---|---|---| +| 1/8 | 9.0e-3 | 9.3e-3 | 1.1e-3 | +| 1/16 | 6.7e-4 | 6.6e-4 | 1.4e-4 | +| 1/32 | 7.8e-5 | 4.9e-5 | 1.7e-5 | +| 1/64 | 1.6e-5 | 4.0e-6 | 2.2e-6 | + +The Galerkin form is spatially limited at every time step in the table: its error is the +same at dt 0.2 as at dt 0.0125 and converges at third order in $h$, three times the +interpolation error. The time integration is not what limits this problem, because the +pattern is steady and only the amplitude decays, and Crank-Nicolson integrates +$e^{-2\nu t}$ with $2\nu\,\Delta t \le 0.004$ almost exactly. What the SUPG column measures +is the stabilisation's consistency error, and it scales with $\tau_s$: halving the time +step raises the transient term in $\tau_s$ and lowers the error by 1.5 to 1.8 until the +advective term takes over, and on the 1/64 mesh the floor is four times the Galerkin +error. The advecting-velocity choices coincide to four digits (1.723e-4 at dt 0.05 for +extrapolated, three Picard passes and Newton). Across the viscosity range at dt 0.025 on +1/32 the SUPG error is 4.2e-4 at $\nu = 1$ (the energy has decayed to 1.8%), 1.9e-5 at +0.1, 1.2e-4 at 0.01 and 5.1e-4 at 0.001 (element Reynolds number 100). The kinetic energy +ratio follows $e^{-4\nu t}$ to 1e-6 with free-slip walls for both forms. + +Imposing the exact velocity on the walls instead (`-uw_bc dirichlet`, the time a runtime +expression in the condition) gives 1.13e-4 at dt 0.025 on 1/32, the free-slip value. It +first gave 4.7e-3 on every mesh and at every time step, with the decay 10% too slow, and +freezing the time deliberately reproduced that number to four digits: the time expression +had been created at the value zero, and sympy's automatic evaluation, reading the +expression's `is_zero` assumption from its value, had evaluated $e^{-2\nu t}$ out of the +boundary formula before the JIT saw it (issue #696, since fixed: a `UWexpression` no longer +reports `is_zero`, `is_positive` or `is_negative` from its current value, so sympy cannot fold on +them; `tests/test_0503` carries the `exp(c)` control). + +### The recovered viscous term: measured and withdrawn + +The SUPG column above is the stabilisation's consistency error: the strong residual the +term weights lacks $\nabla\cdot\boldsymbol{\sigma}$ (second derivatives the kernels do +not see). The Péclet turn-down of $\tau_s$ does not remove it: at low Péclet number +$\tau_s \to h^2/(4\nu)$ while the missing term is $\nu\nabla^2\mathbf{u}$, and the product +is $O(h^2)$ with no $\nu$ in it. Three ways of supplying the term were built and measured +(velocity error at $t = 1$, dt 0.0125; Kovasznay at Re 40; the cylinder on the 1/20 mesh): + +| case | SUPG | Galerkin | projected stress | balance form (Louis) | balance, smoothed L = 0.01 to 0.2 | +|---|---|---|---|---|---| +| vortex 1/32 | 7.8e-5 | 4.9e-5 | 7.8e-5 | 4.9e-5 | 7.8e-5 | +| vortex 1/64 | 1.6e-5 | 4.0e-6 | diverged | 4.1e-6 | | +| Kovasznay 1/16 | 6.6e-4 | 1.1e-4 | | 1.1e-4 | | +| Kovasznay 1/32 | 2.6e-4 | 1.6e-5 | diverged | 1.6e-5 | | +| cylinder $C_D$ / $C_L$ max | 3.046 / 0.897 | 3.098 / 0.909 | | diverged (step 25 to 50) | 3.04 / 0.85 to 0.87 | + +- **Projected stress**: the deviatoric stress of the advecting velocity fitted to a + continuous P2 tensor and differentiated. No change at 1/32, unstable finer: a + differentiated fit to a discontinuous strain rate is not a Laplacian. +- **Balance form**: $\nabla\cdot\boldsymbol{\sigma}^n = \rho(D\mathbf{u}/Dt)^n + \nabla p^n + - \mathbf{f}$ from the stored levels and a stored pressure, so the residual is the + increment of the out-of-balance force between levels. It returns the Galerkin accuracy + to two digits on every resolved case, and it does so because it is a tautology: for any + slowly varying discrete solution the residual it builds is zero, so it does not recover + the viscous term, it switches the stabilisation off. Where the stabilisation is needed it + fails the same way, with the lagged residual fed back as a source (cylinder, drag 7% high + at step 25, linear solve diverged before step 50, on a mesh where plain Galerkin runs). +- **Balance form projected with a smoothing length** (screened Poisson, 0.01 to 0.2 on the + vortex, one to two cylinder cells on the cylinder): the projected term matches the exact + $\nu\nabla^2\mathbf{u}$ to a few per cent and the SUPG error does not move at any length, + while the cylinder stays stable and within 1% of plain SUPG on drag. A continuous + recovery of the viscous term, however accurate, does not touch the error. + +What the three say together: the consistency error on resolved P2 flow is not the smooth +part of the missing viscous term. It is the pointwise, element-wise residual of the +discrete solution (the piecewise-constant P1 pressure gradient and the second derivatives +of the P2 velocity, both O(h) pointwise) that $\tau_s\,\mathbf{a}\cdot\nabla\mathbf{w}$ +integrates; only a term that cancels it pointwise removes it, and that term cancels the +stabilisation with it. The remedy the measurements support is not a recovered Laplacian +but the weight: where the cell Péclet number is small the term is not needed and costs a +fixed multiple of the Galerkin error, second order in $h$ (`supg_weight`, or the Galerkin +form; Kovasznay's recommendation stands). A Péclet-dependent weight is the design +question that remains. The options were removed from the solver after the measurement (a +knob that quietly disables the stabilisation should not ship); the drivers' `-uw_recovered` +switches went with them and the runs are in the study directory (`rec_*`, `bal_*`, `sm_*`). + +### The shape of tau (`tau_shape`) + +The inverse-sum $\tau_s$ is above the optimal 1-D curve at cell Péclet numbers of order +1 to 10 (Louis: the shape of the correction was always the debated trade-off between +accuracy and cost). Two further shapes are selectable, each combined with the same +transient cap $[(C_t c_0/\Delta t)^2 + \tau^{-2}]^{-1/2}$: Brooks-Hughes, +$\tau = (h/2|a|)(\coth Pe - 1/Pe)$, and the doubly asymptotic $(h/2|a|)\min(Pe/3, 1)$, +$Pe = |a|h/2\nu$. Same cases as above (cell Péclet number in brackets): + +| case | inverse sum | Brooks-Hughes | doubly asymptotic | Galerkin | +|---|---|---|---|---| +| vortex 1/32, dt 0.0125 (Pe 5) | 7.8e-5 | 7.8e-5 | 7.8e-5 | 4.9e-5 | +| vortex 1/32, dt 0.1 | 2.1e-4 | 1.7e-4 | 1.8e-4 | 4.8e-5 | +| vortex 1/64, dt 0.0125 (Pe 2.5) | 1.6e-5 | 1.4e-5 | 1.4e-5 | 4.0e-6 | +| Kovasznay 1/16 (Pe 1 to 3) | 6.6e-4 | 4.1e-4 | 4.8e-4 | 1.1e-4 | +| Kovasznay 1/32 | 2.6e-4 | 1.2e-4 | 1.2e-4 | 1.6e-5 | +| cylinder $C_D$ / $C_L$ max (Pe 10 at the wall) | 3.046 / 0.897 | 3.057 / 0.903 | 3.046 / 0.896 | 3.098 / 0.909 | + +The shape matters where the cell Péclet number is near one: on Kovasznay the optimal +form halves the error (1.6 to 2.3 times) and on the 1/64 vortex it takes 15% off; where +the transient term caps $\tau_s$ (the 1/32 vortex at dt 0.0125) or advection dominates +(the cylinder, where all three shapes are $h/2|a|$) nothing moves. What remains after the +optimal shape is still seven times the Galerkin error on Kovasznay at 1/32: the shape +reduces the excess of $\tau_s$ over the 1-D optimum, it cannot remove the $O(h^2)$ product +of $\tau_s$ and the missing viscous term. The two 1-D shapes are exposed as options; the +inverse sum stays the default (smooth, no per-cell Péclet evaluation), and the weight, +by cell Péclet number, remains the lever that reaches the Galerkin value. + +### The weight by cell Péclet number (`peclet_weight`) + +The other side of the same lever: leave $\tau_s$ alone and multiply the term by +$w = Pe^2/(Pe^2 + Pe_c^2)$, $Pe = |a|h/2\nu$, so it is off where the cell is +diffusion-dominated and full where advection dominates. Same cases, three thresholds: + +| case (cell Péclet) | SUPG | $Pe_c = 2$ | $Pe_c = 4$ | $Pe_c = 8$ | Galerkin | +|---|---|---|---|---|---| +| vortex 1/32 (Pe 5) | 7.8e-5 | 6.3e-5 | 5.3e-5 | 4.9e-5 | 4.9e-5 | +| vortex 1/64 (Pe 2.5) | 1.6e-5 | 8.6e-6 | 5.0e-6 | 4.1e-6 | 4.0e-6 | +| Kovasznay 1/16 (Pe 1 to 3) | 6.6e-4 | 3.6e-4 | 1.6e-4 | 1.1e-4 | 1.1e-4 | +| Kovasznay 1/32 | 2.6e-4 | 5.2e-5 | 2.1e-5 | 1.6e-5 | 1.6e-5 | +| cylinder $C_D$ / $C_L$ max (Pe 10 wall, 37 channel) | 3.046 / 0.897 | 3.061 / 0.908 | 3.080 / 0.919 | 3.094 / 0.917 | 3.098 / 0.909 | +| cylinder St | 0.298 | 0.297 | 0.296 | 0.296 | 0.295 | + +This is the measurement that closes the trade-off. At $Pe_c = 8$ every resolved case +sits on the Galerkin value and the cylinder is still stable and within 0.2% of Galerkin +on drag with the weight at 0.6 on the wall cells and 0.95 in the channel; at $Pe_c = 4$ +the resolved cases are within 1.3 times Galerkin and the wall cells keep 86% of the +term. The weight does what neither the recovered viscous term nor the shape of $\tau_s$ +could: it removes the cost of stabilisation where the 1-D analysis says none is needed +and leaves it where it is. The default stays at zero (uniform weight) so that the +recorded benchmarks and the test references do not move; $Pe_c = 4$ is the recommended +setting for resolved or mixed problems. Louis's ruling (2026-09-06, the code being +unreleased): $Pe_c = 4$ is the default of both solvers. The scalar solver on the convection +benchmarks with it (`~/+Simulations/supg_vs_slcn_657/convection_benchmarks/`, runs +`*_pew4`; Blankenbach 1a reference Nu 4.884, Vrms 42.865): + +| case | uniform weight: Vrms / Nu cold / Nu mid | $Pe_c = 4$: Vrms / Nu cold / Nu mid | transport s/step | +|---|---|---|---| +| box, Ra 1e4, 1/32 | 42.790 / 4.913 / 4.872 | 42.868 / 4.920 / 4.884 | 0.066 / 0.070 | +| annulus, Ra 1e4, 0.03 | 38.39 / 2.514 / 2.500 | 38.61 / 2.525 / 2.514 | 0.179 / 0.178 | + +The box lands on the reference to four digits in Vrms and in the mid-plane Nusselt +number (the cells there sit at a Péclet number near one, where the term was costing +accuracy); the annulus moves 0.6% in the same direction; the cost does not move. The +parallel test's serial reference for the Navier-Stokes solver (Kovasznay at 1/8) goes +from 3.83e-3 to 1.42e-3; the pure-advection references are unchanged (weight 1). + +### A defect in the integrals (#695) + +The first error metric of this benchmark, an integral of $|\mathbf{v} - \mathbf{u}(t)|^2$ +with the time as a runtime expression, returned $1 - e^{-2\nu t}$ at every time step: the +exact field inside the integral never left $t = 0$. `uw.maths.Integral`, `BdIntegral` and +`CellWiseIntegral` compile through the same JIT as the solvers, which routes every +`uw.function.expression` to PETSc's constants array, but none of them set the constants on +the DS they integrate with, so the kernels read zeros: any expression in an integrand +integrated to nothing, and a fresh Integral returned the cached zero. The constitutive +viscosity is such an expression, which is where the cylinder drag went (next section). Fixed +on this branch (`petsc_maths.pyx`, the boundary integral sets them on its sandbox DS); +`tests/test_0503_integral_expression_constants.py`. + +## The DDt as the transport plugin + +Louis asked (2026-09-06) whether the history manager could be the object that decides +how transport is done, so that one solver takes SUPG where it is needed and a +semi-Lagrangian history where it is not. It can, and it now is. Every solver that owns an +unknown composes its residual from three contributions of its `DuDt`: + +| contribution | `EulerianSUPG` | `SemiLagrangian`, `Eulerian`, `Lagrangian` | +|---|---|---| +| `time_derivative()` | $(\psi^{n+1}-\psi^n)/\Delta t$ (theta rule) or the BDF stencil over the history | the same, over its own history | +| `advection()` | $\sum_k w_k\,(\mathbf{a}_k\cdot\nabla)\psi^{(k)}$, entry by entry of the unknown | zero (the history carries it) | +| `stabilisation_flux(R)` | $\tau\,R\otimes\mathbf{a}$, one flux row per component of $R$ | zero | +| `states()`, `spatial_weights()` | the levels and the weights $w_k$ of the scheme, for the solver's own flux | the same | + +The scalar solver assembles $f_0 = \dot\phi + \mathbf{u}\cdot\nabla\phi - f$ and +$\mathbf{f}_1 = \sum_k w_k\kappa\nabla\phi^{(k)} + \tau R\mathbf{u}$ from these; the +Navier-Stokes solver multiplies the first two by $\rho$, adds $\nabla p$ to the residual +the flux sees, and keeps the viscous flux of the scheme and the pressure as its own. The +manager owns what the transport needs: the advecting velocity as data (`V_fn`, and +`V_fn_history` for the stored levels, which is the stored velocity itself for momentum), +the timestep as a runtime constant (`delta_t`, written by every flavour's +`update_pre_solve`), the time scheme, the diffusivity that $\tau$ sees (set by the solver +from its constitutive model), and the stabilisation knobs. The nonlinearity of a +self-advected unknown lives in what `V_fn` is: the extrapolated field, the Picard iterate, +or the unknown's own symbol for Newton. + +What this bought, measured: the refactor moved no physics. The Péclet-weight rows of +Kovasznay at 1/16 and 1/32 (1.596e-4, 2.082e-5), the vortex decay at 1/32 (5.272e-5, +energy ratio 0.960789) and the Blankenbach box (42.8675 / 4.9204 / 4.8840) reproduce to +every printed digit, and the two-rank tests keep their serial constants. The cylinder at +1/20 with LU keeps its mean drag, lift extrema, reaction drag and Strouhal number to the +printed digits (3.0532, 0.9193 / -0.9652, 3.1219, 0.2964) while the drag peak moves from +3.0797 to 3.0802 and the pressure difference at peak lift from 2.4134 to 2.4125: the +assembled expressions are the same terms in a different order, and a shedding wake +amplifies the last bits over 1400 steps where a steady state does not. A `SemiLagrangian` manager dropped into `AdvDiffusion` reproduces +`AdvDiffusionSLCN` to the solver tolerance on pure advection (test_1057): the solver's +equation with zero advection and zero stabilisation is the semi-Lagrangian one. A tensor +unknown, flattened to its independent components on a `MATRIX` variable, is transported +through the multi-component solver with a residual that is nothing but the manager's +terms (test_1057, uniform translation of a Gaussian stress to 5%). That is stress +transport without rotation; the rotation of a transported tensor is a constitutive +matter and stays out of the transport. + +Two things the plain `Eulerian` manager keeps: with a velocity it still applies the +explicit splitting correction to the history (its `_advection_mode` is `"split"`), which +is what the Richards and Darcy solvers rely on; `EulerianSUPG` sets the mode to +`"assembled"` and the solver's residual carries the advection instead. And the +semi-Lagrangian Stokes stress history (`DFDt` on a viscoelastic Stokes solve) is +untouched: it advects a stress that is not an unknown of the solve, which is a different +job from the one the contract describes. + +## What the timestep estimate means + +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 + +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 b88e30519..82f67deb2 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -166,6 +166,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/docs/examples/Tutorial_Thermal_Convection_Units.py b/docs/examples/Tutorial_Thermal_Convection_Units.py index 31da1f0d3..8bdc614be 100644 --- a/docs/examples/Tutorial_Thermal_Convection_Units.py +++ b/docs/examples/Tutorial_Thermal_Convection_Units.py @@ -264,7 +264,7 @@ def kelvin_to_celsius(temp): print(" Top/Bottom: No-slip, v = 0") print(" Left/Right: Free-slip, vx = 0") -thermal = uw.systems.AdvDiffusion( +thermal = uw.systems.AdvDiffusionSLCN( mesh, u_Field=temperature, V_fn=velocity, diff --git a/docs/examples/WIP/developer_tools/SOpt.py b/docs/examples/WIP/developer_tools/SOpt.py index e1b06f974..8c713d014 100644 --- a/docs/examples/WIP/developer_tools/SOpt.py +++ b/docs/examples/WIP/developer_tools/SOpt.py @@ -732,7 +732,7 @@ def pipemesh_return_coords_to_bounds(coords): # %% -field_advection = uw.systems.AdvDiffusion(openmesh, u_Field=obstruction, V_fn=v_phi, order=1) +field_advection = uw.systems.AdvDiffusionSLCN(openmesh, u_Field=obstruction, V_fn=v_phi, order=1) field_advection.constitutive_model = uw.constitutive_models.DiffusionModel field_advection.constitutive_model.Parameters.diffusivity = 1.0 field_advection.estimate_dt() diff --git a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSLCN_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSLCN_RotationTest.py index e24a6e80c..907b42a5c 100644 --- a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSLCN_RotationTest.py +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSLCN_RotationTest.py @@ -139,7 +139,7 @@ r_i = params.uw_radius_inner r_o = params.uw_radius_outer -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshball, u_Field=t_soln, V_fn=v_soln, 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..fbe3b378e --- /dev/null +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py @@ -0,0 +1,174 @@ +# --- +# 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.AdvDiffusion`. 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, 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 + +- **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. +- **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_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) +""" + +# %% +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_dt_fraction 0.1 -uw_order 2 +``` +""" + +# %% +params = uw.Params( + uw_res=32, + 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, + 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.AdvDiffusion( + 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) + +# %% [markdown] +""" +## Time loop + +`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_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: + 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 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. +""" + +# %% +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/docs/examples/convection/advanced/Ex_AdvectionDiffusionSwarm_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSwarm_RotationTest.py index 92f00fc10..9f6400590 100644 --- a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSwarm_RotationTest.py +++ b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSwarm_RotationTest.py @@ -91,7 +91,7 @@ # + -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshball, u_Field=t_soln, V_fn = v_soln, diff --git a/docs/examples/convection/advanced/Ex_Convection_Cartesian-Swarm.py b/docs/examples/convection/advanced/Ex_Convection_Cartesian-Swarm.py index a29a59f2a..7c5b093d1 100644 --- a/docs/examples/convection/advanced/Ex_Convection_Cartesian-Swarm.py +++ b/docs/examples/convection/advanced/Ex_Convection_Cartesian-Swarm.py @@ -94,7 +94,7 @@ # + -ad = uw.systems.AdvDiffusion(meshbox, t_soln, T1.sym, order=3) +ad = uw.systems.AdvDiffusionSLCN(meshbox, t_soln, T1.sym, order=3) ad._u_star_projector.smoothing = 0.0 diff --git a/docs/examples/convection/advanced/Ex_Convection_Cylinder.py b/docs/examples/convection/advanced/Ex_Convection_Cylinder.py index a379d19e0..109f339fc 100644 --- a/docs/examples/convection/advanced/Ex_Convection_Cylinder.py +++ b/docs/examples/convection/advanced/Ex_Convection_Cylinder.py @@ -166,7 +166,7 @@ """ # %% -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshball, u_Field=t_soln, V_fn=v_soln, diff --git a/docs/examples/convection/intermediate/Ex_Convection_4_SLCN_Cartesian-NL.py b/docs/examples/convection/intermediate/Ex_Convection_4_SLCN_Cartesian-NL.py index b5d91fb0b..9e037309d 100644 --- a/docs/examples/convection/intermediate/Ex_Convection_4_SLCN_Cartesian-NL.py +++ b/docs/examples/convection/intermediate/Ex_Convection_4_SLCN_Cartesian-NL.py @@ -161,7 +161,7 @@ """ # %% -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshbox, u_Field=t_soln, V_fn=v_soln, diff --git a/docs/examples/convection/intermediate/Ex_Convection_5_SLCN_Cartesian-Yield.py b/docs/examples/convection/intermediate/Ex_Convection_5_SLCN_Cartesian-Yield.py index 8f43a1274..3ce7f4fc8 100644 --- a/docs/examples/convection/intermediate/Ex_Convection_5_SLCN_Cartesian-Yield.py +++ b/docs/examples/convection/intermediate/Ex_Convection_5_SLCN_Cartesian-Yield.py @@ -158,7 +158,7 @@ """ # %% -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshbox, u_Field=t_soln, V_fn=v_soln, diff --git a/docs/examples/convection/intermediate/Ex_Convection_Cartesian_ThermoChem.py b/docs/examples/convection/intermediate/Ex_Convection_Cartesian_ThermoChem.py index d1fc0df60..bcacf2f74 100644 --- a/docs/examples/convection/intermediate/Ex_Convection_Cartesian_ThermoChem.py +++ b/docs/examples/convection/intermediate/Ex_Convection_Cartesian_ThermoChem.py @@ -171,7 +171,7 @@ # %% k = params.uw_diffusivity -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshbox, u_Field=t_soln, V_fn=v_soln, diff --git a/docs/examples/convection/intermediate/Ex_MoresiSolomatov_Convection_Cartesian.py b/docs/examples/convection/intermediate/Ex_MoresiSolomatov_Convection_Cartesian.py index 8ed65b526..a0aadeda0 100644 --- a/docs/examples/convection/intermediate/Ex_MoresiSolomatov_Convection_Cartesian.py +++ b/docs/examples/convection/intermediate/Ex_MoresiSolomatov_Convection_Cartesian.py @@ -173,7 +173,7 @@ """ # %% -adv_diff = uw.systems.AdvDiffusion( +adv_diff = uw.systems.AdvDiffusionSLCN( meshbox, u_Field=t_soln, V_fn=v_soln, diff --git a/docs/examples/fluid_mechanics/README.md b/docs/examples/fluid_mechanics/README.md index 4346c45f7..00c49b7b8 100644 --- a/docs/examples/fluid_mechanics/README.md +++ b/docs/examples/fluid_mechanics/README.md @@ -61,6 +61,16 @@ Fluid mechanics forms the foundation for understanding mantle convection, magma - Interface dynamics and surface tension - Applications: magma-crystal systems, air-water flows +10. **Navier-Stokes on the grid: lid-driven cavity** - `Ex_Navier_Stokes_SUPG_Lid_Driven_Cavity.py` + - The Eulerian SUPG Navier-Stokes solver at Re 100 against Ghia et al. (1982) + - One linear solve per step; Picard or Newton for the fully implicit form + - Introduces: `NavierStokes`, the cell-Peclet weight of the stabilisation + +11. **Navier-Stokes on the grid: Taylor-Green vortex decay** - `Ex_Navier_Stokes_SUPG_Taylor_Green_Vortex.py` + - An exact unsteady solution: the velocity error and the energy decay measured directly + - Free-slip walls as partial Dirichlet conditions + - Introduces: validation against an exact time-dependent solution + ## 🧮 Mathematical Background ### Governing Equations diff --git a/docs/examples/fluid_mechanics/advanced/Ex_NavierStokesRotationTest.py b/docs/examples/fluid_mechanics/advanced/Ex_NavierStokesRotationTest.py index 7171e82dd..cbc0a3579 100644 --- a/docs/examples/fluid_mechanics/advanced/Ex_NavierStokesRotationTest.py +++ b/docs/examples/fluid_mechanics/advanced/Ex_NavierStokesRotationTest.py @@ -148,7 +148,7 @@ """ # %% -navier_stokes = uw.systems.NavierStokes( +navier_stokes = uw.systems.NavierStokesSLCN( meshball, velocityField=v_soln, pressureField=p_soln, diff --git a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d.py b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d.py index fa9489101..b526dfba0 100644 --- a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d.py +++ b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d.py @@ -363,7 +363,7 @@ def pipemesh_return_coords_to_bounds(coords): """ # %% -navier_stokes = uw.systems.NavierStokes( +navier_stokes = uw.systems.NavierStokesSLCN( pipemesh, velocityField=v_soln, pressureField=p_soln, diff --git a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py index 72e6bb64f..dbad10bbb 100644 --- a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py +++ b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py @@ -337,7 +337,7 @@ def pipemesh_return_coords_to_bounds(coords): """ # %% -navier_stokes = uw.systems.NavierStokes( +navier_stokes = uw.systems.NavierStokesSLCN( pipemesh, velocityField=v_soln, pressureField=p_soln, diff --git a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Lid_Driven_Flow_2d.py b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Lid_Driven_Flow_2d.py index 61bf75adf..89b2c9fd7 100644 --- a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Lid_Driven_Flow_2d.py +++ b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_Lid_Driven_Flow_2d.py @@ -41,7 +41,7 @@ import sympy import underworld3 as uw -from underworld3.systems import NavierStokes +from underworld3.systems import NavierStokesSLCN # Ghia et al. (1982) reference data: u-velocity along vertical centreline GHIA_Y = np.array([0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, @@ -68,7 +68,7 @@ def run_cavity(order): p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=True, vtype=uw.VarType.SCALAR) - ns = NavierStokes(mesh, velocityField=v, pressureField=p, rho=1.0, order=order) + ns = NavierStokesSLCN(mesh, velocityField=v, pressureField=p, rho=1.0, order=order) ns.constitutive_model = uw.constitutive_models.ViscousFlowModel ns.constitutive_model.Parameters.viscosity = 1.0 / RE ns.saddle_preconditioner = 1.0 diff --git a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Lid_Driven_Cavity.py b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Lid_Driven_Cavity.py new file mode 100644 index 000000000..27cfcad84 --- /dev/null +++ b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Lid_Driven_Cavity.py @@ -0,0 +1,151 @@ +# %% [markdown] +""" +# Navier-Stokes Lid-Driven Cavity with the Eulerian SUPG solver (Re = 100) + +**PHYSICS:** fluid_mechanics +**DIFFICULTY:** advanced +**RUNTIME:** ~4 minutes + +## Description + +The lid-driven cavity at Re = 100 with `uw.systems.NavierStokes`, the +Navier-Stokes solver that assembles the momentum advection on the grid and +stabilises it with streamline-upwind Petrov-Galerkin weighting. Each step is +one linear Oseen solve with the advecting velocity extrapolated from the two +stored levels. The centreline velocity extrema are compared with Ghia, Ghia & +Shin (1982). + +## Key Concepts + +- Eulerian (grid-based) Navier-Stokes with SUPG stabilisation +- The advecting velocity: extrapolated, Picard-corrected, or implicit +- The cell-Peclet weight of the stabilisation (on by default) +- Marching to a steady state and reading centreline profiles + +## Reference + +Ghia, Ghia & Shin (1982), "High-Re solutions for incompressible flow using +the Navier-Stokes equations and a multigrid method", J. Comp. Physics 48, 387-411. +""" + +# %% [markdown] +""" +## Parameters +""" + +# %% +RE = 100.0 # PARAM: Reynolds number (unit lid speed, unit cavity, viscosity 1/Re) +CELLSIZE = 1 / 32 # PARAM: mesh element size +COURANT = 1.0 # PARAM: time step as a multiple of the cell-crossing time at the lid +NSTEPS = 400 # PARAM: number of time steps +PICARD = 0 # PARAM: extra Picard passes per step (0 = one linear solve per step) + +# %% +import numpy as np +import sympy +from mpi4py import MPI +import underworld3 as uw + +# %% [markdown] +""" +## Mesh and fields + +P2 velocity and P1 pressure (Taylor-Hood) on an unstructured simplex mesh. +""" + +# %% +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=CELLSIZE, qdegree=3) +v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + +# %% [markdown] +""" +## The solver + +`NavierStokes` is a subclass of the Stokes solver: it takes the same +constitutive model and boundary conditions. `rho=1` with viscosity `1/RE` +gives Re on the unit cavity. `advection="extrapolated"` (the default) makes +each step one linear solve; `picard_iterations` re-solves with the latest +iterate for the fully implicit fixed point; `advection="implicit"` lets the +nonlinear solver take Newton steps instead. +""" + +# %% +ns = uw.systems.NavierStokes( + mesh, v, p, rho=1.0, order=1, advection="extrapolated", picard_iterations=PICARD) +ns.constitutive_model = uw.constitutive_models.ViscousFlowModel +ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 / RE +ns.tolerance = 1.0e-6 + +for boundary in ("Left", "Right", "Bottom"): + ns.add_dirichlet_bc((0.0, 0.0), boundary) +ns.add_dirichlet_bc((1.0, 0.0), "Top") # the lid, singular at the corners +ns.bodyforce = sympy.Matrix([[0.0, 0.0]]) + +# %% [markdown] +""" +## Time stepping + +The Courant number is the number of cells the lid crosses in a step. The +implicit scheme has no stability limit on it; Courant 1 is a good accuracy +choice for the transient, and the run is stopped when the velocity stops +changing. +""" + +# %% +dt = COURANT * CELLSIZE / 1.0 +line = np.linspace(0.0, 1.0, 201) +vertical = np.c_[0.5 * np.ones_like(line), line] # x = 0.5: u(y) +horizontal = np.c_[line, 0.5 * np.ones_like(line)] # y = 0.5: v(x) + +def centreline_extrema(): + u_c = uw.function.evaluate(v.sym[0], vertical).reshape(-1) + v_c = uw.function.evaluate(v.sym[1], horizontal).reshape(-1) + comm = uw.mpi.comm + # evaluate() answers for the points this rank owns: reduce the extrema. + return (comm.allreduce(float(u_c.min()), op=MPI.MIN), + comm.allreduce(float(v_c.max()), op=MPI.MAX), + comm.allreduce(float(v_c.min()), op=MPI.MIN)) + +for step in range(NSTEPS): + before = np.array(v.array[...]) + ns.solve(timestep=dt, zero_init_guess=False) + change = float(np.abs(np.asarray(v.array[...]) - before).max()) if before.size else 0.0 + change = uw.mpi.comm.allreduce(change, op=MPI.MAX) + if (step + 1) % 50 == 0 or step == 0: + u_min, v_max, v_min = centreline_extrema() + uw.pprint(f"step {step + 1:4d} t {dt * (step + 1):.3f} " + f"u_min {u_min:.4f} v_max {v_max:.4f} v_min {v_min:.4f} change {change:.2e}") + if change < 1.0e-6: + break + +# %% [markdown] +""" +## Comparison with Ghia et al. (1982) + +On a 1/32 mesh the three extrema come within about 4% of the reference; the +difference is the mesh (the extrema are steady to four digits). +""" + +# %% +GHIA = dict(u_min=-0.2109, v_max=0.1753, v_min=-0.2453) +u_min, v_max, v_min = centreline_extrema() +uw.pprint(f"u_min on x = 0.5: {u_min:.4f} (Ghia {GHIA['u_min']})") +uw.pprint(f"v_max on y = 0.5: {v_max:.4f} (Ghia {GHIA['v_max']})") +uw.pprint(f"v_min on y = 0.5: {v_min:.4f} (Ghia {GHIA['v_min']})") +assert abs(u_min - GHIA["u_min"]) < 0.02 and abs(v_min - GHIA["v_min"]) < 0.02 + +# %% [markdown] +""" +## Notes + +- `ns.supg_weight = 0` gives the plain Galerkin form; on this mesh at Re 100 + it runs as well, because the element Reynolds number is small. Set the + weight back to 1 and raise Re to see where the stabilisation starts to matter. +- `peclet_weight` (default 4) turns the stabilisation off in cells that are + diffusion-dominated, where it is not needed and costs accuracy. +- The design note `docs/developer/design/eulerian-supg-transport.md` records + the benchmarks (Kovasznay flow, this cavity to Re 1000, the DFG cylinder, + Taylor-Green vortex decay). +""" diff --git a/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Taylor_Green_Vortex.py b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Taylor_Green_Vortex.py new file mode 100644 index 000000000..7f3219185 --- /dev/null +++ b/docs/examples/fluid_mechanics/advanced/Ex_Navier_Stokes_SUPG_Taylor_Green_Vortex.py @@ -0,0 +1,138 @@ +# %% [markdown] +""" +# Taylor-Green vortex decay with the Eulerian SUPG Navier-Stokes solver + +**PHYSICS:** fluid_mechanics +**DIFFICULTY:** advanced +**RUNTIME:** ~1 minute + +## Description + +An exact unsteady solution of the Navier-Stokes equations: a lattice of +counter-rotating vortices that decays in place, + + u = (-sin x cos y, cos x sin y) exp(-2 nu t), p = (cos 2x + cos 2y)/4 exp(-4 nu t), + +on the box [0, pi]^2. On that box the walls carry no normal flow and no +tangential stress, so free-slip walls (the normal component fixed) are exact +and nothing on the boundary depends on time. The velocity error against the +exact solution at the end of the run measures the scheme directly. + +## Key Concepts + +- Time-dependent validation against an exact Navier-Stokes solution +- Free-slip walls as partial Dirichlet conditions +- The kinetic energy decay, exp(-4 nu t), as a second check +- Where the SUPG stabilisation costs accuracy and how the Peclet weight removes it +""" + +# %% [markdown] +""" +## Parameters +""" + +# %% +NU = 0.01 # PARAM: viscosity (density 1) +RES = 16 # PARAM: cells across the box +DT = 0.025 # PARAM: time step +T_END = 0.5 # PARAM: end time +PECLET_WEIGHT = 4.0 # PARAM: cell-Peclet weight of the stabilisation (0 = uniform) + +# %% +import numpy as np +import sympy +import underworld3 as uw + +# %% [markdown] +""" +## Mesh, fields and the exact solution +""" + +# %% +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(np.pi, np.pi), cellSize=np.pi / RES, regular=True, qdegree=3) +x, y = mesh.X +v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + +def exact(t): + F = sympy.exp(-2 * NU * t) + U = sympy.Matrix([[-sympy.sin(x) * sympy.cos(y) * F, sympy.cos(x) * sympy.sin(y) * F]]) + P = (sympy.cos(2 * x) + sympy.cos(2 * y)) * F ** 2 / 4 + return U, P + +U0, P0 = exact(0.0) +v.array[:, 0, :] = uw.function.evaluate(U0, v.coords).reshape(-1, 2) +p.array[:, 0, 0] = uw.function.evaluate(P0, p.coords).reshape(-1) + +# %% [markdown] +""" +## The solver with free-slip walls + +A partial Dirichlet condition fixes one component and leaves the other free: +`(0.0, None)` on the vertical walls, `(None, 0.0)` on the horizontal ones. +""" + +# %% +ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=1, peclet_weight=PECLET_WEIGHT) +ns.constitutive_model = uw.constitutive_models.ViscousFlowModel +ns.constitutive_model.Parameters.shear_viscosity_0 = NU +ns.tolerance = 1.0e-8 +ns.add_dirichlet_bc((0.0, None), "Left") +ns.add_dirichlet_bc((0.0, None), "Right") +ns.add_dirichlet_bc((None, 0.0), "Bottom") +ns.add_dirichlet_bc((None, 0.0), "Top") +ns.bodyforce = sympy.Matrix([[0.0, 0.0]]) + +# %% [markdown] +""" +## The error against the exact solution + +The exact velocity is F(t) U0, so the L2 error expands into three integrals +that carry no time dependence: ||v - F U0||^2 = - 2F + F^2 . +""" + +# %% +I_vv = uw.maths.Integral(mesh, v.sym.dot(v.sym)) +I_vU = uw.maths.Integral(mesh, v.sym.dot(U0)) +I_UU = float(uw.maths.Integral(mesh, U0.dot(U0)).evaluate()) + +def velocity_error(t): + F = np.exp(-2 * NU * t) + vv, vU = float(I_vv.evaluate()), float(I_vU.evaluate()) + return np.sqrt(max(vv - 2 * F * vU + F ** 2 * I_UU, 0.0) / (F ** 2 * I_UU)) + +E0 = float(I_vv.evaluate()) +uw.pprint(f"interpolation error of the exact field on this mesh: {velocity_error(0.0):.2e}") + +# %% [markdown] +""" +## March +""" + +# %% +n_steps = int(round(T_END / DT)) +t = 0.0 +for step in range(n_steps): + ns.solve(timestep=DT, zero_init_guess=False) + t += DT + if (step + 1) % 5 == 0 or step + 1 == n_steps: + uw.pprint(f"step {step + 1:3d} t {t:.3f} velocity error {velocity_error(t):.3e} " + f"E/E0 {float(I_vv.evaluate()) / E0:.6f} exact {np.exp(-4 * NU * t):.6f}") + +# %% [markdown] +""" +## What to expect + +On the 1/16 mesh the velocity error at t = 0.5 is a few times 1e-4 (the P2 +interpolation error is 1.4e-4) and the kinetic energy follows exp(-4 nu t) +to six digits. With `PECLET_WEIGHT = 0` the stabilisation acts in every cell +and the error rises by a fixed factor: this flow is resolved and needs no +stabilisation, which is what the weight detects. With `ns.supg_weight = 0` +(plain Galerkin) the error is the same as with the weight. +""" + +# %% +err = velocity_error(t) +uw.pprint(f"final velocity error {err:.3e}, energy ratio {float(I_vv.evaluate()) / E0:.6f} (exact {np.exp(-4 * NU * t):.6f})") +assert err < 2.0e-3 diff --git a/docs/examples/heat_transfer/advanced/Ex_Poisson_v.SLCN.py b/docs/examples/heat_transfer/advanced/Ex_Poisson_v.SLCN.py index c8f2ed7ba..ce57144f6 100644 --- a/docs/examples/heat_transfer/advanced/Ex_Poisson_v.SLCN.py +++ b/docs/examples/heat_transfer/advanced/Ex_Poisson_v.SLCN.py @@ -118,7 +118,7 @@ poisson1 = uw.systems.Poisson(mesh, u_Field=phi) -poisson2 = uw.systems.AdvDiffusion(mesh, +poisson2 = uw.systems.AdvDiffusionSLCN(mesh, u_Field=phi, V_fn = V.sym, order = 1) 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. diff --git a/src/underworld3/cython/petsc_maths.pyx b/src/underworld3/cython/petsc_maths.pyx index 5be87ac74..a4445d780 100644 --- a/src/underworld3/cython/petsc_maths.pyx +++ b/src/underworld3/cython/petsc_maths.pyx @@ -1,5 +1,6 @@ from typing import Union import sympy +import numpy as np import underworld3 import underworld3.timing as timing @@ -15,6 +16,31 @@ cdef extern from "petsc.h" nogil: PetscErrorCode DMPlexComputeCellwiseIntegralFEM( PetscDM, PetscVec, PetscVec, void* ) +def _pack_manifest(manifest): + """The current values of the JIT constants manifest as a contiguous array.""" + from underworld3.utilities._jitextension import _pack_constants + if not manifest: + return None + return np.ascontiguousarray(_pack_constants(manifest), dtype=np.float64) + + +cdef _set_ds_constants(PetscDS ds, manifest): + """Hand the current UWexpression values to the DS the integral kernel reads. + + The JIT routes every ``uw.function.expression`` in the integrand to PETSc's + constants array (the same mechanism the solvers use, so a changed value does + not recompile). A DS that never receives the values hands the kernel zeros: + a viscosity, a time or any other expression in an integrand silently + integrated to nothing (found on the cylinder drag, 2026-09-05). + """ + cdef double[::1] vals + values = _pack_manifest(manifest) + if values is None or len(values) == 0: + return + vals = values + CHKERRQ(PetscDSSetConstants(ds, len(values), &vals[0])) + + def dm_force_coordinate_field(dm): """Force coordinate field creation and strip boundary labels from the coordinate DM. Must be called after createCoordinateSpace and after @@ -122,8 +148,9 @@ class Integral: cdef DS ds = self.dm.getDS() cdef PetscScalar val_array[256] - # Now set callback... + # Now set callback (and the current constant values the kernel reads)... ierr = PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]); CHKERRQ(ierr) + _set_ds_constants(ds.ds, _getext_result.constants_manifest) ierr = DMPlexComputeIntegralFEM(dm.dm, cgvec.vec, &(val_array[0]), NULL); CHKERRQ(ierr) self.dm.restoreGlobalVec(a_global) @@ -290,8 +317,9 @@ class CellWiseIntegral: elif isinstance(self.fn, sympy.vector.Dyadic): raise RuntimeError("Integral evaluation for Dyadic integrands not supported.") - cdef PtrContainer ext = getext(self.mesh, JITCallbackSet(residual=(self.fn,)), - self.mesh.vars.values()).ptrobj + _getext_result = getext(self.mesh, JITCallbackSet(residual=(self.fn,)), + self.mesh.vars.values()) + cdef PtrContainer ext = _getext_result.ptrobj # Pull out vec for variables, and go ahead with the integral self.mesh.update_lvec() @@ -316,6 +344,7 @@ class CellWiseIntegral: cdef DM dm = self.mesh.dm cdef DS ds = self.mesh.dm.getDS() CHKERRQ( PetscDSSetObjective(ds.ds, 0, ext.fns_residual[0]) ) + _set_ds_constants(ds.ds, _getext_result.constants_manifest) # DMPlexComputeCellwiseIntegralFEM writes Nf scalars per cell into a # flat [cell*Nf + field] layout when the output vector carries no @@ -461,6 +490,11 @@ class BdIntegral: cdef PetscDMLabel sandbox_label = NULL CHKERRQ(DMGetLabel(sandbox_dm, boundary_bytes, &sandbox_label)) + # The sandbox has its own DS (DMCreateDS): the constants go there. + cdef PetscDS sandbox_ds = NULL + CHKERRQ(DMGetDS(sandbox_dm, &sandbox_ds)) + _set_ds_constants(sandbox_ds, _getext_result.constants_manifest) + # Output value cdef PetscScalar result = 0.0 diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 56a5222a0..b0844701c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3892,13 +3892,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/function/expressions.py b/src/underworld3/function/expressions.py index 48cfb8ae6..e34c094b1 100644 --- a/src/underworld3/function/expressions.py +++ b/src/underworld3/function/expressions.py @@ -1165,23 +1165,29 @@ def is_extended_real(self): @property def is_positive(self): - """Delegate to wrapped expression.""" - if self._sym is not None and hasattr(self._sym, 'is_positive'): - return self._sym.is_positive + """Unknown, always: a UWexpression is a runtime constant whose value can + change after construction, so sympy must not fold on its current sign or + on it being zero (#696: ``exp(c)`` with ``c`` created at 0 became 1 at + construction, freezing a time ramp). The value is read when the + expression is unwrapped for compilation, not here.""" return None @property def is_negative(self): - """Delegate to wrapped expression.""" - if self._sym is not None and hasattr(self._sym, 'is_negative'): - return self._sym.is_negative + """Unknown, always: a UWexpression is a runtime constant whose value can + change after construction, so sympy must not fold on its current sign or + on it being zero (#696: ``exp(c)`` with ``c`` created at 0 became 1 at + construction, freezing a time ramp). The value is read when the + expression is unwrapped for compilation, not here.""" return None @property def is_zero(self): - """Delegate to wrapped expression.""" - if self._sym is not None and hasattr(self._sym, 'is_zero'): - return self._sym.is_zero + """Unknown, always: a UWexpression is a runtime constant whose value can + change after construction, so sympy must not fold on its current sign or + on it being zero (#696: ``exp(c)`` with ``c`` created at 0 became 1 at + construction, freezing a time ramp). The value is read when the + expression is unwrapped for compilation, not here.""" return None @property diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 4b4b679d3..9788c447a 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -4894,6 +4894,13 @@ def advection( substep). Default ``False`` here; note :meth:`NodalPointSwarm.advection` defaults it to ``True``. """ + if self.local_size < 0: + # DMSwarm reports -1 until particles have been added on some rank (#702); + # an EMPTY rank of a populated swarm is size 0 and is handled below. + raise RuntimeError( + "This swarm has never been populated (no particles were added on any " + "rank): call populate() or add_particles_with_coordinates() before advection." + ) # Convert delta_t to model units if it has units # This ensures consistent arithmetic: velocity is in model units, so time must be too import underworld3 as uw diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index ea6823ab5..17097bbf0 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -18,9 +18,14 @@ Projection : class L2 projection of fields onto mesh variables. AdvDiffusion : class - Advection-diffusion with semi-Lagrangian transport. + Advection-diffusion composed from a DDt transport manager (the default + manager, EulerianSUPG, assembles implicit advection with SUPG). +AdvDiffusionSLCN : class + Advection-diffusion with semi-Lagrangian transport (flux history). NavierStokes : class - Navier-Stokes equations with inertia. + Navier-Stokes composed from a DDt transport manager (EulerianSUPG default). +NavierStokesSLCN : class + Navier-Stokes with semi-Lagrangian transport and a stress history. Diffusion : class Pure diffusion (no advection). TransientDarcy : class @@ -30,8 +35,9 @@ Time Derivative Schemes ----------------------- -Lagrangian_DDt, SemiLagragian_DDt, Eulerian_DDt - Time derivative approximations for transient problems. +Lagrangian_DDt, SemiLagragian_DDt, Eulerian_DDt, EulerianSUPG_DDt + Time derivative approximations for transient problems; EulerianSUPG_DDt + is the transport plugin of the Eulerian solvers (assembled advection, SUPG). See Also -------- @@ -63,7 +69,10 @@ # 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 +# The generic names are the composing solvers: the transport (assembled SUPG +# advection, or a semi-Lagrangian history) is the DDt manager they hold. +from .advection_diffusion_eulerian import SNES_AdvectionDiffusion_Composed as AdvDiffusion +from .navier_stokes_eulerian import SNES_NavierStokes_Composed as NavierStokes # import diffusion-only solver from .solvers import SNES_Diffusion as Diffusion @@ -75,7 +84,6 @@ # These are now implemented the same way using the ddt module from .solvers import SNES_NavierStokes as NavierStokesSwarm from .solvers import SNES_NavierStokes as NavierStokesSLCN -from .solvers import SNES_NavierStokes as NavierStokes from .free_surface import FreeSurface @@ -87,6 +95,7 @@ from .ddt import SemiLagrangian as SemiLagragian_DDt from .ddt import Lagrangian_Swarm as Lagrangian_Swarm_DDt from .ddt import Eulerian as Eulerian_DDt +from .ddt import EulerianSUPG as EulerianSUPG_DDt # δ-continuation driver for hard viscoplastic (Drucker–Prager) yield from .yield_continuation import yield_continuation, YieldHomotopyControl diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py new file mode 100644 index 000000000..06d274943 --- /dev/null +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -0,0 +1,717 @@ +r"""Advection-diffusion composed from a DDt transport manager. + +The solver assembles the diffusive flux and the source on the mesh and takes +its transport (time derivative, advection, stabilisation) from the history +manager it holds. The default manager, :class:`~underworld3.systems.ddt.EulerianSUPG`, +makes it the fully implicit Eulerian scheme with SUPG stabilisation described +below; a :class:`~underworld3.systems.ddt.SemiLagrangian` manager makes it a +semi-Lagrangian scheme on the field history. + +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``. + +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 + +import numpy as np +import sympy +from typing import Callable, Optional, Union + +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 _DDtBase, _as_row_vector +from underworld3.systems.ddt import Eulerian as Eulerian_DDt +from underworld3.systems.ddt import EulerianSUPG as EulerianSUPG_DDt +from underworld3.systems.solvers import ( + _advective_diffusive_dt, + _dimensionalise_dt, + _invalidate_solution_cache, + _nondimensionalise_timestep, +) + + +def _check_supplied_manager(DuDt, order, theta): + """A supplied history manager fixes the scheme: the arguments must agree with it.""" + if DuDt.order != order: + raise ValueError( + f"DuDt supplied is order {DuDt.order} but order={order} was asked for: a " + "supplied manager fixes the scheme, pass the matching order.") + if theta is not None and hasattr(DuDt, "theta") and float(DuDt.theta) != float(theta): + raise ValueError( + f"DuDt supplied has theta={DuDt.theta} but theta={theta} was asked for.") + + +class SNES_AdvectionDiffusion_Composed(SNES_Scalar): + r"""Advection-diffusion solver composed from its DDt transport manager. + + With the default manager (:class:`~underworld3.systems.ddt.EulerianSUPG`): + implicit in time, assembled on the mesh, SUPG in space. + + .. math:: + \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f + + 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.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(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 + ========== ======= ===================================================== + + ``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: + + 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 + + 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 + + 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 + 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, advection, source) 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()``) 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, as :math:`(\mathbf{u}\Delta t)^2` for the second-order schemes. + :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 + serial and needs no departure points in parallel. + + 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 + 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. + verbose : bool, default False + DuDt : DDt history manager, optional + The transport plugin. By default an + :class:`~underworld3.systems.ddt.EulerianSUPG` built from ``V_fn``, + ``order`` and ``theta``. Any history manager that follows the DDt + transport contract (``time_derivative``, ``advection``, + ``stabilisation_flux``) can be supplied instead: a + :class:`~underworld3.systems.ddt.SemiLagrangian` history turns this + solver into a semi-Lagrangian scheme on the field history, with no + assembled advection and no stabilisation. A supplied manager fixes + ``order`` and ``theta``. + 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 every + 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 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 + def __init__( + self, + mesh: uw.discretisation.Mesh, + u_Field: uw.discretisation.MeshVariable, + V_fn, + order: int = 1, + theta: Optional[float] = None, + peclet_weight: float = 4.0, + 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." + ) + 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"AdvDiffusion 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 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) + # 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). + 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 " + "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) + + self.f = sympy.Matrix.zeros(1, 1) + self._last_timestep = None + self._last_change_rate = None + + # The transport plugin: the history manager owns the time scheme, the + # advecting velocity, the assembled advection and the stabilisation. + if DuDt is None: + self.Unknowns.DuDt = EulerianSUPG_DDt( + self.mesh, + u_Field, + V_fn, + vtype=uw.VarType.SCALAR, + degree=u_Field.degree, + continuous=u_Field.continuous, + order=order, + theta=theta, + varsymbol=u_Field.symbol, + verbose=verbose, + bcs=self.essential_bcs, + smoothing=0.0, + peclet_weight=peclet_weight, + ) + else: + if not isinstance(DuDt, _DDtBase): + raise TypeError(f"DuDt must be a DDt history manager, not {type(DuDt).__name__}.") + if sympy.Matrix(DuDt.psi_fn).shape != u_Field.sym.shape: + raise ValueError("DuDt tracks a different unknown from u_Field.") + _check_supplied_manager(DuDt, order, theta) + self.Unknowns.DuDt = DuDt + self._theta = float(getattr(self.DuDt, "theta", theta)) + + # 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 + + # 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 + + super()._object_viewer() + scheme = {("am", 1): f"Adams-Moulton order 1, theta = {self._theta}", + ("bdf", 1): "backward Euler"}.get( + (self.integrator, self.order), + f"{self.integrator.upper()} order {self.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 + # ------------------------------------------------------------------ + + @property + def integrator(self) -> str: + """The multistep family in use: ``"am"`` (the theta rule) at order 1, ``"bdf"`` above.""" + return self.DuDt.integrator + + @property + def order(self) -> int: + """Requested order of the time integration.""" + return self.DuDt.order + + @property + def theta(self) -> float: + """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.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." + ) + if not hasattr(self.DuDt, "theta"): + raise AttributeError(f"{type(self.DuDt).__name__} has no theta to set.") + self._theta = value + self.DuDt.theta = value + + @property + def delta_t(self): + 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.DuDt.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.DuDt.delta_t.sym = dt + self._last_timestep = dt + + @property + def V_fn(self): + """Advecting velocity, ``(1, dim)`` (the history manager's).""" + return self.DuDt.V_fn + + @V_fn.setter + def V_fn(self, value): + self.DuDt.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 + + # The stabilisation knobs live on the history manager; these pass through. + + @property + def peclet_weight(self) -> float: + """The critical cell Péclet number of the weight (constructor choice; 0 = uniform).""" + return self.DuDt.peclet_weight + + @property + def supg_weight(self) -> float: + """Scale of the SUPG term: 1 (default) or 0 for plain Galerkin. No rebuild.""" + return self.DuDt.supg_weight + + @supg_weight.setter + def supg_weight(self, value): + self.DuDt.supg_weight = value + + @property + def tau_weights(self): + r"""The weights :math:`(C_t, C_u, C_\kappa)` of the three terms in :math:`\tau`.""" + return self.DuDt.tau_weights + + @tau_weights.setter + def tau_weights(self, values): + self.DuDt.tau_weights = values + + # ------------------------------------------------------------------ + # The residual, composed from the history manager's contributions + # ------------------------------------------------------------------ + + 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.DuDt.spatial_weights(), self.DuDt.states()): + if w == 0: + continue + grad = self.mesh.vector.gradient(phi[0]) + total = total + w * (grad * c) + return total + + def _strong_residual(self): + """Time derivative, advection and source, as a ``(1, 1)`` matrix.""" + return self.DuDt.time_derivative() + self.DuDt.advection() - self._f + + 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 _stabilisation_flux(self): + if hasattr(self.DuDt, "diffusivity"): + self.DuDt.diffusivity = self._scalar_diffusivity() + return self.DuDt.stabilisation_flux(self._strong_residual()) + + F0 = Template( + r"f_0(\phi)", + lambda self: 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._stabilisation_flux(), + "Diffusive flux of the time scheme plus the SUPG flux tau R u.", + ) + + # ------------------------------------------------------------------ + # Timestep and solve + # ------------------------------------------------------------------ + + @timing.routine_timer_decorator + 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 + ---------- + 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. + """ + 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(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, + zero_init_guess: Optional[bool] = None, + timestep=None, + _force_setup: bool = False, + _evalf: bool = False, + verbose: bool = False, + divergence_retries: int = 0, + ): + r"""Advance :math:`\phi` by one step. + + 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 not None: + self.delta_t = timestep + elif self._last_timestep is None: + raise ValueError( + "solve() needs a timestep: pass timestep=
or set solver.delta_t first." + ) + dt = self._last_timestep + + if _force_setup: + self._needs_function_rewire = True + if not self.constitutive_model._solver_is_setup: + self._needs_function_rewire = True + # 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) + before = np.array(self.u.data).reshape(-1) + 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 (from a copy of the unknown: the manager's + # history may live on a swarm or carry units). + from mpi4py import MPI + change = np.abs(np.asarray(self.u.data).reshape(-1) - before) + 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 + self.constitutive_model._solver_is_setup = True diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e198aad57..f82f58c42 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -169,13 +169,21 @@ class DDtLagrangianSwarmState(_DDtCoreState): def _as_float(value): - """Extract a plain float from various numeric types (Pint, UWQuantity, etc.).""" + """A plain, NON-DIMENSIONAL float from a number, a Pint quantity or a UWQuantity. + + A quantity is scaled by the active model's reference scales (#701: taking + its magnitude gave the kernels a dimensional timestep); without reference + scales the magnitude is what non-dimensionalisation returns. + """ if value is None: return None if isinstance(value, (int, float)): return float(value) - if hasattr(value, "magnitude"): - return float(value.magnitude) + if hasattr(value, "magnitude") or hasattr(value, "dimensionality"): + nd = uw.non_dimensionalise(value) + if hasattr(nd, "magnitude"): + return float(nd.magnitude) + return float(nd) if hasattr(value, "value"): return float(value.value) try: @@ -541,12 +549,44 @@ class _DDtBase(uw_object): def _init_history_tracking(self, order): """Deferred-initialisation and variable-dt bookkeeping attributes.""" + # The timestep as a runtime constant of the compiled kernels: every + # flavour writes it through the ``_dt`` property, so a solver that + # composes its residual from :meth:`time_derivative` never recompiles + # when the step changes. + self._delta_t = _UWexpression( + rf"\Delta t_{{{self.instance_number}}}", 1.0, "DDt timestep", + _unique_name_generation=True) # History tracking: deferred initialization and effective order self._history_initialised = False self._n_solves_completed = 0 self._dt = None # current timestep (set by solver or update_pre_solve) self._dt_history = [None] * order # previous timesteps for variable-dt BDF + @property + def _dt(self): + return self._dt_value + + @_dt.setter + def _dt(self, value): + self._dt_value = value + if value is None: + return + try: + dt = float(_as_float(value)) + except Exception: + return + if dt > 0.0: + self._delta_t.sym = dt + + @property + def delta_t(self): + r"""The timestep :math:`\Delta t` as a UW expression (a runtime constant). + + Written by ``update_pre_solve`` and by a solver's ``delta_t`` setter; + read by :meth:`time_derivative`. + """ + return self._delta_t + def _init_coefficient_expressions(self, order, theta, with_exp): """Create BDF/AM (and optionally ETD-2 exp) coefficient UWexpressions. @@ -666,6 +706,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. @@ -708,6 +771,94 @@ def initiate_history_fn(self): """Deprecated: use ``initialise_history`` instead.""" self.initialise_history() + # ----- The transport contract ----- + # + # A solver that owns an unknown composes its residual from these terms + # and never asks which flavour it holds: + # + # F0 = time_derivative() + advection() - f + # F1 = + # + stabilisation_flux(R) + # + # The history flavours (Symbolic, Eulerian, SemiLagrangian, Lagrangian) + # carry their transport in the history itself, so advection() and the + # stabilisation flux are zero for them; EulerianSUPG assembles both. + + @property + def integrator(self) -> str: + """``"am"`` (the theta rule on the spatial terms) at order 1, ``"bdf"`` above.""" + return "am" if self.order == 1 else "bdf" + + def _unknown_shape(self): + """Shape of the unknown as a matrix (``Symbolic`` stores ``_shape`` as data).""" + psi = self.psi_fn + return psi.shape if isinstance(psi, sympy.MatrixBase) else (1, 1) + + def states(self): + r"""``[psi^{n+1}, psi^{n}, psi^{n-1}, ...]`` as matrices of the unknown's shape.""" + return [sympy.Matrix(self.psi_fn)] + [sympy.Matrix(h) for h in self._history_syms()] + + def spatial_weights(self): + """Weight of a spatial operator at each level of :meth:`states`. + + ``[1, 0, ...]`` for the BDF family (every spatial term at n+1); the + Adams-Moulton weights for the theta rule. + """ + n = len(self.psi_star) + if self.integrator == "bdf": + return [sympy.Integer(1)] + [sympy.Integer(0)] * n + return list(self.am_coefficient_expressions[: n + 1]) + + def time_derivative(self): + r"""The time derivative of the scheme, a matrix of the unknown's shape. + + ``(psi^{n+1} - psi^{n}) / dt`` for the theta rule, the BDF stencil over + the history divided by ``dt`` above order 1, with ``dt`` the runtime + constant :attr:`delta_t`. + """ + if self.integrator == "am": + new, old = self.states()[:2] + return (new - old) / self._delta_t + return sympy.Matrix(self.bdf()) / self._delta_t + + def advection(self): + """The assembled advection term: zero for a history-carrying flavour.""" + return sympy.zeros(*self._unknown_shape()) + + def stabilisation_flux(self, R): + r"""The stabilisation flux for a strong residual ``R``: zero here. + + Shape ``(len(R), dim)``: one flux row per component of ``R``. + """ + mesh = getattr(self, "mesh", None) + if mesh is None: + raise TypeError(f"{type(self).__name__} has no mesh: no flux shape to return.") + return sympy.zeros(len(_as_matrix(R)), mesh.dim) + + +def _as_matrix(R): + """A residual as a sympy Matrix: a bare scalar becomes ``(1, 1)``.""" + return R if isinstance(R, sympy.MatrixBase) else sympy.Matrix([[R]]) + + +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 Symbolic(_DDtBase): r""" @@ -1050,11 +1201,16 @@ def __init__( bcs=[], order=1, smoothing=0.0, + num_components=None, ): super().__init__() self.mesh = mesh self.V_fn = V_fn + # With a velocity, the plain Eulerian flavour applies it as an + # explicit splitting correction of the history ("split"); + # EulerianSUPG assembles it in the solver's residual instead. + self._advection_mode = "split" self.theta = theta self.bcs = bcs self.verbose = verbose @@ -1063,6 +1219,7 @@ def __init__( self.continuous = continuous self.smoothing = smoothing self.evalf = evalf + self.num_components = num_components self._init_history_tracking(order) @@ -1093,6 +1250,7 @@ def __init__( uw.discretisation.MeshVariable( f"psi_star_Eulerian_{self.instance_number}_{i}", self.mesh, + num_components, vtype=vtype, degree=degree, continuous=continuous, @@ -1330,32 +1488,36 @@ def update_pre_solve( _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) _update_am_values(self._am_coeffs, self.effective_order, self.theta) - if self.V_fn is not None and dt is not None: - coords = self.psi_star[0].coords - dim = self.mesh.dim - X = self.mesh.X - - # Build u·∇φ symbolically for each component of psi_fn - # psi_fn is a Matrix; V_fn is also a Matrix. For scalar - # psi_fn the shape is (1,1); for vector it is (1,dim). - psi = self.psi_fn - V = self.V_fn - ncomp = max(psi.shape) # number of tracked components - - for c in range(ncomp): - # ∂φ_c/∂x_i for each spatial dimension - grad_c = sympy.Matrix([psi[c].diff(X[i]) for i in range(dim)]) - # u·∇φ_c = V_i * ∂φ_c/∂x_i - advection_expr = sum(V[i] * grad_c[i] for i in range(dim)) - - advection_vals = uw.function.evaluate( - advection_expr, coords, evalf=evalf, - ).reshape(-1) - - self.psi_star[0].data[:, c] -= dt * advection_vals + if self.V_fn is not None and dt is not None and self._advection_mode == "split": + self._apply_split_advection(dt, evalf) return + def _apply_split_advection(self, dt, evalf=False): + """Explicit operator-splitting correction: ``psi_star[0] -= dt (V . grad) psi``.""" + coords = self.psi_star[0].coords + dim = self.mesh.dim + X = self.mesh.X + + # Build u·∇φ symbolically for each component of psi_fn + # psi_fn is a Matrix; V_fn is also a Matrix. For scalar + # psi_fn the shape is (1,1); for vector it is (1,dim). + psi = self.psi_fn + V = self.V_fn + ncomp = max(psi.shape) # number of tracked components + + for c in range(ncomp): + # ∂φ_c/∂x_i for each spatial dimension + grad_c = sympy.Matrix([psi[c].diff(X[i]) for i in range(dim)]) + # u·∇φ_c = V_i * ∂φ_c/∂x_i + advection_expr = sum(V[i] * grad_c[i] for i in range(dim)) + + advection_vals = uw.function.evaluate( + advection_expr, coords, evalf=evalf, + ).reshape(-1) + + self.psi_star[0].data[:, c] -= dt * advection_vals + def update_post_solve( self, dt, @@ -1390,6 +1552,250 @@ def update_exp_coefficients(self, dt, tau_eff): _update_exp_values(self._exp_coeffs, dt, tau_eff) +class EulerianSUPG(Eulerian): + r"""Eulerian history manager that assembles its transport: implicit advection with SUPG. + + The transport plugin of the Eulerian solvers. It holds the history of + one unknown on the mesh, as :class:`Eulerian` does, and contributes the + three terms a solver composes its residual from: the time derivative of + the multistep scheme, the implicit advection + :math:`\sum_k w_k\,(\mathbf{a}_k\cdot\nabla)\psi^{(k)}` applied + component-wise to a scalar, a vector or a tensor unknown, and the + streamline-upwind Petrov-Galerkin flux :math:`\tau\,R\otimes\mathbf{a}` + of the solver's strong residual :math:`R`. The same solver takes a + :class:`SemiLagrangian` history in its place: that flavour answers zero + for the advection and the flux because its history is already traced + back along the characteristics. + + ``V_fn`` is data: the velocity the transport uses at the new level. The + nonlinearity of a self-advected unknown lives in what ``V_fn`` is (the + unknown's own symbol for Newton, an extrapolated or Picard field for a + linear step), and ``V_fn_history`` names the velocity at the stored + levels when it is not ``V_fn`` (the stored velocity itself for momentum). + + The stabilisation parameter is + + .. math:: + \tau = \frac{w}{\sqrt{(C_t c_0/\Delta t)^2 + (C_u|\mathbf{a}|/h)^2 + (C_\kappa\kappa/h^2)^2}} + + (``tau_shape="inverse_sum"``) with :math:`h` the local cell size, + :math:`c_0` the leading multistep coefficient, :math:`\kappa` the + :attr:`diffusivity` the solver declares (the diffusivity of a scalar, + :math:`\eta/\rho` for momentum, zero for a transported stress) and + :math:`w` the product of ``supg_weight`` and the cell-Péclet weight + :math:`Pe^2/(Pe^2 + Pe_c^2)`, :math:`Pe = |\mathbf{a}|h/2\kappa`, which + switches the term off where diffusion dominates. ``"brooks_hughes"`` and + ``"doubly_asymptotic"`` are the optimal 1-D shapes, each capped by the + transient term. Every weight is a runtime constant of the kernels. + + Parameters + ---------- + mesh, psi_fn, vtype, degree, continuous, varsymbol, verbose, bcs, smoothing + As for :class:`Eulerian`; ``psi_fn`` is the unknown's MeshVariable. + V_fn : MeshVariable or sympy row Matrix + The advecting velocity, ``(1, dim)``. + order : int, default 1 + 1 is the theta rule (Crank-Nicolson at ``theta=0.5``), 2 and 3 BDF. + theta : float, optional + Crank-Nicolson blend at order 1 (0.5 default; 1.0 backward Euler). + Orders 2 and 3 take ``theta=1.0`` and refuse anything else. + diffusivity : expression, default 0 + What :math:`\tau` sees as the diffusive rate; a solver sets it from + its constitutive model when it builds its flux. + supg_weight, tau_weights, tau_shape, peclet_weight + The stabilisation knobs described above. + num_components : tuple, optional + The history variable shape when ``vtype`` is ``MATRIX``. + """ + + _TAU_SHAPES = ("inverse_sum", "brooks_hughes", "doubly_asymptotic") + + @timing.routine_timer_decorator + def __init__( + self, + mesh: uw.discretisation.Mesh, + psi_fn, + V_fn, + vtype: uw.VarType, + degree: int, + continuous: bool, + order: int = 1, + theta: Optional[float] = None, + varsymbol: Optional[str] = r"u", + verbose: Optional[bool] = False, + bcs=None, + smoothing: float = 0.0, + diffusivity=0, + supg_weight: float = 1.0, + tau_weights=(2.0, 2.0, 4.0), + tau_shape: str = "inverse_sum", + peclet_weight: float = 4.0, + num_components=None, + ): + 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 is not None else (0.5 if order == 1 else 1.0) + 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 (a BDF stencil " + "pairs with terms at n+1, not with a centred flux)." + ) + if tau_shape not in self._TAU_SHAPES: + raise ValueError(f"tau_shape must be one of {self._TAU_SHAPES}, got {tau_shape!r}") + + # A caller's list is kept BY REFERENCE on purpose: a solver passes its + # live essential_bcs so conditions added later reach the projections. + # Only the default gets a fresh list, never a shared one. + super().__init__( + mesh, psi_fn, vtype, degree, continuous, V_fn=None, theta=theta, + varsymbol=varsymbol, verbose=verbose, bcs=[] if bcs is None else bcs, + order=order, smoothing=smoothing, num_components=num_components, + ) + self._advection_mode = "assembled" + self._integrator = "am" if order == 1 else "bdf" + self.V_fn = V_fn + self.V_fn_history = None + self.diffusivity = diffusivity + self._tau_shape = str(tau_shape) + self._peclet_weight = float(peclet_weight) + + # The stabilisation knobs are runtime constants. + tag = self.instance_number + unique = dict(_unique_name_generation=True) + self._supg_weight = _UWexpression( + rf"w^{{\mathrm{{SUPG}}}}_{{{tag}}}", 1.0, "SUPG term weight (0 = Galerkin)", **unique) + self._tau_weights = [ + _UWexpression(rf"C^{{\tau}}_{{t,{tag}}}", 2.0, "tau transient weight", **unique), + _UWexpression(rf"C^{{\tau}}_{{u,{tag}}}", 2.0, "tau advective weight", **unique), + _UWexpression(rf"C^{{\tau}}_{{\kappa,{tag}}}", 4.0, "tau diffusive weight", **unique), + ] + self.supg_weight = supg_weight + self.tau_weights = tau_weights + + # ----- data ----- + + @property + def V_fn(self): + """The advecting velocity at the new level, ``(1, dim)``.""" + return self._V_fn + + @V_fn.setter + def V_fn(self, value): + self._V_fn = None if value is None else _as_row_vector(value, self.mesh.dim) + + @property + def integrator(self) -> str: + return self._integrator + + def advecting_velocity(self, level: int = 0): + """The velocity carrying the unknown at ``states()[level]``.""" + if level == 0 or not self.V_fn_history: + return self.V_fn + return _as_row_vector(self.V_fn_history[level - 1], self.mesh.dim) + + @property + def tau_shape(self) -> str: + return self._tau_shape + + @property + def peclet_weight(self) -> float: + return self._peclet_weight + + @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): + for w, v in zip(self._tau_weights, values): + w.sym = float(v) + + # ----- the contract ----- + + def _convective(self, a, psi): + r"""``(a . grad) psi`` entry by entry, a matrix of ``psi``'s shape.""" + dim = self.mesh.dim + grad = self.mesh.vector.gradient + + def entry(r, c): + g = grad(psi[r, c]) + return sum(a[0, i] * g[0, i] for i in range(dim)) + + return sympy.Matrix(*psi.shape, entry) + + def advection(self): + r""":math:`\sum_k w_k\,(\mathbf{a}_k\cdot\nabla)\psi^{(k)}` over the levels of the scheme.""" + total = sympy.zeros(*self._unknown_shape()) + for k, (w, psi_k) in enumerate(zip(self.spatial_weights(), self.states())): + if w == 0: + continue + total = total + w * self._convective(self.advecting_velocity(k), psi_k) + return total + + def tau(self): + r"""The stabilisation parameter :math:`\tau` (times the weights).""" + dim = self.mesh.dim + a = self.advecting_velocity(0) + a_mag2 = sum(a[0, i] ** 2 for i in range(dim)) + h = self.mesh.cell_size() + nu = self.diffusivity + if self.integrator == "bdf": + c0 = self.bdf_coefficient_expressions[0] + else: + c0 = sympy.Integer(1) + ct, cu, cv = self._tau_weights + transient = (ct * c0 / self._delta_t) ** 2 + weight = self._supg_weight + if self._peclet_weight > 0.0: + # Pe^2 / (Pe^2 + Pe_c^2) written without dividing by nu (1 for nu = 0). + ah2 = a_mag2 * h ** 2 + weight = weight * ah2 / (ah2 + 4 * self._peclet_weight ** 2 * nu ** 2 + 1.0e-30) + if self._tau_shape == "inverse_sum": + advective = (cu * sympy.sqrt(a_mag2) / h) ** 2 + viscous = (cv * nu / h ** 2) ** 2 + return weight / sympy.sqrt(transient + advective + viscous + 1.0e-30) + # The 1-D optimal shapes: tau = (h / 2|a|) xi(Pe), Pe = |a| h / (2 nu). + a_mag = sympy.sqrt(a_mag2 + 1.0e-30) + Pe = a_mag * h / (2 * nu + 1.0e-30) # finite at zero diffusivity (the default) + if self._tau_shape == "brooks_hughes": + xi = 1 / sympy.tanh(Pe) - 1 / Pe # coth is not C99: the printer would rewrite it through exp + else: + xi = sympy.Min(Pe / 3, 1) + tau_steady = h / (2 * a_mag) * xi + return weight / sympy.sqrt(transient + 1 / (tau_steady ** 2 + 1.0e-30)) + + def stabilisation_flux(self, R): + r"""The SUPG flux :math:`\tau\,R\otimes\mathbf{a}`, one row per component of ``R``. + + ``R`` is the solver's strong residual of the unknown's shape (first + derivatives only). The result has shape ``(len(R), dim)``: for a + scalar the row :math:`\tau R\mathbf{a}`, for a vector + :math:`F_{ij} = \tau R_i a_j`. + """ + R = _as_matrix(R) + column = R.reshape(len(R), 1) + return self.tau() * (column * self.advecting_velocity(0)) + + def _object_viewer(self): + from IPython.display import Latex, display + + super()._object_viewer() + display(Latex(r"$\quad\mathbf{a} = $ " + self.V_fn._repr_latex_())) + display(Latex(rf"$\quad$ integrator: {self.integrator}, tau shape: {self.tau_shape}")) + + class SemiLagrangian(_DDtBase): r""" Semi-Lagrangian history manager using nodal swarm. diff --git a/src/underworld3/systems/navier_stokes_eulerian.py b/src/underworld3/systems/navier_stokes_eulerian.py new file mode 100644 index 000000000..58f92d85e --- /dev/null +++ b/src/underworld3/systems/navier_stokes_eulerian.py @@ -0,0 +1,568 @@ +r"""Navier-Stokes composed from a DDt transport manager. + +The solver assembles the viscous flux, the pressure and the body force and +takes the momentum transport from the history manager it holds. With the +default manager, :class:`~underworld3.systems.ddt.EulerianSUPG`, it is the +Eulerian scheme with SUPG momentum transport described below. + +The incompressible Navier-Stokes equations solved on the mesh with the +momentum advection assembled implicitly in the saddle-point residual and +stabilised by the streamline-upwind Petrov-Galerkin term, the vector +counterpart of :class:`~underworld3.systems.AdvDiffusion`. The time +scheme is the same multistep family: Crank-Nicolson (the theta rule) at +order 1, BDF2 at order 2, with the history held on the mesh by the +Eulerian history manager. No stress history is carried: the viscous stress +at an earlier level is rebuilt from the stored velocity level through the +constitutive model. + +The advecting velocity :math:`\mathbf{a}` in :math:`(\mathbf{a}\cdot\nabla)\mathbf{u}^{n+1}` +is a choice (``advection=``): ``"extrapolated"`` (default) uses +:math:`2\mathbf{u}^n - \mathbf{u}^{n-1}`, a second-order lag that makes +each step one linear (Oseen) solve; ``picard_iterations`` re-solves with +the latest iterate for the fully implicit fixed point; ``"implicit"`` puts +:math:`\mathbf{u}^{n+1}` itself in the advection and lets the SNES take +Newton steps on the quadratic term. + +Design note: ``docs/developer/design/eulerian-supg-transport.md``. +""" + +import warnings + +import numpy as np +import sympy +from typing import Optional, Union + +import underworld3 as uw +import underworld3.timing as timing +from underworld3.function import expression as public_expression +from underworld3.systems.ddt import _DDtBase +from underworld3.systems.ddt import EulerianSUPG as EulerianSUPG_DDt +from underworld3.systems.advection_diffusion_eulerian import _check_supplied_manager +from underworld3.systems.solvers import SNES_Stokes, _dimensionalise_dt + +_ADVECTION_MODES = ("extrapolated", "implicit") + + +class SNES_NavierStokes_Composed(SNES_Stokes): + r"""Navier-Stokes solver composed from its DDt transport manager + (Eulerian SUPG momentum transport by default). + + Solves + + .. math:: + \rho\left(\frac{\partial \mathbf{u}}{\partial t} + + (\mathbf{u}\cdot\nabla)\mathbf{u}\right) + - \nabla\cdot\left[\boldsymbol{\tau}(\mathbf{u}) - p\mathbf{I}\right] + = \mathbf{f}, \qquad \nabla\cdot\mathbf{u} = 0, + + with :math:`\boldsymbol{\tau}` the deviatoric stress of the constitutive + model. In the pointwise form the momentum residual is + + .. math:: + \mathbf{f}_0 = \mathbf{R}, \qquad + \mathbf{F}_1 = \sum_k w_k\,\boldsymbol{\tau}(\mathbf{u}^{(k)}) + - p_\mathrm{mech}\mathbf{I} + \tau_\mathrm{s}\,\mathbf{R}\otimes\mathbf{a}, + + where :math:`\mathbf{R} = \rho\,(\dot{\mathbf{u}} + \sum_k w_k (\mathbf{a}_k\cdot\nabla)\mathbf{u}^{(k)}) + + \nabla p - \mathbf{f}` is the strong residual of the time scheme (first + derivatives only, so without the viscous term; :math:`\mathbf{f}_0` carries + it without :math:`\nabla p`, which enters through the flux), :math:`w_k` the weights of the spatial operator at each time level + (:math:`w_0 = 1` for BDF, the Adams-Moulton weights for the theta rule), + :math:`\mathbf{a}_0 = \mathbf{a}` the advecting velocity at the new level + and :math:`\mathbf{a}_k = \mathbf{u}^{(k)}` at the stored ones. The last + term of :math:`\mathbf{F}_1` is the Petrov-Galerkin perturbation + :math:`\tau_\mathrm{s}(\mathbf{a}\cdot\nabla)\mathbf{w}` applied to + :math:`\mathbf{R}`, with + + .. math:: + \tau_\mathrm{s} = \left[\left(\frac{C_t c_0}{\Delta t}\right)^2 + + \left(\frac{C_u |\mathbf{a}|}{h}\right)^2 + + \left(\frac{C_\nu\, \nu}{h^2}\right)^2\right]^{-1/2}, + \qquad \nu = \eta / \rho, + + :math:`h` the local cell size and the three weights runtime constants + (``tau_weights``). The pressure equation is the incompressibility + constraint, unchanged from the Stokes solver; Taylor-Hood elements need + no pressure stabilisation. + + Parameters + ---------- + mesh, velocityField, pressureField + As for :class:`~underworld3.systems.Stokes`. + rho : float or expression, default 1.0 + Density. + order : int, default 1 + Time scheme: 1 is the theta rule (Crank-Nicolson at ``theta=0.5``), + 2 is BDF2. + theta : float, optional + Crank-Nicolson blend at order 1 (0.5 default; 1.0 backward Euler). + Order 2 takes ``theta=1.0`` and refuses anything else. + advection : {"extrapolated", "implicit"}, default "extrapolated" + The advecting velocity at the new time level: the second-order + extrapolation :math:`2\mathbf{u}^n - \mathbf{u}^{n-1}` (a linear + step) or the unknown itself (Newton on the quadratic term). + picard_iterations : int, default 0 + With ``"extrapolated"``, the number of further passes per step that + re-solve with the latest iterate as the advecting velocity, stopping + early when the velocity stops changing (``picard_tolerance``). The + fixed point is the fully implicit scheme without a tangent. + picard_tolerance : float, default 1e-4 + Relative change of the velocity (max norm) below which the Picard + passes stop. + tau_shape : {"inverse_sum", "brooks_hughes", "doubly_asymptotic"} + The shape of the stabilisation parameter. ``"inverse_sum"`` (default) + is the Shakib-Tezduyar form above, smooth and cheap but above the + optimal 1-D curve at cell Péclet numbers of order 1 to 10. + ``"brooks_hughes"`` is the optimal 1-D form :math:`\tau = (h/2|a|)\, + (\coth Pe - 1/Pe)`, ``"doubly_asymptotic"`` its two-limit + approximation :math:`(h/2|a|)\min(Pe/3, 1)`, with :math:`Pe = |a| h / + (2\nu)`; both are combined with the transient term as + :math:`[(C_t c_0/\Delta t)^2 + \tau^{-2}]^{-1/2}` so the time step still + caps them. The advective and viscous weights are not used by these two. + peclet_weight : float, default 4 + A critical cell Péclet number. The SUPG term is multiplied by + :math:`Pe^2 / (Pe^2 + Pe_c^2)`, :math:`Pe = |a| h / 2\nu`, so the + stabilisation is off where the cell is diffusion-dominated (where it + is not needed and costs a fixed multiple of the Galerkin error) and + full where advection dominates. Measured on the vortex decay, + Kovasznay and the cylinder: at 4 the resolved cases are within 1.3 + times the Galerkin error and the cylinder wall cells keep 86% of the + term; at 8 the resolved cases sit on Galerkin and the cylinder is still + stable. Zero gives the uniform weight ``supg_weight`` everywhere. + degree, p_continuous, verbose + As for :class:`~underworld3.systems.Stokes`. + + Notes + ----- + - ``DFDt`` (a stress history) is refused: the theta rule forms the + viscous stress at level n from the stored velocity as + :math:`2\eta\,\dot\varepsilon(\mathbf{u}^n)` with the current effective + viscosity, which is exact for a constant viscosity; use ``order=2`` + (all spatial terms at n+1) with a strain-rate dependent viscosity. + - The linear solver is the Stokes fieldsplit: the velocity block is + nonsymmetric, which its smoother and flexible outer solver already + allow for, while the pressure Schur approximation is the viscous-limit + one and costs outer iterations as :math:`\rho|\mathbf{a}|\Delta t/\eta` + grows. + - The velocity history levels, the advecting-velocity field and the + extrapolation level are mesh variables the solver owns. + """ + + @timing.routine_timer_decorator + def __init__( + self, + mesh: uw.discretisation.Mesh, + velocityField: uw.discretisation.MeshVariable, + pressureField: uw.discretisation.MeshVariable, + rho=1.0, + order: int = 1, + theta: Optional[float] = None, + advection: str = "extrapolated", + picard_iterations: int = 0, + picard_tolerance: float = 1.0e-4, + tau_shape: str = "inverse_sum", + peclet_weight: float = 4.0, + degree: Optional[int] = 2, + p_continuous: Optional[bool] = True, + verbose: bool = False, + DuDt: Optional[_DDtBase] = None, + DFDt=None, + restore_points_func=None, + ): + if DFDt is not None: + raise ValueError( + "AdvDiffusion-style Navier-Stokes carries no stress history: " + "the viscous stress at earlier levels is rebuilt from the stored " + "velocity. Do not pass DFDt." + ) + if restore_points_func is not None: + warnings.warn( + "NavierStokes ignores restore_points_func: it configures the " + "semi-Lagrangian trace-back and the Eulerian scheme has none.", + stacklevel=2, + ) + order = int(order) + if order not in (1, 2): + raise ValueError(f"order must be 1 or 2, not {order}.") + theta = float(theta) if theta is not None else (0.5 if order == 1 else 1.0) + 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 takes theta=1.0." + ) + advection = str(advection) + if advection not in _ADVECTION_MODES: + raise ValueError(f"advection must be one of {_ADVECTION_MODES}, not {advection!r}.") + + super().__init__( + mesh, velocityField, pressureField, degree, p_continuous, verbose, + DuDt=None, DFDt=None, + ) + + self._theta = theta + self._advection_mode = advection + self._picard_iterations = int(picard_iterations) + self._picard_tolerance = float(picard_tolerance) + self._picard_count = 0 + self._last_timestep = None + self._last_change_rate = None + + tag = self.instance_number + self._rho = public_expression(rf"\rho_{{{tag}}}", rho, "Density") + + # The advecting velocity at the new level (values set before each + # solve: the extrapolation, or the latest Picard iterate) and the + # level n-1 the extrapolation needs beyond what the history holds. + u = self.Unknowns.u + self._a_var = uw.discretisation.MeshVariable( + f"a_NSSUPG_{tag}", self.mesh, self.mesh.dim, degree=u.degree, + continuous=u.continuous, varsymbol=rf"\mathbf{{a}}_{{{tag}}}") + self._u_prev = uw.discretisation.MeshVariable( + f"u_prev_NSSUPG_{tag}", self.mesh, self.mesh.dim, degree=u.degree, + continuous=u.continuous, varsymbol=rf"\mathbf{{u}}^{{n-1}}_{{{tag}}}") + self._history_primed = False + + # The transport plugin: the history manager owns the time scheme, the + # advecting velocity, the assembled advection and the stabilisation. + # At the stored levels the momentum is carried by the stored velocity. + if DuDt is None: + self.Unknowns.DuDt = EulerianSUPG_DDt( + self.mesh, + u, + self._advecting_velocity(), + vtype=uw.VarType.VECTOR, + degree=u.degree, + continuous=u.continuous, + order=order, + theta=theta, + varsymbol=u.symbol, + verbose=verbose, + bcs=self.essential_bcs, + smoothing=0.0, + tau_shape=tau_shape, + peclet_weight=peclet_weight, + ) + else: + if not isinstance(DuDt, _DDtBase): + raise TypeError(f"DuDt must be a DDt history manager, not {type(DuDt).__name__}.") + if sympy.Matrix(DuDt.psi_fn).shape != u.sym.shape: + raise ValueError("DuDt tracks a different unknown from the velocity.") + _check_supplied_manager(DuDt, order, theta) + self.Unknowns.DuDt = DuDt + self._theta = float(getattr(DuDt, "theta", theta)) + # This solver decides the advecting velocity (per step: the extrapolation, + # the Picard iterate, or the unknown) and names the stored velocity as the + # carrier of the stored levels, whoever built the manager. + if hasattr(self.DuDt, "V_fn_history"): + self.DuDt.V_fn = self._advecting_velocity() + self.DuDt.V_fn_history = [ps.sym for ps in self.DuDt.psi_star] + + # ------------------------------------------------------------------ + # Scheme description and knobs + # ------------------------------------------------------------------ + + @property + def integrator(self) -> str: + """``"am"`` (the theta rule) at order 1, ``"bdf"`` at order 2.""" + return self.DuDt.integrator + + @property + def order(self) -> int: + """Time scheme order.""" + return self.DuDt.order + + @property + def theta(self) -> float: + """Adams-Moulton blend at order 1 (1.0 backward Euler, 0.5 Crank-Nicolson).""" + return self._theta + + @theta.setter + def theta(self, value): + value = float(value) + if value != 1.0 and self.order != 1: + raise ValueError("theta applies at order 1 only; order 2 takes theta=1.0.") + if not hasattr(self.DuDt, "theta"): + raise AttributeError(f"{type(self.DuDt).__name__} has no theta to set.") + self._theta = value + self.DuDt.theta = value + + @property + def advection(self) -> str: + """``"extrapolated"`` (linear Oseen step) or ``"implicit"`` (Newton).""" + return self._advection_mode + + @advection.setter + def advection(self, value): + value = str(value) + if value not in _ADVECTION_MODES: + raise ValueError(f"advection must be one of {_ADVECTION_MODES}, not {value!r}.") + if value != self._advection_mode: + self._advection_mode = value + if hasattr(self.DuDt, "V_fn_history"): + self.DuDt.V_fn = self._advecting_velocity() + self.is_setup = False + + @property + def picard_iterations(self) -> int: + return self._picard_iterations + + @picard_iterations.setter + def picard_iterations(self, value): + self._picard_iterations = int(value) + + @property + def peclet_weight(self) -> float: + """The critical cell Péclet number of the weight (0 = no Péclet weighting).""" + return self.DuDt.peclet_weight + + @property + def tau_shape(self) -> str: + """The shape of the stabilisation parameter (constructor choice).""" + return self.DuDt.tau_shape + + @property + def picard_count(self) -> int: + """Picard passes the last step took beyond the first solve.""" + return self._picard_count + + @property + def rho(self): + """Density (a UW expression).""" + return self._rho + + @rho.setter + def rho(self, value): + self._rho.sym = value + + @property + def delta_t(self): + r"""The timestep :math:`\Delta t` as a UW expression (the history manager's runtime constant).""" + return self.DuDt.delta_t + + @delta_t.setter + def delta_t(self, value): + value = self._nondimensional_time(value) + self.DuDt.delta_t.sym = value + self._last_timestep = value + + # The stabilisation knobs live on the history manager; these pass through. + + @property + def supg_weight(self) -> float: + """Weight of the SUPG term; 0 gives the plain Galerkin scheme.""" + return self.DuDt.supg_weight + + @supg_weight.setter + def supg_weight(self, value): + self.DuDt.supg_weight = value + + @property + def tau_weights(self): + """The three weights of tau: transient, advective, viscous.""" + return self.DuDt.tau_weights + + @tau_weights.setter + def tau_weights(self, values): + self.DuDt.tau_weights = values + + # ------------------------------------------------------------------ + # The residual, composed from the history manager's contributions + # ------------------------------------------------------------------ + + def _advecting_velocity(self): + """The advecting velocity at the new level, as a ``(1, dim)`` row.""" + if self._advection_mode == "implicit": + return self.u.sym + return self._a_var.sym + + def _strong_residual(self, with_pressure=False): + r"""The strong momentum residual of the time scheme, first derivatives only. + + ``with_pressure=False`` gives the terms that live in :math:`\mathbf{f}_0`: + density times the time derivative and the advection, less the body + force. ``with_pressure=True`` adds :math:`\nabla p`, the residual the + SUPG term must see: the pressure is applied through the flux + :math:`-p\mathbf{I}` in :math:`\mathbf{F}_1`, so it must not appear in + :math:`\mathbf{f}_0`, but a strong residual without it is O(1) at the + exact solution and the stabilisation then injects an O(tau) error + (measured on Kovasznay flow: 50 times the Galerkin error). The + viscous term needs second derivatives the kernels do not see; it is + the remaining inconsistency for P2 velocity. + """ + # The body-force setter may store a column; the residual is a row. + dim = self.mesh.dim + f = sympy.Matrix(self.bodyforce.sym).reshape(1, dim) + R = self._rho * (self.DuDt.time_derivative() + self.DuDt.advection()) - f + if with_pressure: + X = self.mesh.X + R = R + sympy.Matrix([[self.p.sym[0].diff(X[i]) for i in range(dim)]]) + return R + + def _viscous_stress(self, u_row): + r"""Deviatoric stress ``2 eta strain(u)`` for a velocity row, with the + current effective viscosity of the constitutive model.""" + eta = self.constitutive_model.K + return 2 * eta * sympy.Matrix(self.mesh.vector.strain_tensor(u_row)) + + def _viscous_flux(self): + states = self.DuDt.states() + weights = self.DuDt.spatial_weights() + total = weights[0] * self.stress_deviator + for w, u_k in zip(weights[1:], states[1:]): + if w == 0: + continue + total = total + w * self._viscous_stress(u_k) + return total + + def _stabilisation_flux(self): + if hasattr(self.DuDt, "diffusivity"): + self.DuDt.diffusivity = self.constitutive_model.K / self._rho + return self.DuDt.stabilisation_flux(self._strong_residual(with_pressure=True)) + + @property + def F0(self): + """Pointwise momentum residual: strong residual of the time scheme.""" + f0 = public_expression( + r"\mathbf{f}_0\left( \mathbf{u} \right)", + self._strong_residual(), + "Navier-Stokes SUPG: strong residual (time derivative, advection, body force)", + ) + self._u_f0 = f0 + return f0 + + @property + def F1(self): + """Pointwise flux: weighted viscous stress, mechanical pressure, SUPG term.""" + dim = self.mesh.dim + mechanical_pressure = ( + self.p.sym[0] - self.penalty * self.constitutive_model.K * self.div_u) + F1 = public_expression( + r"\mathbf{F}_1\left( \mathbf{u} \right)", + self._viscous_flux() - sympy.eye(dim) * mechanical_pressure + + self._stabilisation_flux(), + "Navier-Stokes SUPG: viscous flux of the time scheme, pressure, tau R (x) a", + ) + self._u_f1 = F1 + return F1 + + # ------------------------------------------------------------------ + # Timestep and solve + # ------------------------------------------------------------------ + + def _set_advecting_velocity(self, values): + self._a_var.array[...] = values + + def _prime_history(self): + """First solve: the extrapolation level equals the current velocity.""" + if not self._history_primed: + self._u_prev.array[...] = self.u.array[...] + self._history_primed = True + + @timing.routine_timer_decorator + def estimate_dt(self, fraction: float = 0.02, basis: str = "accuracy"): + r"""A timestep for this scheme. + + ``basis="accuracy"`` (default): the step at which the velocity changes + by ``fraction`` of its range, from the realised rate of the last step; + before the first solve, and whenever nothing has changed yet, the + cell-crossing time of the Stokes solver (``basis="resolution"``). + """ + if basis == "resolution" or self._last_change_rate is None: + return SNES_Stokes.estimate_dt(self) + if basis != "accuracy": + raise ValueError(f"basis must be 'accuracy' or 'resolution', not {basis!r}.") + from mpi4py import MPI + comm = uw.mpi.comm + speed = np.linalg.norm(np.asarray(self.u.array).reshape(-1, self.mesh.dim), axis=1) + hi = comm.allreduce(float(speed.max()) if speed.size else 0.0, op=MPI.MAX) + rate = self._last_change_rate + dt = fraction * hi / rate if rate > 0.0 else np.inf + if np.isinf(dt) or hi <= 0.0: + return SNES_Stokes.estimate_dt(self) + # The same units as the resolution estimate: a quantity when a model + # with reference scales is active, a plain number otherwise. + return _dimensionalise_dt(dt) + + @timing.routine_timer_decorator + def solve( + self, + zero_init_guess: Optional[bool] = None, + timestep=None, + _force_setup: bool = False, + verbose: bool = False, + picard_iterations: Optional[int] = None, + divergence_retries: int = 0, + **kwargs, + ): + r"""Advance the velocity and pressure by one step. + + ``timestep`` sets :attr:`delta_t`; omit it to reuse the last value. + With ``advection="extrapolated"`` the step is one linear solve, plus + up to ``picard_iterations`` further solves with the latest iterate as + the advecting velocity; with ``"implicit"`` the SNES solves the + quadratic term by Newton iteration. + """ + for name in ("time", "order", "evalf", "_evalf", "homotopy"): + kwargs.pop(name, None) + if kwargs: + warnings.warn(f"NavierStokes.solve ignores {sorted(kwargs)}", stacklevel=2) + if timestep is not None: + self.delta_t = timestep + elif self._last_timestep is None: + raise ValueError("solve() needs a timestep: pass timestep=
or set solver.delta_t first.") + dt = self._last_timestep + + if _force_setup: + self._needs_function_rewire = True + if not self.constitutive_model._solver_is_setup: + self._needs_function_rewire = True + # The base _build resolves the preconditioner choice against the mesh + # before the SNES reads its options; the setup stages must not be run + # directly here (they mark the solver set up first, #683). + self._build(verbose) + + self._prime_history() + u_n = np.array(self.u.array[...]) + if self._advection_mode == "extrapolated": + self._set_advecting_velocity(2.0 * u_n - np.asarray(self._u_prev.array[...])) + self.DuDt.update_pre_solve(dt, verbose=verbose) + + passes = 1 + if self._advection_mode == "extrapolated": + n_picard = self._picard_iterations if picard_iterations is None else int(picard_iterations) + passes += max(n_picard, 0) + from mpi4py import MPI + comm = uw.mpi.comm + self._picard_count = 0 + for k in range(passes): + previous = np.array(self.u.array[...]) + if k > 0: + self._set_advecting_velocity(previous) + SNES_Stokes.solve( + self, zero_init_guess if k == 0 else False, + _force_setup=_force_setup if k == 0 else False, + verbose=verbose, picard=0, divergence_retries=divergence_retries, + ) + # The reductions run on every pass, outside any branch: a rank must + # never skip a collective its peers take (tests/test_0052). + change = np.abs(np.asarray(self.u.array[...]) - previous).max() if previous.size else 0.0 + scale = np.abs(np.asarray(self.u.array[...])).max() if previous.size else 0.0 + change = comm.allreduce(float(change), op=MPI.MAX) + scale = comm.allreduce(float(scale), op=MPI.MAX) + if k > 0: + self._picard_count = k + if change <= self._picard_tolerance * max(scale, 1.0e-300): + break + + # Realised rate of change of the velocity, for estimate_dt. + change = np.linalg.norm( + (np.asarray(self.u.array[...]) - u_n).reshape(-1, self.mesh.dim), axis=1) + local = float(change.max()) if change.size else 0.0 + self._last_change_rate = comm.allreduce(local, op=MPI.MAX) / dt + + # Shift the extrapolation level, then the history. + self._u_prev.array[...] = self.DuDt.psi_star[0].array[...] + self.DuDt.update_post_solve(dt, verbose=verbose) + + self.is_setup = True + self.constitutive_model._solver_is_setup = True diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index c89963f99..524a411bf 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. @@ -4362,111 +4460,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( diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 7d63ae183..395cc7231 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -798,6 +798,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, @@ -847,7 +865,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. @@ -873,11 +891,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}]" @@ -923,7 +952,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/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") 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..987b64a1c --- /dev/null +++ b/tests/parallel/test_1077_advdiff_supg_parallel.py @@ -0,0 +1,50 @@ +"""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 (re-recorded with the local cell size, #687; +# np=2 reproduced it to 1.4e-12). +SERIAL_ERROR = 0.030152131513640566 + + +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.AdvDiffusion(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: + # the partition effect this guards against was 5e-4 (#687); platforms differ at 1e-7 + assert abs(err - SERIAL_ERROR) < 1e-6 * SERIAL_ERROR, (err, SERIAL_ERROR) diff --git a/tests/parallel/test_1078_navier_stokes_supg_parallel.py b/tests/parallel/test_1078_navier_stokes_supg_parallel.py new file mode 100644 index 000000000..53aa4e64d --- /dev/null +++ b/tests/parallel/test_1078_navier_stokes_supg_parallel.py @@ -0,0 +1,56 @@ +"""The Eulerian SUPG Navier-Stokes solver gives the serial answer on any number of ranks. + +Kovasznay flow (exact steady Navier-Stokes at Re 40), a few steps from the +exact solution; the integral velocity error must match a serial reference to +solver tolerance, whatever the partition. + +Run: mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1078_navier_stokes_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 1/4 refined once (1/8), geometric multigrid on the velocity +# block, Crank-Nicolson, dt 0.05, 6 steps, peclet_weight 4 (recorded with this file). +# The GAMG fallback without a hierarchy gave a platform-dependent answer (7% on the +# Linux CI), so the mesh carries a refinement hierarchy and the solve is tight. +SERIAL_ERROR = 0.0013227881494769559 + + +def _run(tolerance=1.0e-8): + Re = 40.0 + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-0.5, -0.5), maxCoords=(1.0, 0.5), cellSize=1.0 / 4, qdegree=3, regular=False, + refinement=1) + x, y = mesh.X + lam = Re / 2 - sympy.sqrt(Re ** 2 / 4 + 4 * sympy.pi ** 2) + U_ex = sympy.Matrix([[1 - sympy.exp(lam * x) * sympy.cos(2 * sympy.pi * y), + lam / (2 * sympy.pi) * sympy.exp(lam * x) * sympy.sin(2 * sympy.pi * y)]]) + v = uw.discretisation.MeshVariable("U1078", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("P1078", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0) + ns.constitutive_model = uw.constitutive_models.ViscousFlowModel + ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 / Re + ns.tolerance = tolerance + for b in ("Left", "Right", "Top", "Bottom"): + ns.add_dirichlet_bc(U_ex, b) + v.array[:, 0, :] = uw.function.evaluate(U_ex, v.coords).reshape(-1, 2) + for _ in range(6): + ns.solve(timestep=0.05) + err2 = uw.maths.Integral(mesh, (v.sym - U_ex).dot(v.sym - U_ex)).evaluate() + return float(np.sqrt(err2)) + + +def test_error_is_partition_independent(): + err = _run() + print(f"SERIAL_ERROR_MEASURED {err!r}") + 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: + # the partition effect this guards against was 5e-4 (#687); platforms differ at 1e-7 + assert abs(err - SERIAL_ERROR) < 1e-6 * SERIAL_ERROR, (err, SERIAL_ERROR) diff --git a/tests/test_0006_memory_leak.py b/tests/test_0006_memory_leak.py index 2066546c5..8c9c09bac 100644 --- a/tests/test_0006_memory_leak.py +++ b/tests/test_0006_memory_leak.py @@ -54,7 +54,7 @@ def test_stokes_advdiff_memory_leak(): stokes.add_dirichlet_bc([0.0, sympy.oo], "Right") # AdvDiff - advdiff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v) + advdiff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=v) advdiff.constitutive_model = uw.constitutive_models.DiffusionModel advdiff.constitutive_model.Parameters.diffusivity = 1.0 advdiff.add_dirichlet_bc([0.0], "Top") diff --git a/tests/test_0008_snapshot_realsolver.py b/tests/test_0008_snapshot_realsolver.py index d37f15de0..28c6a055c 100644 --- a/tests/test_0008_snapshot_realsolver.py +++ b/tests/test_0008_snapshot_realsolver.py @@ -78,7 +78,7 @@ def _build(): v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=1) T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) - adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=v) + adv_diff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=v) adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel adv_diff.constitutive_model.Parameters.diffusivity = 1.0 adv_diff.add_dirichlet_bc(0.0, "Left") diff --git a/tests/test_0052_collective_guard_scan.py b/tests/test_0052_collective_guard_scan.py index 5f3aaa6c3..74cf6668e 100644 --- a/tests/test_0052_collective_guard_scan.py +++ b/tests/test_0052_collective_guard_scan.py @@ -90,6 +90,12 @@ "`_domain_boundary_facets`, which is allgathered and deduplicated by " "exact coordinate identity. Both are the same bytes everywhere, so " "every rank computes the same flag.", + ("systems/navier_stokes_eulerian.py", "solve"): + "the Picard loop breaks on `change <= tol * scale`, and both `change` " + "and `scale` are `comm.allreduce(..., MAX)` results taken on every " + "pass by every rank before the test, so every rank leaves the loop on " + "the same pass and the reductions of the next pass are reached by all " + "or by none.", ("utilities/rotated_bc.py", "solve_rotated_freeslip"): "the cache is created and destroyed on all ranks together, so " "`cache is not None` is uniform; the allgather inside exists " diff --git a/tests/test_0116_swarm_never_populated.py b/tests/test_0116_swarm_never_populated.py new file mode 100644 index 000000000..a0e6376ee --- /dev/null +++ b/tests/test_0116_swarm_never_populated.py @@ -0,0 +1,29 @@ +"""A swarm that was never populated says so when advected (#702). + +Before any particles are added on any rank the DMSwarm size is -1 and the +advection died inside numpy ("negative dimensions are not allowed"). An empty +rank of a populated swarm (size 0) is a different, valid case. + +Run: pixi run python -m pytest tests/test_0116_swarm_never_populated.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def test_advecting_a_never_populated_swarm_is_a_clear_error(): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + x, y = mesh.X + swarm = uw.swarm.Swarm(mesh) + assert swarm.local_size < 0 + with pytest.raises(RuntimeError, match="never been populated"): + swarm.advection(sympy.Matrix([[-y, x]]), 0.01) + # control: the same swarm, populated, advects + swarm.populate(fill_param=1) + before = np.array(swarm.data) + swarm.advection(sympy.Matrix([[-y, x]]), 0.01) + assert swarm.local_size > 0 and np.abs(np.asarray(swarm.data) - before).max() > 0 diff --git a/tests/test_0200_solver_smoke.py b/tests/test_0200_solver_smoke.py index 9879b79a7..f8ae9a576 100644 --- a/tests/test_0200_solver_smoke.py +++ b/tests/test_0200_solver_smoke.py @@ -97,7 +97,7 @@ def test_advection_diffusion_solver_runs(self): T.array[:, 0, 0] = 0.5 v.array[:, 0, :] = 0.0 - adv_diff = uw.systems.AdvDiffusion( + adv_diff = uw.systems.AdvDiffusionSLCN( mesh, u_Field=T, V_fn=v.sym, diff --git a/tests/test_0503_integral_expression_constants.py b/tests/test_0503_integral_expression_constants.py new file mode 100644 index 000000000..1e90dbbc3 --- /dev/null +++ b/tests/test_0503_integral_expression_constants.py @@ -0,0 +1,79 @@ +"""A `uw.function.expression` inside an integrand must reach the integral kernel. + +The JIT routes every expression constant to PETSc's constants array (so that a +changed value does not recompile). The integral classes compiled through that +path but never set the values on the DS they integrate with, so the kernel read +zeros: any integrand carrying a viscosity, a time, or any other expression +integrated to nothing, and a fresh Integral returned the same zero from the +cache. Found on the cylinder drag (the viscous traction vanished), 2026-09-05. +""" +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 setup(): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.25, regular=True, qdegree=3) + x, y = mesh.X + T = uw.discretisation.MeshVariable("T_c", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(y ** 2, T.coords).reshape(-1) # dT/dy = 2y + return mesh, x, y, T + + +def test_volume_integral_carries_the_expression_value(setup): + mesh, x, y, T = setup + c = uw.function.expression(r"c_{v}", 2.0, "probe constant") + integral = uw.maths.Integral(mesh, c * T.sym[0].diff(y)) # 2 * int 2y = 2 + assert np.isclose(integral.evaluate(), 2.0, rtol=1e-8) + c.sym = 3.0 # a changed value, no recompile + assert np.isclose(integral.evaluate(), 3.0, rtol=1e-8) + assert np.isclose(uw.maths.Integral(mesh, c * T.sym[0].diff(y)).evaluate(), 3.0, rtol=1e-8) + + +def test_boundary_integral_carries_the_expression_value(setup): + mesh, x, y, T = setup + c = uw.function.expression(r"c_{b}", 2.0, "probe constant") + integral = uw.maths.BdIntegral(mesh, c * T.sym[0].diff(y), "Top") # 2 * 2 * length 1 + assert np.isclose(integral.evaluate(), 4.0, rtol=1e-8) + c.sym = 0.5 + assert np.isclose(integral.evaluate(), 1.0, rtol=1e-8) + + +def test_cellwise_integral_carries_the_expression_value(setup): + mesh, x, y, T = setup + c = uw.function.expression(r"c_{c}", 2.0, "probe constant") + cells = uw.maths.CellWiseIntegral(mesh, c * T.sym[0].diff(y)).evaluate() + assert np.isclose(np.asarray(cells).sum(), 2.0, rtol=1e-8) + + +def test_constitutive_flux_in_a_boundary_integral(setup): + """The case that found it: the viscous traction on a wall.""" + mesh, x, y, T = setup + v = uw.discretisation.MeshVariable("U_c", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("P_c", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, v, p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 2.0 + v.array[:, 0, :] = uw.function.evaluate(sympy.Matrix([[y ** 2, 0.0]]), v.coords).reshape(-1, 2) + sigma_xy = stokes.constitutive_model.flux[0, 1] # 2 eta (du/dy)/2 = 2y * 2 / ... = eta * 2y + assert np.isclose(uw.maths.BdIntegral(mesh, sigma_xy, "Top").evaluate(), 4.0, rtol=1e-8) + + +def test_a_constant_created_at_zero_still_reaches_the_kernel(setup): + """#696: a runtime constant whose value is zero at construction must not be + folded away by sympy (exp(c) with c.is_zero became 1 at construction, so a + ramp that started at t = 0 stayed frozen); setting it later must change the + value.""" + mesh, x, y, T = setup + c0 = uw.function.expression(r"c_{0}", 0.0, "starts at zero") + integrand = sympy.exp(c0) * T.sym[0].diff(y) + assert c0 in integrand.atoms(sympy.Symbol), "sympy folded exp(c) at construction" + integral = uw.maths.Integral(mesh, integrand) + assert abs(float(integral.evaluate()) - 1.0) < 1e-10 # int dT/dy = 1 on this box + c0.sym = 1.0 + assert abs(float(uw.maths.Integral(mesh, integrand).evaluate()) - np.e) < 1e-9 diff --git a/tests/test_0506_tensor_evaluate.py b/tests/test_0506_tensor_evaluate.py index b297e0192..6b543cc30 100644 --- a/tests/test_0506_tensor_evaluate.py +++ b/tests/test_0506_tensor_evaluate.py @@ -89,7 +89,7 @@ def test_navier_stokes_solve_does_not_trigger_ddt_fallback(): v = uw.discretisation.MeshVariable("u", mesh, mesh.dim, degree=2) p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1, continuous=True) - ns = uw.systems.NavierStokes( + ns = uw.systems.NavierStokesSLCN( mesh, velocityField=v, pressureField=p, rho=1.0, order=2 ) ns.constitutive_model = uw.constitutive_models.ViscousFlowModel diff --git a/tests/test_0610_navier_stokes_slcn_projection.py b/tests/test_0610_navier_stokes_slcn_projection.py index a849ebe63..f5892a780 100644 --- a/tests/test_0610_navier_stokes_slcn_projection.py +++ b/tests/test_0610_navier_stokes_slcn_projection.py @@ -33,7 +33,7 @@ def test_navier_stokes_slcn_solve_does_not_raise_shape_error(): v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=True) - ns = uw.systems.NavierStokes( + ns = uw.systems.NavierStokesSLCN( mesh, velocityField=v, pressureField=p, diff --git a/tests/test_0650_recursion_prevention_regression.py b/tests/test_0650_recursion_prevention_regression.py index aa698b289..82f47878c 100644 --- a/tests/test_0650_recursion_prevention_regression.py +++ b/tests/test_0650_recursion_prevention_regression.py @@ -139,7 +139,7 @@ def test_advection_diffusion_parameter_evaluation(self): temperature = uw.discretisation.MeshVariable("T", mesh, 1, degree=1) # Create advection-diffusion solver - adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=temperature, V_fn=velocity) + adv_diff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=temperature, V_fn=velocity) # Set constitutive model with UWexpression diffusivity (this was failing) adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel @@ -254,7 +254,7 @@ def test_estimate_dt_no_recursion(self): temperature = uw.discretisation.MeshVariable("T", mesh, 1, degree=1) # Create solver - adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=temperature, V_fn=velocity) + adv_diff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=temperature, V_fn=velocity) adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel adv_diff.constitutive_model.Parameters.diffusivity = uw.function.expression( r"\kappa", sym=1e-6 diff --git a/tests/test_0820_template_parameter_propagation.py b/tests/test_0820_template_parameter_propagation.py index 333bd85b3..d385e846c 100644 --- a/tests/test_0820_template_parameter_propagation.py +++ b/tests/test_0820_template_parameter_propagation.py @@ -188,7 +188,7 @@ def test_advdiff_diffusivity_parameter_propagation(self): with uw.synchronised_array_update(): v_soln.array[...] = 0.0 - adv_diff = uw.systems.AdvDiffusion( + adv_diff = uw.systems.AdvDiffusionSLCN( self.mesh, u_Field=phi, V_fn=v_soln, diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py new file mode 100644 index 000000000..b919504ac --- /dev/null +++ b/tests/test_1055_advdiff_supg_api.py @@ -0,0 +1,266 @@ +"""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.AdvDiffusion(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_with_the_slcn_defaults(mesh): + adv, _T = _solver(mesh, "a") + assert type(adv).__name__ == "SNES_AdvectionDiffusion_Composed" + # 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.EulerianSUPG) + # V_fn is data on the history manager: the velocity the transport uses + assert adv.DuDt.V_fn == adv.V_fn and adv.DuDt._advection_mode == "assembled" + + +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_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) + 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(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): + _solver(mesh, tag, **kwargs) + + +def test_timestep_is_required(mesh): + adv, _T = _solver(mesh, "b") + with pytest.raises(ValueError, match="needs a timestep"): + adv.solve() + + +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_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.DuDt._integrator = "am" + weights = adv.DuDt.spatial_weights() + assert len(weights) == 3 + states = adv.DuDt.states() + assert len(states) == 3 + # every history state appears (through its derivatives) in the advection operator + names = {str(atom.func) for atom in adv.DuDt.advection().atoms(sympy.Function)} + for s in states[1:]: + assert any(str(s[0].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_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 + 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.AdvDiffusion(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 + + +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_1056_navier_stokes_supg_api.py b/tests/test_1056_navier_stokes_supg_api.py new file mode 100644 index 000000000..2bc4dd32c --- /dev/null +++ b/tests/test_1056_navier_stokes_supg_api.py @@ -0,0 +1,130 @@ +"""API contract of the Eulerian SUPG Navier-Stokes solver. + +Structural checks that run in seconds: the export, argument validation, the +scheme assembled from the velocity history, the advecting-velocity switch and +the Picard passes, and the Stokes limit. + +Run: pixi run python -m pytest tests/test_1056_navier_stokes_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=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.2, qdegree=3) + + +def _cavity(mesh, tag, **kwargs): + """Lid-driven cavity: no-slip walls, a unit lid, unit viscosity.""" + v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, **kwargs) + ns.constitutive_model = uw.constitutive_models.ViscousFlowModel + ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + for b in ("Left", "Right", "Bottom"): + ns.add_dirichlet_bc((0.0, 0.0), b) + ns.add_dirichlet_bc((1.0, 0.0), "Top") + return ns, v, p + + +def test_exported_and_constructs_with_the_scalar_solver_rules(mesh): + ns, _v, _p = _cavity(mesh, "a", rho=1.0) + assert type(ns).__name__ == "SNES_NavierStokes_Composed" + assert ns.integrator == "am" and ns.order == 1 and ns.theta == 0.5 + assert isinstance(ns.DuDt, uw.systems.ddt.EulerianSUPG) + assert ns.DuDt.V_fn == ns._a_var.sym and ns.DuDt.V_fn_history[0] == ns.DuDt.psi_star[0].sym + assert ns.DFDt is None + assert _cavity(mesh, "b", order=2)[0].integrator == "bdf" + with pytest.raises(ValueError, match="theta applies"): + _cavity(mesh, "c", order=2, theta=0.5) + with pytest.raises(ValueError, match="stress history"): + _cavity(mesh, "d", DFDt=object()) + with pytest.raises(ValueError, match="advection must be"): + _cavity(mesh, "e", advection="upwind") + + +def test_a_step_runs_and_the_scheme_is_one_linear_solve(mesh): + ns, v, _p = _cavity(mesh, "s", rho=1.0) + ns.solve(timestep=0.05) + assert ns.snes.getIterationNumber() == 1 + assert np.isfinite(np.asarray(v.array)).all() + assert ns.picard_count == 0 + ns.solve(timestep=0.05, picard_iterations=3) + assert 1 <= ns.picard_count <= 3 + + +def test_stokes_limit_reproduces_the_stokes_solver(mesh): + """With rho -> 0 the momentum equation is the Stokes equation.""" + ns, v, p = _cavity(mesh, "z", rho=0.0) + ns.supg_weight = 0.0 + ns.tolerance = 1.0e-7 + ns.solve(timestep=1.0) + vs = uw.discretisation.MeshVariable("U_stokes", mesh, 2, degree=2) + ps = uw.discretisation.MeshVariable("P_stokes", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, vs, ps) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + for b in ("Left", "Right", "Bottom"): + stokes.add_dirichlet_bc((0.0, 0.0), b) + stokes.add_dirichlet_bc((1.0, 0.0), "Top") + stokes.tolerance = ns.tolerance + stokes.solve() + a, b = np.asarray(v.array).reshape(-1), np.asarray(vs.array).reshape(-1) + assert np.abs(a - b).max() < 1e-5 * np.abs(b).max() + + +def test_timestep_is_a_runtime_constant_and_theta_is_settable(mesh): + ns, _v, _p = _cavity(mesh, "t", rho=1.0) + ns.solve(timestep=0.05) + key = ns._current_jit_cache_key + ns.solve(timestep=0.02) + assert ns._current_jit_cache_key == key + ns.theta = 1.0 + ns.solve(timestep=0.02) + assert ns.theta == 1.0 and ns.DuDt.theta == 1.0 + assert ns._current_jit_cache_key == key + + +def test_tau_shapes_construct_and_step(mesh): + for shape in ("brooks_hughes", "doubly_asymptotic"): + ns, v, _p = _cavity(mesh, f"s_{shape}", rho=1.0, tau_shape=shape) + assert ns.tau_shape == shape + ns.solve(timestep=0.05) + assert np.isfinite(np.asarray(v.array)).all() and ns.snes.getIterationNumber() == 1 + with pytest.raises(ValueError, match="tau_shape"): + _cavity(mesh, "s_bad", tau_shape="optimal") + + +def test_peclet_weight_constructs_and_steps(mesh): + ns, v, _p = _cavity(mesh, "pe", rho=1.0, peclet_weight=2.0) + assert ns.peclet_weight == 2.0 + ns.solve(timestep=0.05) + assert np.isfinite(np.asarray(v.array)).all() and ns.snes.getIterationNumber() == 1 + + +def test_estimate_dt_carries_time_units_on_both_bases(): + """Under a model with reference scales both estimates come back as time + quantities (Copilot on #688: the accuracy basis returned a bare number + while the resolution fallback returned a quantity).""" + orchestration_model = uw.get_default_model() + orchestration_model.set_reference_quantities( + length=uw.quantity(1.0, "m"), time=uw.quantity(1.0, "s")) + try: + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + ns, v, _p = _cavity(mesh, "units", rho=1.0) + ns.solve(timestep=0.05) + before = ns.estimate_dt(basis="resolution") + after = ns.estimate_dt() + for dt in (before, after): + assert hasattr(dt, "dimensionality") and "[time]" in str(dt.dimensionality), dt + assert float(after.magnitude) > 0 + finally: + uw.reset_default_model() diff --git a/tests/test_1057_ddt_transport_plugin.py b/tests/test_1057_ddt_transport_plugin.py new file mode 100644 index 000000000..a3953a1a7 --- /dev/null +++ b/tests/test_1057_ddt_transport_plugin.py @@ -0,0 +1,221 @@ +"""The DDt history manager as the transport plugin of the Eulerian solvers. + +A solver composes its residual from three contributions of its history +manager (``time_derivative``, ``advection``, ``stabilisation_flux``) and +never asks which flavour it holds: the ``EulerianSUPG`` manager assembles +implicit advection with streamline-upwind stabilisation, the history-carrying +flavours answer zero for both. These checks cover the contract on each +flavour, the shapes for scalar, vector and tensor unknowns, the +semi-Lagrangian manager dropped into the SUPG solver, and a tensor unknown +transported through the multi-component solver. + +Run: pixi run python -m pytest tests/test_1057_ddt_transport_plugin.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.utilities._api_tools import Template + +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=1 / 16, qdegree=3) + + +def _gaussian(x, y, x0=0.5, y0=0.0, width=0.03): + return sympy.exp(-((x - x0) ** 2 + (y - y0) ** 2) / width) + + +def _is_zero(M): + return all(e == 0 for e in sympy.Matrix(M)) + + +def test_history_flavours_answer_zero_for_advection_and_stabilisation(mesh): + x, y = mesh.X + T = uw.discretisation.MeshVariable("T_c", mesh, 1, degree=2) + V = sympy.Matrix([[-y, x]]) + sl = uw.systems.ddt.SemiLagrangian( + mesh, T.sym, V, vtype=uw.VarType.SCALAR, degree=2, continuous=True, order=1) + assert sl.integrator == "am" + assert _is_zero(sl.advection()) and sl.advection().shape == (1, 1) + flux = sl.stabilisation_flux(sympy.Matrix([[7]])) + assert flux.shape == (1, 2) and _is_zero(flux) + td = sl.time_derivative() + assert td.shape == (1, 1) + assert sl.states()[1] == sl.psi_star[0].sym and len(sl.spatial_weights()) == 2 + # the timestep is a runtime constant the manager writes on every pre-solve + T.array[:, 0, 0] = uw.function.evaluate(_gaussian(x, y), T.coords).reshape(-1) + sl.update_pre_solve(0.02) + assert float(sl.delta_t.sym) == 0.02 + assert T.sym[0] in td.atoms(sympy.Function) and sl.psi_star[0].sym[0] in td.atoms(sympy.Function) + + eulerian = uw.systems.ddt.Eulerian( + mesh, T, vtype=uw.VarType.SCALAR, degree=2, continuous=True, V_fn=V) + assert eulerian._advection_mode == "split" # the velocity corrects the history + assert _is_zero(eulerian.advection()) and _is_zero(eulerian.stabilisation_flux(sympy.ones(1, 1))) + + +def test_supg_manager_contributions_have_the_unknowns_shape(mesh): + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + SUPG = uw.systems.ddt.EulerianSUPG + + T = uw.discretisation.MeshVariable("T_s", mesh, 1, degree=2) + scalar = SUPG(mesh, T, V, vtype=uw.VarType.SCALAR, degree=2, continuous=True) + assert scalar._advection_mode == "assembled" and scalar.V_fn == V + assert scalar.time_derivative().shape == (1, 1) and scalar.advection().shape == (1, 1) + R = scalar.time_derivative() + scalar.advection() + assert scalar.stabilisation_flux(R).shape == (1, 2) + # the advection is the velocity dotted with the gradient of each level, weighted + w0, w1 = scalar.spatial_weights() + expected = sum(w * V.dot(mesh.vector.gradient(level[0])) + for w, level in zip((w0, w1), scalar.states())) + assert sympy.simplify(scalar.advection()[0] - expected) == 0 + + U = uw.discretisation.MeshVariable("U_s", mesh, 2, degree=2) + vector = SUPG(mesh, U, V, vtype=uw.VarType.VECTOR, degree=2, continuous=True, order=2) + assert vector.integrator == "bdf" and vector.advection().shape == (1, 2) + R = sympy.Matrix([[sympy.Symbol("R_0"), sympy.Symbol("R_1")]]) + F = vector.stabilisation_flux(R) + assert F.shape == (2, 2) + assert sympy.simplify(F - vector.tau() * (R.T * V)) == sympy.zeros(2, 2) # F_ij = tau R_i a_j + # a self-advected unknown names the stored velocity at the stored levels + vector.V_fn_history = [ps.sym for ps in vector.psi_star] + assert vector.advecting_velocity(1) == vector.psi_star[0].sym + + S = uw.discretisation.MeshVariable("S_s", mesh, vtype=uw.VarType.SYM_TENSOR, degree=1) + tensor = SUPG(mesh, S, V, vtype=uw.VarType.SYM_TENSOR, degree=1, continuous=True) + assert tensor.advection().shape == (2, 2) + assert tensor.advection()[0, 1] == tensor.advection()[1, 0] + assert tensor.stabilisation_flux(tensor.advection()).shape == (4, 2) + + with pytest.raises(ValueError, match="tau_shape"): + SUPG(mesh, T, V, vtype=uw.VarType.SCALAR, degree=2, continuous=True, tau_shape="optimal") + with pytest.raises(ValueError, match="theta applies"): + SUPG(mesh, T, V, vtype=uw.VarType.SCALAR, degree=2, continuous=True, order=2, theta=0.5) + + +def test_semi_lagrangian_manager_drops_into_the_supg_solver(): + """With a semi-Lagrangian history the SUPG solver assembles no advection + and no stabilisation: on pure advection its equation is the one the + semi-Lagrangian solver solves, and the two fields agree to the solver + tolerances after a quarter revolution of a Gaussian.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 16, qdegree=3) + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + dt, steps = 0.1, 16 + + def field(tag): + T = uw.discretisation.MeshVariable(f"T_{tag}", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(_gaussian(x, y), T.coords).reshape(-1) + return T + + T_plug = field("plug") + history = uw.systems.ddt.SemiLagrangian( + mesh, T_plug.sym, V, vtype=uw.VarType.SCALAR, degree=2, continuous=True, order=1) + plug = uw.systems.AdvDiffusion(mesh, T_plug, V, DuDt=history) + assert plug.DuDt is history and plug.integrator == "am" and plug.order == 1 + assert _is_zero(plug.DuDt.advection()) and _is_zero(plug._stabilisation_flux()) + with pytest.raises(AttributeError): + plug.supg_weight # no stabilisation knobs on this manager + + T_slcn = field("slcn") + slcn = uw.systems.AdvDiffusionSLCN(mesh, T_slcn, V) + slcn.constitutive_model = uw.constitutive_models.DiffusionModel + slcn.constitutive_model.Parameters.diffusivity = 0.0 + + T_supg = field("supg") + supg = uw.systems.AdvDiffusion(mesh, T_supg, V) + + for solver in (plug, slcn, supg): + for b in ("Left", "Right", "Top", "Bottom"): + solver.add_dirichlet_bc(0.0, b) + for _ in range(steps): + plug.solve(timestep=dt) + slcn.solve(timestep=dt) + supg.solve(timestep=dt) + + a, b, c = (np.asarray(T.array[:, 0, 0]) for T in (T_plug, T_slcn, T_supg)) + assert np.abs(a - b).max() < 1e-5 * np.abs(b).max() # measured 5e-8 per step + # negative control: the assembled scheme is a different discretisation + assert np.abs(a - c).max() > 1e-3 + # and every scheme moved the Gaussian + T0 = uw.function.evaluate(_gaussian(x, y), T_plug.coords).reshape(-1) + assert np.abs(a - T0).max() > 0.3 + + +def test_tensor_unknown_is_transported_through_the_multicomponent_solver(): + """A flattened symmetric tensor (xx, xy, yy) carried by a uniform velocity + with the SUPG manager as the transport of a multi-component solver whose + residual is just the manager's terms.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 16, qdegree=3) + x, y = mesh.X + V = sympy.Matrix([[0.5, 0.0]]) + S = uw.discretisation.MeshVariable("S_t", mesh, (1, 3), vtype=uw.VarType.MATRIX, degree=2) + amplitudes = (1.0, 0.5, -1.0) + g0 = uw.function.evaluate(_gaussian(x, y, x0=-0.25), S.coords).reshape(-1) + for k, amp in enumerate(amplitudes): + S.array[:, 0, k] = amp * g0 + + transport = uw.systems.ddt.EulerianSUPG( + mesh, S, V, vtype=uw.VarType.MATRIX, degree=2, continuous=True, + num_components=(1, 3)) + + class TensorTransport(uw.systems.SNES_MultiComponent): + F0 = Template(r"f_0", lambda self: self.DuDt.time_derivative() + self.DuDt.advection(), + "time derivative and advection of every component") + F1 = Template(r"F_1", lambda self: self.DuDt.stabilisation_flux( + self.DuDt.time_derivative() + self.DuDt.advection()), "the SUPG flux per component") + + solver = TensorTransport(mesh, u_Field=S, DuDt=transport) + solver.constitutive_model = uw.constitutive_models.Constitutive_Model + solver.petsc_options["snes_rtol"] = 1e-8 + solver.petsc_options["ksp_rtol"] = 1e-9 + dt, steps = 0.05, 10 + for _ in range(steps): + transport.update_pre_solve(dt) + solver.solve() + transport.update_post_solve(dt) + assert float(transport.delta_t.sym) == dt + + exact = uw.function.evaluate(_gaussian(x, y, x0=-0.25 + 0.5 * dt * steps), S.coords).reshape(-1) + data = np.asarray(S.array) + for k, amp in enumerate(amplitudes): + err = np.linalg.norm(data[:, 0, k] - amp * exact) / np.linalg.norm(amp * exact) + assert err < 0.05, (k, err) + # negative control: the field moved away from where it started + assert np.linalg.norm(data[:, 0, k] - amp * g0) / np.linalg.norm(amp * g0) > 0.3 + + +def test_dimensional_timestep_reaches_the_manager_non_dimensional(): + """A quantity timestep handed to a manager's pre-solve is scaled by the model's + reference time before it becomes the kernels' runtime constant (#701): the + semi-Lagrangian solver passes its Pint step straight through.""" + from underworld3.systems.ddt import _as_float + q = uw.quantity(100.0, "kyr") + assert _as_float(q) == 100.0 # no reference scales: the magnitude + orchestration_model = uw.get_default_model() + orchestration_model.set_reference_quantities( + length=uw.quantity(1000.0, "km"), time=uw.quantity(1.0, "Myr")) + try: + assert abs(_as_float(q) - 0.1) < 1e-12 + assert abs(_as_float(q._pint_qty) - 0.1) < 1e-12 # a raw Pint quantity too + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + x, y = mesh.X + T = uw.discretisation.MeshVariable("T_dim", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(_gaussian(x, y), T.coords).reshape(-1) + slcn = uw.systems.AdvDiffusionSLCN(mesh, T, sympy.Matrix([[-y, x]])) + slcn.constitutive_model = uw.constitutive_models.DiffusionModel + slcn.constitutive_model.Parameters.diffusivity = 0.0 + slcn.solve(timestep=q) + assert abs(float(slcn.DuDt.delta_t.sym) - 0.1) < 1e-12, slcn.DuDt.delta_t.sym + finally: + uw.reset_default_model() 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) diff --git a/tests/test_1100_AdvDiffCartesian.py b/tests/test_1100_AdvDiffCartesian.py index 7be331e16..28cba14aa 100644 --- a/tests/test_1100_AdvDiffCartesian.py +++ b/tests/test_1100_AdvDiffCartesian.py @@ -128,7 +128,7 @@ def test_advDiff_boxmesh(mesh_type): # #### Create the advDiff solver - adv_diff = uw.systems.AdvDiffusion( + adv_diff = uw.systems.AdvDiffusionSLCN( mesh, u_Field=T, V_fn=v, 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..377e9ab3d --- /dev/null +++ b/tests/test_1100_advdiff_supg_rotating_gaussian.py @@ -0,0 +1,130 @@ +"""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: 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 +""" +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, 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.AdvDiffusion(mesh, T, sympy.Matrix([[-y, x]]), + order=order, 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, 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}") + 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(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_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", 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()) diff --git a/tests/test_1110_advDiffAnnulus.py b/tests/test_1110_advDiffAnnulus.py index 70d9d1d1b..0dd5a68a7 100644 --- a/tests/test_1110_advDiffAnnulus.py +++ b/tests/test_1110_advDiffAnnulus.py @@ -41,7 +41,7 @@ def test_adv_diff_annulus(): r_o = 1.0 delta_t = 0.05 ## 1/20 rotation in one step - adv_diff = uw.systems.AdvDiffusion( + adv_diff = uw.systems.AdvDiffusionSLCN( mesh, u_Field=t_soln, V_fn=v_soln,