diff --git a/docs/advanced/eulerian-advection-diffusion.md b/docs/advanced/eulerian-advection-diffusion.md new file mode 100644 index 000000000..757455cb2 --- /dev/null +++ b/docs/advanced/eulerian-advection-diffusion.md @@ -0,0 +1,123 @@ +# Eulerian advection-diffusion (SUPG): a drop-in for SLCN + +`uw.systems.AdvDiffusionSUPG` solves the same scalar transport equation as the +semi-Lagrangian solver `uw.systems.AdvDiffusionSLCN`, + +$$ +\frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f , +$$ + +but assembles every term on the mesh, implicit in time, with streamline-upwind +(SUPG) stabilisation. There is no trace-back and no departure point. The two +classes share their interface, so switching is one line: + +```python +adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=1) # was AdvDiffusionSLCN +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = 1.0e-3 +adv.add_dirichlet_bc(1.0, "Bottom") +adv.add_dirichlet_bc(0.0, "Top") + +dt = adv.estimate_dt() # accuracy-based: 2% of the field's range per step +adv.solve(timestep=dt) +``` + +The one deliberate difference is the timestep estimate. The semi-Lagrangian +`estimate_dt` reports the cell-crossing time, which for this solver is neither +a stability limit nor an accuracy one. The Eulerian solver's `estimate_dt` +instead returns the step at which the field changes by a given fraction of its +range (0.02 by default), from the advective rate before the first solve and +from the rate the last step actually produced after it. It does not depend on +the mesh, so cells refined for the Stokes problem do not shrink it. A script +that sizes its step in Courant numbers can still ask for +`estimate_dt(basis="resolution")`. + +## What carries over + +| SLCN | SUPG | note | +|---|---|---| +| `order=1, theta=0.5` | same | Crank-Nicolson, the default for both | +| `order=1, theta=1.0` | same | backward Euler | +| `order=2, theta=1.0` | same | SL-BDF2 becomes BDF2 | +| `order=2, theta=0.5` | refused | refused for the same reason: a BDF stencil does not pair with a centred flux | +| `f`, `V_fn`, `constitutive_model`, `delta_t` | same | | +| `estimate_dt()` | accuracy-based by default | the field may change by `fraction` (0.02) of its range per step; `basis="resolution"` returns the cell-crossing time SLCN reports | +| `solve(zero_init_guess, timestep, ...)` | same | | +| `DuDt.set_initial_history(values, dt)` | same | plant an exact history to start at full order | +| `restore_points_func`, `monotone_mode`, `old_frame_traceback`, `DFDt` | ignored, with a warning | they configure the trace-back | + +`order=3` (BDF3) is available; see below for when it is safe. + +## When to use which + +Both solvers are free of any stability limit on the timestep, so cells refined +for the Stokes problem never dictate the transport step. They differ in what +bounds their accuracy and in what a step costs. + +**Eulerian SUPG.** The error is set by how far the transported feature moves per +step relative to its own width, as $(\mathbf{u}\Delta t)^2$ for the second-order +schemes. It does not depend on the cell size at all: on a rotating Gaussian a band +refined to $h/9$, with its cells at a local Courant number of 13, changes the error +in the third digit only. A step costs one nonsymmetric solve, four to six times +less than a semi-Lagrangian step in serial, and it needs no departure points in +parallel. On a moving mesh the field and its history are re-interpolated by the +ordinary remesh transfer, so no special staging is needed. + +**Semi-Lagrangian.** The error is nearly independent of the timestep but +accumulates one interpolation per step, so at small Courant numbers it is the +worse scheme (21% against 0.6% after one revolution at Courant 0.5 on the same +mesh). Its limit is the arc a characteristic turns per step, about 10 degrees for +the RK2 trace-back, a property of the flow rather than the mesh. Above roughly +Courant 2 on the feature's own scale it keeps its accuracy where the Eulerian +scheme loses it. + +In parallel the Eulerian step stays seven times cheaper in serial and about five +times at eight ranks (the departure-point search parallelises perfectly, the +ILU preconditioner a little less), and its answer is identical to ten digits +at every rank count, where the semi-Lagrangian answer moves with the partition. + +A practical rule: if the timestep is chosen so that the temperature field itself +is resolved in time (a fraction of a feature width per step), the Eulerian solver +is cheaper and more accurate; if the step is deliberately long relative to the +transported features, the semi-Lagrangian solver is the one that survives it. + +## Choosing the time scheme + +Measured on a rotating Gaussian, one revolution, relative $L_2$ error; the full +tables are in the design note. + +| scheme | behaviour | +|---|---| +| Crank-Nicolson (`order=1`) | three to four times more accurate than BDF2 at the same timestep below Courant 2; rings once the feature is under-resolved in time | +| BDF2 (`order=2`) | damped and stable at every Courant number; the choice for sharp or under-resolved fields | +| BDF3 (`order=3`) | the most accurate scheme below Courant 1 when diffusion is present; on pure advection it grows slowly at any Courant number, so use it only with diffusion | +| backward Euler (`order=1, theta=1.0`) | 20 to 40% error at any practical timestep; not for transport | +| Adams-Moulton 2, 3 (not offered) | third and fourth order below Courant 1 but blow up on advection from about Courant 1, which is why there is no knob for them | + +All schemes cost the same per step: the history terms are extra kernel inputs, +not extra solves. Changing the timestep between steps changes a runtime constant +of the compiled kernels; nothing is recompiled. + +## Details that differ from SLCN + +- The strong residual used in the SUPG term carries the time derivative and the + advection but no diffusion term, because PETSc's pointwise kernels see first + derivatives only. For linear elements the missing term is identically zero. +- The stabilisation parameter uses the local cell size (`mesh.cell_size()`) and + three weights that are runtime constants (`solver.tau_weights`); + `solver.supg_weight = 0` gives the plain Galerkin scheme for comparison. +- The linear system is nonsymmetric, so the solver uses GMRES with an + additive-Schwarz ILU preconditioner, with the Krylov tolerance matched to the + SNES tolerance so that a step is one Newton iteration. Measured, this is the + cheaper solve at every Courant number up to eight ranks and its iteration + count does not grow with the rank count. `solver.preconditioner = "fmg"` + switches to geometric multigrid over the mesh's refinement hierarchy + (`refinement >= 1`) for very large rank counts. Every option can be + overridden through `solver.petsc_options`. + +## Further reading + +- Design note and measurements: `docs/developer/design/eulerian-supg-transport.md` +- The semi-Lagrangian schemes: {doc}`semi-lagrangian-time-integration` +- Example: `docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py` diff --git a/docs/advanced/index.md b/docs/advanced/index.md index cb47f8dea..be59f652e 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 +level-set-transport porous-flow snapshot-restore troubleshooting diff --git a/docs/advanced/level-set-transport.md b/docs/advanced/level-set-transport.md new file mode 100644 index 000000000..9f8cb8b3b --- /dev/null +++ b/docs/advanced/level-set-transport.md @@ -0,0 +1,93 @@ +# Conservative level sets + +`uw.systems.LevelSetSolver` carries a material interface as the 0.5 contour of a +smoothed indicator + +$$ +\psi = \tfrac12\left(1 + \tanh\frac{\varphi}{2\varepsilon}\right), +$$ + +where $\varphi$ is the signed distance to the interface (positive inside) and +$\varepsilon$ the interface thickness, a fraction of the local cell size. The +field is transported by an ordinary scalar solver; what makes it a level set is +what happens after each step: + +- **reinitialisation** restores the $\tanh$ profile without moving the 0.5 + contour (Parameswaran and Mandal 2023, integrated in pseudo-time with SSP-RK3); +- **mass correction** restores the enclosed volume by a uniform, clipped shift + found by bisection (Zhang, Zou and Greaves 2010). + +Neither depends on the transport scheme, so the solver takes either the Eulerian +SUPG solver (the default) or the semi-Lagrangian one. + +```python +from underworld3.systems import level_set + +psi = uw.discretisation.MeshVariable("psi", mesh, 1, degree=2) +eps = level_set.interface_thickness(mesh, psi, scale=0.35) +level_set.initialise_psi(psi, eps, interface_geometry="polygon", + interface_coordinates=circle_points) # or signed_distance=... + +ls = uw.systems.LevelSetSolver(psi, velocity=v.sym, epsilon=eps) # advection="slcn" to compare +for step in range(n_steps): + ls.solve(dt) # advect, reinitialise when due, restore the volume + +viscosity = level_set.material_property_field(psi.sym[0], [eta_outside, eta_inside], "geometric") +``` + +## Choices + +| argument | meaning | +|---|---| +| `advection` | `"supg"` (default) or `"slcn"`; both run pure advection | +| `order`, `theta` | the transport solver's time scheme; Crank-Nicolson by default, which preserves the profile's amplitude between reinitialisations | +| `reini_frequency`, `reini_steps`, `reini_dt` | how often, how many pseudo-time steps, and how long each is (half the smallest $\varepsilon$ by default) | +| `far_field` | the value of $\psi$ imposed on the domain boundary; set it whenever the flow crosses the boundary (an inflow boundary with no value lets mass in) | +| `conserve_mass` | `"auto"` (default): the global correction is on for `"slcn"`, which loses volume by interpolation, and off for `"supg"`, which conserves it to solver tolerance on its own; the clip to [0, 1] of the ringing at a one-cell band then costs about 0.2% per revolution, which `volume_drift` reports | +| `adv_solver_bc` | box wall labels on which a zero normal gradient is imposed by copying the neighbouring interior nodes | + +**Band thickness.** `interface_thickness(scale=0.35)`, the g-adopt default, +gives a band well under one cell, which a continuous-Galerkin transport rings +at. Measured on a rotating circle at 32 cells across, one revolution, SUPG with +no mass correction: + +| `scale` | $\varepsilon / h$ | volume drift | +|---|---|---| +| 0.35 | 0.12 | +0.84% (clipped ringing) | +| 1.0 | 0.36 | +0.28% | +| 2.0 | 0.71 | -0.18% (ringing gone) | +| 3.0 | 1.07 | -0.85% (reinitialisation curvature error) | + +For the SUPG transport a `scale` of 1.5 to 2, a band of two to three cells, is +the sensible setting; the thickness trades interface resolution for a clean +transport. + +`initialise_psi` accepts a precomputed signed distance, or a polygon, curve or +`shapely` geometry (the latter three need the optional `shapely` package). +`material_property_field` blends a property across one or more level sets with +a sharp, arithmetic, geometric or harmonic transition. + +## Cost + +Per step at 64 by 64 (LeVeque flow, Courant 0.5, reinitialisation every fifth +step): the SUPG advection takes 0.13 s and the SLCN advection 1.24 s; the +reinitialisation 0.06 to 0.11 s averaged; the mass correction 0.13 to 0.24 s. +Since the Eulerian transport does not need the correction, its level-set step +costs about 0.19 s against 1.59 s for the semi-Lagrangian one. + +## Which transport solver + +On the LeVeque swirling flow at 64 by 64 and Courant 0.5 (period 2), the SUPG +level set returns with a shape error of 0.028 against 0.051 for the +semi-Lagrangian one, at half the wall time; the mass correction pins both to +the same volume. The Eulerian solver's advantage is the same as for any scalar: +no interpolation loss per step, and cells refined for the Stokes problem cost +nothing. See {doc}`eulerian-advection-diffusion`. The example +`docs/examples/convection/advanced/Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py` runs the +comparison. + +## Credit + +The level-set pipeline, its SUPG transport and the LeVeque comparison are +NengLu's contribution (issue #657); this module unifies the two variants of +that work on the shared solver interface. diff --git a/docs/advanced/semi-lagrangian-time-integration.md b/docs/advanced/semi-lagrangian-time-integration.md index 23a935db2..f48c36ed7 100644 --- a/docs/advanced/semi-lagrangian-time-integration.md +++ b/docs/advanced/semi-lagrangian-time-integration.md @@ -112,6 +112,19 @@ $[\theta,\,1-\theta]$: `theta` is settable after construction: `adv_diff.DFDt.theta = 1.0`. +## The Eulerian alternative + +`uw.systems.AdvDiffusionSUPG` solves the same equation without a trace-back: +all terms are assembled on the mesh, implicit in time, with SUPG +stabilisation. Its `order=` and `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/developer/design/eulerian-supg-transport.md b/docs/developer/design/eulerian-supg-transport.md new file mode 100644 index 000000000..81dce0202 --- /dev/null +++ b/docs/developer/design/eulerian-supg-transport.md @@ -0,0 +1,343 @@ +# 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). + +## Parallel + +LeVeque flow, conservative level set, 20 steps at 128 by 128 (16,384 cells) and 10 at +256 by 256 (65,536 cells), wall time per advection step max-reduced over ranks; +reinitialisation every fifth step timed separately (`~/+Simulations/supg_vs_slcn_657/parallel/`). + +| ranks | SUPG 128 | SLCN 128 | SUPG 256 | SLCN 256 | +|---|---|---|---|---| +| 1 | 0.55 s | 3.98 s | 2.22 s | 15.7 s | +| 2 | 0.28 s | 1.77 s | | | +| 4 | 0.145 s | 1.17 s | 0.80 s | 4.12 s | +| 8 | 0.126 s | 1.06 s | 0.67 s | 3.50 s | + +Three observations. + +- Per step the Eulerian solve is seven times cheaper in serial at both sizes; the + gap narrows to about five times at eight ranks, because the departure-point + work of the semi-Lagrangian scheme parallelises perfectly while the + additive-Schwarz ILU preconditioner needs more GMRES iterations as its + subdomains shrink (SUPG speed-up 3.3 at eight ranks on 256 by 256 against 4.5 + for SLCN). These SUPG rows were taken with the Krylov tolerance at its default + and so carry a second Newton step (see "Preconditioner" above); with the + tolerance matched, the Schwarz iteration count on the two-cell band is the + same at one and eight ranks and geometric multigrid does not beat it. The + assembly itself scales. +- The Eulerian answer is partition-independent: the enclosed volume agrees to + ten digits at every rank count. The semi-Lagrangian answer is not: it moves in + the sixth digit at two and four ranks and by 1.6% at eight ranks on the + 128 by 128 mesh (0.06957 against 0.07068), which points at departure points + near partition boundaries being sampled wrongly at higher rank counts. That + is a defect in the semi-Lagrangian trace-back to chase separately; the + Eulerian scheme has no such path. +- The 128 by 128 problem is too small for eight ranks (about 2,000 cells each); + the 256 by 256 rows are the ones to read for scaling. + +These are timings at equal step. The fair cost comparison is per unit of simulated +time at equal error, each scheme at its own accuracy-limited step (the field-change +fraction for SUPG, the trace-back arc for SLCN), which is the next measurement. + +## 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 c32976e1f..26a102965 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -165,6 +165,7 @@ design/TURBULENCE_MODEL_DESIGN design/declined-coord-units-proposal design/nonlinear-solver-homotopy-warmstart design/fault-zone-hybrid-architecture +design/eulerian-supg-transport ``` ```{toctree} diff --git a/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py b/docs/examples/convection/advanced/Ex_AdvectionDiffusionSUPG_RotationTest.py new file mode 100644 index 000000000..61da93c87 --- /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.AdvDiffusionSUPG`. The exact +solution is known at every time (`uw.analytic.RotatingGaussian`), so the +error is measured directly rather than inferred from a picture. + +The scheme is stable at any cell Courant number; what limits the timestep +is how far the anomaly moves per step relative to its own width, 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.AdvDiffusionSUPG( + 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_LevelSet_LeVeque_SUPG_vs_SLCN.py b/docs/examples/convection/advanced/Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py new file mode 100644 index 000000000..9eb482238 --- /dev/null +++ b/docs/examples/convection/advanced/Ex_LevelSet_LeVeque_SUPG_vs_SLCN.py @@ -0,0 +1,178 @@ +# --- +# 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] +""" +# Level set in the LeVeque swirling flow: SUPG against SLCN + +**PHYSICS:** convection +**DIFFICULTY:** advanced + +## Description + +The swirling deformation flow of LeVeque (1996), the standard stress test +for interface transport: a circle is stretched into a thin spiral filament +for half a period, then the flow reverses exactly and the circle should +come back. Any irreversible error, whether interpolation loss in a +trace-back or stabilisation diffusion, shows up as a failure to recover +the initial shape. + +The same conservative level set is carried by the two transport solvers, +Eulerian SUPG and semi-Lagrangian, under the same velocity and timestep, +each with its own reinitialisation and mass correction. The script reports +the shape error against the frozen initial field, the enclosed volume, and +the wall time of each. + +The stream function is + +$$\psi(x, y, t) = \frac{1}{\pi}\sin^2(\pi x)\,\sin^2(\pi y)\,\cos(\pi t / T)$$ + +with period `T`: 2 (LeVeque's own value, a gentle round trip) or 8 +(Enright et al. 2002, filaments thinner than the mesh). + +Contributed by NengLu (issue #657); converted to the repository's script +conventions. + +## Parameters + +- `uw_res`: cells across the unit square +- `uw_period`: reversal period `T` +- `uw_courant`: timestep as a multiple of the cell-crossing time +""" + +# %% +import os +import time + +import numpy as np +import sympy + +import underworld3 as uw +from underworld3.systems import level_set + +# %% +params = uw.Params( + uw_res=64, + uw_period=2.0, + uw_courant=0.5, + uw_reini_frequency=5, + uw_outdir="output/levelset_leveque", +) + +# %% [markdown] +""" +## Mesh and the time-dependent velocity +""" + +# %% +mesh = uw.meshing.StructuredQuadBox( + elementRes=(params.uw_res, params.uw_res), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)) +x, y = mesh.X + +v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2) + +stream = (1 / sympy.pi) * sympy.sin(sympy.pi * x) ** 2 * sympy.sin(sympy.pi * y) ** 2 +u_x, u_y = -sympy.diff(stream, y), sympy.diff(stream, x) + + +def set_velocity(t): + modulation = float(np.cos(np.pi * t / params.uw_period)) + v.array[:, 0, 0] = modulation * uw.function.evaluate(u_x, v.coords).reshape(-1) + v.array[:, 0, 1] = modulation * uw.function.evaluate(u_y, v.coords).reshape(-1) + + +# %% [markdown] +""" +## Two level sets, one initial circle, one per solver +""" + +# %% +radius, centre = 0.15, (0.5, 0.75) +angles = np.linspace(0.0, 2.0 * np.pi, 91) +circle = np.column_stack((centre[0] + radius * np.cos(angles), centre[1] + radius * np.sin(angles))) + +solvers = {} +for name in ("supg", "slcn"): + psi = uw.discretisation.MeshVariable(f"psi_{name}", mesh, 1, degree=2) + eps = level_set.interface_thickness(mesh, psi, scale=0.35) + level_set.initialise_psi(psi, eps, interface_geometry="polygon", interface_coordinates=circle) + psi0 = uw.discretisation.MeshVariable(f"psi0_{name}", mesh, 1, degree=2) + psi0.array[...] = psi.array[...] + solver = uw.systems.LevelSetSolver( + psi, velocity=v.sym, epsilon=eps, advection=name, + reini_steps=1, reini_frequency=params.uw_reini_frequency) + solvers[name] = dict(psi=psi, psi0=psi0, solver=solver, wall=0.0) + + +def shape_error(psi, psi0): + return float(np.sqrt(max(uw.maths.Integral(mesh, (psi.sym[0] - psi0.sym[0]) ** 2).evaluate(), 0.0))) + + +# %% [markdown] +""" +## Time loop + +Both solvers take the same step, chosen as a multiple of the cell-crossing +time so the comparison is at equal Courant number. +""" + +# %% +dt = params.uw_courant / params.uw_res +n_steps = int(np.round(params.uw_period / dt)) +dt = params.uw_period / n_steps +report_every = max(1, n_steps // 16) +initial_area = np.pi * radius ** 2 + +t = 0.0 +for step in range(n_steps): + set_velocity(t) + for name, s in solvers.items(): + t0 = time.perf_counter() + s["solver"].solve(dt) + s["wall"] += time.perf_counter() - t0 + t += dt + if step % report_every == 0 or step == n_steps - 1: + for name, s in solvers.items(): + volume = s["solver"].interface_volume() + uw.pprint(f"t = {t:6.3f} {name}: volume drift {100 * (volume - initial_area) / initial_area:+.3f}% " + f"shape error {shape_error(s['psi'], s['psi0']):.3e} wall {s['wall']:.1f} s") + +# %% [markdown] +""" +## Round trip + +At `t = T` the flow has returned the fluid to where it started; the shape +error measures what the transport did not undo. +""" + +# %% +for name, s in solvers.items(): + uw.pprint(f"{name}: round-trip shape error {shape_error(s['psi'], s['psi0']):.3e}, " + f"total wall {s['wall']:.1f} s") + +# %% +if uw.mpi.size == 1: + import pyvista as pv + import underworld3.visualisation as vis + + pl = pv.Plotter(window_size=(900, 450), shape=(1, 2)) + for i, (name, s) in enumerate(solvers.items()): + pvmesh = vis.mesh_to_pv_mesh(mesh) + pvmesh.point_data["psi"] = vis.scalar_fn_to_pv_points(pvmesh, s["psi"].sym) + pl.subplot(0, i) + pl.add_mesh(pvmesh, scalars="psi", cmap="RdBu_r", clim=(0, 1), show_edges=False) + pl.add_mesh(pvmesh.contour([0.5], scalars="psi"), color="black", line_width=2) + pl.add_text(name, font_size=10) + os.makedirs(params.uw_outdir, exist_ok=True) + pl.show(cpos="xy", screenshot=os.path.join(params.uw_outdir, "leveque_round_trip.png")) 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/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d3c67611b..c0e65cf82 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3896,13 +3896,21 @@ def update_lvec(self, swarm_sync=True): # The field decomposition seems to fail if coarse DMs are present names, isets, dms = self.dm.createFieldDecomposition() - # traverse subdms, taking user generated data in the subdm - # local vec, pushing it into a global sub vec - for var, subiset, subdm in zip(self.vars.values(), isets, dms): - # var.vec lazily creates the PETSc local vector on first access - lvec = var.vec + # Traverse the DM's fields BY NAME. `self.vars` holds its + # variables weakly, so a dropped-and-collected variable leaves + # a field behind in the DM; a positional zip would then pack + # every later variable into the wrong field (measured: the + # cell-size field landing in a P2 slot as garbage, NaN + # residuals in a solver that reads it). An orphaned field is + # zeroed so nothing stale can reach a kernel. + for name, subiset, subdm in zip(names, isets, dms): + var = self.vars.get(name) subvec = a_global.getSubVector(subiset) - subdm.localToGlobal(lvec, subvec, addv=False) + if var is None: + subvec.set(0.0) + else: + # var.vec lazily creates the PETSc local vector on first access + subdm.localToGlobal(var.vec, subvec, addv=False) a_global.restoreSubVector(subiset, subvec) for iset in isets: diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index ea6823ab5..3841861b1 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -19,6 +19,10 @@ L2 projection of fields onto mesh variables. AdvDiffusion : class Advection-diffusion with semi-Lagrangian transport. +AdvDiffusionSUPG : class + Advection-diffusion, implicit Eulerian with SUPG stabilisation. +LevelSetSolver : class + Conservative level-set transport (advection, reinitialisation, mass correction). NavierStokes : class Navier-Stokes equations with inertia. Diffusion : class @@ -64,6 +68,7 @@ # These are now implemented the same way using the ddt module from .solvers import SNES_AdvectionDiffusion as AdvDiffusionSLCN from .solvers import SNES_AdvectionDiffusion as AdvDiffusion +from .advection_diffusion_eulerian import SNES_AdvectionDiffusion_SUPG as AdvDiffusionSUPG # import diffusion-only solver from .solvers import SNES_Diffusion as Diffusion @@ -91,3 +96,6 @@ # δ-continuation driver for hard viscoplastic (Drucker–Prager) yield from .yield_continuation import yield_continuation, YieldHomotopyControl from .solve_report import SolveReport + +from . import level_set +from .level_set import LevelSetSolver diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py new file mode 100644 index 000000000..d294bdb35 --- /dev/null +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -0,0 +1,754 @@ +r"""Fully implicit Eulerian advection-diffusion with SUPG stabilisation. + +The scalar transport equation + +.. math:: + \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f + +discretised on the mesh with a linear multistep rule in time and a +streamline-upwind Petrov-Galerkin (SUPG) term in space. Every time level +is a mesh variable held by an :class:`~underworld3.systems.ddt.Eulerian` +history manager, so the scheme's order is a construction argument and the +timestep and multistep coefficients are runtime constants of the compiled +kernels: neither changes the generated code. + +The companion of :class:`~underworld3.systems.solvers.SNES_AdvectionDiffusion` +(semi-Lagrangian). The Eulerian scheme is stable at any cell Courant number +and its accuracy is set by how far the transported feature moves in one +step; the semi-Lagrangian scheme's accuracy is set by how far a characteristic +turns in one step. See ``docs/developer/design/eulerian-supg-transport.md``. + +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 Eulerian as Eulerian_DDt +from underworld3.systems.solvers import ( + _advective_diffusive_dt, + _dimensionalise_dt, + _invalidate_solution_cache, + _nondimensionalise_timestep, +) + + +def _as_row_vector(V_fn, dim): + """Coerce a velocity expression to a ``(1, dim)`` sympy row Matrix.""" + if isinstance(V_fn, uw.discretisation.MeshVariable): + V_fn = V_fn.sym + if isinstance(V_fn, sympy.MatrixBase): + if V_fn.shape == (1, dim): + return V_fn + if V_fn.shape == (dim, 1): + return V_fn.T + raise ValueError( + f"V_fn has shape {V_fn.shape} but the mesh is {dim}-D; expected a " + f"(1, {dim}) row vector such as `v.sym` of a vector MeshVariable." + ) + raise ValueError( + f"V_fn must be a (1, {dim}) sympy Matrix or a vector MeshVariable, " + f"not {type(V_fn).__name__}." + ) + + +class SNES_AdvectionDiffusion_SUPG(SNES_Scalar): + r"""Eulerian advection-diffusion solver, implicit in time, SUPG in space. + + .. math:: + \frac{\partial \phi}{\partial t} + \mathbf{u}\cdot\nabla\phi + - \nabla\cdot(\kappa\nabla\phi) = f + + A drop-in replacement for :class:`~underworld3.systems.solvers.SNES_AdvectionDiffusion` + (``uw.systems.AdvDiffusionSLCN``): the constructor, ``order``, ``theta``, + ``f``, ``V_fn``, ``constitutive_model``, ``delta_t``, ``estimate_dt`` and + ``solve`` all keep the semi-Lagrangian solver's meaning, so a script changes + the class name and nothing else:: + + adv = uw.systems.AdvDiffusionSUPG(mesh, T, v.sym, order=1) # was AdvDiffusionSLCN + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = 1.0e-3 + adv.add_dirichlet_bc(0.0, "Left") + adv.solve(timestep=dt) + + The arguments that only make sense for a trace-back + (``restore_points_func``, ``monotone_mode``, ``old_frame_traceback``, + ``DFDt``) are accepted and ignored with a warning. + + **Time schemes.** ``order`` and ``theta`` select the same schemes as for + the semi-Lagrangian solver: + + ========== ======= ===================================================== + ``order`` ``theta`` scheme + ========== ======= ===================================================== + 1 0.5 Crank-Nicolson (default; the SLCN convention) + 1 1.0 backward Euler + 2 1.0 BDF2, all spatial terms at :math:`n+1` (the SL-BDF2 convention) + 3 1.0 BDF3 + ========== ======= ===================================================== + + ``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 : Eulerian, optional + A pre-built history manager (order at least ``order``, no ``V_fn``). + restore_points_func, monotone_mode, old_frame_traceback, DFDt + Semi-Lagrangian arguments, accepted for drop-in compatibility and + ignored with a warning: there is no trace-back here. + + Notes + ----- + The diffusivity is set through the constitutive model, as for 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, + 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"AdvDiffusionSUPG ignores {', '.join(ignored)}: these configure " + "the semi-Lagrangian trace-back and the Eulerian scheme has none.", + stacklevel=2, + ) + order = int(order) + if order not in (1, 2, 3): + raise ValueError(f"order must be 1, 2 or 3, not {order}.") + # theta 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). + integrator = "am" if order == 1 else "bdf" + if theta != 1.0 and order != 1: + raise ValueError( + "theta applies at order 1 only (0.5 is Crank-Nicolson, 1.0 is " + "backward Euler); order 2 and 3 take theta=1.0, the same rule as " + "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._integrator = integrator + self._time_order = order + self._theta = theta + self._V_fn = _as_row_vector(V_fn, mesh.dim) + + tag = self.instance_number + self._delta_t = public_expression( + rf"\Delta t_{{{tag}}}", 1.0, "Eulerian advection-diffusion timestep") + self._last_timestep = None + self._last_change_rate = None + + # SUPG on/off and the three tau weights are runtime constants: the + # compiled kernels read them from PETSc's constants[] array. + self._supg_weight = public_expression( + rf"w^{{\mathrm{{SUPG}}}}_{{{tag}}}", 1.0, "SUPG term weight (0 = Galerkin)") + self._tau_weights = [ + public_expression(rf"C^{{\tau}}_{{t,{tag}}}", 2.0, "tau transient weight"), + public_expression(rf"C^{{\tau}}_{{u,{tag}}}", 2.0, "tau advective weight"), + public_expression(rf"C^{{\tau}}_{{\kappa,{tag}}}", 4.0, "tau diffusive weight"), + ] + + if DuDt is None: + self.Unknowns.DuDt = Eulerian_DDt( + self.mesh, + u_Field, + vtype=uw.VarType.SCALAR, + degree=u_Field.degree, + continuous=u_Field.continuous, + V_fn=None, + theta=theta, + varsymbol=u_Field.symbol, + verbose=verbose, + bcs=self.essential_bcs, + order=order, + smoothing=0.0, + ) + else: + if DuDt.order < order: + raise ValueError( + f"DuDt supplied is order {DuDt.order} but order {order} was requested." + ) + if getattr(DuDt, "V_fn", None) is not None: + raise ValueError( + "DuDt must be built with V_fn=None: advection is assembled " + "implicitly by this solver, not as an explicit history correction." + ) + self.Unknowns.DuDt = DuDt + + # Diffusivity lives on the constitutive model, as for every scalar + # solver; kappa = 0 until the user sets it. + self.constitutive_model = uw.constitutive_models.DiffusionModel + self.constitutive_model.Parameters.diffusivity = 0.0 + + # 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._time_order), + f"{self._integrator.upper()} order {self._time_order}") + display(Latex(r"$\quad\mathrm{u} = $ " + self.u.sym._repr_latex_())) + display(Latex(r"$\quad\mathbf{v} = $ " + self._V_fn._repr_latex_())) + display(Latex(r"$\quad\Delta t = $ " + self._delta_t._repr_latex_())) + display(Latex(rf"$\quad$ time scheme: {scheme}")) + + # ------------------------------------------------------------------ + # Scheme description + # ------------------------------------------------------------------ + + @property + def integrator(self) -> str: + """The multistep family in use: ``"am"`` (the theta rule) at order 1, ``"bdf"`` above.""" + return self._integrator + + @property + def order(self) -> int: + """Requested order of the time integration.""" + return self._time_order + + @property + def theta(self) -> float: + """Adams-Moulton blend at order 1 (1.0 backward Euler, 0.5 Crank-Nicolson). + + Settable after construction, as on the semi-Lagrangian solver: the + blend is a runtime constant of the compiled kernels, refreshed from + the history manager before every solve, so nothing is recompiled. + """ + return self._theta + + @theta.setter + def theta(self, value): + value = float(value) + if value != 1.0 and self._time_order != 1: + raise ValueError( + "theta applies at order 1 only (0.5 is Crank-Nicolson, 1.0 is " + "backward Euler); order 2 and 3 take theta=1.0." + ) + self._theta = value + self.DuDt.theta = value + + @property + def delta_t(self): + r"""The timestep :math:`\Delta t` as a UW expression. + + Set by :meth:`solve`, or assign it directly (a number or a quantity + with time units) and call ``solve()`` without ``timestep``, as with + the semi-Lagrangian solver. A new value updates a runtime constant of + the compiled kernels; nothing is recompiled. + """ + return self._delta_t + + @delta_t.setter + def delta_t(self, value): + dt = float(_nondimensionalise_timestep(value)) + if dt <= 0.0: + raise ValueError(f"timestep must be positive, not {dt}.") + if dt != self._last_timestep: + self._delta_t.sym = dt + self._last_timestep = dt + + @property + def V_fn(self): + """Advecting velocity, ``(1, dim)``.""" + return self._V_fn + + @V_fn.setter + def V_fn(self, value): + self._V_fn = _as_row_vector(value, self.mesh.dim) + self.is_setup = False + + @property + def f(self): + """Volumetric source term.""" + return self._f + + @f.setter + def f(self, value): + self._f = sympy.Matrix((value,)) + self._needs_function_rewire = True + + @property + def supg_weight(self) -> float: + """Scale of the SUPG term: 1 (default) or 0 for plain Galerkin. No rebuild.""" + return float(self._supg_weight.sym) + + @supg_weight.setter + def supg_weight(self, value): + self._supg_weight.sym = float(value) + + @property + def tau_weights(self): + r"""The weights :math:`(C_t, C_u, C_\kappa)` of the three terms in :math:`\tau`.""" + return tuple(float(w.sym) for w in self._tau_weights) + + @tau_weights.setter + def tau_weights(self, values): + ct, cu, ck = (float(v) for v in values) + self._tau_weights[0].sym = ct + self._tau_weights[1].sym = cu + self._tau_weights[2].sym = ck + + # ------------------------------------------------------------------ + # Residual pieces (raw field symbols only, so the Jacobian sees them) + # ------------------------------------------------------------------ + + def _states(self): + r"""``[phi^{n+1}, phi^{n}, phi^{n-1}, ...]`` as scalar field symbols.""" + return [self.u.sym[0]] + [ps.sym[0] for ps in self.DuDt.psi_star] + + def _spatial_weights(self): + """Weight of the spatial operator at each time level of ``_states``.""" + n = len(self.DuDt.psi_star) + if self._integrator == "bdf": + return [sympy.Integer(1)] + [sympy.Integer(0)] * n + return self.DuDt.am_coefficient_expressions[: n + 1] + + def _time_derivative(self): + if self._integrator == "bdf": + return self.DuDt.bdf()[0] / self._delta_t + phi_new, phi_old = self._states()[:2] + return (phi_new - phi_old) / self._delta_t + + def _advection(self): + dim = self.mesh.dim + u = self._V_fn + total = sympy.Integer(0) + for w, phi in zip(self._spatial_weights(), self._states()): + if w == 0: + continue + grad = self.mesh.vector.gradient(phi) + total = total + w * sum(u[0, i] * grad[0, i] for i in range(dim)) + return total + + def _diffusive_flux(self): + r"""``(1, dim)`` flux :math:`\sum_k w_k\,\nabla\phi^{(k)}\cdot\kappa` from the constitutive tensor.""" + dim = self.mesh.dim + c = self.constitutive_model.c + total = sympy.zeros(1, dim) + for w, phi in zip(self._spatial_weights(), self._states()): + if w == 0: + continue + grad = self.mesh.vector.gradient(phi) + total = total + w * (grad * c) + return total + + def _strong_residual(self): + return self._time_derivative() + self._advection() - self._f[0] + + def _scalar_diffusivity(self): + kappa = self.constitutive_model.Parameters.diffusivity + if isinstance(kappa, sympy.MatrixBase): + raise ValueError( + "The SUPG parameter needs a scalar diffusivity; anisotropic " + "diffusion is not supported by this solver." + ) + return kappa + + def _tau(self): + dim = self.mesh.dim + u = self._V_fn + u_mag2 = sum(u[0, i] ** 2 for i in range(dim)) + h = self.mesh.cell_size() + kappa = self._scalar_diffusivity() + if self._integrator == "bdf": + c0 = self.DuDt.bdf_coefficient_expressions[0] + else: + c0 = sympy.Integer(1) + ct, cu, ck = self._tau_weights + transient = (ct * c0 / self._delta_t) ** 2 + advective = (cu * sympy.sqrt(u_mag2) / h) ** 2 + diffusive = (ck * kappa / h ** 2) ** 2 + return self._supg_weight / sympy.sqrt(transient + advective + diffusive + 1.0e-30) + + F0 = Template( + r"f_0(\phi)", + lambda self: sympy.Matrix([[self._strong_residual()]]), + "Strong residual of the time scheme: time derivative, advection and source.", + ) + F1 = Template( + r"\mathbf{F}_1(\phi)", + lambda self: self._diffusive_flux() + self._tau() * self._strong_residual() * self._V_fn, + "Diffusive flux of the time scheme plus the SUPG flux tau R u.", + ) + + # ------------------------------------------------------------------ + # Timestep and solve + # ------------------------------------------------------------------ + + @timing.routine_timer_decorator + def estimate_dt(self, 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) + super().solve(zero_init_guess, _force_setup, divergence_retries=divergence_retries) + _invalidate_solution_cache(self.u) + # The realised rate of change of the field over this step feeds the + # accuracy-based estimate_dt; psi_star[0] still holds phi^n here. + from mpi4py import MPI + change = np.abs(np.asarray(self.u.array).reshape(-1) + - np.asarray(self.DuDt.psi_star[0].array).reshape(-1)) + local = float(change.max()) if change.size else 0.0 + self._last_change_rate = uw.mpi.comm.allreduce(local, op=MPI.MAX) / dt + self.DuDt.update_post_solve(dt, verbose=verbose) + + self.is_setup = True + self.constitutive_model._solver_is_setup = True diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e198aad57..a37533cbe 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -666,6 +666,29 @@ def bdf_coefficients(self): """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + @property + def bdf_coefficient_expressions(self): + r"""The BDF coefficient symbols :math:`[c_0, c_1, \dots]` as UWexpressions. + + For a solver that assembles its own weighted sum of history terms + (an Eulerian scheme applying the multistep rule to a spatial + operator, say). The symbols are routed through PETSc's + ``constants[]`` array, so their values follow ``effective_order`` + and the timestep without a recompile; ``bdf_coefficients`` gives + the current values. + """ + return list(self._bdf_coeffs) + + @property + def am_coefficient_expressions(self): + r"""The Adams-Moulton coefficient symbols :math:`[a_0, a_1, \dots]` as UWexpressions. + + :math:`a_0` weights the new state, :math:`a_k` the history slot + ``psi_star[k-1]``. Same constants-routing as + :attr:`bdf_coefficient_expressions`. + """ + return list(self._am_coeffs) + def _history_syms(self): """History terms as sympy expressions for the weighted sums. diff --git a/src/underworld3/systems/level_set.py b/src/underworld3/systems/level_set.py new file mode 100644 index 000000000..7d4a5e5da --- /dev/null +++ b/src/underworld3/systems/level_set.py @@ -0,0 +1,672 @@ +r"""Conservative level-set transport: advection, reinitialisation, mass correction. + +A conservative level set represents an interface by the 0.5 contour of a +smoothed indicator + +.. math:: + \psi = \tfrac12\left(1 + \tanh\frac{\varphi}{2\varepsilon}\right), + +with :math:`\varphi` the signed distance to the interface and +:math:`\varepsilon` the interface thickness. Transport of :math:`\psi` is +ordinary scalar advection; what makes it a level set is what happens between +steps: a reinitialisation that restores the :math:`\tanh` profile without +moving the 0.5 contour, and a global correction that restores the enclosed +volume. Both are post-step operations on the field, so any scalar transport +solver can carry it. + +:class:`LevelSetSolver` takes the Eulerian SUPG solver +(:class:`~underworld3.systems.AdvDiffusionSUPG`) by default and the +semi-Lagrangian solver on request. The reinitialisation equation is that of +Parameswaran and Mandal (2023), integrated in pseudo-time with SSP-RK3 +(Gottlieb and Shu 1998); the mass correction is the uniform shift of Zhang, +Zou and Greaves (2010). The signed-distance helpers accept a polygon or a +curve (they use ``shapely``, an optional dependency) or a precomputed +distance array. :func:`material_property_field` blends material properties +across the interface, in the manner of g-adopt's ``field_interface``. + +The level-set pipeline, the SUPG transport it drove and the LeVeque +swirling-flow comparison are NengLu's (issue #657, branch ``levelset``); +this module unifies the two variants of that work on one solver interface. + +References +---------- +Parameswaran, S. and Mandal, J. C. (2023). A stable interface-preserving +reinitialization equation for conservative level set method. European +Journal of Mechanics B/Fluids, 98, 40-63. + +Zhang, Y., Zou, Q. and Greaves, D. (2010). Numerical simulation of +free-surface flow using the level-set method with global mass correction. +International Journal for Numerical Methods in Fluids, 63, 651-680. +""" + +import warnings +from typing import Optional, Union + +import numpy as np +import sympy + +import underworld3 as uw +from underworld3 import discretisation, systems + + +# --------------------------------------------------------------------------- +# Initial condition +# --------------------------------------------------------------------------- + +def _tanh_profile(distance, epsilon): + r"""The conservative level-set profile :math:`\tfrac12(1 + \tanh(\varphi/2\varepsilon))`.""" + return (1.0 + np.tanh(np.asarray(distance) / (2.0 * np.asarray(epsilon)))) / 2.0 + + +def _shapely(): + """The optional ``shapely`` dependency, or a clear error.""" + try: + from shapely import geometry + from shapely import prepare + except ImportError as exc: + raise ImportError( + "The level-set geometry helpers need the optional package " + "'shapely' (pip install shapely). Pass signed_distance= to " + "initialise_psi to avoid it." + ) from exc + return geometry, prepare + + +def initialise_psi( + psi: discretisation.MeshVariable, + epsilon, + *, + signed_distance=None, + interface_geometry: Optional[str] = None, + interface=None, + interface_coordinates=None, + boundary_coordinates=None, +) -> None: + r"""Fill ``psi`` with the conservative level-set profile of an interface. + + Parameters + ---------- + psi : MeshVariable + Scalar field to fill; 1 inside the interface, 0 outside, 0.5 on it. + epsilon : MeshVariable or float + Interface thickness (see :func:`interface_thickness`). + signed_distance : ndarray, optional + Precomputed signed distance at ``psi``'s nodes, positive inside. + When given, the geometry arguments are ignored and ``shapely`` is + not needed. + interface_geometry : {"curve", "polygon", "shapely"}, optional + How the interface is described. + interface : shapely LineString or Polygon, optional + For ``interface_geometry="shapely"``. + interface_coordinates : sequence of (x, y), optional + Vertices of the curve or polygon. + boundary_coordinates : sequence of (x, y), optional + Extra vertices that close an open curve into the polygon that + defines the inside. + """ + eps = epsilon.array[:, 0, 0] if hasattr(epsilon, "array") else float(epsilon) + if signed_distance is not None: + psi.array[:, 0, 0] = _tanh_profile(signed_distance, eps) + return + if interface_geometry is None: + raise ValueError("Provide either signed_distance or interface_geometry.") + if interface_coordinates is None and interface_geometry != "shapely": + raise ValueError( + f"interface_coordinates is required for interface_geometry={interface_geometry!r}.") + distance = _signed_distance_from_geometry( + interface_geometry, interface, interface_coordinates, boundary_coordinates, + np.asarray(psi.coords)) + psi.array[:, 0, 0] = _tanh_profile(distance, eps) + + +def interface_thickness( + mesh: discretisation.Mesh, + phi: discretisation.MeshVariable, + *, + scale: float = 0.35, + use_min_edge_length: bool = False, +) -> discretisation.MeshVariable: + r"""Interface thickness :math:`\varepsilon` from the local cell size. + + :math:`\varepsilon = \mathrm{scale}\cdot V^{1/d}/\sqrt{d}` per cell + (or ``scale`` times the shortest edge), carried to ``phi``'s nodes from + the nearest cell centroid. Returned as a scalar MeshVariable of the same + degree as ``phi``. + + ``scale=0.35`` (the default, following the discontinuous-Galerkin + setting of g-adopt) gives :math:`\varepsilon \approx h/8` on triangles, + a band well under one cell. A continuous-Galerkin transport rings at + that: on a rotating circle at 32 cells across, the clip of the ringing + changed the volume by 0.84% per revolution at 0.35, 0.28% at 1.0 and + 0.18% at 2.0 (:math:`\varepsilon \approx 0.7h`, a band of two to three + cells, ringing gone); at 3.0 the reinitialisation's own curvature error + takes over (0.85%). For the SUPG transport, ``scale`` between 1.5 and 2 + is the sensible setting. + """ + from scipy.spatial import cKDTree + + dm = mesh.dm + dim = mesh.dim + c_start, c_end = dm.getHeightStratum(0) + n_cells = c_end - c_start + cell_epsilon = np.empty(n_cells) + cell_centroids = np.empty((n_cells, dim)) + + if use_min_edge_length: + if mesh.qdegree > 1: + raise ValueError("use_min_edge_length=True needs a straight-edged mesh (qdegree=1).") + v_start, v_end = dm.getDepthStratum(0) + coords = np.asarray(mesh.X.coords) + for i, cell in enumerate(range(c_start, c_end)): + closure, _ = dm.getTransitiveClosure(cell) + verts = [p - v_start for p in closure if v_start <= p < v_end] + v_coords = coords[verts] + edges = [np.linalg.norm(v_coords[a] - v_coords[b]) + for a in range(len(verts)) for b in range(a + 1, len(verts))] + cell_epsilon[i] = scale * min(edges) + cell_centroids[i, :] = v_coords.mean(axis=0) + else: + factor = scale / np.sqrt(dim) + for i, cell in enumerate(range(c_start, c_end)): + vol, centroid, _ = dm.computeCellGeometryFVM(cell) + cell_epsilon[i] = factor * float(np.asarray(vol).ravel()[0]) ** (1.0 / dim) + cell_centroids[i, :] = np.asarray(centroid).ravel()[:dim] + + epsilon = discretisation.MeshVariable( + r"\epsilon", mesh, 1, degree=phi.degree, continuous=phi.continuous) + _, nearest = cKDTree(cell_centroids).query(np.asarray(phi.coords)) + epsilon.array[:, 0, 0] = cell_epsilon[nearest] + return epsilon + + +def _signed_distance_closed(polygon, points): + """Positive inside a closed polygon, negative outside.""" + geometry, prepare = _shapely() + prepare(polygon) + boundary = polygon.boundary + inside = np.array([polygon.contains(geometry.Point(p)) for p in points]) + distance = np.array([boundary.distance(geometry.Point(p)) for p in points]) + return np.where(inside, distance, -distance) + + +def _signed_distance_open(curve, enclosed, points): + """Distance to an open curve, signed by the polygon that defines the inside.""" + geometry, prepare = _shapely() + prepare(enclosed) + inside = np.array([enclosed.intersects(geometry.Point(p)) for p in points]) + distance = np.array([curve.distance(geometry.Point(p)) for p in points]) + return np.where(inside, distance, -distance) + + +def _signed_distance_from_geometry(interface_geometry, interface, interface_coordinates, + boundary_coordinates, points): + geometry, _prepare = _shapely() + + def closed_from(coords): + if boundary_coordinates is not None: + raise ValueError("boundary_coordinates is only for an open interface.") + return _signed_distance_closed(geometry.Polygon(coords), points) + + def open_from(curve, coords): + if boundary_coordinates is None: + raise ValueError("boundary_coordinates must close an open interface.") + enclosed = geometry.Polygon(np.vstack((coords, boundary_coordinates))) + return _signed_distance_open(curve, enclosed, points) + + if interface_geometry == "curve": + curve = geometry.LineString(interface_coordinates) + return closed_from(interface_coordinates) if curve.is_closed else open_from(curve, interface_coordinates) + if interface_geometry == "polygon": + if boundary_coordinates is None: + return _signed_distance_closed(geometry.Polygon(interface_coordinates), points) + return open_from(geometry.LineString(interface_coordinates), interface_coordinates) + if interface_geometry == "shapely": + if interface is None: + raise ValueError("interface must be given for interface_geometry='shapely'.") + if isinstance(interface, geometry.Polygon): + return _signed_distance_closed(interface, points) + return open_from(interface, np.asarray(interface.coords)) + raise ValueError( + f"Unknown interface_geometry={interface_geometry!r}; choose 'curve', 'polygon' or 'shapely'.") + + +# --------------------------------------------------------------------------- +# The solver +# --------------------------------------------------------------------------- + +class LevelSetSolver: + r"""Conservative level-set transport on a scalar transport solver. + + Each call to :meth:`solve` advects the level set by one step, applies + the optional wall correction, reinitialises when due, and restores the + enclosed volume. + + **Reinitialisation** integrates, in pseudo-time :math:`\tau`, + + .. math:: + \frac{\partial\psi}{\partial\tau} + = -\psi(1-\psi)(1-2\psi) + \varepsilon(1-2\psi)|\nabla\psi| + + (Parameswaran and Mandal 2023, eq. 17) with SSP-RK3. Both terms carry the + factor :math:`(1-2\psi)`, so :math:`\psi = 0.5` is a fixed point: the + profile sharpens to width :math:`\varepsilon` without moving the + interface. :math:`|\nabla\psi|` is an :math:`L_2` projection onto the + mesh at each stage. + + **Mass correction** finds the uniform shift :math:`\delta` with + :math:`\int \mathrm{clip}(\psi + \delta, 0, 1)\,d\Omega` equal to the + initial enclosed volume and leaves the field in that clipped, shifted + state (Zhang, Zou and Greaves 2010). The map is monotone and its slope is + the area of the transition band, so a bracketed secant iteration started + from that slope reaches the target in a few integrals (five against about + thirty for the bisection it replaces, to the same 1e-10). + + Parameters + ---------- + level_set : MeshVariable + Continuous scalar field holding :math:`\psi`. + velocity : MeshVariable or sympy Matrix + Advecting velocity. + epsilon : MeshVariable + Interface thickness (:func:`interface_thickness`). + advection : {"supg", "slcn"}, default "supg" + The transport solver: the Eulerian SUPG solver or the semi-Lagrangian + one. Both are pure advection here. + order, theta : int, float + Time scheme of the transport solver (Crank-Nicolson by default, the + same meaning for both solvers). + reini_dt : float, optional + Pseudo-time step of the reinitialisation (default half the smallest + :math:`\varepsilon`). + reini_steps : int, default 5 + Pseudo-time steps per reinitialisation. + reini_frequency : int, optional + Advection steps between reinitialisations; by default from the + domain size and :math:`\varepsilon`. + far_field : float, optional + Value of :math:`\psi` imposed on every mesh boundary (0 outside the + interface, 1 inside). Set it whenever the flow crosses the domain + boundary: a continuous-Galerkin scheme with no value on an inflow + boundary lets mass in, measured as a 4% volume drift in twenty steps + of a rotating circle against 8e-5 with the value imposed. Leave it + unset only when the boundary is a streamline. + adv_solver_opts : dict, optional + PETSc options forwarded to the transport solver. + adv_solver_bc : sequence of str, optional + Box wall labels on which a zero normal gradient is imposed after + each step by copying the neighbouring interior nodes (a box-mesh + convenience). + conserve_mass : bool or "auto", default "auto" + Apply the global mass correction after each step. ``"auto"`` turns + it on for the semi-Lagrangian transport, which loses volume through + interpolation, and off for the Eulerian one, which conserves the + enclosed volume to solver tolerance on its own once ``far_field`` is + set where the flow crosses the boundary (measured: 8e-5 over twenty + steps; the reinitialisation changes it at second order only). What + does change it is the clip to [0, 1] of the transport's ringing at a + thin band: 0.84% per revolution of a circle at the default thickness + (``scale=0.35``), 0.18% at ``scale=2.0`` (see + :func:`interface_thickness`). Turn the correction on if that + matters; it costs about as much as the Eulerian advection step. + mass_correction_tol, mass_correction_max_iter + Bisection tolerance on the volume and iteration cap. + + Examples + -------- + >>> mesh = uw.meshing.UnstructuredSimplexBox(cellSize=1 / 32) + >>> psi = uw.discretisation.MeshVariable("psi", mesh, 1, degree=2) + >>> eps = uw.systems.level_set.interface_thickness(mesh, psi) + >>> uw.systems.level_set.initialise_psi(psi, eps, interface_geometry="polygon", + ... interface_coordinates=circle_points) + >>> ls = uw.systems.LevelSetSolver(psi, velocity=v.sym, epsilon=eps) + >>> for step in range(100): + ... ls.solve(dt) + """ + + def __init__( + self, + level_set: discretisation.MeshVariable, + *, + velocity, + epsilon: discretisation.MeshVariable, + advection: str = "supg", + order: int = 1, + theta: float = 0.5, + reini_dt: Optional[float] = None, + reini_steps: int = 5, + reini_frequency: Optional[int] = None, + far_field: Optional[float] = None, + adv_solver_opts: Optional[dict] = None, + adv_solver_bc=None, + conserve_mass: Union[bool, str] = "auto", + mass_correction_tol: float = 1.0e-10, + mass_correction_max_iter: int = 40, + ) -> None: + if level_set.num_components != 1: + raise ValueError("level_set must be a scalar MeshVariable.") + if not level_set.continuous: + raise ValueError("level_set must be a continuous MeshVariable.") + if advection not in ("supg", "slcn"): + raise ValueError(f"advection must be 'supg' or 'slcn', not {advection!r}.") + + self.phi = level_set + self.mesh = level_set.mesh + self.velocity = velocity.sym if isinstance(velocity, discretisation.MeshVariable) else velocity + self.epsilon = epsilon + self.advection = advection + self.reini_dt = float(reini_dt) if reini_dt is not None else 0.5 * self._global_min_epsilon() + self.reini_steps = int(reini_steps) + self.step = 0 + + if advection == "supg": + self._adv_solver = systems.AdvDiffusionSUPG( + self.mesh, self.phi, self.velocity, order=order, theta=theta) + else: + history = systems.ddt.SemiLagrangian( + self.mesh, self.phi.sym, self.velocity, + vtype=uw.VarType.SCALAR, degree=self.phi.degree, continuous=self.phi.continuous, + varsymbol="cphi", bcs=[], order=order, smoothing=0.0, + monotone_mode="clamp", theta=theta) + self._adv_solver = systems.AdvDiffusionSLCN( + self.mesh, u_Field=self.phi, V_fn=self.velocity, order=order, + DuDt=history, theta=theta) + self._adv_solver.constitutive_model = uw.constitutive_models.DiffusionModel + self._adv_solver.constitutive_model.Parameters.diffusivity = 0.0 + if far_field is not None: + for boundary in self.mesh.boundaries: + self._adv_solver.add_dirichlet_bc(float(far_field), boundary.name) + self._adv_solver_bc = adv_solver_bc + for key, value in (adv_solver_opts or {}).items(): + self._adv_solver.petsc_options[key] = value + + # |grad psi| for the reinitialisation, projected onto the mesh + self._grad_magnitude = sympy.sqrt(sum(g ** 2 for g in self.mesh.vector.gradient(self.phi.sym[0]))) + self.phi_grad = discretisation.MeshVariable( + r"|\nabla\psi|", self.mesh, 1, degree=self.phi.degree, continuous=self.phi.continuous) + self._grad_projector = systems.Projection(self.mesh, self.phi_grad, degree=self.phi.degree) + self._grad_projector.uw_function = self._grad_magnitude + + self._reini_frequency = int(reini_frequency) if reini_frequency is not None else self._default_frequency() + + if conserve_mass == "auto": + conserve_mass = advection == "slcn" + self.conserve_mass = bool(conserve_mass) + self._mass_correction_tol = float(mass_correction_tol) + self._mass_correction_max_iter = int(mass_correction_max_iter) + self._clip_volume_change = 0.0 + self._target_volume = self.interface_volume() + + # ------------------------------------------------------------------ + # Public interface + # ------------------------------------------------------------------ + + @property + def advection_solver(self): + """The transport solver carrying the level set.""" + return self._adv_solver + + @property + def reini_frequency(self) -> int: + """Advection steps between reinitialisations.""" + return self._reini_frequency + + def estimate_dt(self, **kwargs): + """The transport solver's timestep estimate (see its ``estimate_dt``).""" + return self._adv_solver.estimate_dt(**kwargs) + + def solve(self, dt: float, *, reinitialise: bool = True) -> None: + """Advance the level set by one step of size ``dt``.""" + self._adv_solver.solve(timestep=dt) + # The reinitialisation equation assumes psi in [0, 1]; the transport + # can overshoot at a band a cell wide, so clip before anything reads + # the field. Records what the clip removed, for the volume budget. + self._clip_volume_change += self._clip_in_place() + if self._adv_solver_bc: + self._apply_boundary_neumann(labels=self._adv_solver_bc) + self.step += 1 + + if reinitialise and self.step % self._reini_frequency == 0: + self.reinitialise() + self._clip_volume_change += self._clip_in_place() # RK stages can leave 1e-6 undershoots + if self._adv_solver_bc: + self._apply_boundary_neumann(labels=self._adv_solver_bc) + + if self.conserve_mass: + self._correct_mass(self._target_volume) + + def reinitialise(self) -> None: + """Run ``reini_steps`` SSP-RK3 pseudo-time steps of the reinitialisation equation.""" + for _ in range(self.reini_steps): + self._reini_ssprk3_step(self.reini_dt) + + @property + def volume_drift(self) -> float: + """Relative change of the enclosed volume since construction.""" + return (self.interface_volume() - self._target_volume) / self._target_volume + + def interface_volume(self) -> float: + r"""The enclosed volume :math:`\int\psi\,d\Omega`.""" + return uw.maths.Integral(self.mesh, self.phi.sym[0]).evaluate() + + def clamp(self, lo: float = 0.0, hi: float = 1.0) -> None: + """Clip the field to ``[lo, hi]`` in place. Not mass-conserving on its own.""" + self.phi.array[:, 0, 0] = np.clip(self.phi.array[:, 0, 0], lo, hi) + + def _clip_in_place(self) -> float: + """Clip to [0, 1]; return the volume the clip changed (a nodal estimate).""" + values = np.asarray(self.phi.array[:, 0, 0]) + clipped = np.clip(values, 0.0, 1.0) + if np.array_equal(values, clipped): + return 0.0 + before = self.interface_volume() + self.phi.array[:, 0, 0] = clipped + return self.interface_volume() - before + + # ------------------------------------------------------------------ + # Reinitialisation + # ------------------------------------------------------------------ + + def _rhs(self, values: np.ndarray) -> np.ndarray: + """The right-hand side of the reinitialisation equation at nodal values.""" + self.phi.array[:, 0, 0] = values + self._grad_projector.uw_function = self._grad_magnitude + self._grad_projector.solve() + grad = np.asarray(self.phi_grad.array[:, 0, 0]) + eps = np.asarray(self.epsilon.array[:, 0, 0]) + sharpening = -values * (1.0 - values) * (1.0 - 2.0 * values) + balance = eps * (1.0 - 2.0 * values) * grad + return sharpening + balance + + def _reini_ssprk3_step(self, dtau: float) -> None: + psi0 = np.array(self.phi.array[:, 0, 0]) + psi1 = psi0 + dtau * self._rhs(psi0) + psi2 = 0.75 * psi0 + 0.25 * psi1 + 0.25 * dtau * self._rhs(psi1) + self.phi.array[:, 0, 0] = psi0 / 3.0 + 2.0 * psi2 / 3.0 + 2.0 * dtau * self._rhs(psi2) / 3.0 + + def _global_min_epsilon(self) -> float: + from mpi4py import MPI + values = np.asarray(self.epsilon.array[:, 0, 0]) + local = float(values.min()) if values.size else np.inf + return uw.mpi.comm.allreduce(local, op=MPI.MIN) + + def _default_frequency(self) -> int: + """Reinitialise every step on a coarse mesh, less often as it refines.""" + from mpi4py import MPI + coords = np.asarray(self.mesh.X.coords) + dim = coords.shape[1] + hi = np.array([uw.mpi.comm.allreduce(float(coords[:, i].max()) if len(coords) else -np.inf, op=MPI.MAX) + for i in range(dim)]) + lo = np.array([uw.mpi.comm.allreduce(float(coords[:, i].min()) if len(coords) else np.inf, op=MPI.MIN) + for i in range(dim)]) + domain_size = float(np.sqrt(np.sum((hi - lo) ** 2))) + return max(1, round(4.9e-3 * domain_size / self._global_min_epsilon() - 0.25)) + + # ------------------------------------------------------------------ + # Wall correction (box meshes) + # ------------------------------------------------------------------ + + def _apply_boundary_neumann(self, labels=("Left", "Right", "Top", "Bottom")) -> None: + """Zero normal gradient on box walls: copy the neighbouring interior row or column.""" + from mpi4py import MPI + comm = uw.mpi.comm + + coords = np.asarray(self.phi.coords) + n_local = coords.shape[0] + axis_for_label = {"Left": 0, "Right": 0, "Top": 1, "Bottom": 1} + is_min_side = {"Left": True, "Right": False, "Top": False, "Bottom": True} + + for label in labels: + axis = axis_for_label[label] + tang = 1 - axis + op = MPI.MIN if is_min_side[label] else MPI.MAX + if n_local: + local_extreme = coords[:, axis].min() if is_min_side[label] else coords[:, axis].max() + else: + local_extreme = np.inf if is_min_side[label] else -np.inf + wall_val = comm.allreduce(float(local_extreme), op=op) + + local_axis_vals = np.unique(coords[:, axis]) if n_local else np.empty(0) + all_axis_vals = np.unique(np.concatenate(comm.allgather(local_axis_vals))) + if all_axis_vals.size < 2: + continue + inner_val = all_axis_vals[np.argsort(np.abs(all_axis_vals - wall_val))][1] + + inner_idx = np.where(np.isclose(coords[:, axis], inner_val, atol=1e-8))[0] if n_local else np.empty(0, dtype=int) + local_pairs = (np.column_stack((coords[inner_idx, tang], np.asarray(self.phi.array[inner_idx, 0, 0]))) + if len(inner_idx) else np.empty((0, 2))) + gathered = [p for p in comm.allgather(local_pairs) if p.shape[0] > 0] + if not gathered: + continue + table = np.vstack(gathered) + table = table[np.argsort(table[:, 0])] + table = table[np.concatenate(([True], np.diff(table[:, 0]) > 1e-10))] + + wall_idx = np.where(np.isclose(coords[:, axis], wall_val, atol=1e-8))[0] if n_local else np.empty(0, dtype=int) + if len(wall_idx) == 0: + continue + wall_tang = coords[wall_idx, tang] + pos = np.clip(np.searchsorted(table[:, 0], wall_tang), 1, len(table) - 1) + left_err = np.abs(wall_tang - table[pos - 1, 0]) + right_err = np.abs(table[pos, 0] - wall_tang) + nearest = np.where(right_err < left_err, pos, pos - 1) + good = np.minimum(left_err, right_err) <= 1e-6 + if not np.all(good): + warnings.warn( + f"Wall correction on {label!r}: {np.count_nonzero(~good)} wall node(s) " + "have no interior counterpart; left unchanged.", stacklevel=2) + self.phi.array[wall_idx[good], 0, 0] = table[nearest[good], 1] + + # ------------------------------------------------------------------ + # Mass correction + # ------------------------------------------------------------------ + + def _correct_mass(self, target: float, lo: float = 0.0, hi: float = 1.0) -> None: + r"""Uniform shift, clipped, restoring the enclosed volume to ``target``. + + Finds :math:`\delta` with :math:`\int\mathrm{clip}(\psi+\delta, lo, hi)\,d\Omega + = V_{\rm target}` and leaves the field in that state. The clip makes + the map :math:`\delta \mapsto V` nonlinear but monotone: it only moves + the transition band, so its slope is the band's area. A secant + iteration started from that slope converges in a few evaluations; + every evaluation is one integral over the mesh, which is what made + the bisection this replaces cost most of a level-set step. The + bracket is kept as a safeguard: an iterate that leaves it falls back + to its midpoint. + """ + data0 = np.array(self.phi.array[:, 0, 0]) + + def volume_for_shift(delta: float) -> float: + self.phi.array[:, 0, 0] = np.clip(data0 + delta, lo, hi) + return self.interface_volume() + + v0 = volume_for_shift(0.0) + residual = target - v0 + if abs(residual) < self._mass_correction_tol: + return + + # First guess: only the band moves, so dV/d(delta) is about its area. + # The nodal fraction of the field inside (lo, hi) times the domain area + # is a fair estimate of that on a reasonably uniform mesh. + span = hi - lo + band = float(np.mean((data0 > lo + 1e-6 * span) & (data0 < hi - 1e-6 * span))) + domain = self._domain_volume() + slope = max(band * domain, 1e-12) + d_prev, v_prev = 0.0, v0 + d_cur = residual / slope + lo_d, hi_d = (0.0, np.inf) if residual > 0 else (-np.inf, 0.0) + + for _ in range(self._mass_correction_max_iter): + v_cur = volume_for_shift(d_cur) + r_cur = v_cur - target + if abs(r_cur) < self._mass_correction_tol: + return + # keep the bracket [lo_d, hi_d] around the root (V is monotone) + if r_cur < 0: + lo_d = max(lo_d, d_cur) + else: + hi_d = min(hi_d, d_cur) + dv = v_cur - v_prev + if abs(dv) > 0: + d_next = d_cur - r_cur * (d_cur - d_prev) / dv + else: + d_next = d_cur + 2.0 * (d_cur - d_prev) + if not (lo_d < d_next < hi_d) or not np.isfinite(d_next): + if np.isfinite(lo_d) and np.isfinite(hi_d): + d_next = 0.5 * (lo_d + hi_d) + else: + d_next = 2.0 * d_cur + d_prev, v_prev, d_cur = d_cur, v_cur, d_next + + warnings.warn( + f"Mass correction did not reach the target volume {target:.6g} in " + f"{self._mass_correction_max_iter} iterations; leaving the field at the last shift.", + stacklevel=2) + + def _domain_volume(self) -> float: + if not hasattr(self, "_domain_volume_value"): + self._domain_volume_value = float(uw.maths.Integral(self.mesh, sympy.Integer(1)).evaluate()) + return self._domain_volume_value + + +# --------------------------------------------------------------------------- +# Material properties across the interface +# --------------------------------------------------------------------------- + +def material_property_field(level_set, field_values, interface: str): + r"""A material property blended across one or more level sets. + + Parameters + ---------- + level_set : sympy expression or list of them + The level-set field(s), e.g. ``psi.sym[0]``; with several, the last + is the innermost. + field_values : list of float + One value per material, innermost last. + interface : {"sharp", "sharp_adjoint", "arithmetic", "geometric", "harmonic"} + How the property crosses the interface. + """ + kinds = ("sharp", "sharp_adjoint", "arithmetic", "geometric", "harmonic") + if interface not in kinds: + raise ValueError(f"interface must be one of {kinds}, not {interface!r}.") + + level_sets = list(level_set) if isinstance(level_set, (list, tuple)) else [level_set] + values = list(field_values) + + result = None + while level_sets: + ls = sympy.Max(sympy.Min(level_sets.pop(), 1), 0) + value = values.pop() + other = values.pop() if not level_sets else result + if interface == "sharp": + result = sympy.Piecewise((value, ls > sympy.Rational(1, 2)), (other, True)) + elif interface == "sharp_adjoint": + shifted = ls - sympy.Rational(1, 2) + heaviside = (shifted + sympy.Abs(shifted)) / 2 / shifted + result = value * heaviside + other * (1 - heaviside) + elif interface == "arithmetic": + result = value * ls + other * (1 - ls) + elif interface == "geometric": + result = value ** ls * other ** (1 - ls) + else: + result = 1 / (ls / value + (1 - ls) / other) + return result diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index e8b82d845..2044c8571 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -321,6 +321,104 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True): return vel +def _advective_diffusive_dt(constitutive_K, V_fn, mesh, direction_aware=False, + percentile=0.0): + r"""Per-element resolution timestep, reduced to one global value. + + The minimum over cells of the advective crossing time :math:`h/|v|` and + the diffusive time :math:`h^2/\kappa`, nondimensional. Shared by the + semi-Lagrangian and the Eulerian advection-diffusion solvers: for both + it is a *resolution* estimate, not a stability limit. The semi-Lagrangian + scheme is unconditionally stable and the implicit Eulerian scheme is + stable at any cell Courant number; what bounds either one is accuracy + on the feature being transported, which the mesh cannot know. + + Parameters + ---------- + constitutive_K : sympy expression or number + Diffusivity (the constitutive model's unified ``K``). + V_fn : sympy Matrix + Advecting velocity, evaluated at cell centroids. + mesh : Mesh + direction_aware : bool, default False + Use the per-cell extent along the local velocity instead of the + isotropic radius (triangles only; falls back otherwise). + percentile : float, default 0.0 + ``0`` takes the strict global minimum; ``> 0`` takes that global + percentile of the per-element timesteps, so a few sliver cells + cannot collapse the estimate. + + Returns + ------- + (dt, dt_adv, dt_diff) : floats + The estimate and its two components; ``inf`` where a component does + not apply (zero velocity, zero diffusivity). + """ + from mpi4py import MPI + + comm = uw.mpi.comm + + diffusivity_glob = _global_max_diffusivity(constitutive_K, mesh) + vel = _centroid_velocities_nd(V_fn, mesh) + vel_magnitudes = np.linalg.norm(vel, axis=1) + element_radii = mesh._radii + + def _reduce_dt(per_elem): + fin = per_elem[np.isfinite(per_elem)] if len(per_elem) else per_elem + if percentile and percentile > 0: + gathered = comm.allgather(np.ascontiguousarray(fin, dtype=float)) + allv = (np.concatenate([a for a in gathered if a.size]) + if any(a.size for a in gathered) else np.empty(0)) + return float(np.percentile(allv, percentile)) if allv.size else np.inf + loc = float(np.min(fin)) if len(fin) else np.inf + return comm.allreduce(loc, op=MPI.MIN) + + if diffusivity_glob > 0: + dt_diff_per_element = (element_radii ** 2) / diffusivity_glob + else: + dt_diff_per_element = np.array([np.inf]) + + if direction_aware: + from underworld3.meshing.smoothing import _tri_cells + tris = _tri_cells(mesh.dm) + if tris is None: + h_per_element = element_radii + else: + coords = np.asarray(mesh.X.coords) + centroids = coords[tris].mean(axis=1) + vhat = np.where( + vel_magnitudes[:, None] > 0, + vel / np.maximum(vel_magnitudes[:, None], 1.0e-30), + 0.0) + D = coords[tris] - centroids[:, None, :] + # Signed projections of the cell vertices along v-hat: the + # extent material actually traverses through the cell. + s = np.einsum('cvd,cd->cv', D, vhat) + h_per_element = np.maximum(s.max(axis=1) - s.min(axis=1), 0.0) + else: + h_per_element = element_radii + + with np.errstate(divide='ignore', invalid='ignore'): + dt_adv_per_element = np.where( + vel_magnitudes > 0, h_per_element / vel_magnitudes, np.inf) + + dt_diff = _reduce_dt(dt_diff_per_element) + dt_adv = _reduce_dt(dt_adv_per_element) + return min(dt_diff, dt_adv), dt_adv, dt_diff + + +def _dimensionalise_dt(dt_estimate): + """Return a timestep estimate with physical time units when a model with + reference scales is active, otherwise as a plain nondimensional scalar.""" + try: + return uw.dimensionalise(np.squeeze(dt_estimate), {'[time]': 1}) + except Exception: + # Sanctioned fallback: no active scaling model. _as_scalar because + # np.squeeze promotes a Python float to a 0-d array, which is not a + # number any caller expects (see _apply_unit_aware_scaling). + return _as_scalar(np.squeeze(dt_estimate)) + + def _invalidate_solution_cache(u): """Drop the cached data view of a just-solved variable. @@ -4302,111 +4400,19 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): with reference scales is available, otherwise nondimensional. """ - ### required modules - from mpi4py import MPI - - comm = uw.mpi.comm - - ## global max diffusivity (unified .K property: diffusivity for - ## diffusion models) - diffusivity_glob = _global_max_diffusivity( - self.constitutive_model.K, self.mesh) - - ### velocity values at element centroids (nondimensional) - vel = _centroid_velocities_nd(self.V_fn, self.mesh) - - # Get per-element velocity magnitudes - vel_magnitudes = np.linalg.norm(vel, axis=1) - - # Get per-element radii (characteristic element size) - element_radii = self.mesh._radii - - ## estimate dt of adv and diff components using per-element approach - ## dt_adv_i = h_i / |v_i| for advection - ## dt_diff_i = h_i^2 / κ for diffusion (using global κ for now) - - # Reduce per-element dt to one global value. Default (percentile=0) = - # strict global MINIMUM — one cell sets the limit. percentile>0 takes the - # Nth global percentile (50 = median) of the per-element dt instead, so a - # few anisotropic SLIVER cells (velocity ACROSS a thin cell) don't collapse - # dt. SLCN is unconditionally stable, and ``direction_aware`` already - # credits cells stretched ALONG the flow — together they give the - # orientation-aware + sliver-robust timestep. - def _reduce_dt(per_elem): - fin = per_elem[np.isfinite(per_elem)] if len(per_elem) else per_elem - if percentile and percentile > 0: - gathered = comm.allgather(np.ascontiguousarray(fin, dtype=float)) - allv = (np.concatenate([a for a in gathered if a.size]) - if any(a.size for a in gathered) else np.empty(0)) - return float(np.percentile(allv, percentile)) if allv.size else np.inf - loc = float(np.min(fin)) if len(fin) else np.inf - return comm.allreduce(loc, op=MPI.MIN) - - # Per-element diffusive timestep (all elements use same diffusivity) - if diffusivity_glob > 0: - dt_diff_per_element = (element_radii ** 2) / diffusivity_glob - else: - dt_diff_per_element = np.array([np.inf]) - - # Per-element advective timestep — either isotropic - # (mesh._radii / |v|) or direction-aware (v-aligned cell - # extent / |v|). - if direction_aware: - # Per-cell vertex indices (triangle / tet). - from underworld3.meshing.smoothing import _tri_cells - tris = _tri_cells(self.mesh.dm) - if tris is None: - # Fall back to isotropic for non-triangle meshes. - h_per_element = element_radii - else: - coords = np.asarray(self.mesh.X.coords) - centroids = coords[tris].mean(axis=1) - # v-hat per cell (use centroid v we already have) - vhat = np.where( - vel_magnitudes[:, None] > 0, - vel / np.maximum(vel_magnitudes[:, None], - 1.0e-30), - 0.0) - D = coords[tris] - centroids[:, None, :] - # Signed projections along v̂ per cell vertex - s = np.einsum('cvd,cd->cv', D, vhat) - h_per_element = s.max(axis=1) - s.min(axis=1) - # Sanity-floor — for zero-velocity cells s=0 - # ⇒ h_eff=0 ⇒ dt_adv=inf via the where below - h_per_element = np.maximum( - h_per_element, 0.0) - else: - h_per_element = element_radii - - with np.errstate(divide='ignore', invalid='ignore'): - dt_adv_per_element = np.where( - vel_magnitudes > 0, - h_per_element / vel_magnitudes, - np.inf - ) - # Global reduction — strict min (percentile=0) or Nth percentile (median). - min_dt_diff_glob = _reduce_dt(dt_diff_per_element) - min_dt_adv_glob = _reduce_dt(dt_adv_per_element) + dt_estimate, dt_adv, dt_diff = _advective_diffusive_dt( + self.constitutive_model.K, self.V_fn, self.mesh, + direction_aware=direction_aware, percentile=percentile) # Store for user inspection - self.dt_adv = min_dt_adv_glob if not np.isinf(min_dt_adv_glob) else 0.0 - self.dt_diff = min_dt_diff_glob if not np.isinf(min_dt_diff_glob) else 0.0 + self.dt_adv = dt_adv if not np.isinf(dt_adv) else 0.0 + self.dt_diff = dt_diff if not np.isinf(dt_diff) else 0.0 - # Take overall minimum (respecting infinity for zero velocity/diffusivity cases) - dt_estimate = min(min_dt_diff_glob, min_dt_adv_glob) - - # If both are infinite (no velocity and no diffusivity), return infinity + # Both infinite (no velocity and no diffusivity): nothing to bound if np.isinf(dt_estimate): return np.inf - # Dimensionalise the result to physical time - try: - return uw.dimensionalise(np.squeeze(dt_estimate), {'[time]': 1}) - except Exception: - # Fallback: return plain nondimensional number. _as_scalar because - # np.squeeze promotes a Python float to a 0-d array, which is not - # a number any caller expects (see _apply_unit_aware_scaling). - return _as_scalar(np.squeeze(dt_estimate)) + return _dimensionalise_dt(dt_estimate) @timing.routine_timer_decorator def solve( diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 01e1a0582..21d190348 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -773,6 +773,24 @@ def getext( @timing.routine_timer_decorator +def _aux_component_offsets(mesh): + """Component offset of every field of the mesh DM, keyed by field id. + + Read from the DM itself, not from ``mesh.vars``: a MeshVariable that + was dropped and collected leaves its PETSc field in the DM (a DMPlex + cannot shed a field), and PETSc lays the auxiliary arrays out over + ALL fields in field order. The offsets therefore have to count the + orphaned fields too. + """ + offsets = {} + total = 0 + for field_id in range(mesh.dm.getNumFields()): + fe, _label = mesh.dm.getField(field_id) + offsets[field_id] = total + total += fe.getNumComponents() + return offsets + + def generate_c_source( name, mesh: underworld3.discretisation.Mesh, @@ -822,7 +840,7 @@ def generate_c_source( count_bd_residual_sig, count_bd_jacobian_sig = callbacks.counts # `_ccode` patching - def ccode_patch_fns(varlist, prefix_str): + def ccode_patch_fns(varlist, prefix_str, component_offsets=None): """ This function patches uw functions with the necessary ccode routines for the code printing. @@ -848,11 +866,22 @@ def ccode_patch_fns(varlist, prefix_str): ordered according to their `field_id`. prefix_str: str The string prefix to write. + component_offsets: dict, optional + Component offset of every field in the DM, by ``field_id`` + (see ``_aux_component_offsets``). When given, each variable + is patched from ITS OWN field's offset instead of a running + count over ``varlist``: a field whose Python variable has + been dropped stays in the DM and still occupies its slots, + so a running count would shift every later variable onto + the wrong data. """ u_i = 0 # variable increment u_x_i = 0 # variable gradient increment lambdafunc = lambda self, printer: self._ccodestr for var in varlist: + if component_offsets is not None: + u_i = component_offsets[var.field_id] + u_x_i = u_i * mesh.cdim if var.vtype == VarType.SCALAR: # monkey patch this guy into the function type(var.fn)._ccodestr = f"{prefix_str}[{u_i}]" @@ -898,7 +927,8 @@ def ccode_patch_fns(varlist, prefix_str): # is important, as the secondary call will overwrite # those patched in the first call. - ccode_patch_fns(_stable_sorted(mesh.vars.values()), "petsc_a") + ccode_patch_fns(_stable_sorted(mesh.vars.values()), "petsc_a", + component_offsets=_aux_component_offsets(mesh)) ccode_patch_fns(primary_field_list, "petsc_u") # Also patch `BaseScalar` types. Nothing fancy - patch the overall type, diff --git a/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..ae9938fd9 --- /dev/null +++ b/tests/parallel/test_1077_advdiff_supg_parallel.py @@ -0,0 +1,49 @@ +"""The Eulerian SUPG solver gives the serial answer on any number of ranks. + +The scheme has no rank-local step: history is a mesh variable, the residual +is assembled by PETSc, the timestep is a runtime constant. So the integral +error against the rotating-Gaussian oracle after a few steps must match a +serial reference to solver tolerance, whatever the partition. + +Run: mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1077_advdiff_supg_parallel.py +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a, pytest.mark.mpi] + +# Serial reference, res 16, BDF2, dt 0.05, 8 steps (recorded with this file; +# np=2 reproduced it to 1.4e-12). +SERIAL_ERROR = 0.0301522514 + + +def _run(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, + qdegree=3, regular=False) + x, y = mesh.X + sol = uw.analytic.RotatingGaussian(mesh, sigma=0.12, centre_radius=0.5, omega=1.0) + T = uw.discretisation.MeshVariable("T1077", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), order=2) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + dt = 0.05 + adv.DuDt.set_initial_history( + [uw.function.evaluate(sol.at(-k * dt), T.coords).reshape(-1, 1, 1) for k in range(2)], + dt=dt) + for _ in range(8): + adv.solve(timestep=dt) + return sol.error(sol.at(8 * dt), T, norm="integral") + + +def test_error_is_partition_independent(): + err = _run() + assert np.isfinite(err) and err < 0.05, err + gathered = uw.mpi.comm.allgather(err) + assert max(gathered) - min(gathered) < 1e-12, gathered + if SERIAL_ERROR is not None: + assert abs(err - SERIAL_ERROR) < 1e-8, (err, SERIAL_ERROR) diff --git a/tests/test_1055_advdiff_supg_api.py b/tests/test_1055_advdiff_supg_api.py new file mode 100644 index 000000000..9b39dc983 --- /dev/null +++ b/tests/test_1055_advdiff_supg_api.py @@ -0,0 +1,265 @@ +"""API contract of the Eulerian SUPG advection-diffusion solver. + +Structural checks that run in seconds: the export, argument validation, the +scheme assembled from the history manager, and the rule that a change of +timestep is a change of a runtime constant, never a recompile. + +Run: pixi run python -m pytest tests/test_1055_advdiff_supg_api.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +@pytest.fixture(scope="module") +def mesh(): + return uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.25, qdegree=3) + + +def _solver(mesh, tag, **kwargs): + x, y = mesh.X + T = uw.discretisation.MeshVariable(f"T_{tag}", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate( + sympy.exp(-((x - 0.5) ** 2 + y ** 2) / 0.03), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSUPG(mesh, T, sympy.Matrix([[-y, x]]), **kwargs) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + return adv, T + + +def test_exported_and_constructs_with_the_slcn_defaults(mesh): + adv, _T = _solver(mesh, "a") + assert type(adv).__name__ == "SNES_AdvectionDiffusion_SUPG" + # order 1, theta 0.5: Crank-Nicolson, the semi-Lagrangian solver's default + assert adv.integrator == "am" and adv.order == 1 and adv.theta == 0.5 + assert isinstance(adv.DuDt, uw.systems.ddt.Eulerian) + assert adv.DuDt.V_fn is None, "advection is implicit, not a history correction" + + +def test_slcn_order_theta_pairs_select_the_documented_schemes(mesh): + assert _solver(mesh, "p1", order=1, theta=1.0)[0].integrator == "am" # backward Euler + assert _solver(mesh, "p2", order=2, theta=1.0)[0].integrator == "bdf" # SL-BDF2's counterpart + assert _solver(mesh, "p3", order=2)[0].integrator == "bdf" # theta 0.5 only bites at order 1 + with pytest.raises(ValueError, match="theta applies"): + _solver(mesh, "p4", order=2, theta=0.5) + + +def test_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._integrator = "am" + weights = adv._spatial_weights() + assert len(weights) == 3 + states = adv._states() + assert len(states) == 3 + # every history state appears (through its derivatives) in the advection operator + names = {str(atom.func) for atom in adv._advection().atoms(sympy.Function)} + for s in states[1:]: + assert any(str(s.func) in n for n in names), (s, names) + + +def test_timestep_change_is_a_constant_update_not_a_recompile(mesh): + adv, _T = _solver(mesh, "e", order=2) + adv.solve(timestep=0.01) + key = adv._current_jit_cache_key + names = [getattr(c, "name", str(c)) for c in adv.constants_manifest] + assert any(r"\Delta t" in n for n in names), names + assert any("BDF" in n for n in names), names + adv.solve(timestep=0.013) + assert adv._current_jit_cache_key == key + assert float(adv.delta_t.sym) == 0.013 + + +def test_timestep_change_reaches_the_kernels(mesh): + """A solver stepped 0.01 then 0.02 gives the same field as a fresh solver + stepped 0.02 from the same state: the constant is really updated.""" + adv1, T1 = _solver(mesh, "f1") + adv1.solve(timestep=0.01) + state = np.array(T1.array) + adv1.solve(timestep=0.02) + + adv2, T2 = _solver(mesh, "f2") + T2.array[...] = state + adv2.DuDt.initialise_history() + adv2.solve(timestep=0.02) + # to the linear-solver tolerance (measured 2e-11 against a 2e-2 control) + assert np.allclose(np.asarray(T1.array), np.asarray(T2.array), rtol=0, atol=1e-8) + + # negative control: a different timestep gives a visibly different field + adv3, T3 = _solver(mesh, "f3") + T3.array[...] = state + adv3.DuDt.initialise_history() + adv3.solve(timestep=0.01) + assert np.abs(np.asarray(T2.array) - np.asarray(T3.array)).max() > 1e-3 + + +def test_order_ramps_from_one_unless_history_is_planted(mesh): + adv, T = _solver(mesh, "g", order=2) + adv.solve(timestep=0.01) + assert adv.DuDt.effective_order == 1 + adv.solve(timestep=0.01) + assert adv.DuDt.effective_order == 2 + + adv2, T2 = _solver(mesh, "h", order=2) + adv2.DuDt.set_initial_history([np.array(T2.array), np.array(T2.array)], dt=0.01) + adv2.solve(timestep=0.01) + assert adv2.DuDt.effective_order == 2 + + +def test_galerkin_baseline_needs_no_rebuild(mesh): + adv, _T = _solver(mesh, "i") + adv.solve(timestep=0.01) + key = adv._current_jit_cache_key + adv.supg_weight = 0.0 + adv.solve(timestep=0.01) + assert adv._current_jit_cache_key == key + assert adv.supg_weight == 0.0 + + +def test_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.AdvDiffusionSUPG(child, T, sympy.Matrix([[-yc, xc]])) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + adv.solve(timestep=0.02) + assert adv.snes.getKSP().getPC().getType() == "asm" + assert adv._custom_mg is None + data = np.asarray(T.array[:, 0, 0]) + assert np.isfinite(data).all() and 0.9 < data.max() < 1.01 + + +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_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_advdiff_supg_rotating_gaussian.py b/tests/test_1100_advdiff_supg_rotating_gaussian.py new file mode 100644 index 000000000..4415bdc03 --- /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.AdvDiffusionSUPG(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_1100_levelset_rotation.py b/tests/test_1100_levelset_rotation.py new file mode 100644 index 000000000..c609499d1 --- /dev/null +++ b/tests/test_1100_levelset_rotation.py @@ -0,0 +1,92 @@ +"""Conservative level set under rigid rotation, on both transport solvers. + +A circle carried once round the box centre must come back: the enclosed +volume is held throughout (by the scheme itself for the Eulerian transport, +by the mass correction for the semi-Lagrangian one), the 0.5 contour returns +to its initial position, and the reinitialisation keeps the profile at its +thickness rather than letting it smear. + +Run: pixi run python -m pytest tests/test_1100_levelset_rotation.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.systems import level_set + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +RADIUS = 0.15 +CENTRE = (0.5, 0.75) + + +def _setup(tag, advection): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 32, qdegree=3) + x, y = mesh.X + psi = uw.discretisation.MeshVariable(f"psi_{tag}", mesh, 1, degree=2) + # a band of two to three cells; the g-adopt default 0.35 is under one + # cell and a continuous-Galerkin transport rings at it + eps = level_set.interface_thickness(mesh, psi, scale=2.0) + distance = RADIUS - np.sqrt((psi.coords[:, 0] - CENTRE[0]) ** 2 + (psi.coords[:, 1] - CENTRE[1]) ** 2) + level_set.initialise_psi(psi, eps, signed_distance=distance) + psi0 = np.array(psi.array) + # rigid rotation about the box centre, one revolution in 2 pi + velocity = sympy.Matrix([[-(y - 0.5), x - 0.5]]) + # the rotating flow crosses the walls: impose the far-field value there + solver = uw.systems.LevelSetSolver(psi, velocity=velocity, epsilon=eps, advection=advection, + far_field=0.0, reini_steps=1, reini_frequency=5) + return mesh, psi, psi0, solver + + +@pytest.mark.parametrize("advection", ["supg", "slcn"]) +def test_circle_returns_after_one_revolution(advection): + mesh, psi, psi0, solver = _setup(advection, advection) + area0 = solver.interface_volume() + # int psi exceeds pi r^2 by O(eps^2 / r^2) for a curved band: 12% here + assert area0 == pytest.approx(np.pi * RADIUS ** 2, rel=0.15) + + # the corrector is on by default for SLCN only; the Eulerian transport + # holds the volume by itself to solver tolerance, and what remains is + # the clip of residual ringing and the reinitialisation's curvature + # error, 0.18% per revolution at this band thickness + assert solver.conserve_mass == (advection == "slcn") + volume_tol = 1e-6 if advection == "slcn" else 5e-3 + n_steps = 200 + dt = 2.0 * np.pi / n_steps + for _ in range(n_steps): + solver.solve(dt) + assert abs(solver.volume_drift) < volume_tol, solver.volume_drift + + data = np.asarray(psi.array).reshape(-1) + assert data.min() >= 0.0 and data.max() <= 1.0 + # profile stays sharp: the transition band holds a small share of nodes + in_band = np.mean((data > 0.05) & (data < 0.95)) + assert in_band < 0.2, in_band + # the 0.5 contour is back: nodal mismatch of the indicator is small + mismatch = np.mean(np.abs((data > 0.5).astype(float) - (psi0.reshape(-1) > 0.5).astype(float))) + assert mismatch < 0.01, mismatch + + +def test_initialise_from_polygon_needs_shapely_or_a_distance(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, qdegree=2) + psi = uw.discretisation.MeshVariable("psi_poly", mesh, 1, degree=1) + eps = level_set.interface_thickness(mesh, psi) + pytest.importorskip("shapely") + angles = np.linspace(0.0, 2.0 * np.pi, 33) + circle = np.column_stack((CENTRE[0] + RADIUS * np.cos(angles), CENTRE[1] + RADIUS * np.sin(angles))) + level_set.initialise_psi(psi, eps, interface_geometry="polygon", interface_coordinates=circle) + data = np.asarray(psi.array).reshape(-1) + assert 0.0 <= data.min() and data.max() <= 1.0 and 0.0 < data.mean() < 0.2 + + +def test_material_property_field_blends(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, qdegree=2) + psi = uw.discretisation.MeshVariable("psi_mat", mesh, 1, degree=1) + field = level_set.material_property_field(psi.sym[0], [1.0, 100.0], "arithmetic") + assert field.subs(psi.sym[0], 1) == 100.0 and field.subs(psi.sym[0], 0) == 1.0 + with pytest.raises(ValueError, match="interface must be one of"): + level_set.material_property_field(psi.sym[0], [1.0, 100.0], "cubic")