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

docs/advanced/index.md

@@ -150,6 +150,7 @@ particle-population-and-materials
 semi-lagrangian-time-integration
 eulerian-advection-diffusion
 eulerian-navier-stokes
+level-set-transport
 porous-flow
 snapshot-restore
 troubleshooting

docs/advanced/level-set-transport.md

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

docs/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"))