diff --git a/docs/advanced/index.md b/docs/advanced/index.md index e30e123e6..1c58fce67 100644 --- a/docs/advanced/index.md +++ b/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 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/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/docs/examples/free_surface/advanced/Ex_Crameri2012Case1_uw3_FreeSurf.py b/docs/examples/free_surface/advanced/Ex_Crameri2012Case1_uw3_FreeSurf.py new file mode 100644 index 000000000..f62f6a2e5 --- /dev/null +++ b/docs/examples/free_surface/advanced/Ex_Crameri2012Case1_uw3_FreeSurf.py @@ -0,0 +1,335 @@ +#!/usr/bin/env python +# coding: utf-8 +""" +Crameri et al. (2012, GJI 189, 38-54) -- CASE 1: relaxation of a 2-D cosine +surface perturbation on a two-layer (lid + mantle) box. + + Ex_Crameri2012Case1_uw3_FreeSurf.py + +Setup, from Fig. 1(a) and Section 2.4 of the paper +-------------------------------------------------- + box 2800 km wide x 700 km high + lid 100 km thick, rho_L = 3300 kg/m3, eta_L = 1e23 Pa s + mantle rho_M = 3300 kg/m3, eta_M = 1e21 Pa s + g 10 m/s2 + surface cosine, initial maximum topography h_init = 7 km, + wavelength 2800 km (= the box width) + lid base FLAT at 600 km -- so the lid thickness varies laterally by + +/-7 per cent, as the paper states + sides symmetric / free slip + bottom no slip + +Analytical solution (paper eq. 23) +---------------------------------- + h(t) = h_init * exp(gamma t), gamma = -0.2139e-11 1/s, tau = 14.825 ka + +Two features of Case 1 make it the sharpest test available here: + +1. rho_L = rho_M, so the DRIVING buoyancy is identically zero: the only force + is the topographic self-load, the equilibrium surface is flat (h_inf = 0), + and the relaxation is a SINGLE mode. That is precisely the regime in which + FreeSurface's exponential update -- and the global least-squares gamma -- + are exact, so any error in the measured relaxation time is attributable to + the spatial discretisation and to the material representation, not to the + time integrator. + +2. Density is uniform, so the material model affects ONLY the viscosity. + Case 1 is therefore a clean, analytically-backed test of the viscosity + mixing rule across the lid/mantle interface -- the participating codes in + the paper do not agree on one (MILAMIN and SULEC use harmonic, FDCON and + I2VIS arithmetic, STAGYY geometric; Table 1). +""" + +import os +import sys +import numpy as np +import sympy +from mpi4py import MPI + +import underworld3 as uw +from underworld3.cython.petsc_discretisation import petsc_dm_find_labeled_points_local + + +# ===================================================================== +# switches +# ===================================================================== +MATERIAL = "swarm" # "swarm" (MaterialSwarm) or "cls" (MaterialCLS) +ETA_MIX = "arithmetic" # "arithmetic" | "harmonic" | "geometric" +res = 128 # cell ~21.9 km; 256 -> ~10.9 km, the paper's Case 1 resolution +h_init_km = 7.0 # paper's Case 1; 0.7 is the small-amplitude check it also reports + +if len(sys.argv) > 1: # allow: python Ex_Crameri2012Case1_uw3_FreeSurf.py cls + MATERIAL = sys.argv[1] +if len(sys.argv) > 2: + ETA_MIX = sys.argv[2] + + +# ===================================================================== +# dimensional setup and the non-dimensionalisation +# ===================================================================== +# Same scaling as the Case 2 scripts: +# length L* = 2800 km (box width -> 1.0) +# viscosity eta* = 1e21 Pa s (mantle -> 1.0) +# bodyforce (rho g)* = 33 Pa/m (rho_M g = 33000 Pa/m -> rho_g = 1000) +L_km, H_km, lid_km = 2800.0, 700.0, 100.0 + +L, H = 1.0, H_km / L_km # 1.0, 0.25 +H_mantle = (H_km - lid_km) / L_km # 0.21429 -- FLAT lid base +h_init = h_init_km / L_km # 0.0025 for 7 km +lam = 1.0 # wavelength = box width = 2800 km + +lid_eta, mantle_eta = 100.0, 1.0 # 1e23, 1e21 Pa s +lid_rho, mantle_rho = 3300.0, 3300.0 # equal: no driving buoyancy +rho_g = 1000.0 # rho_M g in ND units + +# time scale: u* = (rho g)* L*^2 / eta*, t* = L* / u* +_L = L_km * 1e3 +_rhog = mantle_rho * 10.0 / rho_g # 33 Pa/m per ND unit +_u = _rhog * _L**2 / 1.0e21 +t_star = _L / _u # 1.0823e13 s per ND time unit +_yr = 365.0 * 24 * 3600 +ka_per_nd = t_star / _yr / 1.0e3 # 343.18 ka per ND time unit + +# paper eq. 23 +gamma_dim = 0.2139e-11 # 1/s +tau_ka = 1.0 / gamma_dim / _yr / 1.0e3 # 14.825 ka +tau_nd = tau_ka / ka_per_nd # 0.043198 + +max_time = 100.0 / ka_per_nd # run to 100 ka, as in Fig. 2 +dt_set = tau_nd / 20.0 # 0.741 ka; the integrator is exact for + # one mode, this bounds mesh distortion +save_every = 5 + +outputPath = "op_Ex_Crameri2012Case1_uw3_FreeSurf_%s_%s_res%d/" % (MATERIAL, ETA_MIX, res) +if uw.mpi.rank == 0: + if os.path.exists(outputPath): + for i in os.listdir(outputPath): + os.remove(outputPath + i) + else: + os.makedirs(outputPath) +uw.mpi.barrier() + +if uw.mpi.rank == 0: + print("Crameri (2012) Case 1 -- %s materials, eta mixing = %s, res = %d" + % (MATERIAL, ETA_MIX, res)) + print(" 1 ND time unit = %.2f ka tau = %.4f ka = %.6f ND" + % (ka_per_nd, tau_ka, tau_nd)) + print(" h_init = %.1f km = %.5f ND run to 100 ka = %.4f ND (%.2f tau)" + % (h_init_km, h_init, max_time, max_time / tau_nd)) + print(" cell = %.2f km, lid = %.1f cells\n" % (L_km / res, lid_km / (L_km / res))) + + +# ===================================================================== +# mesh and unknowns +# ===================================================================== +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(L, H), cellSize=L / res, qdegree=3 +) + +vdegree, predegree, psidegree = 2, 1, 2 +v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=vdegree, continuous=True) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=predegree, continuous=True) +timeField = uw.discretisation.MeshVariable("time", mesh, 1, degree=1) + + +# ===================================================================== +# initial cosine topography +# ===================================================================== +# -cos so the maximum sits at the box CENTRE rather than on the side walls. +# With reflective/free-slip sides both signs are the same eigenmode at the same +# rate, but reading the extremum at an interior node avoids the corner where +# the rotated free-slip wall meets the stress-free top. +deform_fn = -h_init * sympy.cos(2.0 * sympy.pi * mesh.X[0] / lam) + +Dz = uw.discretisation.MeshVariable("Dz", mesh, 1, degree=1) +diffuser = uw.systems.Poisson(mesh, Dz) +diffuser.constitutive_model = uw.constitutive_models.DiffusionModel +diffuser.constitutive_model.Parameters.diffusivity = 1.0 +diffuser.add_essential_bc((deform_fn,), "Top") +diffuser.add_essential_bc((0.0,), "Bottom") +diffuser.solve() + +displacement = np.zeros((mesh.X.coords.shape[0], mesh.dim)) +displacement[:, -1] = uw.function.evaluate(Dz.sym[0], mesh.X.coords)[:, 0, 0] +mesh.deform(mesh.X.coords + displacement) + + +# ===================================================================== +# materials -- built AFTER the deform, and painted by PHYSICAL y +# ===================================================================== +# Order matters. The Laplacian carrier above moves every interior node, not +# just the surface, so a lid painted before the deform would be carried down +# with it and its base would follow the cosine. Painting afterwards, against +# the deformed coordinates, puts the lid base FLAT at 600 km and leaves the +# lid thickness varying by +/-7 per cent -- which is the configuration the +# paper specifies. +x, y = mesh.X +yhat = sympy.Matrix([[0, 1]]) + +if MATERIAL == "swarm": + params = uw.Params( + uw_proxy_location="integration_points", + uw_proxy_sampling="nearest", + uw_cell_size=L / res, # the ACTUAL cell size + uw_fill_param=3, + ) + materials = uw.swarm.MaterialSwarm( + mesh, + fill_param=params.uw_fill_param, + proxy_location=params.uw_proxy_location, + proxy_sampling=params.uw_proxy_sampling, + ) + rho_swarm = materials.add_variable("rho", size=1, dtype=float) + eta_swarm = materials.add_variable("eta", size=1, dtype=float) + + materials.add("mantle", rho=mantle_rho, eta=mantle_eta) + materials.add("lid", rho=lid_rho, eta=lid_eta) + # a swarm region is a BOOLEAN mask over particle coordinates + materials["lid"] = y > H_mantle + + v_mat = None + +elif MATERIAL == "cls": + # A dedicated velocity MeshVariable, NOT v.sym. The level-set solvers are + # compiled at construction and cannot be retargeted afterwards; given a + # genuine MeshVariable, FreeSurface keeps its DATA in sync with the + # consistent surface velocity (minus the mesh velocity) every step. Handed + # a bare sympy expression it cannot, and the materials would silently keep + # advecting with the free velocity and no ALE correction. + v_mat = uw.discretisation.MeshVariable( + "Vmat", mesh, mesh.dim, degree=vdegree, continuous=True + ) + materials = uw.level_set.MaterialCLS( + mesh, velocity=v_mat, degree=psidegree, epsilon_scale=1.75, + advection="supg", reini_steps=1, reini_frequency=5, + ) + materials.add("mantle", rho=mantle_rho, eta=mantle_eta) + materials.add("lid", rho=lid_rho, eta=lid_eta) + # a level-set region is a SIGNED DISTANCE; for a horizontal plane + # (y - H_mantle) is the exact one, which is what reinitialisation wants + materials["lid"] = y - H_mantle + +else: + raise ValueError("MATERIAL must be 'swarm' or 'cls', not %r" % MATERIAL) + +materials.mixing(eta=ETA_MIX, rho="arithmetic") +rho = materials.rho +eta = materials.eta + + +# ===================================================================== +# Stokes + free surface +# ===================================================================== +bodyforce = (-1.0 * rho / 3300.0 * rho_g) * yhat + +# rho_L == rho_M, so (-rho/3300*rho_g + rho_g) is identically zero: there is +# no driving buoyancy anywhere and the equilibrium surface is flat. Passing the +# zero vector explicitly says so, and spares the held lid a solve whose answer +# is known. background_buoyancy stays None because what we hand over here is +# already the driving-only part. +driving_buoyancy = sympy.Matrix([[0.0, 0.0]]) + +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta +stokes.bodyforce = bodyforce +stokes.materials = materials +stokes.add_essential_bc((0.0, 0.0), "Bottom") # no slip (paper, both cases) +stokes.add_rotated_freeslip_bc(0.0, "Left") # symmetric / free slip +stokes.add_rotated_freeslip_bc(0.0, "Right") +stokes.petsc_options["snes_type"] = "ksponly" +stokes.tolerance = 1.0e-5 + +fs = uw.systems.FreeSurface( + stokes, "Top", + buoyancy_scale=rho_g, + composition=materials, + driving_buoyancy=driving_buoyancy, +) + + +# ===================================================================== +# topography diagnostic +# ===================================================================== +topwall = petsc_dm_find_labeled_points_local(mesh.dm, "Top") + + +def max_topography(): + """Global maximum surface elevation above the reference level H.""" + local = -1.0e300 + if topwall is not None and np.size(topwall) > 0: + local = float((mesh.X.coords[topwall][:, -1] - H).max()) + return uw.mpi.comm.allreduce(local, op=MPI.MAX) + + +def surface_profile(): + """(x, h) along the top boundary on this rank, sorted by x.""" + if topwall is None or np.size(topwall) == 0: + return np.empty(0), np.empty(0) + c = mesh.X.coords[topwall] + o = np.argsort(c[:, 0]) + return c[o, 0], c[o, -1] - H + + +# ===================================================================== +# time loop +# ===================================================================== +step, time = 0, 0.0 +history = [(0.0, max_topography())] +if uw.mpi.rank == 0: + print(" step t [ka] h_max [km] analytic [km] rel.err") + print(" %4d %10.3f %11.4f %14.4f %9s" + % (0, 0.0, history[0][1] * L_km, h_init_km, "-")) + +while time < max_time: + fs.solve() + + dt = min(fs.estimate_dt(), dt_set, max_time - time) + fs.advance(dt) + time += dt + step += 1 + + h = max_topography() + history.append((time, h)) + + if step % save_every == 0 or time >= max_time: + t_ka = time * ka_per_nd + ana = h_init * np.exp(-t_ka / tau_ka) + if uw.mpi.rank == 0: + print(" %4d %10.3f %11.4f %14.4f %8.2f%%" + % (step, t_ka, h * L_km, ana * L_km, + 100.0 * (h - ana) / ana if ana != 0 else float("nan"))) + with mesh.access(timeField): + timeField.data[:, 0] = time * ka_per_nd + mesh.write_timestep( + "mesh", meshUpdates=True, meshVars=[v, p, timeField], + outputPath=outputPath, index=step, + ) + +# ===================================================================== +# output +# ===================================================================== +if uw.mpi.rank == 0: + hist = np.asarray(history) + t_ka = hist[:, 0] * ka_per_nd + h_km = hist[:, 1] * L_km + ana = h_init_km * np.exp(-t_ka / tau_ka) + np.savetxt( + outputPath + "topography.csv", + np.column_stack([t_ka, h_km, ana]), + delimiter=",", header="t_ka,h_max_km,h_analytic_km", comments="", + ) + + # relaxation time by log-linear fit over the first 2 tau, where the signal + # is well above the discretisation floor + m = (t_ka > 0) & (t_ka < 2 * tau_ka) & (h_km > 0) + slope = np.polyfit(t_ka[m], np.log(h_km[m]), 1)[0] + tau_num = -1.0 / slope + print("\n fitted relaxation time = %.3f ka (analytic %.3f ka, error %+.2f%%)" + % (tau_num, tau_ka, 100.0 * (tau_num - tau_ka) / tau_ka)) + print(" wrote " + outputPath + "topography.csv") + + xs, hs = surface_profile() + np.savetxt(outputPath + "surface_final.csv", + np.column_stack([xs * L_km, hs * L_km]), + delimiter=",", header="x_km,h_km", comments="") diff --git a/docs/examples/free_surface/advanced/Ex_Crameri2012Case2_uw3_FreeSurf.py b/docs/examples/free_surface/advanced/Ex_Crameri2012Case2_uw3_FreeSurf.py new file mode 100644 index 000000000..173da3a97 --- /dev/null +++ b/docs/examples/free_surface/advanced/Ex_Crameri2012Case2_uw3_FreeSurf.py @@ -0,0 +1,353 @@ +#!/usr/bin/env python +# coding: utf-8 +""" +Crameri et al. (2012, GJI 189, 38-54) -- CASE 2: a buoyant plume rises and +loads the base of the lithosphere, deflecting the free surface. + + Ex_Crameri2012Case2_uw3_FreeSurf.py + +Setup, from Section 2.3 and Fig. 1(b) of the paper +-------------------------------------------------- + box 2800 km wide x 700 km high + lithosphere 100 km thick, rho_L = 3300 kg/m3, eta_L = 1e23 Pa s + mantle 600 km thick, rho_M = 3300 kg/m3, eta_M = 1e21 Pa s + plume r_P = 50 km, rho_P = 3200 kg/m3, eta_P = 1e20 Pa s, + centred 1400 km from the sides and 300 km above the bottom + of the mantle layer + g 10 m/s2 + sides free slip; bottom no slip; top free surface + +Reference results (true free-surface codes, Figs 6a and 7; DIGITISED off the +figures, so treat them as +/- 20 m) +----------------------------------------------------------------------------- + ~0.2 Ma ~250 m end of the fast isostatic adjustment (Fig. 6b) + 4 Ma ~460 m + 8 Ma ~660 m + 16 Ma ~810 m maximum is reached around 16-17 Ma + 20 Ma ~830 m +The profile (Fig. 7) is a bump ~700 km wide at x = 1400 km with NEGATIVE side +lobes near -100 m -- the flexural moat -- which is a useful shape check. + +Timestep: the paper notes the isostatic relaxation time is ~15 ka and that +constant steps of 250-4000 yr were tested, with FDCON stable up to 16-32 ka. +Note that `fs.advance(1.0)` in ND units here would be 343 ka -- more than +twenty isostatic relaxation times in a single step. +""" + +import os +import sys +import numpy as np +import sympy +from mpi4py import MPI + +import underworld3 as uw +from underworld3.cython.petsc_discretisation import petsc_dm_find_labeled_points_local + + +# ===================================================================== +# switches +# ===================================================================== +MATERIAL = "swarm" # "swarm" (MaterialSwarm) or "cls" (MaterialCLS) +ETA_MIX = "arithmetic" # "arithmetic" | "harmonic" | "geometric" -- one switch, + # so both material models necessarily use the same rule +res = 128 # cell 21.9 km. The paper's UNDERWORLD Case 2 free-surface + # run used ~10 km cells; see the plume-resolution line + # printed at startup before trusting any amplitude. +MAX_TIME_MA = 20.0 # Fig. 6 runs to 20 Ma; use 1.0 for a smoke test +DT_KA = 20.0 # upper bound on the step; fs.estimate_dt() also applies +EPS_SCALE = 1.75 # MaterialCLS interface half-width, in cells + +# python Ex_Crameri2012Case2_uw3_FreeSurf.py [material] [eta_mix] [res] [max_time_Ma] +if len(sys.argv) > 1: + MATERIAL = sys.argv[1] +if len(sys.argv) > 2: + ETA_MIX = sys.argv[2] +if len(sys.argv) > 3: + res = int(sys.argv[3]) +if len(sys.argv) > 4: + MAX_TIME_MA = float(sys.argv[4]) + + +# ===================================================================== +# dimensional setup and non-dimensionalisation +# ===================================================================== +# length L* = 2800 km (box width -> 1.0) +# viscosity eta* = 1e21 Pa s (mantle -> 1.0) +# bodyforce (rho g)* = 33 Pa/m (rho_M g = 33000 Pa/m -> rho_g = 1000) +L_km, H_km, lid_km = 2800.0, 700.0, 100.0 + +L, H = 1.0, H_km / L_km # 1.0, 0.25 +H_mantle = (H_km - lid_km) / L_km # 0.21429 + +plume_r = 50.0 / L_km # 0.017857 +plume_x0 = 1400.0 / L_km # 0.5 -- 1400 km from the sides +plume_y0 = 300.0 / L_km # 0.10714 -- 300 km above the base + +lid_eta, mantle_eta, plume_eta = 100.0, 1.0, 0.1 # 1e23, 1e21, 1e20 Pa s +lid_rho, mantle_rho, plume_rho = 3300.0, 3300.0, 3200.0 +rho_g = 1000.0 # rho_M g in ND units + +_L = L_km * 1e3 +_rhog = mantle_rho * 10.0 / rho_g # 33 Pa/m per ND unit +_u = _rhog * _L**2 / 1.0e21 +t_star = _L / _u +_yr = 365.0 * 24 * 3600 +ka_per_nd = t_star / _yr / 1.0e3 # 343.18 ka per ND time unit + +max_time = MAX_TIME_MA * 1.0e3 / ka_per_nd +dt_set = DT_KA / ka_per_nd +SNAP_MA = [4.0, 8.0, 12.0, 16.0, 20.0] # Figs 5 and 7 snapshots + +# digitised from Figs 6(a)/7 -- see the docstring; +/- 20 m +REFERENCE_M = {4.0: 460.0, 8.0: 660.0, 12.0: 760.0, 16.0: 810.0, 20.0: 830.0} + +outputPath = "op_Ex_Crameri2012Case2_uw3_FreeSurf_%s_%s_res%d/" % (MATERIAL, ETA_MIX, res) +if uw.mpi.rank == 0: + if os.path.exists(outputPath): + for i in os.listdir(outputPath): + os.remove(outputPath + i) + else: + os.makedirs(outputPath) +uw.mpi.barrier() + + +# ===================================================================== +# mesh and unknowns +# ===================================================================== +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(L, H), cellSize=L / res, qdegree=3 +) + +vdegree, predegree, psidegree = 2, 1, 2 +v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=vdegree, continuous=True) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=predegree, continuous=True) +timeField = uw.discretisation.MeshVariable("time", mesh, 1, degree=1) + +x, y = mesh.X +yhat = sympy.Matrix([[0, 1]]) + +if uw.mpi.rank == 0: + cell_km = L_km / res + print("Crameri (2012) Case 2 -- %s materials, eta mixing = %s, res = %d" + % (MATERIAL, ETA_MIX, res)) + print(" 1 ND time unit = %.2f ka run to %.1f Ma = %.3f ND, dt <= %.1f ka" + % (ka_per_nd, MAX_TIME_MA, max_time, DT_KA)) + print(" cell = %.2f km : lid = %.1f cells, PLUME RADIUS = %.2f cells" + % (cell_km, lid_km / cell_km, 50.0 / cell_km)) + if 50.0 / cell_km < 4.0: + print(" ** the plume is under-resolved (the paper used ~10 km cells, ~5 cells") + print(" per plume radius). Amplitudes here are resolution-limited. **") + print("") + + +# ===================================================================== +# materials -- THE ONLY BLOCK THAT DIFFERS BETWEEN THE TWO TWINS +# ===================================================================== +if MATERIAL == "swarm": + params = uw.Params( + uw_proxy_location="integration_points", + uw_proxy_sampling="nearest", + uw_cell_size=L / res, # the ACTUAL cell size + uw_fill_param=3, + ) + materials = uw.swarm.MaterialSwarm( + mesh, + fill_param=params.uw_fill_param, + proxy_location=params.uw_proxy_location, + proxy_sampling=params.uw_proxy_sampling, + ) + rho_swarm = materials.add_variable("rho", size=1, dtype=float) + eta_swarm = materials.add_variable("eta", size=1, dtype=float) + + materials.add("mantle", rho=mantle_rho, eta=mantle_eta) + materials.add("lid", rho=lid_rho, eta=lid_eta) + materials.add("plume", rho=plume_rho, eta=plume_eta) + + # A swarm region is a BOOLEAN mask over particle coordinates. (A bare signed- + # distance expression is refused here -- there is no convention telling the + # swarm which sign is "inside".) + materials["lid"] = y > H_mantle + materials["plume"] = sympy.sqrt((x - plume_x0)**2 + (y - plume_y0)**2) < plume_r + + + def write_materials(index): + """Particle fields, for the viscosity snapshots of Fig. 5.""" + rho_swarm.data[:, 0] = uw.function.evaluate(materials.rho, materials.coords)[:, 0, 0] + eta_swarm.data[:, 0] = uw.function.evaluate(materials.eta, materials.coords)[:, 0, 0] + materials.write_timestep( + filename="swarm", swarmname="materials", outputPath=outputPath, index=index, + swarmVars=[materials.index, rho_swarm, eta_swarm], + ) + +elif MATERIAL == "cls": + # A dedicated velocity MeshVariable, NOT v.sym. Each material's LevelSetSolver + # compiles its velocity into its weak form at construction and cannot be + # retargeted afterwards; given a genuine MeshVariable, FreeSurface keeps that + # variable's DATA in sync with its consistent surface velocity (minus the mesh + # velocity) every step. Handed a bare sympy expression it cannot, and the + # materials would silently keep advecting with the free velocity and no ALE + # correction -- on a mesh that deforms every step, with SUPG transport. + v_mat = uw.discretisation.MeshVariable( + "Vmat", mesh, mesh.dim, degree=vdegree, continuous=True + ) + materials = uw.level_set.MaterialCLS( + mesh, velocity=v_mat, degree=psidegree, epsilon_scale=EPS_SCALE, + advection="supg", reini_steps=1, reini_frequency=5, + ) + eta_mesh = uw.discretisation.MeshVariable( + "etaM", mesh, vtype=uw.VarType.SCALAR, degree=psidegree, continuous=True) + rho_mesh = uw.discretisation.MeshVariable( + "rhoM", mesh, vtype=uw.VarType.SCALAR, degree=psidegree, continuous=True) + + materials.add("mantle", rho=mantle_rho, eta=mantle_eta) + materials.add("lid", rho=lid_rho, eta=lid_eta) + materials.add("plume", rho=plume_rho, eta=plume_eta) + + # A level-set region is a SIGNED DISTANCE, positive inside. (y - H_mantle) is + # the exact distance to a horizontal plane, and (r_P - |x - x0|) the exact + # distance to a circle -- which is what reinitialisation wants to see. + materials["lid"] = y - H_mantle + materials["plume"] = plume_r - sympy.sqrt((x - plume_x0)**2 + (y - plume_y0)**2) + + + def write_materials(index): + """Blended fields on the mesh, for the viscosity snapshots of Fig. 5.""" + eta_mesh.data[:, 0] = uw.function.evaluate(materials.eta, eta_mesh.coords)[:, 0, 0] + rho_mesh.data[:, 0] = uw.function.evaluate(materials.rho, rho_mesh.coords)[:, 0, 0] + mesh.write_timestep( + "materials", meshUpdates=True, meshVars=[eta_mesh, rho_mesh], + outputPath=outputPath, index=index, + ) + +else: + raise ValueError("MATERIAL must be 'swarm' or 'cls', not %r" % MATERIAL) + +materials.mixing(eta=ETA_MIX, rho="arithmetic") +rho = materials.rho +eta = materials.eta + + +# ===================================================================== +# Stokes + free surface +# ===================================================================== +# rho_P != rho_M, so the driving buoyancy is non-zero -- but ONLY inside the +# plume, since rho_L == rho_M makes the lithosphere neutrally buoyant. The +# held lid therefore recovers an h_inf driven purely by the plume load, and +# background_buoyancy stays None because what is handed over below is already +# the driving-only part of the body force. +bodyforce = (-1.0 * rho / 3300.0 * rho_g) * yhat +driving_buoyancy = (-1.0 * rho / 3300.0 * rho_g + rho_g) * yhat + +# Integrated driving buoyancy: the decisive number for comparing the two +# material models. It is non-zero only in the plume, so it measures how much +# buoyant material each representation actually carries. If the twins disagree +# here, they will disagree on topography amplitude for that reason alone, and +# no amount of surface-side debugging will explain it. +if uw.mpi.rank == 0: + print(" integrated driving buoyancy = %.6e <-- compare between the twins" + % float(uw.maths.Integral(mesh, driving_buoyancy[1]).evaluate())) + print(" (analytic for a sharp disc = %.6e)" + % (np.pi * plume_r**2 * (1.0 - plume_rho / 3300.0) * rho_g)) + print("") + +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta +stokes.bodyforce = bodyforce +stokes.materials = materials +stokes.add_essential_bc((0.0, 0.0), "Bottom") # no slip +stokes.add_rotated_freeslip_bc(0.0, "Left") # free slip +stokes.add_rotated_freeslip_bc(0.0, "Right") +stokes.petsc_options["snes_type"] = "ksponly" +stokes.tolerance = 1.0e-5 + +fs = uw.systems.FreeSurface( + stokes, "Top", + buoyancy_scale=rho_g, + composition=materials, + driving_buoyancy=driving_buoyancy, +) + + +# ===================================================================== +# topography diagnostics +# ===================================================================== +topwall = petsc_dm_find_labeled_points_local(mesh.dm, "Top") + + +def max_topography(): + """Global maximum surface elevation above the reference level H, in ND.""" + local = -1.0e300 + if topwall is not None and np.size(topwall) > 0: + local = float((mesh.X.coords[topwall][:, -1] - H).max()) + return uw.mpi.comm.allreduce(local, op=MPI.MAX) + + +def save_profile(tag): + """Surface profile in km/m, for comparison with Fig. 7.""" + if topwall is None or np.size(topwall) == 0: + return + c = mesh.X.coords[topwall] + o = np.argsort(c[:, 0]) + np.savetxt( + outputPath + "surface_%s.csv" % tag, + np.column_stack([c[o, 0] * L_km, (c[o, -1] - H) * L_km * 1000.0]), + delimiter=",", header="x_km,h_m", comments="", + ) + + +# ===================================================================== +# time loop +# ===================================================================== +step, time = 0, 0.0 +history = [(0.0, 0.0)] +snaps = list(SNAP_MA) + +if uw.mpi.rank == 0: + print(" step t [Ma] h_max [m] paper [m]") + print(" %4d %9.4f %10.2f %11s" % (0, 0.0, 0.0, "0")) + +while time < max_time: + fs.solve() + + dt = min(fs.estimate_dt(), dt_set, max_time - time) + fs.advance(dt) + time += dt + step += 1 + + t_Ma = time * ka_per_nd / 1.0e3 + h_m = max_topography() * L_km * 1000.0 + history.append((t_Ma, h_m)) + + crossed = snaps and t_Ma >= snaps[0] + if crossed or step % 25 == 0: + ref = "" + if crossed: + key = snaps.pop(0) + ref = "%.0f" % REFERENCE_M.get(key, float("nan")) + save_profile("%02dMa" % int(round(key))) + with mesh.access(timeField): + timeField.data[:, 0] = t_Ma + mesh.write_timestep( + "mesh", meshUpdates=True, meshVars=[v, p, timeField], + outputPath=outputPath, index=step, + ) + write_materials(step) + if uw.mpi.rank == 0: + print(" %4d %9.4f %10.2f %11s" % (step, t_Ma, h_m, ref)) + +# ===================================================================== +# output +# ===================================================================== +save_profile("final") +if uw.mpi.rank == 0: + hist = np.asarray(history) + np.savetxt(outputPath + "topography.csv", hist, + delimiter=",", header="t_Ma,h_max_m", comments="") + print("\n %-10s %10s %10s %8s" % ("t [Ma]", "model [m]", "paper [m]", "diff")) + for t_ref, h_ref in sorted(REFERENCE_M.items()): + if t_ref <= hist[-1, 0]: + h = float(np.interp(t_ref, hist[:, 0], hist[:, 1])) + print(" %-10.1f %10.1f %10.1f %7.1f%%" + % (t_ref, h, h_ref, 100.0 * (h - h_ref) / h_ref)) + print("\n wrote " + outputPath + "topography.csv") diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index 3680040ad..9f6ed1f88 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -209,6 +209,7 @@ def view(): import underworld3.parameters import underworld3.materials import underworld3.checkpoint +import underworld3.level_set # After underworld3.function: the analytic suite still sources SolCx from it. import underworld3.analytic diff --git a/src/underworld3/level_set.py b/src/underworld3/level_set.py new file mode 100644 index 000000000..1f9cedcc8 --- /dev/null +++ b/src/underworld3/level_set.py @@ -0,0 +1,1308 @@ +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.AdvDiffusion`, composed with the default +:class:`~underworld3.systems.ddt.EulerianSUPG` transport manager) 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 +from underworld3.materials import MaterialDistribution, _property_name + + +# --------------------------------------------------------------------------- +# 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.AdvDiffusion( + 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. + # Named from self.phi.name (not a fixed literal) -- a hardcoded + # name here means every SECOND (and later) LevelSetSolver on the + # same mesh collides on this auxiliary variable: PETSc/UW3 then + # prints "Variable ... already exists - Skipping" and SILENTLY + # hands back the FIRST instance's phi_grad, so the second + # material's reinitialisation reads/writes the first material's + # gradient-magnitude field with no exception anywhere -- precisely + # the "silently shared level sets" hazard + # MaterialDistribution._check_level_sets_are_new's own docstring + # warns about, one level deeper (inside LevelSetSolver itself, + # not just at the distribution's own level-set naming). Every + # multi-material script with 2+ materials on one mesh (both + # MultiMaterialLevelSet and MaterialCLS) hits this. + self._grad_magnitude = sympy.sqrt(sum(g ** 2 for g in self.mesh.vector.gradient(self.phi.sym[0]))) + self.phi_grad = discretisation.MeshVariable( + rf"|\nabla {self.phi.name}|", 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 + + +# ============================================================================= +# MaterialCLS: materials on conservative level sets, in MaterialSwarm's own +# style +# ============================================================================= +# +# materials = uw.level_set.MaterialCLS(mesh, velocity=v.sym) +# materials.add("mantle", shear_viscosity_0=1.0, density=3300) +# materials.add("slab", shear_viscosity_0=1.0e3, density=3400) +# materials["slab"] = mesh.X[1] > 0.53 +# stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +# stokes.materials = materials +# stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1 +# +# reads exactly like the uw.swarm.MaterialSwarm example it is modelled on -- +# `.add()`, `materials["name"] = region`, `materials.`, +# `stokes.materials = materials` are all inherited, unmodified, from the +# SAME shared base class MaterialSwarm and MaterialRegions already use +# (underworld3.materials.MaterialDistribution) -- see that module's +# docstring: "Both present the same face to a solver". MaterialCLS is a +# third _distribution_ (where a material is), sharing the same _definition_ +# machinery (what a material is) as the other two. `stokes.materials = ...` +# therefore works identically for all three, with no special-casing. +# +# What's genuinely different from MaterialSwarm here: +# * No particles: each material (other than the implicit default -- see +# below) is its OWN conservative level set (LevelSetSolver), advected +# directly by `velocity` -- no swarm, no population control, no +# proxy-sampling choice. +# * `materials["name"] = region` needs a DISTANCE, not just a boolean +# mask: a level set's tanh profile is built from how FAR a point is +# from the interface, not merely which side it's on. A bare boolean +# half-space condition (`mesh.X[1] > 0.53`, exactly MaterialSwarm's own +# example) is auto-converted to its EXACT signed distance when it is +# linear/affine in the mesh coordinates -- the common case (a lid, a +# slab dipping at a fixed angle). For curved interfaces, pass a signed- +# distance expression directly (positive inside), e.g. +# `r - sympy.sqrt((x-x0)**2+(y-y0)**2)` for a circular inclusion -- see +# `_normalise_region`'s docstring for exactly what is and is not +# accepted, and why arbitrary nonlinear conditions are refused rather +# than silently guessed at. +# * Exactly ONE material may be left unpainted -- the FIRST one declared +# (materials.materials[0], index 0) -- and it becomes the IMPLICIT +# background, `max(1 - sum(others), 0)`, with no level set or solver of +# its own: this mirrors MaterialSwarm/MaterialRegions' own "unpainted +# particles/cells default to index 0" convention, translated to a +# continuous field. A SECOND unpainted material is refused (ambiguous +# for a continuous partition, unlike a discrete index). +# * `.solve(dt)` is the CLS-specific transport step (advects every +# explicit material's level set, with reinitialisation/mass +# correction) -- the direct analogue of MaterialSwarm's own +# `.advection(v, dt)`, called explicitly by the model script each step +# rather than implied by a solver. +# * `mixing()`/`blend()`/`materials.` add a third rule, +# "geometric" (`exp(sum(w_i log v_i))`), alongside the base class's +# "arithmetic" and "harmonic" -- the log-mean is the standard choice +# for a property spanning orders of magnitude (viscosity), where +# arithmetic under-weights a thin weak/strong inclusion and harmonic +# assumes an iso-stress (series) geometry that does not generally hold +# for an arbitrary, evolving multi-material arrangement. See +# `level_set.material_property_field`, which already offers +# "geometric" for the older two-material nested-list form this class +# supersedes for N independent, arbitrarily-adjacent materials. + + +class _DerivedLevelSet: + """The implicit default (index-0) material's "level set": never + advected or reinitialised on its own -- it simply IS the complement of + every explicitly-tracked material, by construction, recomputed fresh + every time it's read. Duck-types just enough of a MeshVariable + (`.sym[0]`) for MaterialDistribution.blend()'s `masks[i].sym[0]` to work + unmodified.""" + + def __init__(self, others): + self._others = others # explicit sibling level-set MeshVariables + + @property + def sym(self): + total = sum((sympy.Max(o.sym[0], 0) for o in self._others), sympy.Integer(0)) + return sympy.Matrix([[sympy.Max(1 - total, 0)]]) + + +class MaterialCLS(MaterialDistribution): + r"""Materials on conservative level sets, read the same way + :class:`~underworld3.swarm.MaterialSwarm` is. + + Parameters + ---------- + mesh : Mesh + The computational mesh. + velocity : MeshVariable or sympy expression + Shared advecting velocity for every explicit material's level set. + registry : MaterialRegistry, optional + Where the materials are defined. A fresh one is made if omitted; + pass one to share definitions with another distribution (e.g. a + :class:`~underworld3.swarm.MaterialSwarm` on a different mesh). + name : str, optional + Base name for the level-set variables. Defaults to ``material``, + then ``material_1`` and so on (shared counter with the other + distributions), so two on one mesh do not collide. + epsilon : MeshVariable or float, optional + Interface thickness, shared by every explicit material. Defaults + to :func:`interface_thickness` computed once the first explicit + material's level set exists, at ``scale=epsilon_scale``. Pass this + explicitly (and keep it consistent with whatever value governs + ongoing advection) rather than relying on two independent defaults + drifting apart. + epsilon_scale : float, default 0.35 + Scale passed to :func:`interface_thickness` when ``epsilon`` is + left as ``None`` -- since this class builds its own level-set + fields internally (lazily, at first build), there is no ``phi`` + of your own to call :func:`interface_thickness` on ahead of time + with a custom scale, hence this passthrough. + :func:`interface_thickness`'s own docstring recommends + ``1.5-2.0`` specifically for SUPG transport (the default here, + ``advection="supg"``) rather than its module default of ``0.35``, + which targets SLCN. + degree : int, default 2 + Polynomial degree of each level-set field. + advection, reini_steps, reini_frequency, theta, conserve_mass, + **levelset_kwargs + Forwarded to every explicit material's own + :class:`LevelSetSolver` (see its docstring — pure advection only; + no ``diffusivity`` parameter exists there). Pass a single value to + share it across all materials, or a dict keyed by material name to + vary it per material. + + Examples + -------- + >>> materials = uw.level_set.MaterialCLS(mesh, velocity=v.sym) + >>> materials.add("mantle", shear_viscosity_0=1.0, density=3300) + >>> materials.add("slab", shear_viscosity_0=1.0e3, density=3400) + >>> materials["slab"] = mesh.X[1] > 0.53 # affine -> exact distance + >>> stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + >>> stokes.materials = materials + >>> stokes.bodyforce = -materials.density * mesh.CoordinateSystem.unit_e_1 + + A curved interface needs an explicit signed distance (positive inside), + not a boolean condition -- there is no general closed-form distance to + an arbitrary nonlinear condition's boundary: + + >>> x, y = mesh.X + >>> materials.add("plume", shear_viscosity_0=1.0e-1, density=3200) + >>> materials["plume"] = 0.05 - sympy.sqrt((x - 0.5) ** 2 + (y - 0.1) ** 2) + + Advected once per step, alongside the Stokes/FreeSurface solve: + + >>> materials.solve(dt) + + See Also + -------- + underworld3.swarm.MaterialSwarm : the same interface, on particles. + underworld3.materials.MaterialRegions : the same interface, on mesh regions. + underworld3.materials.MaterialRegistry : where the definitions live. + """ + + def __init__( + self, + mesh, + *, + velocity, + registry=None, + name: Optional[str] = None, + epsilon=None, + epsilon_scale: float = 0.35, + degree: int = 2, + advection: str = "supg", + reini_steps: int = 1, + reini_frequency: Optional[int] = None, + theta: float = 0.5, + conserve_mass: Union[bool, str] = "auto", + **levelset_kwargs, + ): + # Bookkeeping BEFORE _init_distribution touches __getattr__-visible + # state, matching MaterialSwarm's own ordering rationale. + self._painted = set() + self._pending = {} + self._phi = {} + self._cls_solvers = {} # material name -> LevelSetSolver + # NB: named to not collide with MaterialDistribution's OWN + # `self._solvers`, which is an unrelated list of ATTACHED PDE + # solvers (stokes.materials = ... appends to that one). + self._built = False + self._default_name = None + self._degree = degree + self._epsilon = epsilon + self._epsilon_scale = epsilon_scale + self.mesh = mesh + # Keep the ORIGINAL MeshVariable too (not just its .sym), if one + # was given -- see sync_velocity_from()'s docstring for why: a + # MeshVariable's identity never changes, only its data, so writing + # into it updates every already-built LevelSetSolver's advection + # in place, with no rebuild -- unlike reassigning `self.velocity` + # itself, which only ever takes effect for solvers built AFTER the + # reassignment. + self._velocity_variable = velocity if isinstance(velocity, discretisation.MeshVariable) else None + self.velocity = velocity.sym if hasattr(velocity, "sym") else velocity + self._levelset_kwargs = dict( + advection=advection, reini_steps=reini_steps, + reini_frequency=reini_frequency, theta=theta, + conserve_mass=conserve_mass, + **levelset_kwargs, + ) + self._init_distribution(registry=registry, name=name) + + # -- level sets (MaterialDistribution's required hooks) --------------- + + @property + def built(self): + """Whether the level-set fields/solvers have been created yet (on + the first ``add()``-then-read, or an explicit ``build()``/ + ``populate()``). Public so callers taking ownership of a + not-yet-built ``MaterialCLS`` — e.g. ``FreeSurface``, which needs + to set its own velocity on it before that happens — can check + without reaching into a private attribute.""" + return self._built + + def _check_can_declare(self, name): + if self._built: + raise RuntimeError( + f"cannot add material {name!r}: already built. Declare every " + "material before the first region assignment or property read." + ) + + def _per_material(self, key, value, name): + if isinstance(value, dict): + if name not in value: + raise KeyError(f"{key!r} dict has no entry for material {name!r}.") + return value[name] + return value + + def _ensure_built(self): + if self._built: + return + if len(self._registry) == 0: + raise RuntimeError( + "no materials have been declared — call add() before reading " + "a material property" + ) + + materials_in_order = list(self._registry.materials) + default = materials_in_order[0] # index 0 -- same convention as + explicit = materials_in_order[1:] # MaterialSwarm/MaterialRegions + + unpainted = [m.name for m in explicit if m.name not in self._painted] + if unpainted: + raise RuntimeError( + f"material(s) {unpainted} were declared but never given a " + f"region (materials[name] = ...). Only the FIRST-declared " + f"material ({default.name!r}) may be left unpainted -- it " + "becomes the implicit background (whatever the others don't " + "cover), the same index-0-is-default convention " + "MaterialSwarm/MaterialRegions use." + ) + if default.name in self._painted: + # Explicitly painted too -- treat as a real, independent + # material rather than the derived background. + default = None + + self._phi = { + m.name: discretisation.MeshVariable( + rf"\phi_{{{m.name}}}", self.mesh, 1, degree=self._degree, continuous=True) + for m in materials_in_order + if default is None or m.name != default.name + } + if not self._phi: + raise RuntimeError( + "every material was left unpainted, including the implicit " + "default — at least one explicit region is required." + ) + + eps = self._epsilon + if eps is None: + eps = interface_thickness( + self.mesh, next(iter(self._phi.values())), scale=self._epsilon_scale) + + for name, region in self._pending.items(): + phi = self._phi[name] + sd = (region if isinstance(region, np.ndarray) + else uw.function.evaluate(region, phi.coords).flatten()) + initialise_psi(phi, eps, signed_distance=sd) + + for name, phi in self._phi.items(): + kwargs = { + k: self._per_material(k, v, name) + for k, v in self._levelset_kwargs.items() + } + kwargs.setdefault("epsilon", eps) + self._cls_solvers[name] = LevelSetSolver(phi, velocity=self.velocity, **kwargs) + + self._default_name = default.name if default is not None else None + self._pending = {} + self._built = True + + def _level_sets(self): + self._ensure_built() + order = [m.name for m in self._registry.materials] + others = [self._phi[n] for n in order if n != self._default_name] + return [ + _DerivedLevelSet(others) if n == self._default_name else self._phi[n] + for n in order + ] + + # -- painting ----------------------------------------------------------- + + def __setitem__(self, name, region): + """Give ``name`` a region -- see the class docstring for what + ``region`` may be (an affine half-space condition, a signed- + distance expression, or a precomputed array).""" + if name not in {m.name for m in self._registry.materials}: + raise KeyError(f"material {name!r} was not declared — call add({name!r}, ...) first") + if self._built: + raise RuntimeError( + "cannot paint a region: already built. Region painting must " + "happen before the first read (same rule as add())." + ) + self._painted.add(name) + self._pending[name] = self._normalise_region(region) + + def _normalise_region(self, region): + r"""Accept, in order of preference: + + - a precomputed numpy signed-distance array (used as-is; must match + the eventual ``phi.coords`` ordering for this material); + - a sympy ``Relational`` half-space condition + (``mesh.X[1] > 0.53``-style) that is AFFINE (linear) in the mesh + coordinates -- converted automatically to its EXACT signed + distance: :math:`(\ell - r)/\|\nabla(\ell-r)\|` for + :math:`\ell > r` (sign flipped for :math:`\ell < r`), which is + exact for any linear :math:`\ell - r` since a hyperplane's + distance field is itself linear. This is the direct level-set + analogue of ``materials["slab"] = mesh.X[1] > 0.53`` from + ``MaterialSwarm`` -- same syntax, exact result for this common + case; + - any other sympy expression, taken to ALREADY be a signed + distance (positive inside), e.g. ``r - sympy.sqrt((x-x0)**2 + + (y-y0)**2)`` for a circular inclusion. + + A NON-affine boolean condition (curved or compound, e.g. + ``(mesh.X[1] > 0.53) & (mesh.X[0] < 0.2)``) is refused rather than + silently approximated: there is no general closed-form distance to + an arbitrary condition's boundary, and a crude "inside=+1, + outside=-1" substitute would give the tanh profile no graded + transition zone at all (see the module docstring's warning about + exactly this trap). Pass a signed-distance expression directly + instead. + """ + if isinstance(region, np.ndarray): + return region + if isinstance(region, sympy.core.relational.Relational): + return self._affine_signed_distance(region) + if isinstance(region, sympy.Basic): + return region + raise TypeError( + "a region must be a signed-distance sympy expression (positive " + "inside), an affine half-space condition on the mesh coordinates " + "(e.g. mesh.X[1] > 0.53), or a precomputed numpy signed-distance " + f"array — not {type(region).__name__}" + ) + + def _affine_signed_distance(self, condition): + lhs, rhs = condition.lhs, condition.rhs + if isinstance(condition, (sympy.StrictGreaterThan, sympy.GreaterThan)): + signed = lhs - rhs + elif isinstance(condition, (sympy.StrictLessThan, sympy.LessThan)): + signed = rhs - lhs + else: + raise TypeError( + f"unsupported condition type {type(condition).__name__} for " + "automatic signed-distance conversion" + ) + coords = list(self.mesh.X) + expanded = sympy.expand(signed) + if not expanded.is_polynomial(*coords) or sympy.Poly(expanded, *coords).total_degree() > 1: + raise TypeError( + "automatic signed-distance conversion only supports a LINEAR " + "(half-space) condition in the mesh coordinates, e.g. " + "mesh.X[1] > 0.53 — for a curved or compound interface, pass " + "a signed-distance expression directly (positive inside), " + "e.g. `r - sympy.sqrt((x-x0)**2+(y-y0)**2)`." + ) + grad = [sympy.diff(signed, c) for c in coords] + norm = sympy.sqrt(sum(g ** 2 for g in grad)) + return signed / norm + + # -- transport (the CLS-specific step; MaterialSwarm's analogue is + # its own .advection()) ------------------------------------------------- + + def solve(self, dt: float, *, reinitialise: bool = True) -> None: + """Advance every explicit material's level set by one step of size + ``dt``. The implicit default material needs no step of its own — + it is always exactly the complement of the others.""" + self._ensure_built() + for solver in self._cls_solvers.values(): + solver.solve(dt, reinitialise=reinitialise) + + def estimate_dt(self, **kwargs) -> float: + """The most restrictive of the explicit materials' own transport + solvers' timestep estimates. + + ``basis=`` (only meaningful for a material using + ``advection="supg"`` — its ``LevelSetSolver`` forwards to + ``AdvDiffusion.estimate_dt``, which has a ``basis`` concept; + ``"slcn"``'s ``estimate_dt`` does not and would raise ``TypeError`` + if handed one) is applied ONLY to materials that are actually + SUPG-advected — this class allows mixed per-material advection + schemes (see ``__init__``'s per-material dict kwargs), so a + uniform forward would break the moment any material uses "slcn". + Every other kwarg is forwarded to all materials unchanged. + """ + self._ensure_built() + basis = kwargs.pop("basis", None) + estimates = [] + for solver in self._cls_solvers.values(): + call_kwargs = dict(kwargs) + if basis is not None and solver.advection == "supg": + call_kwargs["basis"] = basis + estimates.append(solver.estimate_dt(**call_kwargs)) + return min(estimates) + + def sync_velocity_from(self, other) -> None: + r"""Copy ``other``'s CURRENT DATA into this distribution's own + dedicated velocity ``MeshVariable`` (whatever was passed as + ``velocity=`` at construction) — for updating the ADVECTING + velocity's VALUES in place, without needing to rebuild any + already-built :class:`LevelSetSolver`. + + Why this works when reassigning :attr:`velocity` itself doesn't, + once :attr:`built` is True: each ``LevelSetSolver``'s compiled + weak form references the velocity variable's SYMBOL — a symbol + whose IDENTITY never changes here, only its underlying data does. + Reassigning ``self.velocity`` to a NEW symbolic expression only + takes effect for solvers built AFTER the reassignment (there is no + such thing as "already built" solvers retroactively noticing a + Python attribute changed); writing new VALUES into the SAME + variable those solvers already reference has no such problem — + PETSc reads a variable's current data at solve time regardless of + when it was last written. + + Requires ``velocity=`` to have been given as a genuine + ``MeshVariable`` at construction (not a bare sympy expression, + e.g. ``v.sym``) — otherwise there is no fixed-identity field to + write updated values into, and this raises. This is exactly the + fix for the case documented on ``FreeSurface``'s ``composition`` + parameter: construct ``MaterialCLS`` with a dedicated + ``MeshVariable`` (not ``v.sym`` directly) as ``velocity=`` so + ``FreeSurface`` can keep it in sync with its own consistent + surface velocity every step, even though the materials were + necessarily already built by the time ``FreeSurface`` sees them. + + Parameters + ---------- + other : MeshVariable + Must share the same discretisation (mesh, degree, number of + components) as this distribution's own velocity variable — + typically the case when syncing from another velocity field on + the SAME mesh (e.g. ``FreeSurface``'s + ``self.consistent.u``). A mismatch surfaces as numpy's own + broadcasting error on the array assignment below. + """ + if self._velocity_variable is None: + raise RuntimeError( + "sync_velocity_from() needs `velocity=` to have been a " + "genuine MeshVariable at construction (not a bare sympy " + "expression, e.g. v.sym) -- there is no dedicated field to " + "write updated values into. Construct this MaterialCLS with " + "a MeshVariable instead, e.g. `velocity=v_adv` (not " + "`velocity=v_adv.sym`)." + ) + # .array is the sanctioned write path and manages its own vector + # state, so no access context is needed (Charter S7). + self._velocity_variable.array[...] = other.array + + def level_set(self, name): + """The underlying level-set MeshVariable for an explicit material — + the machinery, exposed for cases this interface does not cover + (writing it to a checkpoint, plotting it directly), the same escape + hatch ``MaterialSwarm.index`` is for the particle case. The implicit + default material has none of its own (it's the complement of the + others, recomputed on read — see the class docstring) and raises.""" + self._ensure_built() + if name == self._default_name: + raise KeyError( + f"{name!r} is the implicit default material and has no " + "level-set field of its own — it's always the complement of " + "the explicit materials, not a stored field. See " + "MaterialCLS's class docstring." + ) + try: + return self._phi[name] + except KeyError: + raise KeyError(f"no explicit material named {name!r}") from None + + def interface_volumes(self): + """Per-explicit-material enclosed volume/area (see + ``LevelSetSolver.interface_volume``), as a ``{name: volume}`` dict. + Excludes the implicit default material.""" + self._ensure_built() + return {name: s.interface_volume() for name, s in self._cls_solvers.items()} + + def total_volume_error(self) -> float: + """``max_x |sum_i max(phi_i(x), 0) - 1|`` sampled at each explicit + material's own level-set nodes — a diagnostic for how far the + independently advected/reinitialised fields have drifted from + partition-of-unity before renormalisation (`blend()`/attribute + access always sums to exactly 1 by construction regardless).""" + self._ensure_built() + raw = [np.clip(np.asarray(phi.array[:, 0, 0]), 0, None) + for phi in self._phi.values()] + total = sum(raw) + if self._default_name is not None: + total = total + np.clip(1 - total, 0, None) + return float(np.max(np.abs(total - 1))) + + # -- mixing / blend: extend the base class's arithmetic/harmonic with + # "geometric" ------------------------------------------------------------ + + _MIXING_RULES_CLS = ("arithmetic", "harmonic", "geometric") + + def mixing(self, **rules): + """Same as ``MaterialDistribution.mixing()``, extended to also + accept ``"geometric"`` — see the module docstring for why that + matters for a viscosity-like property here.""" + declared = self._registry.declared_properties() + for name, rule in rules.items(): + name = _property_name(name) + if rule not in self._MIXING_RULES_CLS: + raise ValueError( + f"mixing rule for {name!r} must be one of " + f"{self._MIXING_RULES_CLS}, not {rule!r}" + ) + if name not in declared: + raise KeyError( + f"no material declares {name!r}, so a mixing rule for it " + f"would have no effect; declared: {sorted(declared)}" + ) + self._mixing_rules[name] = rule + self._push_all() + return self + + def _weights(self, clip: bool = True): + r"""Renormalised partition-of-unity weights for every material, in + registry order — :math:`w_i = \max(\varphi_i,0) / + (\sum_j \max(\varphi_j,0) + \delta)`. + + NEEDED, not just a nicety: ``MaterialDistribution.blend()`` (the + base class) does a raw weighted sum straight off ``masks[i].sym[0]`` + with no renormalisation, which is fine for ``MaterialSwarm``/ + ``MaterialRegions`` (their masks are already exactly one-hot or a + well-formed partition) but WRONG here — independently advected/ + reinitialised level sets genuinely drift apart from summing to 1 + (see :meth:`total_volume_error`), and an un-renormalised blend + would silently stop being a true weighted average as soon as they + do. :meth:`blend` below is therefore a full override, not a + geometric-only addition delegating elsewhere for the other two + rules — it must renormalise for "arithmetic"/"harmonic" too, not + just "geometric". + """ + self._ensure_built() + masks = self._level_sets() + raw = [sympy.Max(m.sym[0], 0) if clip else m.sym[0] for m in masks] + total = sum(raw) + 1.0e-12 + return [r / total for r in raw] + + def weight(self, name): + """The renormalised partition-of-unity weight for ONE material + (see :meth:`_weights`) — for when a formula needs one material's + own fraction directly (e.g. a body-force term keyed to a single + material), not a property blended across all of them.""" + self._ensure_built() + order = [m.name for m in self._registry.materials] + if name not in order: + raise KeyError(f"no material named {name!r}") + return self._weights()[order.index(name)] + + def blend(self, name, mixing=None): + r"""Same as ``MaterialDistribution.blend()``, but (a) always + renormalises the masks first (see :meth:`_weights` — REQUIRED + here, not optional, unlike the base class) and (b) additionally + supports ``mixing="geometric"``: + :math:`\exp(\sum_i w_i \log v_i)` — the log-mean, the usual choice + for a property spanning orders of magnitude. + """ + name = _property_name(name) + missing = [m.name for m in self._registry.materials + if name not in m.properties] + if missing: + raise KeyError( + f"property {name!r} is not declared by {missing}. Every " + "material must declare a property that is blended." + ) + self._ensure_built() + self._read_properties.add(name) + rule = mixing or self._mixing_rules.get(name, "arithmetic") + if rule not in self._MIXING_RULES_CLS: + raise ValueError( + f"mixing rule for {name!r} must be one of " + f"{self._MIXING_RULES_CLS}, not {rule!r}" + ) + + values = [m.resolved(name) for m in self._registry.materials] + w = self._weights() + + if rule == "arithmetic": + return sum(wi * vi for wi, vi in zip(w, values)) + if rule == "geometric": + return sympy.exp(sum(wi * sympy.log(vi) for wi, vi in zip(w, values))) + + # harmonic — same zero-value guard as the base class, for parity: + # a material whose value is zero poisons the symbol at declaration + # time (1/0 -> ComplexInfinity, a bare C-printer traceback later). + zeros = [m.name for m, value in zip(self._registry.materials, values) + if getattr(value, "is_zero", value == 0)] + if zeros: + raise ValueError( + f"harmonic mixing of {name!r} divides by its value in " + f"{zeros}, which is zero. Use arithmetic mixing, or give it " + "a small non-zero value." + ) + return 1 / sum(wi / vi for wi, vi in zip(w, values)) diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index 4d8522ef2..baef1b927 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -22,6 +22,8 @@ manager, EulerianSUPG, assembles implicit advection with SUPG). AdvDiffusionSLCN : class Advection-diffusion with semi-Lagrangian transport (flux history). +LevelSetSolver : class + Conservative level-set transport (advection, reinitialisation, mass correction). NavierStokes : class Navier-Stokes composed from a DDt transport manager (EulerianSUPG default). NavierStokesSLCN : class @@ -104,3 +106,10 @@ from .yield_continuation import (yield_continuation, YieldHomotopyControl, rate_strengthening_continuation, RateStrengtheningControl) from .solve_report import SolveReport + +# level_set lives at the package top level, not under systems: the Crameri +# examples reach uw.level_set.MaterialCLS and free_surface.py imports +# underworld3.level_set. It imports `systems` itself, so this import sits +# last, after the names it needs are bound. +from .. import level_set +from ..level_set import LevelSetSolver diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index ededaae64..df2a8b571 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -16,6 +16,7 @@ import numpy as np import sympy +import warnings from mpi4py import MPI import underworld3 as uw @@ -23,7 +24,6 @@ from underworld3.function._function import diff_wrt_field from underworld3.coordinates import CoordinateSystemType - class FreeSurface: r"""Evolve a free surface by exponential relaxation toward stress equilibrium. @@ -38,9 +38,17 @@ class FreeSurface: .. math:: - \gamma = \frac{\dot h}{h_\infty - h}, \qquad + \gamma = \frac{\oint w\,\dot h\,(h_\infty - h)\,\mathrm{d}s} + {\oint w\,(h_\infty - h)^2\,\mathrm{d}s}, \qquad h \leftarrow h_\infty + (h - h_\infty)\,e^{-\gamma\,\Delta t}. + :math:`\gamma` is ONE rate for the whole surface, recovered by least squares + against the FE trace mass. For a single relaxation mode — what this update is + derived for — that is exactly the pointwise ratio :math:`\dot h/(h_\infty-h)`, + which is then the same number at every node. Under a multi-wavelength load it + is not, and the pointwise form manufactures node-scale structure that is in + neither field; see :meth:`_relaxation_rate`. + Parameters ---------- stokes : uw.systems.Stokes @@ -58,16 +66,82 @@ class FreeSurface: Analytic outward surface normal for the rotated constraint and the normal-velocity read (e.g. ``X/|X|`` on a spherical cap). ``None`` uses the per-node geometric facet normal — correct as the surface deforms. - composition : uw.discretisation.MeshVariable or None, optional - A material field advected with the consistent surface velocity (e.g. a mesh - level set holding a signed distance, or temperature). ``None`` evolves the - surface alone (e.g. a topographic-relaxation benchmark). + composition : uw.discretisation.MeshVariable, underworld3.level_set.MaterialCLS, + underworld3.swarm_materials.MaterialSwarm, or None, optional + A material field advected with the consistent surface velocity. A + single ``MeshVariable`` (e.g. a mesh level set holding a signed + distance, or temperature) uses ``advection``/``conserve`` below, + transported by ONE solver FreeSurface builds and owns. A + :class:`~underworld3.level_set.MaterialCLS` instead already owns + its OWN per-material :class:`~underworld3.level_set.LevelSetSolver`\\ s + (one each, with their own advection scheme and mass conservation) — + FreeSurface takes it over by calling its own ``.solve(dt)`` each + step, rather than building a second, redundant transport layer on + top; ``advection``/``conserve`` below do not apply and must be left + at their defaults (see each). If the ``MaterialCLS`` is NOT yet + built (no property read, no ``.solve()``/``.build()`` call) when + passed in, FreeSurface also points its velocity at its own + consistent surface velocity before anything gets built — the ideal + case. If it's already built (the common, often unavoidable case + when the same materials also supply ``stokes.bodyforce``/ + ``shear_viscosity_0``, which must exist before ``stokes`` — and + hence ``FreeSurface`` — can be constructed), a warning says so: + those materials keep advecting with whatever velocity they were + originally built with, not the surface-consistent correction. See + :meth:`_build_composition_transport`. + + A :class:`~underworld3.swarm_materials.MaterialSwarm` is the third + supported form, and the simplest of the three to reason about: it + carries the materials on LAGRANGIAN PARTICLES holding absolute + physical coordinates, so FreeSurface builds no transport solver for + it at all and instead calls the swarm's own + ``advection(V_fn, dt, ...)`` each step, with its own consistent + surface velocity IN FULL — no ALE mesh-velocity correction, because + a particle has no moving frame to correct for (see + :meth:`_advect_swarm`). ``advection`` resolves to ``"swarm"`` on its + own here, and ``conserve`` is not valid (a particle population + conserves its material by construction). + + ``None`` evolves the surface alone (e.g. a topographic-relaxation + benchmark). + advection : {None, "slcn", "supg", "swarm"}, default None + Transport scheme. ``None`` (the default) selects the only scheme that + makes sense for the ``composition``'s own type: ``"swarm"`` for a + ``MaterialSwarm``, ``"slcn"`` otherwise — so the previous default + behaviour of a ``MeshVariable`` composition is unchanged. ``"swarm"`` + is resolved automatically and requires a ``MaterialSwarm``: passing it + with a mesh field raises, and passing a mesh scheme with a + ``MaterialSwarm`` warns and is ignored. + + For a single-``MeshVariable`` ``composition`` the choice is between the + two mesh schemes below. It is + ignored (with a warning if explicitly set) when ``composition`` is + a ``MaterialCLS`` — each material there already chose its own + scheme via its own ``LevelSetSolver``. ``"slcn"`` (default, + unchanged from before this option existed) is the semi-Lagrangian + solver (``uw.systems.AdvDiffusionSLCN``) with + ``old_frame_traceback=False`` — see + :meth:`_build_composition_transport` for why that specific setting + matters on a per-step-deforming mesh + (underworldcode/underworld3#423). ``"supg"`` uses the Eulerian + SUPG solver (``uw.systems.AdvDiffusion``) instead: no trace-back at + all, so the deforming-mesh hazard that motivated + ``old_frame_traceback=False`` doesn't apply to it in the first + place, at the cost of the usual SUPG behaviour on an + under-resolved front (see ``underworld3.systems.AdvDiffusion``'s + own docstring). Both call ``.solve(timestep=dt)`` identically, so + switching this does not change anything else about + ``FreeSurface``'s per-step behaviour. conserve : sympy expression or None, optional - An integrand whose domain integral is held fixed by a uniform shift of - ``composition`` after each advection — the enclosed-area conservation of a - mesh level set. Requires ``composition``; the shift derivative is taken - symbolically from this expression, so pass the same smoothed-Heaviside used - in the body force. + An integrand whose domain integral is held fixed by a uniform shift + of ``composition`` after each advection — the enclosed-area + conservation of a mesh level set. Requires a single-``MeshVariable`` + ``composition``; the shift derivative is taken symbolically from + this expression, so pass the same smoothed-Heaviside used in the + body force. Not valid with a ``MaterialCLS`` ``composition`` + (raises) — each material there already conserves its own volume via + its own ``LevelSetSolver``'s ``conserve_mass``, so a single blanket + shift would not even apply to a well-defined single field. driving_buoyancy : sympy expression or None, optional Body force for the held solve. ``None`` reuses ``stokes.bodyforce``; because :math:`\sigma_{nn}` is mean-removed and a constant body force gives a uniform @@ -107,6 +181,16 @@ class FreeSurface: (the background body force, e.g. ``-rho_g * rhat``) instead runs a twin held solve and subtracts its recovery — the exact reference mode. ``None`` (default) is the reduced/driving-only formulation, unchanged. + swarm_advection_order : int, optional + Runge-Kutta order handed to ``MaterialSwarm.advection``. Ignored unless + ``composition`` is a ``MaterialSwarm``, and dropped silently if that + build's ``advection`` does not take an ``order`` keyword. + swarm_restore_fn : callable or None, optional + Hook passed as ``restore_points_to_domain_func`` to + ``MaterialSwarm.advection``, for particles whose RK feet land outside the + domain. ``None`` (default) uses the mesh's own deform-aware + ``return_coords_to_bounds`` when it exists — the same primitive the + semi-Lagrangian trace-back relies on — and otherwise passes nothing. verbose : bool, optional Report per-step surface diagnostics through :func:`uw.mpi.rank`-safe output. """ @@ -118,6 +202,7 @@ def __init__( buoyancy_scale=1.0, normal=None, composition=None, + advection=None, conserve=None, driving_buoyancy=None, smooth_length=0.0, @@ -130,14 +215,99 @@ def __init__( consistent_constraint="strong", consistent_penalty=1.0e5, background_buoyancy=None, + swarm_advection_order=2, + swarm_restore_fn=None, verbose=False, ): + if advection not in (None, "slcn", "supg", "swarm"): + raise ValueError( + f"advection must be None, 'slcn', 'supg' or 'swarm', not {advection!r}." + ) + + # Lazy imports, and they MUST stay inside __init__: this module lives in + # underworld3.systems, while level_set.py and swarm_materials.py sit + # outside it and themselves import FROM underworld3.systems. Hoisting + # either to module level is a circular import. By the time user code + # constructs a FreeSurface, both modules are already loaded, so the cost + # here is a dict lookup. + from underworld3.level_set import MaterialCLS + try: + from underworld3.swarm_materials import MaterialSwarm + except ImportError: # build without swarm materials + MaterialSwarm = () + + # ONE discriminator, resolved once here; every per-step dispatch below + # reads _comp_kind rather than re-running isinstance checks: + # "none" no material field — the surface evolves alone + # "cls" MaterialCLS — owns per-material LevelSetSolvers; .solve(dt) + # "swarm" MaterialSwarm — Lagrangian particles; .advection(V, dt) + # "field" a single MeshVariable — ONE solver FreeSurface builds/owns + if composition is None: + self._comp_kind = "none" + elif isinstance(composition, MaterialCLS): + self._comp_kind = "cls" + elif MaterialSwarm and isinstance(composition, MaterialSwarm): + self._comp_kind = "swarm" + else: + self._comp_kind = "field" + # Back-compatible alias; existing internals and user code read this. + self._is_multi_material = self._comp_kind == "cls" + self._composition_velocity_sync = False # set for real in _build_composition_transport + + # Resolve the transport scheme against the composition's own type. None + # picks the only scheme that can apply; an explicit value is honoured + # where it is meaningful and refused where it is not -- a mesh field + # cannot be particle-advected, and particles are not moved by a PDE. + # Without this, `advection` can reach self._advection still set to None, + # where the "supg" test below silently falls through to the slcn branch. + if self._comp_kind == "swarm": + if advection not in (None, "swarm"): + warnings.warn( + f"FreeSurface: advection={advection!r} is ignored with a " + "MaterialSwarm `composition` -- particles are advected by " + "the swarm's own Lagrangian integrator, not by a mesh PDE " + "solver.", + stacklevel=2, + ) + advection = "swarm" + elif advection == "swarm": + raise ValueError( + "FreeSurface: advection='swarm' requires a MaterialSwarm " + f"`composition`; got {type(composition).__name__}. A mesh " + "field is transported with 'slcn' or 'supg'." + ) + elif advection is None: + advection = "slcn" + + if self._comp_kind == "swarm" and conserve is not None: + raise ValueError( + "FreeSurface: `conserve` is not valid with a MaterialSwarm " + "`composition` -- a Lagrangian particle population carries its " + "material with it and conserves it by construction; there is " + "no field for a blanket `conserve` shift to apply to." + ) + if self._comp_kind == "cls" and conserve is not None: + if conserve is not None: + raise ValueError( + "FreeSurface: `conserve` is not valid with a MaterialCLS " + "`composition` -- each material there already conserves " + "its own volume via its own LevelSetSolver's " + "conserve_mass; there is no single field for a blanket " + "`conserve` shift to apply to." + ) + self.free = stokes self.mesh = stokes.mesh self.surface = surface self.buoyancy_scale = buoyancy_scale self.normal = normal self.composition = composition + self._advection = advection + self._swarm_advection_order = int(swarm_advection_order) + self._swarm_restore_fn = ( + swarm_restore_fn if swarm_restore_fn is not None + else getattr(stokes.mesh, "return_coords_to_bounds", None)) + self._swarm_adv_kwargs = {} self._conserve_integrand = conserve self._smooth_length = smooth_length # sigma_nn de-smear: "lumped" is the monotone 2D default; on a 3D P2 @@ -242,7 +412,14 @@ def __init__( self._build_consistent() self._build_interior_diffuser() self._filter_iters = int(surface_filter) + # The ALE mesh-velocity correction exists for EULERIAN transport of a + # field whose nodes move with the mesh. A MaterialSwarm carries absolute + # particle coordinates and needs no such correction (see _advect_swarm), + # so its tracking variables and projection solver are never built. + self._track_mesh_velocity = self._comp_kind in ("cls", "field") if composition is not None: + if self._track_mesh_velocity: + self._build_mesh_velocity_tracking() self._build_composition_transport() self._h_inf = None # recovered in solve(), consumed in advance() @@ -789,35 +966,264 @@ def _opposite_boundary(self): others = [n for n in self._walls if n != self.surface] return others[0] if others else None + def _build_mesh_velocity_tracking(self): + r"""Mesh velocity, for the ALE advection correction that SUPG's + Eulerian transport needs but does not do on its own. + + SemiLagrangian (the ``"slcn"`` path) already handles a deforming + mesh internally — it temporarily deforms BACK to the old geometry + to evaluate history values there before restoring (confirmed from + ``underworld3.systems.ddt``'s own source: ``mesh._deform_mesh`` + inside ``SemiLagrangian``'s trace-back — see also + ``old_frame_traceback`` and underworldcode/underworld3#423, in + :meth:`_build_composition_transport`'s SLCN branch). The Eulerian + SUPG history manager (``EulerianSUPG``, what ``AdvDiffusion``/ + ``LevelSetSolver`` use for ``advection="supg"``) has NO such + handling anywhere in its implementation — confirmed by inspection, + not inferred. So on a mesh that moves every step (as FreeSurface's + does), an uncorrected SUPG advection term implicitly measures the + GRID time derivative :math:`D\psi/Dt|_{mesh}` (comparing values at + what are now different physical locations across a step), not the + physical :math:`\partial\psi/\partial t|_{fixed}` the transport + equation actually wants. The standard ALE fix (e.g. Donea et al. + 2004; see also ASPECT's own ALE documentation) is to advect with + the velocity RELATIVE to the mesh, :math:`\mathbf u - \mathbf + u_{mesh}`, not the raw material velocity. + + :math:`\mathbf u_{mesh}` is built at the mesh's own (P1) node + resolution — the same resolution :meth:`_carry_and_deform`'s + ``displacement``/``_normal_direction`` are already computed at — + then projected up to the velocity space (P2, matching + ``consistent.u``) with a ``Vector_Projection``, mirroring + :meth:`_build_surface_velocity_projection`'s existing P2→P1 + projection in the opposite direction. + + ORDERING CAVEAT: :meth:`_carry_and_deform` (which knows THIS + step's actual mesh motion) runs AFTER composition transport in + :meth:`advance`. So the correction applied to THIS step's + advection necessarily uses LAST step's mesh velocity — a + one-step lag, the same standard Picard-style trade-off already + used elsewhere in this project for an analogous ordering + conflict (the discontinuity-capturing coefficient lagged at + ``phi_old`` in ``solver_supg.py``, so its own nonlinearity + doesn't stall Newton). Zero on the very first step, since the + mesh hasn't moved yet — the correct initial value. + """ + self._mesh_disp_p1 = uw.discretisation.MeshVariable( + "mesh_disp_p1_fs", self.mesh, vtype=uw.VarType.VECTOR, degree=1, continuous=True) + self._mesh_velocity = uw.discretisation.MeshVariable( + "mesh_vel_fs", self.mesh, vtype=uw.VarType.VECTOR, degree=2, continuous=True) + self._mesh_vel_proj = uw.systems.Vector_Projection(self.mesh, self._mesh_velocity) + self._mesh_vel_proj.uw_function = self._mesh_disp_p1.sym + # Zero on the very first step, since the mesh has not moved yet. + # .data writes directly (Charter S7); no access context needed. + self._mesh_disp_p1.data[...] = 0.0 + self._mesh_velocity.data[...] = 0.0 + # Scratch space for consistent.u - mesh_velocity, at the velocity + # space's own degree. Only the MaterialCLS sync path needs actual + # DATA written here; the single-MeshVariable SUPG path instead + # uses self._mesh_velocity.sym directly in its V_fn expression + # (see _build_composition_transport), which stays current + # automatically — no explicit sync needed there. + self._adv_velocity_relative = uw.discretisation.MeshVariable( + "V_adv_rel_fs", self.mesh, vtype=uw.VarType.VECTOR, degree=2, continuous=True) + def _build_composition_transport(self): - """Semi-Lagrangian transport of the material field by the consistent surface - velocity. Diffusion is negligible; a level-set distance field carries no - monotone clamp (it is not a bounded [0,1] field).""" - self._comp_ddt = uw.systems.ddt.SemiLagrangian( - self.mesh, self.composition.sym, self._adv_velocity.sym, - vtype=uw.VarType.SCALAR, degree=self.composition.degree, continuous=True, - varsymbol="phi", bcs=[], order=1, smoothing=0.0, - # old_frame_traceback must be FALSE on a per-step-deforming mesh - # (underworldcode/underworld3#423). The old-frame reach-back — introduced as - # the fix for the earlier high-Ra blow-up — is itself an exponential - # amplifier once the surface deformation squeezes the near-boundary cells - # (onset ~5% of radius): the record→trace→solve loop then grows both T - # extremes ~10% per CYCLE (worse at smaller dt), mesh-locked, until T is - # unbounded. A minimal reproducer with no free surface at all (prescribed - # velocity + a ±0.1%/step mesh wobble) shows the same runaway with old-frame - # ON at any theta and is bounded with it OFF. The standard ALE path is safe - # here because the hazards that motivated old-frame are covered by fixes - # landed since: departure feet are restored by the deform-aware - # return_coords_to_bounds, and the monotone clamp bounds the sample. - # Measured on this problem (rho_g 2e5): old-frame ON is unusable beyond - # ~5.4% deformation; OFF holds T in [0,1] to 1e-3 through 17% deformation. - monotone_mode="clamp", theta=0.5, old_frame_traceback=False, - ) - self._comp_adv = uw.systems.AdvDiffusionSLCN( - self.mesh, u_Field=self.composition, V_fn=self._adv_velocity.sym, - order=1, DuDt=self._comp_ddt, - ) - self._comp_adv.constitutive_model = uw.constitutive_models.DiffusionModel + """Transport of the material field by the consistent surface velocity. + + MaterialCLS branch: no separate solver is built at all — the + MaterialCLS already owns one LevelSetSolver per explicit material. + If the materials are NOT yet built, FreeSurface points their + velocity at its own consistent surface velocity before anything + gets built (the ideal case — construct FreeSurface before reading + any material property). If they're ALREADY built — the common, + often unavoidable case when the same MaterialCLS also supplies + stokes.bodyforce/shear_viscosity_0, which must exist before `stokes` + (and hence FreeSurface, which reads stokes.bodyforce immediately) + can be constructed — reassigning `.velocity` now would be a no-op: + each material's own LevelSetSolver already baked its (then-current) + velocity into its compiled weak form at construction time. In that + case, IF the materials were given a genuine MeshVariable as + `velocity=` (not a bare sympy expression like `v.sym`), FreeSurface + instead keeps that variable's DATA in sync with its own consistent + surface velocity every step (:meth:`advance`, via + ``MaterialCLS.sync_velocity_from`` — see its docstring for why this + works even though reassigning `.velocity` itself does not). + Otherwise there is no fixed-identity field to update, and the + materials keep advecting with whatever velocity they were built + with — a warning says so every time it happens. + + Either way ``self._comp_adv`` is left at None, a marker + :meth:`advance` uses to call ``self.composition.solve(dt)`` instead + of a single solve. + + MaterialSwarm branch: no solver is built either, for a different + reason — particle advection is a METHOD CALL on the swarm, not a PDE + solve, so there is no weak form to compile, no history manager to + configure and no diffusivity or tolerance to set. All this branch + does is fix the velocity the particles will be integrated along and + work out, ONCE, which optional keywords this build's + ``advection`` actually accepts (checking the signature here rather + than catching TypeError at the call site, where a TypeError raised + *inside* advection would otherwise be mistaken for a bad keyword and + the particles advected a second time). See :meth:`_advect_swarm`. + + Single-MeshVariable branch (unchanged from before MaterialCLS + support existed): diffusion is negligible; a level-set distance + field carries no monotone clamp (it is not a bounded [0,1] field). + See ``advection`` in the class docstring for the "slcn" vs "supg" + choice. + """ + if self._comp_kind == "swarm": + # The velocity is the CONSISTENT surface velocity, in FULL -- not + # the ALE-relative (u - u_mesh) form the mesh-field branches use. + # Particles store absolute physical coordinates, so they must move + # at the material velocity; subtracting the mesh velocity would + # drag them backwards by exactly the mesh motion. Using + # `consistent` rather than `free` is what keeps a particle sitting + # ON the surface coincident with it -- the consistent solve is + # constrained so u.n = increment/dt there, which is precisely how + # far the mesh surface moves this step (see _solve_consistent). + self._swarm_V_fn = self._adv_velocity.sym + + advect = getattr(self.composition, "advection", None) + if not callable(advect): + raise AttributeError( + "FreeSurface: the MaterialSwarm `composition` has no " + "callable `advection` method -- this build predates " + "particle advection on material swarms, or the object " + f"passed ({type(self.composition).__name__}) is not a " + "MaterialSwarm." + ) + # Probe the signature ONCE so the per-step call is a single clean + # invocation with only the keywords this build understands. + import inspect + try: + params = inspect.signature(advect).parameters + except (TypeError, ValueError): # C-implemented / wrapped + params = {} + takes_kwargs = any( + p.kind is inspect.Parameter.VAR_KEYWORD for p in params.values() + ) + if "order" in params or takes_kwargs: + self._swarm_adv_kwargs["order"] = self._swarm_advection_order + if self._swarm_restore_fn is not None and ( + "restore_points_to_domain_func" in params or takes_kwargs + ): + self._swarm_adv_kwargs["restore_points_to_domain_func"] = ( + self._swarm_restore_fn + ) + self._comp_adv = None + return + + if self._is_multi_material: + self._composition_velocity_sync = False + if self.composition.built: + schemes = {s.advection for s in self.composition._cls_solvers.values()} + if "slcn" in schemes: + warnings.warn( + "FreeSurface: this MaterialCLS includes material(s) " + f"using advection='slcn' (schemes present: " + f"{sorted(schemes)}). The ALE mesh-velocity " + "correction below (needed for SUPG-advected " + "materials on this deforming mesh) is applied " + "through MaterialCLS's single SHARED velocity " + "field -- there is no way to correct only the " + "SUPG-advected materials and leave the SLCN-" + "advected ones alone. SLCN already handles " + "deformation internally (its own coordinate-remap " + "trace-back), so this OVER-corrects those specific " + "materials. Prefer a uniform advection scheme " + "across one MaterialCLS used under FreeSurface.", + stacklevel=2, + ) + if getattr(self.composition, "_velocity_variable", None) is not None: + # The fix: can't rebind the SYMBOL any already-built + # LevelSetSolver compiled against, but writing new DATA + # into that same variable (same identity) every step + # works regardless of when it was last written. The + # data written (advance()) is consistent.u MINUS the + # (one-step-lagged) mesh velocity -- see + # _build_mesh_velocity_tracking. + self._composition_velocity_sync = True + else: + warnings.warn( + "FreeSurface: the MaterialCLS `composition` was " + "already built before FreeSurface could point it at " + "its own consistent surface velocity, AND it was " + "constructed with a bare sympy expression for " + "`velocity=` (e.g. v.sym) rather than a genuine " + "MeshVariable -- there is no dedicated field " + "FreeSurface can update in place instead. Its " + "materials will keep advecting with whatever " + "velocity they were originally built with, NOT " + "FreeSurface's surface-consistent correction (and " + "not the ALE mesh-velocity correction either). Pass " + "a MeshVariable (not its .sym) as `velocity=` when " + "constructing MaterialCLS to enable the in-place " + "fix -- see MaterialCLS.sync_velocity_from.", + stacklevel=2, + ) + else: + # Not yet built: a clean symbolic assignment, no per-step + # sync needed -- see the single-MeshVariable SUPG branch + # below for why this stays current automatically. Same + # mixed-scheme caveat as above applies if materials are + # later declared with per-material advection="slcn"; it + # isn't checkable yet here (nothing is built). + self.composition.velocity = self._adv_velocity.sym - self._mesh_velocity.sym + self._comp_adv = None + return + + if self._advection == "supg": + # No trace-back here, so the deforming-mesh hazard documented below + # for old_frame_traceback on the SLCN path doesn't apply -- there is + # no "old frame" to reach back into. AdvDiffusion is a drop-in for + # AdvDiffusionSLCN (same .solve(timestep=dt), same + # .constitutive_model.Parameters.diffusivity API), so everything + # from self._comp_adv.tolerance down is unchanged for either branch. + # + # V_fn is the ALE-corrected RELATIVE velocity (consistent.u - + # mesh_velocity), not consistent.u alone -- see + # _build_mesh_velocity_tracking for why EulerianSUPG needs this + # and SemiLagrangian (the "slcn" branch, below) doesn't. Both + # are live symbols (self._adv_velocity.sym, + # self._mesh_velocity.sym), so this expression stays current + # automatically every solve -- no explicit per-step sync + # needed here (unlike the MaterialCLS branch above, whose + # solvers were already compiled before FreeSurface could get + # at them). + self._comp_adv = uw.systems.AdvDiffusion( + self.mesh, u_Field=self.composition, + V_fn=self._adv_velocity.sym - self._mesh_velocity.sym, + order=1, theta=0.5, + ) + else: + self._comp_ddt = uw.systems.ddt.SemiLagrangian( + self.mesh, self.composition.sym, self._adv_velocity.sym, + vtype=uw.VarType.SCALAR, degree=self.composition.degree, continuous=True, + varsymbol="phi", bcs=[], order=1, smoothing=0.0, + # old_frame_traceback must be FALSE on a per-step-deforming mesh + # (underworldcode/underworld3#423). The old-frame reach-back — introduced as + # the fix for the earlier high-Ra blow-up — is itself an exponential + # amplifier once the surface deformation squeezes the near-boundary cells + # (onset ~5% of radius): the record→trace→solve loop then grows both T + # extremes ~10% per CYCLE (worse at smaller dt), mesh-locked, until T is + # unbounded. A minimal reproducer with no free surface at all (prescribed + # velocity + a ±0.1%/step mesh wobble) shows the same runaway with old-frame + # ON at any theta and is bounded with it OFF. The standard ALE path is safe + # here because the hazards that motivated old-frame are covered by fixes + # landed since: departure feet are restored by the deform-aware + # return_coords_to_bounds, and the monotone clamp bounds the sample. + # Measured on this problem (rho_g 2e5): old-frame ON is unusable beyond + # ~5.4% deformation; OFF holds T in [0,1] to 1e-3 through 17% deformation. + monotone_mode="clamp", theta=0.5, old_frame_traceback=False, + ) + self._comp_adv = uw.systems.AdvDiffusionSLCN( + self.mesh, u_Field=self.composition, V_fn=self._adv_velocity.sym, + order=1, DuDt=self._comp_ddt, + ) + self._comp_adv.constitutive_model = uw.constitutive_models.DiffusionModel self._comp_adv.constitutive_model.Parameters.diffusivity = 1.0e-7 self._comp_adv.tolerance = 1.0e-4 if self._conserve_integrand is not None: @@ -999,9 +1405,25 @@ def estimate_dt(self, advect_scale=1.0): otherwise keep the surface near-flat and the solves ill-conditioned) — while the surface-motion cap still binds, so the mesh cannot be over-deformed and the surface-flow feedback cannot run away. Scale the advective step, never the - surface safety.""" + surface safety. + + With ``advection="supg"``, ``_comp_adv.estimate_dt()`` is called with + ``basis="resolution"`` explicitly: that class's default + (``basis="accuracy"``) returns a fundamentally different quantity (the + timestep for a target per-step CHANGE in the field, not a stability/ + cell-crossing limit), which would silently corrupt the + ``min(dt_advect, dt_surface)`` combination below if left at its default. + ``basis="resolution"`` is the one that means the same thing as the + semi-Lagrangian solver's ``estimate_dt()``.""" if self.composition is not None: - dt_advect = float(self._comp_adv.estimate_dt()) + if self._comp_kind == "cls": + dt_advect = float(self.composition.estimate_dt(basis="resolution")) + elif self._comp_kind == "swarm": + dt_advect = self._velocity_cfl() + elif self._advection == "supg": + dt_advect = float(self._comp_adv.estimate_dt(basis="resolution")) + else: + dt_advect = float(self._comp_adv.estimate_dt()) else: dt_advect = self._velocity_cfl() dt_advect *= float(advect_scale) @@ -1043,9 +1465,7 @@ def advance(self, dt): shape, h_inf, u_n = self._apply_tangential_transport(shape, h_inf, u_n, dt) displacement = h_inf - shape - with np.errstate(divide="ignore", invalid="ignore"): - gamma = np.where(np.abs(displacement) > 1.0e-9, u_n / displacement, 0.0) - gamma = np.abs(gamma) # relax toward equilibrium; L-stable for any rate + gamma = self._relaxation_rate(u_n, displacement) shape_new = h_inf + (shape - h_inf) * np.exp(-gamma * dt) # Total height change relative to the CURRENT mesh: the tangential transport # (shape - shape0) plus the normal relaxation (shape_new - shape). @@ -1053,10 +1473,52 @@ def advance(self, dt): if self.composition is not None: self._solve_consistent(increment, dt) - self._comp_adv.solve(timestep=dt, zero_init_guess=False) - self._conserve_composition() + if self._comp_kind == "cls": + # Each material already conserves its own volume via its + # own LevelSetSolver's conserve_mass -- no separate + # _conserve_composition() step (that's validated at + # construction: conserve=None is required here). + if self._composition_velocity_sync: + # Materials were already built when handed over (see + # _build_composition_transport): keep their DEDICATED + # velocity variable's data current every step, since + # their compiled solvers can't be retargeted any other + # way once built. ALE-corrected: consistent.u minus + # the mesh velocity from the LAST _carry_and_deform + # call (one-step lag -- see + # _build_mesh_velocity_tracking's ordering caveat; + # this step's own mesh motion isn't known until this + # same advance() call reaches _carry_and_deform, + # below). .array returns a TensorMeshArrayView, which + # supports indexing/assignment (the pattern used + # throughout this file, e.g. _solve_consistent's + # `.array[...] = other.array`) but NOT direct binary + # arithmetic between two whole views -- np.asarray() + # first to get plain ndarrays for the subtraction. + self._adv_velocity_relative.array[...] = ( + np.asarray(self._adv_velocity.array) + - np.asarray(self._mesh_velocity.array) + ) + self.composition.sync_velocity_from(self._adv_velocity_relative) + self.composition.solve(dt) + elif self._comp_kind == "swarm": + # Lagrangian: move the PARTICLES themselves, in absolute + # space, and do it BEFORE the mesh deforms -- the velocity + # they integrate is defined on the geometry the Stokes solves + # have just run on. Nothing is subtracted for mesh motion. + self._advect_swarm(dt) + else: + self._comp_adv.solve(timestep=dt, zero_init_guess=False) + self._conserve_composition() self._carry_and_deform(increment, dt) + if self._comp_kind == "swarm": + # The mesh has now moved underneath a particle population that did + # not move with it, so anything the swarm caches about where its + # particles sit is stale: cell ownership first, and then whatever + # it projects onto the mesh (the proxy fields the body force and + # the viscosity are built from). + self._sync_swarm_after_deform() if self.verbose: uw.pprint( f"FreeSurface: |increment|={np.abs(increment).max():.3e} " @@ -1065,6 +1527,63 @@ def advance(self, dt): # -- step internals ------------------------------------------------------- + def _relaxation_rate(self, u_n, displacement): + r"""ONE relaxation rate for the whole surface, by FE-trace-mass-weighted least + squares: :math:`\gamma = \oint w\,\dot h\,(h_\infty-h) / \oint w\,(h_\infty-h)^2`. + + The exponential update is derived for a surface relaxing as a SINGLE mode, and + with one mode the pointwise ratio :math:`\dot h/(h_\infty-h)` is the same number + at every node — so this least-squares rate recovers it EXACTLY and the update + stays exact at any :math:`\Delta t` (the cosine topographic-relaxation benchmark + is unchanged by this function). + + Under a multi-wavelength load it is not. A compact load beneath a stiff lid + excites a spectrum, each wavelength relaxing at its own rate, and + :math:`\dot h` then comes out spatially BROADER than :math:`h_\infty`, with + different zero crossings. Evaluating the ratio node by node there does not + approximate a rate; it manufactures structure that is in neither field: + + * where :math:`\dot h = 0` but :math:`h_\infty \neq h`, the pointwise + :math:`\gamma` is zero, so :math:`e^{-\gamma\Delta t}=1` and that node is + FROZEN at its current height while both its neighbours relax past it — a + one-node spike, sign-flipped relative to the local trend. This is the + damaging case, and it is invisible to a guard placed on the denominator. + * where :math:`h_\infty \to h` the ratio has a pole, but that one is benign: + the node is driven to :math:`h_\infty`, which is where it already is. + * where the two fields disagree in sign — normal near a load's flanks, since + neighbouring surface points are mechanically coupled and a point may + legitimately move away from its own pointwise equilibrium — an ``abs()`` + silently reverses the motion. + + A single global rate cannot do any of these: no node can be frozen, no pole + exists, and the sign is set by the surface as a whole. It is a genuine + restriction — a surface whose parts really do relax at different rates (a lid + that thickens across the domain, say) is served by one average rate — but that + is the approximation the exponential update already assumes, made honestly and + in one place, rather than a per-node ratio that is unbounded where it is least + meaningful. Measured on Crameri (2012) Case 2: node-to-node curvature drops by + a factor of ten, 2.98e-6 to 2.97e-7 against a 1.1e-4 signal. + + Reduction follows :meth:`_surface_mean`: owned-facet partial weights, counted + exactly once globally, so two scalar allreduces give the exact weighted rate + with no gather and no ordering, identical on every rank. + """ + w = self._surface_weights_local() + if w.size: + num = float(np.dot(w * np.asarray(u_n, dtype=float), displacement)) + den = float(np.dot(w * np.asarray(displacement, dtype=float), displacement)) + else: + num = den = 0.0 + comm = uw.mpi.comm + num = comm.allreduce(num, op=MPI.SUM) + den = comm.allreduce(den, op=MPI.SUM) + # den is a weighted sum of squares: zero only if the surface sits exactly at + # equilibrium everywhere, where a zero rate is the right answer. + gamma = abs(num / den) if den > 1.0e-300 else 0.0 + if self.verbose: + uw.pprint(f"FreeSurface: gamma={gamma:.4e} gamma*dt-safe (L-stable)") + return gamma + def _current_shape(self): """The surface height anomaly (topography direction), mean-removed — read from the deformed mesh geometry, which IS the surface state.""" @@ -1224,6 +1743,58 @@ def _conserve_composition(self): relax = 0.3 # full correction injects mesh-motion jitter (design note) self.composition.array[:, 0, 0] += relax * (current - self._conserve_target) / rate + def _advect_swarm(self, dt): + r"""Advect a ``MaterialSwarm``'s particles with the consistent surface velocity. + + Three things set this apart from the mesh-field branches, and all three + follow from a particle holding an ABSOLUTE coordinate rather than a value + at a node: + + *No ALE correction.* The velocity is the material velocity itself, not + :math:`\mathbf u - \mathbf u_{mesh}`. That relative form exists because an + Eulerian field lives at nodes that move underneath it (see + :meth:`_build_mesh_velocity_tracking`); a particle has no such frame, so + subtracting the mesh velocity would displace it backwards by exactly the + mesh motion every step. + + *Consistent, not free.* :attr:`_adv_velocity` is the consistent solve's + velocity, whose surface-normal component is constrained to the realised + relaxed rate :math:`\Delta h/\Delta t`. A particle on the surface therefore + travels exactly as far as the surface does, and the free boundary stays a + material boundary. Advecting with the FREE velocity instead moves + near-surface particles at the un-relaxed rate — opening a gap below the lid, + or pushing particles through it, by a margin that compounds every step. + + *Before the deform.* Called from :meth:`advance` while the mesh still holds + the geometry the Stokes solves ran on, so the velocity is sampled where it + is defined. :meth:`_carry_and_deform` moves the mesh immediately afterwards + and :meth:`_sync_swarm_after_deform` repairs the bookkeeping. + """ + self.composition.advection(self._swarm_V_fn, dt, **self._swarm_adv_kwargs) + + def _sync_swarm_after_deform(self): + """Re-establish particle bookkeeping after the mesh moved underneath it. + + Two things can be stale once :meth:`_carry_and_deform` has run: which cell + owns each particle, and whatever the swarm projects onto the mesh (the proxy + fields ``rho``/``eta`` are built from, and which therefore set the NEXT + step's body force and viscosity). Both are repaired through whichever hook + this build exposes; a build that refreshes them inside ``Mesh.deform`` — + which is handed ``dt`` precisely so it can do ALE bookkeeping — exposes none + of these names and this is then a no-op, which is the correct outcome. + """ + for hook in ("update_particle_owners", "_update_particle_owners", + "update_particle_cells", "migrate"): + fn = getattr(self.composition, hook, None) + if callable(fn): + fn() + break + for hook in ("update_proxies", "_update_proxies", "rebuild_proxies"): + fn = getattr(self.composition, hook, None) + if callable(fn): + fn() + break + def _carry_and_deform(self, increment, dt): """Carry the surface increment inward with the Laplacian diffuser and apply it as a mesh deformation along the topography direction (vertical on a box, radial @@ -1242,5 +1813,15 @@ def _carry_and_deform(self, increment, dt): displacement = np.asarray( function.evaluate(self._carry.sym[0], coords) ).flatten() - new_coords = coords + displacement[:, None] * self._normal_direction(coords) + disp_vec = displacement[:, None] * self._normal_direction(coords) + new_coords = coords + disp_vec + if self._track_mesh_velocity: + # THIS step's actual mesh motion, at the mesh's own (P1) node + # resolution, projected up to the velocity space -- consumed + # NEXT step (see _build_mesh_velocity_tracking's ordering + # caveat: this step's own composition transport, above, + # already ran before this method does). + target_shape = self._mesh_disp_p1.array.shape + self._mesh_disp_p1.array[...] = (disp_vec / dt).reshape(target_shape) + self._mesh_vel_proj.solve() self.mesh.deform(new_coords, dt=dt) diff --git a/tests/test_1105_levelset_rotation.py b/tests/test_1105_levelset_rotation.py new file mode 100644 index 000000000..dd8a66ac7 --- /dev/null +++ b/tests/test_1105_levelset_rotation.py @@ -0,0 +1,159 @@ +"""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_the_geometry_helpers_refuse_clearly_without_shapely(monkeypatch): + """shapely is deliberately NOT a dependency of underworld3, so for most + installs this refusal is what the polygon route actually does. The import + is forced to fail so the message is asserted on every machine, with or + without shapely present; the signed_distance= route it points at is the + one exercised by test_circle_returns_after_one_revolution above.""" + import builtins + + real_import = builtins.__import__ + + def no_shapely(name, *args, **kwargs): + if name == "shapely" or name.startswith("shapely."): + raise ImportError("no shapely") + return real_import(name, *args, **kwargs) + + monkeypatch.setattr(builtins, "__import__", no_shapely) + with pytest.raises(ImportError, match="signed_distance"): + level_set._shapely() + + +def test_initialise_from_polygon(): + 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") + + +def test_sync_velocity_from_reaches_the_solve_time_vector(): + """``sync_velocity_from`` must update what PETSc reads, not just numpy. + + The method's whole purpose is to change the advecting velocity's VALUES + in place so that already-built solvers -- which hold the variable's + symbol, not its data -- pick the new values up at the next solve. That + only works if the write lands in the variable's global vector, which is + what FreeSurface relies on every step (free_surface.py calls this with + its own consistent surface velocity). Comparing the global norms is the + assertion with teeth: a write that stopped at a detached numpy copy + would leave them apart while the array view still looked right. + """ + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, qdegree=2) + v_adv = uw.discretisation.MeshVariable("v_adv_sync", mesh, 2, degree=2) + v_src = uw.discretisation.MeshVariable("v_src_sync", mesh, 2, degree=2) + + x, y = mesh.X + materials = uw.level_set.MaterialCLS(mesh, velocity=v_adv) + materials.add("blob", 0.2 - sympy.sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)) + materials.build() + + coords = np.asarray(v_src.coords) + v_src.array[:, 0, 0] = 3.0 * coords[:, 0] + 1.0 + v_src.array[:, 0, 1] = -2.0 * coords[:, 1] - 0.5 + + before = np.array(v_adv.array[...]) + materials.sync_velocity_from(v_src) + + after = np.array(v_adv.array[...]) + assert not np.allclose(before, after) + assert np.allclose(after, np.array(v_src.array[...])) + assert v_adv._gvec.norm() == pytest.approx(v_src._gvec.norm(), rel=1e-12) + + +def test_sync_velocity_from_refuses_a_bare_expression(): + """No fixed-identity field to write into, so it says so rather than + silently doing nothing -- the failure mode documented on + ``FreeSurface``'s ``composition`` parameter.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=1.0 / 8, qdegree=2) + v = uw.discretisation.MeshVariable("v_sync_sym", mesh, 2, degree=2) + materials = uw.level_set.MaterialCLS(mesh, velocity=v.sym) + with pytest.raises(RuntimeError, match="genuine MeshVariable"): + materials.sync_velocity_from(v)