From 1ed2c1f7ee2564f501bd4ec857ff0acd65e24750 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 26 Sep 2026 13:24:20 -0700 Subject: [PATCH 1/5] Forward semi-Lagrangian history launched from the field's nodes and element interiors Forward from where the values are known best. A field is known exactly at its own nodes and, since its interpolant is a polynomial inside each element, at any interior point. ForwardNodesSemiLagrangian launches the values at the nodes and at the interior lattice of every element (the discontinuous basis one degree up), carries them one step forward on the shared characteristic trace, fits the arrivals in each cell at the field's own degree (the cell polynomial projector, with its thin-cell and conditioning fallbacks) and reads the fit back at the nodes. The fit never crosses an element boundary, where the interpolant has a kink; a neighbourhood fit across cells was tried first and smoothed (disc error 3.6e-2). As the DuDt of the SLCN advection-diffusion solver on the rotating diffusing Gaussian (P2 temperature, dt 0.02, kappa 0.01): half revolution on the disc 2.00e-2 (backward nodal history 2.52e-2, swarm 2.04e-2, SUPG 1.72e-2), peak 0.1873 against 0.1865 exact; the square box, where the flow crosses all four walls, 6.17e-2 (backward 6.43e-2); a full revolution 7.63e-2 (backward 8.94e-2, SUPG and swarm 8.2e-2). About twice the backward cost per step. First order, serial (the fit near a partition seam needs the other rank's arrivals). Test test_1102: disc and box, hard baselines. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- src/underworld3/systems/ddt.py | 186 ++++++++++++++++++ ...st_1102_forward_nodes_rotating_gaussian.py | 54 +++++ 2 files changed, 240 insertions(+) create mode 100644 tests/test_1102_forward_nodes_rotating_gaussian.py diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index b21f6a02..148da5f3 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -5766,3 +5766,189 @@ def commit_flux_to_history(self, flux, verbose=False): self._launch_values = self._evaluate_at_launch(flux) self._fit_arrivals(self._launch, self._launch_values, cell=self._launch_cell) self._history_committed = True + + +@dataclass +class DDtForwardNodesState(_DDtCoreState): + """Snapshot of a :class:`ForwardNodesSemiLagrangian` instance.""" + psi_star_var_names: list[str] = field(default_factory=list) + + +class ForwardNodesSemiLagrangian(_DDtBase): + r"""Forward semi-Lagrangian history launched from where the field is known. + + A field is known exactly at its own nodes (they are its unknowns) and, since + its interpolant is a polynomial inside each element, at any point of an + element's interior. Each step launches the values at the nodes and at a + lattice inside every element (the points of the discontinuous basis one + degree up), carries each one step forward along the characteristic, fits + the arrivals in each cell to a polynomial of the field's own degree, and + reads that fit back at the nodes. The fit never reaches across an element + boundary, where the interpolant has a kink. A cell with too few arrivals, or + arrivals on a line, falls back to a linear fit over the nearest arrivals; a + cell nothing reached keeps its previous fit (see + :class:`~underworld3.utilities.cell_polynomial_projection.CellPolynomialProjector`). + Arrivals that leave the domain are dropped; a node the flow reached from + outside takes :attr:`inflow_value` when one is set. + + The integration-point counterpart is :class:`ForwardSemiLagrangian`, which + launches from the quadrature points, where a flux (a stress) is formed. + First order; serial. + + Parameters + ---------- + mesh : Mesh + psi_fn : sympy expression or matrix + The carried quantity, for example ``T.sym``. + V_fn : sympy matrix + The velocity that carries it. + vtype : VarType + degree : int + Degree of the store and of the per-cell fit; use the field's own degree. + units : optional + Units of the carried quantity (see :func:`_history_units`). + """ + + applies_inflow_value = True + instances = 0 + + def __init__(self, mesh, psi_fn, V_fn, vtype=VarType.SCALAR, degree=1, varsymbol=None, + order=1, theta=0.5, units=None, **_unsupported): + super().__init__() + if order != 1: + raise NotImplementedError("ForwardNodesSemiLagrangian carries one level; order must be 1") + if mesh.cdim != mesh.dim: + raise NotImplementedError("ForwardNodesSemiLagrangian fits in the embedding coordinates; no manifolds") + if uw.mpi.size > 1: + raise NotImplementedError("ForwardNodesSemiLagrangian is serial: the fit at a node near a " + "partition seam needs the arrivals on the other rank") + if _unsupported: + warnings.warn(f"ForwardNodesSemiLagrangian ignores {sorted(_unsupported)}", stacklevel=2) + self.vtype = vtype + self.mesh = mesh + self.degree = int(degree) + self.continuous = True + self.order = 1 + self.theta = float(theta) + self.V_fn = V_fn + self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) + self._init_history_tracking(1) + if varsymbol is None: + varsymbol = rf"u_{{ [{self.instance_number}] }}" + self.psi_star = [ + uw.discretisation.MeshVariable( + f"psi_star_fwn_{self.instance_number}", mesh, vtype=vtype, degree=self.degree, + continuous=True, varsymbol=rf"{{ {varsymbol}^{{ * }} }}", + units=_history_units(self._psi_fn, units)) + ] + self._components = _storage_components(vtype, tuple(self.psi_star[0].sym.shape)) + # the per-cell fit: a discontinuous variable of the field's degree, and + # the interior lattice the values are also launched from + from underworld3.utilities.cell_polynomial_projection import CellPolynomialProjector + self._fit_var = uw.discretisation.MeshVariable( + f"fit_fwn_{self.instance_number}", mesh, vtype=vtype, degree=self.degree, continuous=False) + self._projector = CellPolynomialProjector(self._fit_var) + self._interior = np.asarray(mesh._get_coords_for_basis(self.degree + 1, continuous=False) + ).reshape(-1, mesh.cdim) + self._fit = None + self._n_v = 2 + self._init_coefficient_expressions(1, self.theta, with_exp=True) + self._register_with_default_model() + + @property + def psi_fn(self): + return self._psi_fn + + @psi_fn.setter + def psi_fn(self, new_fn): + self._psi_fn = new_fn if isinstance(new_fn, sympy.Matrix) else sympy.Matrix([[new_fn]]) + + @property + def state(self) -> "DDtForwardNodesState": + return DDtForwardNodesState( + **self._core_state_kwargs(), + psi_star_var_names=[ps.clean_name for ps in self.psi_star], + ) + + @state.setter + def state(self, s: "DDtForwardNodesState") -> None: + self._validate_state_schema(s, DDtForwardNodesState) + self._validate_psi_star_names(s.psi_star_var_names) + self._restore_core_state(s, am_theta=self.theta) + + # ------------------------------------------------------------------ + def _nodes(self): + return np.asarray(_to_nondim_ndarray(self.psi_star[0].coords)).reshape(-1, self.mesh.cdim) + + def _values_at(self, expr, X): + """``expr`` at the points ``X``, one column per stored component.""" + expr = sympy.Matrix(expr) + return np.column_stack([ + np.asarray(_to_nondim_ndarray(uw.function.evaluate(expr[i, j], X))).reshape(-1) + for (i, j) in self._components]) + + def _reconstruct(self, arrivals, values, nodes): + """Fit the arrivals in each cell at the field's degree; the fit at the nodes.""" + inside = np.asarray(self.mesh.points_in_domain(arrivals), dtype=bool) + self._fit = self._projector.fit(arrivals[inside], values[inside], old=self._fit) + return np.nan_to_num(self._projector.interpolate(self._fit, nodes)) + + # ------------------------------------------------------------------ + def initialise_history(self): + """Start from the current field at the nodes. A history already placed + by :meth:`commit_flux_to_history` is the start, and is kept.""" + self.characteristics.initialise_levels(self._n_v) + if not self._history_committed: + self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) + self._history_initialised = True + + def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, **_ignored): + """Carry the field forward one step from the nodes and rebuild it there. + + ``store_result=False`` says the store already holds the values to launch + (placed by :meth:`commit_flux_to_history`); otherwise the tracked field + is read at the nodes first. + """ + self._dt = dt = self._nondim_timestep(dt) + if not self._history_initialised: + self.initialise_history() + _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) + _update_am_values(self._am_coeffs, self.effective_order, self.theta) + trace = self.characteristics + if self._owns_characteristics: + trace.begin_step(dt) + nodes = self._nodes() + launch = np.vstack([nodes, self._interior]) + # the tracked field, or (when a flux was committed) the store, which is + # a polynomial inside each element, so its interior values are exact + source = self._psi_fn if store_result else self.psi_star[0].sym + values = np.vstack([self._values_at(source, nodes) if store_result else np.array(self.psi_star[0].data), + self._values_at(source, self._interior)]) + key = (_basis_key_of(self.psi_star[0]), "launch") + arrivals = np.asarray(trace.departure_points(key, launch, (("first", 0, -float(dt)),), + evalf=evalf, clamp_final=False)) + self.psi_star[0].data[:, :] = self._reconstruct(arrivals, values, nodes) + arrivals = arrivals[:nodes.shape[0]] + if self._inflow_value is not None: + # a node whose back-trace leaves the domain holds fluid that entered this step + back = 2.0 * nodes - arrivals + restored = np.asarray(self.mesh.return_coords_to_bounds(back.copy())).reshape(back.shape) + entered = np.any(restored != back, axis=1) + if entered.any(): + self._write_inflow(self.psi_star[0], nodes, entered) + if self._owns_characteristics: + trace.finish_step() + + def update(self, dt, evalf=False, verbose=False, **kwargs): + self.update_pre_solve(dt, evalf=evalf, verbose=verbose, **kwargs) + + def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): + self._dt = dt = self._nondim_timestep(dt) + self._dt_history[0] = dt + if self._n_solves_completed < self.order: + self._n_solves_completed += 1 + + def commit_flux_to_history(self, flux, verbose=False): + """Read the new flux at the nodes and leave it in the store until the next carry.""" + self.psi_star[0].data[:, :] = self._values_at(flux, self._nodes()) + self._history_committed = True diff --git a/tests/test_1102_forward_nodes_rotating_gaussian.py b/tests/test_1102_forward_nodes_rotating_gaussian.py new file mode 100644 index 00000000..fe232067 --- /dev/null +++ b/tests/test_1102_forward_nodes_rotating_gaussian.py @@ -0,0 +1,54 @@ +"""The forward-from-nodes history on the rotating diffusing Gaussian. + +Launched from the field's own nodes and from a lattice inside every element +(the field's interpolant is a polynomial there, so both are known exactly), +carried one step forward, fitted per cell at the field's degree and read back at +the nodes. Used as the DuDt of the SLCN advection-diffusion solver, half a +revolution, kappa 0.01, dt 0.02, P2 temperature, against the test_1101 fixture: +the backward nodal history gives 2.52e-2 on the disc and 6.43e-2 on the box +(where the flow crosses all four walls). +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from test_1101_advdiff_swarm_rotating_gaussian import SIGMA, KAPPA, T_END, EXACT_PEAK, _disc + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _run(mesh, walls): + x, y = mesh.X + sol = uw.analytic.RotatingGaussian(mesh, sigma=SIGMA, centre_radius=0.5, omega=1.0, diffusivity=KAPPA) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) + V = sympy.Matrix([[-y, x]]) + duDt = uw.systems.ddt.ForwardNodesSemiLagrangian(mesh, T.sym, V, vtype=uw.VarType.SCALAR, degree=T.degree) + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, DuDt=duDt, order=1) + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = KAPPA + for wall in walls: + adv.add_dirichlet_bc(0.0, wall) + dt = 0.02 + nsteps = int(round(T_END / dt)); dt = T_END / nsteps + for _ in range(nsteps): + adv.solve(timestep=dt) + return float(sol.error(sol.at(T_END), T, norm="integral")), float(np.asarray(T.data).max()) + + +def test_forward_from_nodes_on_the_disc(): + uw.reset_default_model() + err, peak = _run(_disc(24), ("Upper",)) + assert abs(err - 0.01996) < 0.002, err # BASELINE (2026-09-26) + assert abs(peak - EXACT_PEAK) < 0.002, (peak, EXACT_PEAK) + + +def test_forward_from_nodes_holds_the_box(): + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), + cellSize=2.0 / 24, qdegree=3, regular=False) + err, peak = _run(mesh, ("Left", "Right", "Top", "Bottom")) + assert abs(err - 0.0617) < 0.006, err # BASELINE (2026-09-26) + assert abs(peak - EXACT_PEAK) < 0.005, (peak, EXACT_PEAK) From c434d8ac293daed57baacac03e71c1dcc1ac3092 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 26 Sep 2026 19:35:25 -0700 Subject: [PATCH 2/5] Semi-Lagrangian histories named by trace and launch; one entry point; forward from nodes as an advection-diffusion option ddt.SemiLagrangian(mesh, psi_fn, V_fn, vtype, trace=, launch=) selects one of four schemes: BackwardNodesSemiLagrangian (was SemiLagrangian), BackwardIntegrationPointsSemiLagrangian (was IntegrationPointSemiLagrangian), ForwardIntegrationPointsSemiLagrangian (was ForwardSemiLagrangian) and ForwardNodesSemiLagrangian. A keyword the chosen scheme does not take is a TypeError. The old class names resolve with a FutureWarning. stress_transport takes "backward_nodes" (default), "backward_integration_points", "forward_integration_points", "lagrangian", "eulerian"; the old strings map with a FutureWarning. AdvDiffusionSLCN(transport=...) chooses the value history among the four semi-Lagrangian schemes. ForwardNodesSemiLagrangian carries a field known at its nodes and refuses a flux: a stress is formed at the integration points. A cell nothing reached now keeps the field it launched (no fit carried between steps); inflow is detected by the back-reflected point leaving the domain; a moved mesh is refused. ForwardIntegrationPointsSemiLagrangian takes degree and refuses anything but 1. BackwardNodesSemiLagrangian's continuous now defaults to True. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- docs/api/systems_ddt.md | 34 ++- .../subsystems/integration-point-variables.md | 8 +- docs/developer/subsystems/stress-transport.md | 17 +- src/underworld3/constitutive_models.py | 8 +- .../cython/petsc_generic_snes_solvers.pyx | 16 +- src/underworld3/systems/__init__.py | 7 +- src/underworld3/systems/ddt.py | 219 ++++++++++++------ src/underworld3/systems/solver_template.py | 6 +- src/underworld3/systems/solvers.py | 213 +++++++++++------ .../test_1062_forward_stress_history_mpi.py | 4 +- tests/test_0066_integration_point_slcn.py | 34 +-- tests/test_1056_units_slcn_traceback.py | 2 +- tests/test_1059_stress_transport.py | 87 +++++-- tests/test_1060_stress_store_smoothing.py | 10 +- tests/test_1061_stress_forward_history.py | 4 +- tests/test_1063_stress_history_restart.py | 2 +- tests/test_1064_stress_history_units.py | 8 +- tests/test_1065_log_conformation.py | 10 +- ...st_1102_forward_nodes_rotating_gaussian.py | 6 +- 19 files changed, 467 insertions(+), 228 deletions(-) diff --git a/docs/api/systems_ddt.md b/docs/api/systems_ddt.md index d9ded7df..e3086044 100644 --- a/docs/api/systems_ddt.md +++ b/docs/api/systems_ddt.md @@ -37,7 +37,37 @@ bypasses the order ramp, so the first solve runs at full BDF order. ### SemiLagrangian ```{eval-rst} -.. autoclass:: underworld3.systems.ddt.SemiLagrangian +.. autofunction:: underworld3.systems.ddt.SemiLagrangian +``` + +### BackwardNodesSemiLagrangian + +```{eval-rst} +.. autoclass:: underworld3.systems.ddt.BackwardNodesSemiLagrangian + :members: + :show-inheritance: +``` + +### BackwardIntegrationPointsSemiLagrangian + +```{eval-rst} +.. autoclass:: underworld3.systems.ddt.BackwardIntegrationPointsSemiLagrangian + :members: + :show-inheritance: +``` + +### ForwardIntegrationPointsSemiLagrangian + +```{eval-rst} +.. autoclass:: underworld3.systems.ddt.ForwardIntegrationPointsSemiLagrangian + :members: + :show-inheritance: +``` + +### ForwardNodesSemiLagrangian + +```{eval-rst} +.. autoclass:: underworld3.systems.ddt.ForwardNodesSemiLagrangian :members: :show-inheritance: ``` @@ -63,6 +93,6 @@ bypasses the order ramp, so the first solve runs at full BDF order. The following aliases are available via ``underworld3.systems``: - ``Lagrangian_DDt`` → {class}`~underworld3.systems.ddt.Lagrangian` -- ``SemiLagragian_DDt`` → {class}`~underworld3.systems.ddt.SemiLagrangian` +- ``SemiLagragian_DDt`` → {class}`~underworld3.systems.ddt.BackwardNodesSemiLagrangian` - ``Lagrangian_Swarm_DDt`` → {class}`~underworld3.systems.ddt.Lagrangian_Swarm` - ``Eulerian_DDt`` → {class}`~underworld3.systems.ddt.Eulerian` diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index a8e106e8..1a78882e 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -145,7 +145,7 @@ reading it, `evaluate`, the guards). ## Semi-Lagrangian history on the integration points -`uw.systems.ddt.IntegrationPointSemiLagrangian` is the SLCN history built on +`uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian` is the SLCN history built on this variable. Its slots `psi_star[k]` are integration-point variables, so the value the weak form sees at each integration point is the solution from `k+1` steps ago evaluated exactly at the departure point of that @@ -158,7 +158,7 @@ evaluating the snapshot from time `n-k` at the foot. Every slot carries one evaluation error rather than one per generation. ```python -DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V_fn, degree=2, order=1) +DuDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V_fn, degree=2, order=1) adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V_fn, DuDt=DuDt, order=1) ``` @@ -182,11 +182,11 @@ component per integration point: ```python # a momentum history for Navier-Stokes -DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian( +DuDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, v, v.sym, vtype=uw.VarType.VECTOR, degree=2, order=2) # a viscoelastic stress history -DFDt = uw.systems.ddt.IntegrationPointSemiLagrangian( +DFDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, stress, v.sym, vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) DFDt.psi_star[0].sym # a 2x2 symbolic matrix diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index a1def6f9..f5652da8 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -8,7 +8,7 @@ what limits it, and how to keep a run inside those limits. ```python stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) -stokes.stress_transport = "integration_point" # or "semi_lagrangian" (the default), "forward", "eulerian" +stokes.stress_transport = "forward_integration_points" # see the table for the others stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( stokes.Unknowns, order=1, integrator="bdf", objective_rate="upper_convected") stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_p @@ -19,11 +19,20 @@ stokes.constitutive_model.Parameters.dt_elastic = dt ## The five histories +The four semi-Lagrangian names are `_`, the arguments of +`uw.systems.ddt.SemiLagrangian(..., trace=, launch=)`: a backward trace follows the +characteristic back from each storage point and samples the old stress at the foot; a +forward trace carries the old stress from where it is known and fits the arrivals in +each cell. The fourth combination, forward from nodes, carries a field known at its +nodes (a temperature, say); a stress is formed at the integration points, so it is not +a stress history. The former names `semi_lagrangian`, `integration_point` and `forward` +are accepted with a warning. + | `stress_transport` | storage | carried by | stable at | fails by | |---|---|---|---|---| -| `semi_lagrangian` (nodal) | continuous P1 at the vertices | vertex trace-back, interpolation at the foot | any Courant number | excess stress in the first cells off a no-slip wall; on the confined cylinder that excess loses the conformation and the solve hangs | -| `integration_point` | continuous P1 store, sampled at the quadrature points | trace-back of every quadrature point | Courant near one, or below one with store smoothing | a cell-scale mode of the stress that grows below Courant one when the solvent viscosity is small | -| `forward` | discontinuous P1 per cell, fitted from the arrivals | fixed launch set of interior points (the integration points), one forward trajectory a step; the flux is read back at the launch points through a continuous P1 projection; an inflow cell's uncovered share is filled with the inflow value | the cylinder walls at dt 0.04; below Courant one with `flux_smoothing` at c = 0.023 (Waters-King 1/16, dt 0.0125: 0.9543 at t 1 and 0.5185 at t 6.5, against nodal 0.9622 and 0.5171) | the same cell-scale mode as the integration-point history without that smoothing (diverges at t 2.4 there); first order only; does not cross a periodic seam or follow a moving mesh | +| `backward_nodes` (the default) | continuous P1 at the vertices | vertex trace-back, interpolation at the foot | any Courant number | excess stress in the first cells off a no-slip wall; on the confined cylinder that excess loses the conformation and the solve hangs | +| `backward_integration_points` | continuous P1 store, sampled at the quadrature points | trace-back of every quadrature point | Courant near one, or below one with store smoothing | a cell-scale mode of the stress that grows below Courant one when the solvent viscosity is small | +| `forward_integration_points` | discontinuous P1 per cell, fitted from the arrivals | fixed launch set of interior points (the integration points), one forward trajectory a step; the flux is read back at the launch points through a continuous P1 projection; an inflow cell's uncovered share is filled with the inflow value | the cylinder walls at dt 0.04; below Courant one with `flux_smoothing` at c = 0.023 (Waters-King 1/16, dt 0.0125: 0.9543 at t 1 and 0.5185 at t 6.5, against nodal 0.9622 and 0.5171) | the same cell-scale mode as the integration-point history without that smoothing (diverges at t 2.4 there); first order only; does not cross a periodic seam or follow a moving mesh | | `lagrangian` (particles) | a swarm the solver owns and advects, one value per particle, read through a discontinuous cells proxy | the material points themselves: the constitutive flux is evaluated at the particles each step and never projected back to the mesh; a particle that entered through an inflow takes the inflow value | any Courant number; no numerical diffusion of the history | the cost and bookkeeping of a swarm, and a proxy that needs its cells kept populated (population control refills them); the conformation check does not read a per-point tensor from it | | `eulerian` (SUPG grid) | continuous P1 | assembled transport equation with streamline upwinding | with DEVSS | without DEVSS the velocity block loses its preconditioner as the stress grows | diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index a6bb2a52..6b65eaef 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -47,7 +47,7 @@ from underworld3.utilities._api_tools import uw_object from underworld3.swarm import IndexSwarmVariable from underworld3.discretisation import MeshVariable -from underworld3.systems.ddt import SemiLagrangian as SemiLagrangian_DDt +from underworld3.systems.ddt import BackwardNodesSemiLagrangian from underworld3.systems.ddt import _bdf_coefficients, _as_float from underworld3.function.quantities import UWQuantity from underworld3.systems.ddt import Lagrangian as Lagrangian_DDt @@ -459,7 +459,7 @@ def DuDt(self): Returns ------- - SemiLagrangian_DDt or Lagrangian_DDt or None + BackwardNodesSemiLagrangian or Lagrangian_DDt or None The material derivative operator, or None if not set. """ return self._DuDt @@ -467,7 +467,7 @@ def DuDt(self): @DuDt.setter def DuDt( self, - DuDt_value: Union[SemiLagrangian_DDt, Lagrangian_DDt], + DuDt_value: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt], ): """Set the material derivative operator for the unknown.""" self._DuDt = DuDt_value @@ -483,7 +483,7 @@ def DFDt(self): @DFDt.setter def DFDt( self, - DFDt_value: Union[SemiLagrangian_DDt, Lagrangian_DDt], + DFDt_value: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt], ): """Set the material derivative operator for flux history.""" self._DFDt = DFDt_value diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c7d53b68..aed254a3 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3352,8 +3352,8 @@ class SNES_Scalar(SolverBaseClass): u_Field : uw.discretisation.MeshVariable = None, degree: int = 2, verbose = False, - DuDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, - DFDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DuDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DFDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, ): super().__init__(mesh) @@ -4265,8 +4265,8 @@ class SNES_Vector(SolverBaseClass): u_Field : uw.discretisation.MeshVariable = None, degree = 2, verbose = False, - DuDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, - DFDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DuDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DFDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, ): @@ -5268,8 +5268,8 @@ class SNES_MultiComponent(SolverBaseClass): n_components : int = None, degree = 2, verbose = False, - DuDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, - DFDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DuDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DFDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, ): super().__init__(mesh) @@ -6026,8 +6026,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): degree : Optional[int] = 2, p_continuous : Optional[bool] = True, verbose : Optional[bool] =False, - DuDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, - DFDt : Union[uw.systems.ddt.SemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DuDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, + DFDt : Union[uw.systems.ddt.BackwardNodesSemiLagrangian, uw.systems.ddt.Lagrangian] = None, ): diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index 5ba37b93..300aea38 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -93,9 +93,10 @@ # are the Lagrangian implementations actually distinct in reality ? from .ddt import Lagrangian as Lagrangian_DDt -from .ddt import SemiLagrangian as SemiLagragian_DDt -from .ddt import IntegrationPointSemiLagrangian as IntegrationPointSemiLagrangian_DDt -from .ddt import ForwardSemiLagrangian as ForwardSemiLagrangian_DDt +# former names of three of the four schemes ddt.SemiLagrangian selects +from .ddt import BackwardNodesSemiLagrangian as SemiLagragian_DDt +from .ddt import BackwardIntegrationPointsSemiLagrangian as IntegrationPointSemiLagrangian_DDt +from .ddt import ForwardIntegrationPointsSemiLagrangian as ForwardSemiLagrangian_DDt from .ddt import Lagrangian_Swarm as Lagrangian_Swarm_DDt from .ddt import Eulerian as Eulerian_DDt from .ddt import EulerianSUPG as EulerianSUPG_DDt diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 148da5f3..96d942ba 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -50,6 +50,7 @@ underworld3.systems.solvers : PDE solvers using these time derivatives. """ +import inspect import math import warnings @@ -132,7 +133,7 @@ class DDtEulerianState(_DDtCoreState): @dataclass class DDtSemiLagrangianState(_DDtCoreState): - """Snapshot of a :class:`SemiLagrangian` DDt instance. + """Snapshot of a :class:`BackwardNodesSemiLagrangian` DDt instance. Like :class:`DDtEulerianState`, plus an optional ``forcing_star`` variable (when ``with_forcing_history=True``) used by ETD-2 @@ -146,7 +147,7 @@ class DDtSemiLagrangianState(_DDtCoreState): @dataclass class DDtIntegrationPointState(_DDtCoreState): - """Snapshot of an :class:`IntegrationPointSemiLagrangian` instance: the + """Snapshot of an :class:`BackwardIntegrationPointsSemiLagrangian` instance: the point-value slots and their nodal snapshots are mesh variables captured by name; this carries the bookkeeping.""" psi_star_var_names: list[str] = field(default_factory=list) @@ -157,7 +158,7 @@ class DDtIntegrationPointState(_DDtCoreState): @dataclass class DDtForwardState(_DDtCoreState): - """Snapshot of a :class:`ForwardSemiLagrangian` instance: the fitted + """Snapshot of a :class:`ForwardIntegrationPointsSemiLagrangian` instance: the fitted field and the launch values are mesh variables captured by name.""" psi_star_var_names: list[str] = field(default_factory=list) launch_var_name: str = "" @@ -569,7 +570,7 @@ class _DDtBase(uw_object): r"""Shared machinery for the DDt history-manager flavors. The five flavors (:class:`Symbolic`, :class:`Eulerian`, - :class:`SemiLagrangian`, :class:`Lagrangian`, + :class:`BackwardNodesSemiLagrangian`, :class:`Lagrangian`, :class:`Lagrangian_Swarm`) share the same BDF/Adams-Moulton coefficient bookkeeping, effective-order startup ramp, fixed-structure ``bdf()`` / ``adams_moulton_flux()`` expressions, @@ -989,7 +990,7 @@ def _refresh_source_snapshot(self): commits_flux_in_post_solve = False #: The forcing (strain-rate) history the second-order exponential - #: integrator reads; only :class:`SemiLagrangian` allocates one, on + #: integrator reads; only :class:`BackwardNodesSemiLagrangian` allocates one, on #: request. ``None`` means the integrator runs at first order (#739). forcing_star = None @@ -1073,8 +1074,8 @@ def inflow_value(self): the relaxed stress of the incoming flow. :class:`EulerianSUPG` compiles the value into a boundary term of its - transport solve; :class:`IntegrationPointSemiLagrangian` gives it to a - departure point restored to the boundary; :class:`ForwardSemiLagrangian` + transport solve; :class:`BackwardIntegrationPointsSemiLagrangian` gives it to a + departure point restored to the boundary; :class:`ForwardIntegrationPointsSemiLagrangian` fills the uncovered share of an inflow cell with it; :class:`Lagrangian` gives it to every particle that entered through an inflow: one whose back-trace over the step, or over one cell for a @@ -1102,7 +1103,7 @@ def inflow_value(self, value): "restores an out-of-bounds departure point to the boundary " "and reads the transported field there, which constrains " "the inflow but is not the value you set. EulerianSUPG, " - "IntegrationPointSemiLagrangian, ForwardSemiLagrangian and " + "BackwardIntegrationPointsSemiLagrangian, ForwardIntegrationPointsSemiLagrangian and " "Lagrangian apply it (#733, #783).", stacklevel=2) self._inflow_value = value @@ -1917,7 +1918,7 @@ class EulerianSUPG(Eulerian): component-wise to a scalar, a vector or a tensor unknown, and the streamline-upwind Petrov-Galerkin flux :math:`\tau\,R\otimes\mathbf{a}` of the solver's strong residual :math:`R`. The same solver takes a - :class:`SemiLagrangian` history in its place: that flavour answers zero + :class:`BackwardNodesSemiLagrangian` history in its place: that flavour answers zero for the advection and the flux because its history is already traced back along the characteristics. @@ -2559,12 +2560,12 @@ def _matrix_of(V): # TODO(BUG): a non-symmetric psi_fn under vtype=SYM_TENSOR is silently # reduced here, and this class keeps the LOWER entry where -# IntegrationPointSemiLagrangian keeps the UPPER one. Measured 2026-09-10 on +# BackwardIntegrationPointsSemiLagrangian keeps the UPPER one. Measured 2026-09-10 on # [[1+x, 2+y], [100.0, 3+x*y]]: nodal psi_star -> [[1.45, 100.0], [100.0, 3.21]], # integration-point -> [[1.46, 2.47], [2.47, 3.21]]. Neither averages and # neither warns. The integration-point path now warns; this one should too, # and the two should agree on which triangle wins. -class SemiLagrangian(_DDtBase): +class BackwardNodesSemiLagrangian(_DDtBase): r""" Semi-Lagrangian history manager. @@ -2604,7 +2605,7 @@ class SemiLagrangian(_DDtBase): sampling discretisation (the former ``swarm_degree`` / ``swarm_continuous`` were never read, issue #704). Denser sampling at the integration points is a separate history manager - (``IntegrationPointSemiLagrangian``, PR #703). + (``BackwardIntegrationPointsSemiLagrangian``, PR #703). varsymbol : str, optional LaTeX symbol for display. verbose : bool, default=False @@ -2711,7 +2712,7 @@ def __init__( V_fn: sympy.Function, vtype: uw.VarType, degree: int, - continuous: bool, + continuous: bool = True, varsymbol: Optional[str] = None, verbose: Optional[bool] = False, bcs=[], @@ -4734,7 +4735,7 @@ def _storage_components(vtype, shape): return [(i, j) for i in range(shape[0]) for j in range(shape[1])] -class IntegrationPointSemiLagrangian(_DDtBase): +class BackwardIntegrationPointsSemiLagrangian(_DDtBase): r"""Semi-Lagrangian history stored at the mesh integration points. The history slots ``psi_star[k]`` are @@ -4743,12 +4744,12 @@ class IntegrationPointSemiLagrangian(_DDtBase): solution from ``k+1`` steps ago evaluated **exactly** at the departure point of that integration point. There is no nodal history field and no second interpolation: only the FE solution's own error remains in the - advected term. Compare :class:`SemiLagrangian`, which samples at the + advected term. Compare :class:`BackwardNodesSemiLagrangian`, which samples at the nodes, stores a nodal ``psi_star`` and lets the assembler interpolate it to the integration points. Because a delta field cannot be sampled off its points, the chain - ``psi_star[k] <- psi_star[k-1]`` of :class:`SemiLagrangian` is replaced + ``psi_star[k] <- psi_star[k-1]`` of :class:`BackwardNodesSemiLagrangian` is replaced by nodal **snapshots** of the solution and of the velocity at the last ``order`` times. Slot ``k`` is filled by tracing ``k+1`` segments back from every integration point (segment ``j`` with the velocity at time @@ -4766,13 +4767,13 @@ class IntegrationPointSemiLagrangian(_DDtBase): What is not here (yet): units-aware velocity reduction, ALE / old-frame trace-back. Use - :class:`SemiLagrangian` for those, or :class:`Lagrangian_Swarm` when the + :class:`BackwardNodesSemiLagrangian` for those, or :class:`Lagrangian_Swarm` when the history should ride on particles rather than on the rule. Parameters ---------- mesh, psi_fn, V_fn, degree, continuous, varsymbol, verbose, bcs, order, theta - As for :class:`SemiLagrangian`. ``psi_fn`` may be a ``MeshVariable`` + As for :class:`BackwardNodesSemiLagrangian`. ``psi_fn`` may be a ``MeshVariable`` (its nodal data is then copied into the snapshot rather than re-evaluated) or an expression, of any ``vtype``. ``V_fn`` may be any expression (``-v``, ``v/2``, ``c(t) v``); the @@ -4862,7 +4863,7 @@ def __init__( self.store_smoothing = store_smoothing self.mesh = mesh if bcs: - raise ValueError("IntegrationPointSemiLagrangian applies no boundary conditions to its " + raise ValueError("BackwardIntegrationPointsSemiLagrangian applies no boundary conditions to its " "store; an inflow is set through inflow_value") self.bcs = [] self.verbose = verbose @@ -4923,7 +4924,7 @@ def __init__( ) if len(self._components) != self.num_components: raise RuntimeError( - f"IntegrationPointSemiLagrangian: {vtype} maps " + f"BackwardIntegrationPointsSemiLagrangian: {vtype} maps " f"{len(self._components)} components onto " f"{self.num_components} stored columns" ) @@ -5070,7 +5071,7 @@ def _check_rule_oversampling(self, degree): need = self._qdegree_with_at_least(local_dofs + 1) want = self._qdegree_with_at_least(2 * local_dofs) raise RuntimeError( - f"IntegrationPointSemiLagrangian: the mesh rule has {Nq} points per cell " + f"BackwardIntegrationPointsSemiLagrangian: the mesh rule has {Nq} points per cell " f"but a degree-{degree} history has {local_dofs} local dofs; the " "least-squares fit is not oversampled and is unstable at small Courant " f"number. For this cell type and history degree the rule needs at least " @@ -5079,7 +5080,7 @@ def _check_rule_oversampling(self, degree): ) if Nq < 2 * local_dofs: warnings.warn( - f"IntegrationPointSemiLagrangian: {Nq} rule points per cell for " + f"BackwardIntegrationPointsSemiLagrangian: {Nq} rule points per cell for " f"{local_dofs} local dofs is under 2x oversampling: weakly unstable under pure " "advection (growth ~1.005/step at 1.5x, Courant 0.25) and stable with physical " "diffusion at cell Peclet <= 100. 2x (qdegree 3 for P2 on triangles) is neutral.", @@ -5121,7 +5122,7 @@ def _check_psi_shape(self, psi_fn): expected = _psi_shape_for(self.vtype, self.mesh.cdim) if expected is not None and tuple(psi_fn.shape) != expected: raise ValueError( - f"IntegrationPointSemiLagrangian: psi_fn has shape " + f"BackwardIntegrationPointsSemiLagrangian: psi_fn has shape " f"{tuple(psi_fn.shape)} but vtype={self.vtype} on a cdim=" f"{self.mesh.cdim} mesh needs {expected}. Pass the vtype that " "matches the field, or reshape psi_fn." @@ -5137,7 +5138,7 @@ def _check_psi_shape(self, psi_fn): if supplied_var is not None and expected is not None: if int(supplied_var.num_components) != wanted: raise ValueError( - f"IntegrationPointSemiLagrangian: psi_fn stores " + f"BackwardIntegrationPointsSemiLagrangian: psi_fn stores " f"{supplied_var.num_components} components but vtype=" f"{self.vtype} stores {wanted}. A full tensor and a " "symmetric tensor share a shape; pass the vtype the field " @@ -5146,7 +5147,7 @@ def _check_psi_shape(self, psi_fn): components = getattr(self, "num_components", None) if components is not None and expected is not None and wanted != components: raise ValueError( - f"IntegrationPointSemiLagrangian: psi_fn needs {wanted} stored " + f"BackwardIntegrationPointsSemiLagrangian: psi_fn needs {wanted} stored " f"components but this history has {components}; vtype=" f"{self.vtype} is probably not the vtype of the field." ) @@ -5164,7 +5165,7 @@ def _check_psi_shape(self, psi_fn): import warnings warnings.warn( - f"IntegrationPointSemiLagrangian: psi_fn is not symmetric at " + f"BackwardIntegrationPointsSemiLagrangian: psi_fn is not symmetric at " f"{dropped} but vtype=SYM_TENSOR stores only the upper " "triangle, so the lower entries are discarded (not averaged). " "Symmetrise psi_fn explicitly, or use VarType.TENSOR.", @@ -5182,7 +5183,7 @@ def _object_viewer(self): def _nudged_node_coords(self, var): """ND node coordinates of ``var`` moved 0.1 % toward their cell centroids so boundary nodes locate unambiguously (see - :meth:`SemiLagrangian._centroid_shifted_node_coords`).""" + :meth:`BackwardNodesSemiLagrangian._centroid_shifted_node_coords`).""" coords = np.asarray(var.coords_nd) cellid = self.mesh.get_closest_cells(coords).reshape(-1) cent = np.asarray(self.mesh._centroids)[cellid] @@ -5351,7 +5352,7 @@ def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): self._n_solves_completed += 1 -class ForwardSemiLagrangian(_DDtBase): +class ForwardIntegrationPointsSemiLagrangian(_DDtBase): r"""Semi-Lagrangian history carried forward from a fixed set of launch points inside the cells, read by the weak form through a per-cell fit. @@ -5403,6 +5404,7 @@ def __init__( psi_fn, V_fn, vtype=VarType.SCALAR, + degree: int = 1, varsymbol: Optional[str] = None, order: int = 1, theta: float = 0.5, @@ -5411,11 +5413,14 @@ def __init__( ): super().__init__() if order != 1: - raise NotImplementedError("ForwardSemiLagrangian carries one level; order must be 1") + raise NotImplementedError("ForwardIntegrationPointsSemiLagrangian carries one level; order must be 1") + if degree != 1: + raise NotImplementedError("ForwardIntegrationPointsSemiLagrangian fits a linear polynomial per " + f"cell; degree must be 1, not {degree}") if mesh.cdim != mesh.dim: - raise NotImplementedError("ForwardSemiLagrangian fits in the embedding coordinates; no manifolds") + raise NotImplementedError("ForwardIntegrationPointsSemiLagrangian fits in the embedding coordinates; no manifolds") if _unsupported: - warnings.warn(f"ForwardSemiLagrangian ignores {sorted(_unsupported)}: it has one level, " + warnings.warn(f"ForwardIntegrationPointsSemiLagrangian ignores {sorted(_unsupported)}: it has one level, " "a linear fit per cell and no smoothing or monotone option", stacklevel=2) self.vtype = vtype self.mesh = mesh @@ -5427,7 +5432,7 @@ def __init__( self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) expected = _psi_shape_for(vtype, mesh.cdim) if expected is not None and tuple(self._psi_fn.shape) != expected: - raise ValueError(f"ForwardSemiLagrangian: psi_fn has shape {tuple(self._psi_fn.shape)} " + raise ValueError(f"ForwardIntegrationPointsSemiLagrangian: psi_fn has shape {tuple(self._psi_fn.shape)} " f"but vtype={vtype} on a cdim={mesh.cdim} mesh needs {expected}") self._init_history_tracking(1) if varsymbol is None: @@ -5452,7 +5457,7 @@ def __init__( self._launch = np.array(np.asarray(launch_var.coords_nd).reshape(-1, mesh.cdim)) self._nq = int(launch_var.num_points_per_cell) if self._nq < mesh.dim + 1: - raise ValueError(f"ForwardSemiLagrangian needs at least {mesh.dim + 1} integration points per " + raise ValueError(f"ForwardIntegrationPointsSemiLagrangian needs at least {mesh.dim + 1} integration points per " f"cell for a linear fit; this mesh's rule has {self._nq} (raise qdegree)") w_ref = np.asarray(mesh.integration_rule.getData()[1]).reshape(-1) self._cell_measure = self._cell_measures() @@ -5530,7 +5535,7 @@ def _boundary_faces(self): # the mesh's own boundary label tells them apart. label = dm.getLabel("All_Boundaries") if dm.hasLabel("All_Boundaries") else None if label is None and uw.mpi.size > 1: - raise RuntimeError("ForwardSemiLagrangian: the mesh has no All_Boundaries label, so " + raise RuntimeError("ForwardIntegrationPointsSemiLagrangian: the mesh has no All_Boundaries label, so " "partition faces cannot be told from domain faces in parallel") for f in range(f0, f1): if dm.getSupportSize(f) != 1: @@ -5573,7 +5578,7 @@ def psi_fn(self, new_fn): new_fn = new_fn if isinstance(new_fn, sympy.Matrix) else sympy.Matrix([[new_fn]]) expected = _psi_shape_for(self.vtype, self.mesh.cdim) if expected is not None and tuple(new_fn.shape) != expected: - raise ValueError(f"ForwardSemiLagrangian: psi_fn has shape {tuple(new_fn.shape)}, needs {expected}") + raise ValueError(f"ForwardIntegrationPointsSemiLagrangian: psi_fn has shape {tuple(new_fn.shape)}, needs {expected}") self._psi_fn = new_fn def _object_viewer(self): @@ -5729,7 +5734,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** self._dt = dt = self._nondim_timestep(dt) if self._geometry_stamp() != self._launch_geometry: raise NotImplementedError( - "ForwardSemiLagrangian: the launch set, cell measures and boundary faces were " + "ForwardIntegrationPointsSemiLagrangian: the launch set, cell measures and boundary faces were " "built for the mesh as it was, and the mesh has moved or been re-meshed since; " "this flavour does not follow a changing mesh") if not self._history_initialised: @@ -5786,14 +5791,16 @@ class ForwardNodesSemiLagrangian(_DDtBase): reads that fit back at the nodes. The fit never reaches across an element boundary, where the interpolant has a kink. A cell with too few arrivals, or arrivals on a line, falls back to a linear fit over the nearest arrivals; a - cell nothing reached keeps its previous fit (see + cell nothing reached keeps the field it launched (see :class:`~underworld3.utilities.cell_polynomial_projection.CellPolynomialProjector`). Arrivals that leave the domain are dropped; a node the flow reached from outside takes :attr:`inflow_value` when one is set. - The integration-point counterpart is :class:`ForwardSemiLagrangian`, which - launches from the quadrature points, where a flux (a stress) is formed. - First order; serial. + The integration-point counterpart is :class:`ForwardIntegrationPointsSemiLagrangian`, which + launches from the quadrature points, where a flux (a stress) is formed, and + is the forward scheme for one: this class carries a field and refuses a flux. + First order; serial (a node near a partition seam needs arrivals from the + other rank); a fixed mesh. Parameters ---------- @@ -5850,7 +5857,6 @@ def __init__(self, mesh, psi_fn, V_fn, vtype=VarType.SCALAR, degree=1, varsymbol self._projector = CellPolynomialProjector(self._fit_var) self._interior = np.asarray(mesh._get_coords_for_basis(self.degree + 1, continuous=False) ).reshape(-1, mesh.cdim) - self._fit = None self._n_v = 2 self._init_coefficient_expressions(1, self.theta, with_exp=True) self._register_with_default_model() @@ -5878,7 +5884,7 @@ def state(self, s: "DDtForwardNodesState") -> None: # ------------------------------------------------------------------ def _nodes(self): - return np.asarray(_to_nondim_ndarray(self.psi_star[0].coords)).reshape(-1, self.mesh.cdim) + return np.asarray(self.psi_star[0].coords_nd).reshape(-1, self.mesh.cdim) def _values_at(self, expr, X): """``expr`` at the points ``X``, one column per stored component.""" @@ -5888,28 +5894,30 @@ def _values_at(self, expr, X): for (i, j) in self._components]) def _reconstruct(self, arrivals, values, nodes): - """Fit the arrivals in each cell at the field's degree; the fit at the nodes.""" + """Fit the arrivals in each cell at the field's degree; the fit at the nodes. + A cell nothing reached keeps the field it launched.""" inside = np.asarray(self.mesh.points_in_domain(arrivals), dtype=bool) - self._fit = self._projector.fit(arrivals[inside], values[inside], old=self._fit) - return np.nan_to_num(self._projector.interpolate(self._fit, nodes)) + launched = self._values_at(self._psi_fn, np.asarray(self._fit_var.coords_nd)) + fit = self._projector.fit(arrivals[inside], values[inside], old=launched) + at_nodes = self._projector.interpolate(fit, nodes) + if np.isnan(at_nodes).any(): + raise RuntimeError("ForwardNodesSemiLagrangian: a node lies in no cell of the fit") + return at_nodes # ------------------------------------------------------------------ def initialise_history(self): - """Start from the current field at the nodes. A history already placed - by :meth:`commit_flux_to_history` is the start, and is kept.""" + """Start from the current field at the nodes.""" self.characteristics.initialise_levels(self._n_v) - if not self._history_committed: - self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) + self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) self._history_initialised = True - def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, **_ignored): - """Carry the field forward one step from the nodes and rebuild it there. - - ``store_result=False`` says the store already holds the values to launch - (placed by :meth:`commit_flux_to_history`); otherwise the tracked field - is read at the nodes first. - """ + def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): + """Carry the field forward one step from the nodes and rebuild it there.""" self._dt = dt = self._nondim_timestep(dt) + if self._projector.mesh_version != self.mesh._mesh_version: + raise NotImplementedError( + "ForwardNodesSemiLagrangian: the launch lattice and the per-cell fit were built " + "for the mesh as it was, and the mesh has moved or been re-meshed since") if not self._history_initialised: self.initialise_history() _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) @@ -5919,11 +5927,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** trace.begin_step(dt) nodes = self._nodes() launch = np.vstack([nodes, self._interior]) - # the tracked field, or (when a flux was committed) the store, which is - # a polynomial inside each element, so its interior values are exact - source = self._psi_fn if store_result else self.psi_star[0].sym - values = np.vstack([self._values_at(source, nodes) if store_result else np.array(self.psi_star[0].data), - self._values_at(source, self._interior)]) + values = self._values_at(self._psi_fn, launch) key = (_basis_key_of(self.psi_star[0]), "launch") arrivals = np.asarray(trace.departure_points(key, launch, (("first", 0, -float(dt)),), evalf=evalf, clamp_final=False)) @@ -5932,8 +5936,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** if self._inflow_value is not None: # a node whose back-trace leaves the domain holds fluid that entered this step back = 2.0 * nodes - arrivals - restored = np.asarray(self.mesh.return_coords_to_bounds(back.copy())).reshape(back.shape) - entered = np.any(restored != back, axis=1) + entered = ~np.asarray(self.mesh.points_in_domain(back), dtype=bool) if entered.any(): self._write_inflow(self.psi_star[0], nodes, entered) if self._owns_characteristics: @@ -5949,6 +5952,90 @@ def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): self._n_solves_completed += 1 def commit_flux_to_history(self, flux, verbose=False): - """Read the new flux at the nodes and leave it in the store until the next carry.""" - self.psi_star[0].data[:, :] = self._values_at(flux, self._nodes()) - self._history_committed = True + """Refused: a flux is formed at the integration points, not known at the nodes.""" + raise NotImplementedError( + "ForwardNodesSemiLagrangian carries a field known at its nodes; a flux (a stress) " + "is formed at the integration points, so carry it with " + "ForwardIntegrationPointsSemiLagrangian") + + +_SEMI_LAGRANGIAN_SCHEMES = { + ("backward", "nodes"): BackwardNodesSemiLagrangian, + ("backward", "integration_points"): BackwardIntegrationPointsSemiLagrangian, + ("forward", "integration_points"): ForwardIntegrationPointsSemiLagrangian, + ("forward", "nodes"): ForwardNodesSemiLagrangian, +} + + +def SemiLagrangian(mesh, psi_fn, V_fn, vtype=VarType.SCALAR, *, trace="backward", launch="nodes", + **kwargs): + r"""Semi-Lagrangian history of ``psi_fn`` carried by ``V_fn``. + + A semi-Lagrangian history holds the carried quantity at the points where + the weak form reads it, one step along the characteristic. The schemes + differ in two choices: + + ============ ====================== ============================================== + ``trace`` ``launch`` scheme + ============ ====================== ============================================== + backward nodes :class:`BackwardNodesSemiLagrangian` + backward integration_points :class:`BackwardIntegrationPointsSemiLagrangian` + forward integration_points :class:`ForwardIntegrationPointsSemiLagrangian` + forward nodes :class:`ForwardNodesSemiLagrangian` + ============ ====================== ============================================== + + ``trace="backward"`` follows the characteristic back from each storage + point and samples the old field at the departure point. + ``trace="forward"`` launches the old field from where it is known, carries + it one step forward, and fits the arrivals in each cell. ``launch`` names + the storage points: the field's ``nodes``, or the mesh's + ``integration_points``, where a flux such as a stress is formed. + + The remaining arguments are keywords for the scheme; see its class for + them. A keyword the chosen scheme does not take is an error. + + Parameters + ---------- + mesh : Mesh + psi_fn : sympy expression or matrix + The carried quantity, for example ``T.sym``. + V_fn : sympy matrix + The velocity that carries it. + vtype : VarType + trace : {"backward", "forward"} + launch : {"nodes", "integration_points"} + + Examples + -------- + .. code-block:: python + + DuDt = uw.systems.ddt.SemiLagrangian( + mesh, T.sym, V.sym, uw.VarType.SCALAR, degree=T.degree, + trace="forward", launch="nodes") + """ + try: + scheme = _SEMI_LAGRANGIAN_SCHEMES[(trace, launch)] + except KeyError: + raise ValueError( + f"no semi-Lagrangian scheme traces {trace!r} from {launch!r}: trace is " + "'backward' or 'forward', launch is 'nodes' or 'integration_points'") from None + taken = inspect.signature(scheme).parameters + refused = sorted(k for k in kwargs if k not in taken or taken[k].kind is taken[k].VAR_KEYWORD) + if refused: + raise TypeError(f"{scheme.__name__} takes no {refused}") + return scheme(mesh, psi_fn, V_fn, vtype, **kwargs) + + +_RENAMED = { + "IntegrationPointSemiLagrangian": "BackwardIntegrationPointsSemiLagrangian", + "ForwardSemiLagrangian": "ForwardIntegrationPointsSemiLagrangian", +} + + +def __getattr__(name): + if name in _RENAMED: + warnings.warn( + f"ddt.{name} is now ddt.{_RENAMED[name]}, or " + f"ddt.SemiLagrangian(..., trace=..., launch=...)", FutureWarning, stacklevel=2) + return globals()[_RENAMED[name]] + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/src/underworld3/systems/solver_template.py b/src/underworld3/systems/solver_template.py index 873d298b..5164e5d1 100644 --- a/src/underworld3/systems/solver_template.py +++ b/src/underworld3/systems/solver_template.py @@ -9,7 +9,7 @@ from typing import Optional, Union, Callable import underworld3 as uw from underworld3.systems import SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt -from underworld3.systems.ddt import SemiLagrangian, Lagrangian, Eulerian +from underworld3.systems.ddt import BackwardNodesSemiLagrangian, Lagrangian, Eulerian from underworld3 import timing from underworld3.systems.solvers import expression @@ -89,8 +89,8 @@ def __init__( u_Field: Optional[uw.discretisation.MeshVariable] = None, degree: int = 2, verbose: bool = False, - DuDt: Optional[Union[SemiLagrangian, Lagrangian, Eulerian]] = None, - DFDt: Optional[Union[SemiLagrangian, Lagrangian, Eulerian]] = None, + DuDt: Optional[Union[BackwardNodesSemiLagrangian, Lagrangian, Eulerian]] = None, + DFDt: Optional[Union[BackwardNodesSemiLagrangian, Lagrangian, Eulerian]] = None, ): """ Initialize the solver. diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 1c6614f3..d8b80618 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -51,6 +51,8 @@ >>> stokes.solve() """ +import warnings + import sympy from sympy import sympify import numpy as np @@ -453,12 +455,29 @@ def _invalidate_solution_cache(u): target_var._canonical_data = None -from .ddt import SemiLagrangian as SemiLagrangian_DDt +from .ddt import BackwardNodesSemiLagrangian from .ddt import Lagrangian as Lagrangian_DDt from .ddt import Lagrangian_Swarm as Lagrangian_Swarm_DDt from .ddt import Eulerian as Eulerian_DDt from .ddt import Symbolic as Symbolic_DDt +# The semi-Lagrangian schemes a solver can build its history with, named +# "_" after the arguments of ddt.SemiLagrangian. A stress is +# formed at the integration points, so it is not carried forward from nodes. +_SEMI_LAGRANGIAN_TRANSPORTS = ( + "backward_nodes", "backward_integration_points", + "forward_integration_points", "forward_nodes", +) +_STRESS_TRANSPORTS = ( + "backward_nodes", "backward_integration_points", "forward_integration_points", + "lagrangian", "eulerian", +) +_RENAMED_TRANSPORTS = { + "semi_lagrangian": "backward_nodes", + "integration_point": "backward_integration_points", + "forward": "forward_integration_points", +} + class _ConstitutiveModelStateMixin: """Single definition of the constitutive-model readiness flag. @@ -504,9 +523,9 @@ class SNES_Poisson(_ConstitutiveModelStateMixin, SNES_Scalar): Polynomial degree for the solution field (default: 2). verbose : bool, optional Enable verbose output during solve. - DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DuDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for time-dependent problems. - DFDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DFDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for the flux. Notes @@ -527,8 +546,8 @@ def __init__( u_Field: uw.discretisation.MeshVariable = None, degree=2, verbose=False, - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): if type(degree) is bool: # Legacy positional order (mesh, u_Field, verbose, degree): the @@ -1394,9 +1413,9 @@ class SNES_Stokes(_ConstitutiveModelStateMixin, SNES_Stokes_SaddlePt): If True (default), pressure is continuous. Set False for discontinuous pressure. verbose : bool, optional Enable verbose output during solving. Default is False. - DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DuDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Material derivative operator for velocity (used in derived classes). - DFDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DFDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Material derivative operator for flux (used in viscoelastic models). Notes @@ -1450,8 +1469,8 @@ def __init__( p_continuous: Optional[bool] = True, verbose: Optional[bool] = False, # Not used in Stokes, but may be used in NS, VE etc - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): super().__init__( mesh, @@ -1544,43 +1563,56 @@ def set_jacobian_F1_source(self, F1_source, linesearch="cp"): @property def stress_transport(self) -> str: - """How a viscoelastic stress history is carried: ``"semi_lagrangian"`` - (default), ``"integration_point"``, ``"forward"``, ``"lagrangian"`` or - ``"eulerian"``. - - ``"forward"`` carries the stress from a fixed set of launch points inside - the cells (the integration points), one forward trajectory a step, and - fits the arrivals per cell; the constitutive flux is read at the launch - points through a continuous P1 projection. It holds the Maxwell - start-up below Courant one where the integration-point history rings - (see :class:`~underworld3.systems.ddt.ForwardSemiLagrangian`). - - The semi-Lagrangian history traces the stress back along characteristics - and stores it on a nodal field, which the assembler then interpolates to - the integration points: two interpolations a step. ``"integration_point"`` - traces back to the integration points themselves and holds the history - there, so it carries one evaluation error and needs no projection. The - Eulerian one transports the stress on the grid with the same + """How a viscoelastic stress history is carried: ``"backward_nodes"`` + (default), ``"backward_integration_points"``, + ``"forward_integration_points"``, ``"lagrangian"`` or ``"eulerian"``. + + The first three are semi-Lagrangian schemes of + :func:`~underworld3.systems.ddt.SemiLagrangian`, named by the direction + of the trace and the points the history is held at. A backward trace + follows the characteristic back from each storage point and samples the + old stress at the departure point. ``"backward_nodes"`` stores the + history on a nodal field, which the assembler then interpolates to the + integration points: two interpolations a step. + ``"backward_integration_points"`` traces back to the integration points + themselves and holds the history there: one evaluation error and no + projection. A forward trace launches the old stress from where it is + known, carries it one step forward and fits the arrivals in each cell. + ``"forward_integration_points"`` launches from the integration points, + where the stress is formed, and reads the constitutive flux there + through a continuous P1 projection; it holds the Maxwell start-up below + Courant one where the backward integration-point history rings (see + :class:`~underworld3.systems.ddt.ForwardIntegrationPointsSemiLagrangian`). + The fourth semi-Lagrangian scheme, forward from nodes, carries a field + known at its nodes; a stress is formed at the integration points, so it + is not offered here. + + ``"eulerian"`` transports the stress on the grid with the same streamline-upwind stabilisation the Eulerian solvers use, and gives the same answer on any partition. ``"lagrangian"`` carries the stress on a swarm of material points the solver creates and advects, reading the constitutive flux at the particles each step and never projecting it back to the mesh: no numerical diffusion of the history, at the cost of - the swarm (see :class:`~underworld3.systems.ddt.Lagrangian`). The default - is ``"semi_lagrangian"``. Set - it before the constitutive model is assigned: assigning the model - creates the history, and the choice cannot change after that. + the swarm (see :class:`~underworld3.systems.ddt.Lagrangian`). + + Set it before the constitutive model is assigned: assigning the model + creates the history, and the choice cannot change after that. The + former names ``"semi_lagrangian"``, ``"integration_point"`` and + ``"forward"`` are accepted, with a warning. """ - return getattr(self, "_stress_transport", "semi_lagrangian") + return getattr(self, "_stress_transport", "backward_nodes") @stress_transport.setter def stress_transport(self, value): value = str(value) - if value not in ("semi_lagrangian", "integration_point", "forward", - "lagrangian", "eulerian"): + if value in _RENAMED_TRANSPORTS: + warnings.warn( + f"stress_transport={value!r} is now {_RENAMED_TRANSPORTS[value]!r}", + FutureWarning, stacklevel=2) + value = _RENAMED_TRANSPORTS[value] + if value not in _STRESS_TRANSPORTS: raise ValueError( - "stress_transport must be 'semi_lagrangian', 'integration_point', " - f"'forward', 'lagrangian' or 'eulerian', not {value!r}.") + f"stress_transport must be one of {_STRESS_TRANSPORTS}, not {value!r}.") if self.Unknowns.DFDt is not None: raise RuntimeError( "the stress history already exists: set stress_transport before the " @@ -1790,51 +1822,51 @@ def _create_stress_history_ddt(self, order=2): # dimensionless log-conformation of one units=uw.units.Pa if getattr(cm, "_stress_history", "stress") == "stress" else None, ) - if self.stress_transport == "integration_point": + if self.stress_transport == "backward_integration_points": unsupported = set(ddt_kwargs) - {"with_forcing_history"} if unsupported: raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " f"{sorted(unsupported)}, which the integration-point flavour " - "does not provide; use stress_transport='semi_lagrangian'.") - self.Unknowns.DFDt = uw.systems.ddt.IntegrationPointSemiLagrangian( + "does not provide; use stress_transport='backward_nodes'.") + self.Unknowns.DFDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( self.mesh, sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), self.u.sym, **ddt_kwargs, **{k: v for k, v in common.items() if k != "smoothing"}, ) - elif self.stress_transport == "forward": + elif self.stress_transport == "forward_integration_points": if ddt_kwargs: raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " f"{sorted(ddt_kwargs)}, which the forward flavour does not provide; " - "use stress_transport='semi_lagrangian' for it.") - self.Unknowns.DFDt = uw.systems.ddt.ForwardSemiLagrangian( + "use stress_transport='backward_nodes' for it.") + self.Unknowns.DFDt = uw.systems.ddt.ForwardIntegrationPointsSemiLagrangian( self.mesh, sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), self.u.sym, - vtype=common["vtype"], varsymbol=common["varsymbol"], order=order, - units=common["units"], + vtype=common["vtype"], degree=common["degree"], varsymbol=common["varsymbol"], + order=order, units=common["units"], ) elif self.stress_transport == "lagrangian": if ddt_kwargs: raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " f"{sorted(ddt_kwargs)}, which the particle Lagrangian flavour does " - "not provide; use stress_transport='semi_lagrangian' for it.") + "not provide; use stress_transport='backward_nodes' for it.") # Order 1 BDF only for now: the particle flavour has no exponential # coefficients (it is built with_exp=False), and order 2 is not yet # validated. Refuse cleanly rather than crash inside the first solve. if getattr(cm, "_integrator", "bdf") != "bdf": raise NotImplementedError( "the particle Lagrangian stress history supports the BDF " - "integrator only; use stress_transport='semi_lagrangian' for the " + "integrator only; use stress_transport='backward_nodes' for the " "exponential one.") if order > 1: raise NotImplementedError( "the particle Lagrangian stress history is first order for now; " - "use stress_transport='semi_lagrangian' for order 2.") + "use stress_transport='backward_nodes' for order 2.") # The solver owns the swarm: Lagrangian creates and populates it, and # carries the stress on it. Lagrangian_Swarm (a user-supplied swarm) # stays available by passing DFDt= to the constructor. @@ -1856,7 +1888,7 @@ def _create_stress_history_ddt(self, order=2): raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " f"{sorted(ddt_kwargs)}, which only the semi-Lagrangian flavour " - "provides; use stress_transport='semi_lagrangian' for it.") + "provides; use stress_transport='backward_nodes' for it.") self.Unknowns.DFDt = uw.systems.ddt.EulerianSUPG( self.mesh, sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), @@ -2807,8 +2839,8 @@ def __init__( order: Optional[int] = 2, p_continuous: Optional[bool] = True, verbose: Optional[bool] = False, - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): import warnings warnings.warn( @@ -2935,8 +2967,8 @@ def __init__( degree: Optional[int] = 2, p_continuous: Optional[bool] = True, verbose: Optional[bool] = False, - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): super().__init__( mesh, @@ -4382,13 +4414,13 @@ class SNES_AdvectionDiffusion(SNES_Scalar): Function to restore particles to valid domain. verbose : bool, default=False Enable verbose output. - DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DuDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for the unknown. - DFDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DFDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for the flux. monotone_mode : str or None, optional Monotonicity limiter for the semi-Lagrangian trace-back. - Forwarded to the internally-constructed ``SemiLagrangian_DDt`` + Forwarded to the internally-constructed ``BackwardNodesSemiLagrangian`` instances for ``DuDt`` and ``DFDt``. - ``None`` (default): pure FE trace-back. Can overshoot at @@ -4411,7 +4443,7 @@ class SNES_AdvectionDiffusion(SNES_Scalar): theta : float, default=0.5 Adams-Moulton theta for the diffusive flux at order 1. Forwarded to the internally-constructed - ``SemiLagrangian_DDt`` instances (same forwarding rule as + ``BackwardNodesSemiLagrangian`` instances (same forwarding rule as ``monotone_mode``). - ``0.5`` (default): Crank-Nicolson, A-stable but not @@ -4422,6 +4454,17 @@ class SNES_AdvectionDiffusion(SNES_Scalar): SLCN+CN ringing dominates the discretisation error. - ``0.0``: Forward Euler — unstable for stiff diffusion; included for completeness. + transport : str, default="backward_nodes" + The semi-Lagrangian scheme of the internally-constructed ``DuDt``: + ``"backward_nodes"``, ``"backward_integration_points"``, + ``"forward_integration_points"`` or ``"forward_nodes"``, named by the + ``trace`` and ``launch`` arguments of + :func:`~underworld3.systems.ddt.SemiLagrangian`; an option the scheme + does not take (``monotone_mode`` on a forward scheme, say) is refused. + Only the value history is chosen: the diffusive flux history ``DFDt`` + is always backward from the nodes. ``"forward_integration_points"`` + fits a linear polynomial per cell and so needs a degree-1 field; + ``"forward_nodes"`` runs in serial only. old_frame_traceback : bool, default=False Use the old-frame semi-Lagrangian reach-back for the advective ``DuDt`` history on a moving mesh (free surface or interior-node @@ -4482,11 +4525,12 @@ def __init__( order: int = 1, restore_points_func: Callable = None, verbose=False, - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, monotone_mode: Optional[str] = None, theta: float = 0.5, old_frame_traceback: bool = False, + transport: str = "backward_nodes", ): ## Parent class will set up default values etc super().__init__( @@ -4519,8 +4563,14 @@ def __init__( ## NB - Smoothing is generally required for stability. 0.0001 is effective ## at the various resolutions tested. - if DuDt is None: - self.Unknowns.DuDt = SemiLagrangian_DDt( + if transport not in _SEMI_LAGRANGIAN_TRANSPORTS: + raise ValueError(f"transport must be one of {_SEMI_LAGRANGIAN_TRANSPORTS}, " + f"not {transport!r}") + if DuDt is not None and transport != "backward_nodes": + raise ValueError("transport chooses the DuDt the solver builds; it cannot " + "apply to a DuDt that is supplied") + if DuDt is None and transport == "backward_nodes": + self.Unknowns.DuDt = BackwardNodesSemiLagrangian( self.mesh, u_Field.sym, # Symbolic expression - SemiLagrangian evaluates this at each update self._V_fn, @@ -4536,6 +4586,29 @@ def __init__( theta=theta, old_frame_traceback=old_frame_traceback, ) + elif DuDt is None: + if not u_Field.continuous: + raise NotImplementedError( + f"transport={transport!r} holds a continuous history; " + "use transport='backward_nodes' for a discontinuous field") + # options a scheme does not take are refused by ddt.SemiLagrangian, + # so only those that were asked for are passed on + asked = {k: v for k, v in (("monotone_mode", monotone_mode), + ("old_frame_traceback", old_frame_traceback)) if v} + trace, launch = transport.split("_", 1) + self.Unknowns.DuDt = uw.systems.ddt.SemiLagrangian( + self.mesh, + u_Field.sym, + self._V_fn, + uw.VarType.SCALAR, + trace=trace, + launch=launch, + degree=u_Field.degree, + varsymbol=u_Field.symbol, + order=1, + theta=theta, + **asked, + ) else: # validation @@ -4554,7 +4627,7 @@ def __init__( # flux vector (volume meshes have dim==cdim so this is # unchanged; manifold meshes have cdim > dim and need the # extra component). - self.Unknowns.DFDt = SemiLagrangian_DDt( + self.Unknowns.DFDt = BackwardNodesSemiLagrangian( self.mesh, sympy.Matrix([[0] * self.mesh.cdim]), # Actual function is not defined at this point self._V_fn, @@ -4585,7 +4658,7 @@ def __init__( # ALE trace-back path (and REMAP it correctly on an OT opt-out # reset). Only meaningful when DuDt traces back (SemiLagrangian); # Eulerian/Lagrangian fields keep the default policy. - if isinstance(self.Unknowns.DuDt, SemiLagrangian_DDt): + if isinstance(self.Unknowns.DuDt, BackwardNodesSemiLagrangian): from underworld3.discretisation.remesh import RemeshPolicy self.u.remesh_policy = RemeshPolicy.CARRY self.u._remesh_managed_by = self.Unknowns.DuDt @@ -5047,9 +5120,9 @@ class SNES_Diffusion(SNES_Scalar): Numerically evaluate symbolic expressions during setup. verbose : bool, default=False Enable verbose output. - DuDt : Eulerian_DDt, SemiLagrangian_DDt, or Lagrangian_DDt, optional + DuDt : Eulerian_DDt, BackwardNodesSemiLagrangian, or Lagrangian_DDt, optional Time derivative operator for the unknown. - DFDt : Eulerian_DDt, SemiLagrangian_DDt, or Lagrangian_DDt, optional + DFDt : Eulerian_DDt, BackwardNodesSemiLagrangian, or Lagrangian_DDt, optional Time derivative operator for the flux. Notes @@ -5081,8 +5154,8 @@ def __init__( theta: float = 0.0, evalf: Optional[bool] = False, verbose=False, - DuDt: Union[Eulerian_DDt, SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[Eulerian_DDt, SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[Eulerian_DDt, BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[Eulerian_DDt, BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): ## Parent class will set up default values etc super().__init__( @@ -5373,9 +5446,9 @@ class SNES_NavierStokes(SNES_Stokes_SaddlePt): If False, use discontinuous pressure elements. verbose : bool, default=False Enable verbose output. - DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DuDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for velocity. - DFDt : SemiLagrangian_DDt or Lagrangian_DDt, optional + DFDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for stress. Notes @@ -5421,8 +5494,8 @@ def __init__( flux_order: Optional[int] = None, p_continuous: Optional[bool] = False, verbose: Optional[bool] = False, - DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, - DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, + DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, ): ## Parent class will set up default values and load u_Field into the solver super().__init__( @@ -5637,7 +5710,7 @@ def DuDt(self): @DuDt.setter def DuDt( self, - DuDt_value: Union[SemiLagrangian_DDt, Lagrangian_DDt], + DuDt_value: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt], ): """Set the time derivative operator for velocity.""" self.Unknowns.DuDt = DuDt_value diff --git a/tests/parallel/test_1062_forward_stress_history_mpi.py b/tests/parallel/test_1062_forward_stress_history_mpi.py index db421412..8327099b 100644 --- a/tests/parallel/test_1062_forward_stress_history_mpi.py +++ b/tests/parallel/test_1062_forward_stress_history_mpi.py @@ -21,7 +21,7 @@ POINTS = np.array([[0.5, 0.2], [-0.3, -0.35]]) -def turned_over_maxwell_box(transport="forward", steps=10, dt=0.1): +def turned_over_maxwell_box(transport="forward_integration_points", steps=10, dt=0.1): Lx, H = 1.0, 0.5 mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-Lx, -H), maxCoords=(Lx, H), cellSize=0.125, qdegree=3, regular=False) @@ -48,7 +48,7 @@ def turned_over_maxwell_box(transport="forward", steps=10, dt=0.1): def test_the_forward_history_gives_the_serial_stress_on_every_rank(): kind, values, relocated = turned_over_maxwell_box() - assert kind == "ForwardSemiLagrangian" + assert kind == "ForwardIntegrationPointsSemiLagrangian" # The serial values are the hard baseline (they pin the physics; regenerate # them if a default changes). Parallel matches them to 1e-4, not to # round-off: a cell's arrivals are summed into its least-squares fit in an diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 1aefe935..d2b7df65 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -31,7 +31,7 @@ def test_slots_are_exact_departure_point_values(): V = sympy.Matrix([[v[0], v[1]]]) dt = 0.1 - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=2) + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V, degree=2, order=2) assert all(ps.is_integration_point for ps in ddt.psi_star) ddt.update_pre_solve(dt) @@ -70,7 +70,7 @@ def _rotating_gaussian(mesh, kind, dt, nsteps): T = uw.discretisation.MeshVariable(f"T_{kind}", mesh, 1, degree=2) T.data[:, 0] = gauss(np.asarray(T.coords), x0, 0.0) if kind == "ip": - DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=1) + DuDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V, degree=2, order=1) adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, DuDt=DuDt, order=1) else: adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, order=1) @@ -97,10 +97,10 @@ def test_undersampled_rule_is_refused(): T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) V = sympy.Matrix([[1.0, 0.0]]) with pytest.raises(RuntimeError, match="oversampled"): - uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2) + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V, degree=2) # P1 on the same rule is 2x oversampled and accepted. T1 = uw.discretisation.MeshVariable("T1", mesh, 1, degree=1) - uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T1, V, degree=1) + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T1, V, degree=1) @pytest.mark.level_2 @@ -136,7 +136,7 @@ def _unsteady_uniform_flow_check(kind, vform="var"): "ramp": (c * v_var.sym, 1.0)}[vform] if kind == "ip": - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V_fn, degree=2, order=1) + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V_fn, degree=2, order=1) else: ddt = uw.systems.ddt.SemiLagrangian(mesh, T, V_fn, uw.VarType.SCALAR, degree=2, continuous=True, order=1) @@ -200,7 +200,7 @@ def test_composed_advdiffusion_reachability(config): def run(solver_cls, kwargs): T = uw.discretisation.MeshVariable(f"T_{config}_{solver_cls.__name__}", mesh, 1, degree=2) T.data[:, 0] = gauss(np.asarray(T.coords)) - D = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=order, theta=theta) + D = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V, degree=2, order=order, theta=theta) adv = solver_cls(mesh, u_Field=T, V_fn=V, DuDt=D, order=order, **kwargs) adv.constitutive_model = uw.constitutive_models.DiffusionModel adv.constitutive_model.Parameters.diffusivity = 1e-9 @@ -357,7 +357,7 @@ def test_a_vector_history_holds_the_departure_point_values(): with uw.synchronised_array_update(): U.data[...] = _vector_field(np.asarray(U.coords)) - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, U, _velocity(), vtype=uw.VarType.VECTOR, degree=2, order=2) assert ddt.num_components == 2 assert all(ps.is_integration_point for ps in ddt.psi_star) @@ -390,7 +390,7 @@ def test_a_symmetric_tensor_history_transports_every_component(): with uw.synchronised_array_update(): S.data[...] = _pack(_tensor_entries(np.asarray(S.coords)), columns) - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, S, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) assert ddt.num_components == 3 assert ddt._components == columns @@ -425,7 +425,7 @@ def test_a_scalar_history_is_unchanged(): with uw.synchronised_array_update(): T.data[:, 0] = _scalar_field(np.asarray(T.coords)) - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, T, _velocity(), degree=2, order=1) assert ddt.num_components == 1 assert ddt.bdf().shape == (1, 1) @@ -449,7 +449,7 @@ def test_the_history_symbol_participates_in_expressions(vtype): columns = _storage_components(vtype, tuple(var.sym.shape)) with uw.synchronised_array_update(): var.data[...] = 2.0 - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, var, _velocity(), vtype=vtype, degree=2, order=1) ddt.update_pre_solve(DT) @@ -470,7 +470,7 @@ def test_the_refusal_is_gone_but_the_rule_check_is_not(): mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) U = uw.discretisation.MeshVariable("Ur", mesh, vtype=uw.VarType.VECTOR, degree=2) with pytest.raises(RuntimeError, match="qdegree|rule|oversample"): - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, U, _velocity(), vtype=uw.VarType.VECTOR, degree=2, order=1) @@ -492,12 +492,12 @@ def test_a_vtype_that_does_not_match_psi_fn_is_refused(): before = len(mesh.vars) with pytest.raises(ValueError, match="psi_fn has shape"): - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, sympy.Matrix([[1.0]]), _velocity(), vtype=uw.VarType.VECTOR, degree=2, order=1) with pytest.raises(ValueError, match="psi_fn has shape"): - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, sympy.Matrix([[1.0, 2.0]]), _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) @@ -513,7 +513,7 @@ def test_the_shape_guard_is_on_the_setter_not_only_the_constructor(): mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=3) S = uw.discretisation.MeshVariable("Sg", mesh, vtype=uw.VarType.SYM_TENSOR, degree=2) - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian( + ddt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, S, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) with pytest.raises(ValueError, match="psi_fn has shape"): @@ -534,7 +534,7 @@ def test_a_full_tensor_is_not_accepted_as_a_symmetric_one(): full = uw.discretisation.MeshVariable("Tg", mesh, vtype=uw.VarType.TENSOR, degree=2) with pytest.raises(ValueError, match="stores 4 components"): - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, full, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) @@ -550,7 +550,7 @@ def test_an_asymmetric_psi_fn_under_sym_tensor_says_so(): asymmetric = sympy.Matrix([[1 + x, 2 + y], [100.0, 3 + x * y]]) with pytest.warns(UserWarning, match="not symmetric"): - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, asymmetric, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) @@ -558,6 +558,6 @@ def test_an_asymmetric_psi_fn_under_sym_tensor_says_so(): symmetric = sympy.Matrix([[1 + x, 2 + y], [2 + y, 3 + x * y]]) with warnings.catch_warnings(): warnings.simplefilter("error", UserWarning) - uw.systems.ddt.IntegrationPointSemiLagrangian( + uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, symmetric, _velocity(), vtype=uw.VarType.SYM_TENSOR, degree=2, order=1) diff --git a/tests/test_1056_units_slcn_traceback.py b/tests/test_1056_units_slcn_traceback.py index c8ae585d..4c0d827d 100644 --- a/tests/test_1056_units_slcn_traceback.py +++ b/tests/test_1056_units_slcn_traceback.py @@ -53,7 +53,7 @@ def _advect_blob(use_units, vy=20.0, nsteps=5, dt=2.0, scheme="slcn"): T.data[:, 0] = np.exp(-(((c[:, 0] - 500) / 120) ** 2 + ((c[:, 1] - 300) / 120) ** 2)) if scheme == "slcn_ip": - DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V.sym, degree=2, order=1) + DuDt = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian(mesh, T, V.sym, degree=2, order=1) adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym, DuDt=DuDt, order=1) else: adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym) diff --git a/tests/test_1059_stress_transport.py b/tests/test_1059_stress_transport.py index f9d14cc5..099d1409 100644 --- a/tests/test_1059_stress_transport.py +++ b/tests/test_1059_stress_transport.py @@ -223,9 +223,9 @@ def _maxwell_shear(transport, order, steps=20, dt=0.1, integrator="bdf", solver= KINDS = { - "semi_lagrangian": "SemiLagrangian", - "integration_point": "IntegrationPointSemiLagrangian", - "forward": "ForwardSemiLagrangian", + "backward_nodes": "BackwardNodesSemiLagrangian", + "backward_integration_points": "BackwardIntegrationPointsSemiLagrangian", + "forward_integration_points": "ForwardIntegrationPointsSemiLagrangian", "lagrangian": "Lagrangian", "eulerian": "EulerianSUPG", } @@ -247,8 +247,8 @@ def test_every_stress_history_solves_the_maxwell_shear_box(order, integrator, to for transport, expected_kind in KINDS.items(): if integrator == "etd" and order == 2 and transport == "eulerian": continue # the grid flavour has no forcing-history slot yet - if order == 2 and transport == "forward": - continue # the forward flavour carries one level + if order == 2 and transport.startswith("forward_"): + continue # the forward flavours carry one level if transport == "lagrangian" and (order == 2 or integrator == "etd"): continue # order 1 BDF only for now (no exponential coefficients # on the particle flavour; order 2 deferred) @@ -278,7 +278,46 @@ def test_stress_transport_is_validated_and_fixed_once_the_history_exists(): stokes.constitutive_model.Parameters.dt_elastic = 0.1 assert type(stokes.DFDt).__name__ == "EulerianSUPG" with pytest.raises(RuntimeError, match="already exists"): - stokes.stress_transport = "semi_lagrangian" + stokes.stress_transport = "backward_nodes" + + +def test_the_former_stress_transport_names_still_select_their_scheme(): + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("U_old", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_old", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + for old, new in (("semi_lagrangian", "backward_nodes"), + ("integration_point", "backward_integration_points"), + ("forward", "forward_integration_points")): + with pytest.warns(FutureWarning, match=new): + stokes.stress_transport = old + assert stokes.stress_transport == new + # a stress is formed at the integration points: it is not carried from the nodes + with pytest.raises(ValueError, match="stress_transport must be"): + stokes.stress_transport = "forward_nodes" + + +def test_semi_lagrangian_selects_its_scheme_by_trace_and_launch(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25) + T = uw.discretisation.MeshVariable("T_sel", mesh, 1, degree=1) + V = sympy.Matrix([[1.0, 0.0]]) + for (trace, launch), kind in ( + (("backward", "nodes"), "BackwardNodesSemiLagrangian"), + (("backward", "integration_points"), "BackwardIntegrationPointsSemiLagrangian"), + (("forward", "integration_points"), "ForwardIntegrationPointsSemiLagrangian"), + (("forward", "nodes"), "ForwardNodesSemiLagrangian")): + history = uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, uw.VarType.SCALAR, + trace=trace, launch=launch, degree=1) + assert type(history).__name__ == kind + with pytest.raises(ValueError, match="no semi-Lagrangian scheme"): + uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, trace="sideways") + # an option the chosen scheme does not take is refused, not dropped + with pytest.raises(TypeError, match="monotone_mode"): + uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, trace="forward", degree=1, monotone_mode="clamp") + with pytest.raises(NotImplementedError, match="integration points"): + uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, trace="forward", degree=1).commit_flux_to_history(T.sym) + with pytest.warns(FutureWarning, match="ForwardIntegrationPointsSemiLagrangian"): + assert uw.systems.ddt.ForwardSemiLagrangian is uw.systems.ddt.ForwardIntegrationPointsSemiLagrangian def _sheared_varying_modulus(transport, order, steps=10, dt=0.1, res=6): @@ -324,7 +363,7 @@ def _sheared_varying_modulus(transport, order, steps=10, dt=0.1, res=6): def test_the_two_stress_histories_agree_when_the_stress_moves_and_evolves(): """With a stress that is carried as well as relaxed the schemes must agree to within their own time-discretisation error, and more tightly at order 2.""" - traced, traced_slope = _sheared_varying_modulus("semi_lagrangian", 2) + traced, traced_slope = _sheared_varying_modulus("backward_nodes", 2) grid, grid_slope = _sheared_varying_modulus("eulerian", 2) assert traced_slope > 0.1 and grid_slope > 0.1, "the stress must not be uniform" # BASELINES (the L2 norm of the shear stress; see the ledger): each scheme @@ -368,9 +407,9 @@ def run(transport): return type(ns.DFDt).__name__, float(np.asarray(uw.function.evaluate( ns.DFDt.psi_star[0].sym[0, 1], np.array([[0.0, 0.0]]))).reshape(-1)[0]) - traced_kind, traced = run("semi_lagrangian") + traced_kind, traced = run("backward_nodes") grid_kind, grid = run("eulerian") - assert traced_kind == "SemiLagrangian" and grid_kind == "EulerianSUPG" + assert traced_kind == "BackwardNodesSemiLagrangian" and grid_kind == "EulerianSUPG" # BASELINE: the shear stress at the origin after ten steps (see the ledger) assert abs(traced - NS_ORDER2_XY) < 1.0e-3 * NS_ORDER2_XY, traced assert abs(grid - traced) / traced < 1.0e-3, (grid, traced) @@ -429,7 +468,7 @@ def run(transport): return float(np.asarray(uw.function.evaluate( ns.DFDt.psi_star[0].sym[0, 1], np.array([[0.0, 0.0]]))).reshape(-1)[0]) - traced, grid = run("semi_lagrangian"), run("eulerian") + traced, grid = run("backward_nodes"), run("eulerian") # BASELINE: the shear stress at the origin after ten steps (see the ledger) assert abs(traced - CN_ORDER1_XY) < 1.0e-3 * CN_ORDER1_XY, traced assert abs(grid - traced) / traced < 1.0e-3, (grid, traced) @@ -442,7 +481,7 @@ def test_devss_is_off_by_default_and_vanishes_on_a_uniform_strain_rate(monkeypat off unless asked for. The varying-modulus box is the negative control: there D differs from edot by projection error, the term is live, and the answer must move -- by projection error, which is small, but not by nothing.""" - _kind, off, exact = _maxwell_shear("integration_point", 1) + _kind, off, exact = _maxwell_shear("backward_integration_points", 1) assert abs(off - exact) / exact < 0.02 def with_devss(builder, *args, **kw): @@ -456,14 +495,14 @@ def patched(self, *a, **k): m.setattr(uw.systems.Stokes, "__init__", patched) return builder(*args, **kw) - _kind, on, _ = with_devss(_maxwell_shear, "integration_point", 1) + _kind, on, _ = with_devss(_maxwell_shear, "backward_integration_points", 1) assert abs(on - off) < 1e-8 * abs(exact), (on, off) # the pair cancelled # the varying-modulus helper differentiates the history in a weak form, # which an integration-point variable refuses; the term is on the solver # and flavour-independent, so the nodal history serves for this half - norm_off, _ = _sheared_varying_modulus("semi_lagrangian", 1) - norm_on, _ = with_devss(_sheared_varying_modulus, "semi_lagrangian", 1) + norm_off, _ = _sheared_varying_modulus("backward_nodes", 1) + norm_on, _ = with_devss(_sheared_varying_modulus, "backward_nodes", 1) moved = abs(norm_on - norm_off) / norm_off assert 1e-6 < moved < 5e-2, moved # live, and only projection-sized @@ -512,7 +551,7 @@ def test_the_exponential_integrator_runs_on_the_trace_back_navier_stokes(): as the Stokes family does; it did neither, so the memory term was absent and the exponential integrator ran in its viscous limit (#741).""" for integrator, tolerance in (("etd", 1e-3), ("bdf", 0.02)): - _, stress, exact = _maxwell_shear("semi_lagrangian", 1, integrator=integrator, solver="ns_slcn") + _, stress, exact = _maxwell_shear("backward_nodes", 1, integrator=integrator, solver="ns_slcn") assert abs(stress - exact) / exact < tolerance, (integrator, stress, exact) @@ -523,8 +562,8 @@ def test_a_preset_velocity_gives_both_integrators_the_same_first_stress(): exponential one read alpha = phi = 0 and recorded the full viscous stress, seven times the BDF value on the cylinder (#740). One step: the two first-order integrators agree to O(dt / t_r).""" - _, bdf, _ = _maxwell_shear("semi_lagrangian", 1, steps=1, integrator="bdf", initial_velocity=True) - _, etd, _ = _maxwell_shear("semi_lagrangian", 1, steps=1, integrator="etd", initial_velocity=True) + _, bdf, _ = _maxwell_shear("backward_nodes", 1, steps=1, integrator="bdf", initial_velocity=True) + _, etd, _ = _maxwell_shear("backward_nodes", 1, steps=1, integrator="etd", initial_velocity=True) # gammadot = 1, eta = lambda = 1, dt = 0.1. The history's first level is the # constitutive flux of the preset velocity (#740), so one step gives # BDF-1: eta_eff gdot (1 + eta_eff / (mu dt)) = (1/11)(1 + 10/11) = 21/121 @@ -548,7 +587,7 @@ def test_the_integration_point_history_takes_the_inflow_value_at_an_inlet(): incoming = sympy.Matrix([[1.0, 0.25], [0.25, -1.0]]) results = {} for tag, inflow in (("with", incoming), ("without", None)): - manager = uw.systems.ddt.IntegrationPointSemiLagrangian( + manager = uw.systems.ddt.BackwardIntegrationPointsSemiLagrangian( mesh, sympy.Matrix.zeros(2, 2), velocity, vtype=uw.VarType.SYM_TENSOR, degree=1, continuous=True, order=1, varsymbol=rf"S^{{{tag}}}") assert manager.applies_inflow_value @@ -573,7 +612,7 @@ def test_the_integration_point_history_takes_the_inflow_value_at_an_inlet(): -@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points"]) def test_the_upper_convected_element_builds_the_first_normal_stress_in_shear(transport): """Start-up of simple shear for the UCM fluid has the closed form sigma_xy = eta gdot (1 - e^{-t/lambda}) and @@ -626,7 +665,7 @@ def test_a_solvent_viscosity_adds_its_newtonian_stress(): shear cannot tell (any uniform stress satisfies momentum): the channel flow's speed is set by the total viscosity.""" steps, dt, eta_s = 20, 0.1, 0.5 - _, polymer_xy, polymer_exact, _, solvent_xy = _maxwell_shear("semi_lagrangian", 1, steps=steps, dt=dt, solvent=eta_s) + _, polymer_xy, polymer_exact, _, solvent_xy = _maxwell_shear("backward_nodes", 1, steps=steps, dt=dt, solvent=eta_s) gdot = 1.0 assert abs(polymer_xy - polymer_exact) / polymer_exact < 0.02, (polymer_xy, polymer_exact) assert abs(solvent_xy - eta_s * gdot) < 1e-6, solvent_xy @@ -671,7 +710,7 @@ def _one_shear_step(transport, dt=1.0, objective_rate="none"): return stokes -@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points"]) def test_the_conformation_after_one_shear_step_is_one_minus_half_the_step(transport): stokes = _one_shear_step(transport, dt=1.0) health = stokes.constitutive_model.conformation_min_eigenvalue() @@ -693,7 +732,7 @@ def test_the_conformation_check_sees_a_lost_conformation(): v = uw.discretisation.MeshVariable("U_lost", mesh, mesh.dim, degree=2) p = uw.discretisation.MeshVariable("P_lost", mesh, 1, degree=1) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, verbose=False) - stokes.stress_transport = "integration_point" + stokes.stress_transport = "backward_integration_points" stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( stokes.Unknowns, order=1, integrator="bdf") stokes.constitutive_model.Parameters.shear_viscosity_0 = eta @@ -712,7 +751,7 @@ def test_the_conformation_check_sees_a_lost_conformation(): assert health["where"] is not None -@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points"]) def test_the_elastic_timestep_is_the_safety_factor_over_the_shear_rate(transport): stokes = _one_shear_step(transport, dt=1.0, objective_rate="upper_convected") # gammadot = 2 speed / height = 1 everywhere: dt_max = safety / 1. The rate is @@ -724,7 +763,7 @@ def test_the_elastic_timestep_is_the_safety_factor_over_the_shear_rate(transport def test_the_store_smoothing_is_the_coefficient_times_the_local_cell_size_squared(): - stokes = _one_shear_step("integration_point", dt=1.0) + stokes = _one_shear_step("backward_integration_points", dt=1.0) history = stokes.DFDt assert history.store_smoothing == 0.0 assert history._commit_projection.smoothing == 0.0 diff --git a/tests/test_1060_stress_store_smoothing.py b/tests/test_1060_stress_store_smoothing.py index 90b93db0..f5bfc334 100644 --- a/tests/test_1060_stress_store_smoothing.py +++ b/tests/test_1060_stress_store_smoothing.py @@ -30,7 +30,7 @@ def _cell_scale_content(history, mesh): return rms(raw - fit) / max(rms(raw), 1.0e-300) -def waters_king_start_up(store_smoothing, res=16, dt=0.0125, t_end=2.0, transport="integration_point", +def waters_king_start_up(store_smoothing, res=16, dt=0.0125, t_end=2.0, transport="backward_integration_points", return_kind=False): """Waters and King start-up on the integration-point history, pure Maxwell, below Courant one. Returns u at the centre at t 1, the cell-scale content @@ -50,15 +50,15 @@ def waters_king_start_up(store_smoothing, res=16, dt=0.0125, t_end=2.0, transpor ns.add_dirichlet_bc((0.0, 0.0), "Top"); ns.add_dirichlet_bc((0.0, 0.0), "Bottom") ns.add_dirichlet_bc((sympy.oo, 0.0), "Left"); ns.add_dirichlet_bc((sympy.oo, 0.0), "Right") ns.bodyforce = sympy.Matrix([[G, 0.0]]); ns.tolerance = 1e-6 - if transport == "integration_point": + if transport == "backward_integration_points": ns.DFDt.store_smoothing = store_smoothing - elif transport == "forward": + elif transport == "forward_integration_points": ns.DFDt.flux_smoothing = store_smoothing * mesh.cell_size() ** 2 # The content has to be read after the trace-back and before the solve: after # the store the point values are a P1 field sampled at the points and the # cell-scale part is zero by construction, whatever the run is doing. latest = {"content": float("nan")} - if transport == "integration_point": + if transport == "backward_integration_points": carry = ns.DFDt.update_pre_solve def carry_and_measure(*args, **kwargs): out = carry(*args, **kwargs) @@ -95,7 +95,7 @@ def test_the_store_smoothing_holds_the_cell_scale_mode_of_the_integration_point_ """ u_plain, plain, _, kind = waters_king_start_up(0.0, return_kind=True) u_smooth, smooth, _ = waters_king_start_up(0.07) - assert kind == "IntegrationPointSemiLagrangian" + assert kind == "BackwardIntegrationPointsSemiLagrangian" growth = plain[2.0] / plain[1.0] assert 4.0 < growth < 20.0, growth # e^{gamma}, gamma between 1.4 and 3 per unit time assert smooth[2.0] < smooth[1.0], smooth # held: decaying, not growing diff --git a/tests/test_1061_stress_forward_history.py b/tests/test_1061_stress_forward_history.py index 195cacc2..e5af9c33 100644 --- a/tests/test_1061_stress_forward_history.py +++ b/tests/test_1061_stress_forward_history.py @@ -16,8 +16,8 @@ def test_the_forward_history_with_its_read_back_smoothing_holds_the_maxwell_start_up(): - u1, _, u65, kind = waters_king_start_up(0.023, t_end=6.5, transport="forward", return_kind=True) - assert kind == "ForwardSemiLagrangian" + u1, _, u65, kind = waters_king_start_up(0.023, t_end=6.5, transport="forward_integration_points", return_kind=True) + assert kind == "ForwardIntegrationPointsSemiLagrangian" # measured with this dose at 1/16, dt 0.0125: 0.95427 at t 1.00, 0.51851 at t 6.5 # (nodal 0.96215 and 0.51710) assert abs(u1 - 0.9543) < 0.003 diff --git a/tests/test_1063_stress_history_restart.py b/tests/test_1063_stress_history_restart.py index 24db68a2..2c634071 100644 --- a/tests/test_1063_stress_history_restart.py +++ b/tests/test_1063_stress_history_restart.py @@ -42,7 +42,7 @@ def _stress_at_origin(stokes): return float(np.asarray(uw.function.evaluate(stokes.DFDt.psi_star[0].sym[0, 1], np.array([[0.0, 0.0]]))).reshape(-1)[0]) -@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward", "lagrangian"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points", "lagrangian"]) def test_a_restored_history_continues_where_it_left_off(transport): orchestration_model, stokes, dt = _shear_box(transport) for _ in range(6): diff --git a/tests/test_1064_stress_history_units.py b/tests/test_1064_stress_history_units.py index 6d72401f..e825b939 100644 --- a/tests/test_1064_stress_history_units.py +++ b/tests/test_1064_stress_history_units.py @@ -24,8 +24,8 @@ STRESS_SCALE_PA = 1.0e21 / MYR_S STEPS = 12 ORIGIN = np.array([[0.0, 0.0]]) -CASES = [(t, 1) for t in ("semi_lagrangian", "integration_point", "forward", "lagrangian", "eulerian")] \ - + [(t, 2) for t in ("semi_lagrangian", "integration_point", "eulerian")] +CASES = [(t, 1) for t in ("backward_nodes", "backward_integration_points", "forward_integration_points", "lagrangian", "eulerian")] \ + + [(t, 2) for t in ("backward_nodes", "backward_integration_points", "eulerian")] def _shear_box(transport, order, with_units, **model_options): @@ -35,7 +35,7 @@ def _shear_box(transport, order, with_units, **model_options): length=uw.quantity(1.0, "km"), viscosity=uw.quantity(1.0e21, "Pa*s"), time=uw.quantity(1.0, "Myr")) mesh = uw.meshing.StructuredQuadBox(elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)) - tag = f"{transport[:3]}{order}{'u' if with_units else 'n'}" + tag = f"{"".join(w[0] for w in transport.split("_"))}{order}{'u' if with_units else 'n'}" v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2, units="km/Myr" if with_units else None) p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1, units="Pa" if with_units else None) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) @@ -86,7 +86,7 @@ def test_units_model_history_is_the_nondimensional_history(transport, order): assert abs(value_pa - expected_pa) < 1.0e-6 * abs(expected_pa), (transport, order, value_pa, expected_pa) -@pytest.mark.parametrize("transport", ["semi_lagrangian", "forward"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "forward_integration_points"]) def test_units_model_log_conformation_history(transport): """The log-conformation record is dimensionless: the same in both runs; the stress the model reads from it, G (exp(psi) - I), is in Pa.""" diff --git a/tests/test_1065_log_conformation.py b/tests/test_1065_log_conformation.py index ba14f502..47b5b869 100644 --- a/tests/test_1065_log_conformation.py +++ b/tests/test_1065_log_conformation.py @@ -59,7 +59,7 @@ def _extension_step(transport, convected_step, stress_history, integrator="bdf") uw.reset_default_model() mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), minCoords=(-1, -1), maxCoords=(1, 1)) x, y = mesh.X - tag = f"{transport[:3]}{convected_step[:3]}{stress_history[:3]}{integrator}" + tag = f"{"".join(w[0] for w in transport.split("_"))}{convected_step[:3]}{stress_history[:3]}{integrator}" v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2) p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) @@ -89,7 +89,7 @@ def _extension_step(transport, convected_step, stress_history, integrator="bdf") @pytest.mark.parametrize("integrator", ["bdf", "etd"]) @pytest.mark.parametrize("transport, stress_history", [ - ("semi_lagrangian", "stress"), ("semi_lagrangian", "log_conformation"), + ("backward_nodes", "stress"), ("backward_nodes", "log_conformation"), ("eulerian", "stress"), ("eulerian", "log_conformation")]) def test_deformation_step_is_the_closed_form(transport, stress_history, integrator): c = _deformation_step(np.eye(2), np.diag([RATE, -RATE]), DT, 1.0, integrator) @@ -100,14 +100,14 @@ def test_deformation_step_is_the_closed_form(transport, stress_history, integrat def test_linear_step_loses_the_conformation(): """The defect the deformation step removes, measured against its closed form.""" - c_min, _ = _extension_step("semi_lagrangian", "linear", "stress") + c_min, _ = _extension_step("backward_nodes", "linear", "stress") assert abs(c_min - (1 - 2 * DT * RATE + DT) / (1 + DT)) < 1.0e-6, c_min def _shear_startup(transport, integrator, steps=10, dt=0.1): uw.reset_default_model() mesh = uw.meshing.StructuredQuadBox(elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)) - tag = f"s{transport[:3]}{integrator}" + tag = f"s{"".join(w[0] for w in transport.split("_"))}{integrator}" v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2) p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) @@ -139,7 +139,7 @@ def _shear_startup(transport, integrator, steps=10, dt=0.1): @pytest.mark.parametrize("integrator", ["bdf", "etd"]) -@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward", "eulerian"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points", "eulerian"]) def test_log_conformation_history_follows_the_recurrence(transport, integrator): """Uniform stress, so transport is a no-op: this checks the encoding, the decoding and the step together, against the discrete recurrence.""" diff --git a/tests/test_1102_forward_nodes_rotating_gaussian.py b/tests/test_1102_forward_nodes_rotating_gaussian.py index fe232067..9a990ba9 100644 --- a/tests/test_1102_forward_nodes_rotating_gaussian.py +++ b/tests/test_1102_forward_nodes_rotating_gaussian.py @@ -3,7 +3,7 @@ Launched from the field's own nodes and from a lattice inside every element (the field's interpolant is a polynomial there, so both are known exactly), carried one step forward, fitted per cell at the field's degree and read back at -the nodes. Used as the DuDt of the SLCN advection-diffusion solver, half a +the nodes. Selected by transport="forward_nodes" on the SLCN advection-diffusion solver, half a revolution, kappa 0.01, dt 0.02, P2 temperature, against the test_1101 fixture: the backward nodal history gives 2.52e-2 on the disc and 6.43e-2 on the box (where the flow crosses all four walls). @@ -25,8 +25,8 @@ def _run(mesh, walls): T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) V = sympy.Matrix([[-y, x]]) - duDt = uw.systems.ddt.ForwardNodesSemiLagrangian(mesh, T.sym, V, vtype=uw.VarType.SCALAR, degree=T.degree) - adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, DuDt=duDt, order=1) + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, order=1, transport="forward_nodes") + assert type(adv.DuDt).__name__ == "ForwardNodesSemiLagrangian" adv.constitutive_model = uw.constitutive_models.DiffusionModel adv.constitutive_model.Parameters.diffusivity = KAPPA for wall in walls: From b20c992e3cee8ed8de8147ca8450ed952ff9aed7 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 26 Sep 2026 23:28:04 -0700 Subject: [PATCH 3/5] Every stress and advection history gives the serial answer in parallel; forward from nodes carries a stress Forward from nodes as a stress history: the stress is committed by L2 projection onto the continuous store and launched from that field. It is offered by stress_transport="forward_nodes". The per-cell fits are projected onto the store (a shared node weighs every cell), each cell's fit uses only its own arrivals (a linear fit, their mean, or the launched field when too few), the interior lattice is two degrees up, and an arrival on a shared face is fitted in every cell that contains it. In parallel a seam node is launched by its owner, and arrivals that left a rank or sit on a face are offered to every rank. Partition dependence removed from the existing histories: - a field that is a mesh variable (T.sym as well as T) is recorded into the nodal and integration-point histories by copying its nodal data, in serial as in parallel, instead of evaluating it at nudged nodes; - the nodal trace starts from the node itself: the 0.1% nudge toward the closest cell's centroid depended on which cells a rank holds and gave a no-slip wall node a velocity; - every history projection is solved to 1e-10; - global_evaluate: a point that the migration stranded but that a rank's cell contains is evaluated there with the FE interpolant, not by the nearest-centroid rbf extrapolation (3e-2 at a steep stress). tests/parallel/test_1066: six stress and four advection histories equal serial to 1e-6 at np >= 3 (verified np 3, 4, 6); the particle history is a strict xfail at np >= 4 (1.2e-5, cause open). test_1062 tightened to 1e-6. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- docs/developer/subsystems/stress-transport.md | 29 +- src/underworld3/function/_function.pyx | 33 +- src/underworld3/swarm.py | 4 + src/underworld3/systems/ddt.py | 496 ++++++++++-------- src/underworld3/systems/solvers.py | 36 +- .../utilities/cell_polynomial_projection.py | 80 ++- .../test_1062_forward_stress_history_mpi.py | 10 +- .../test_1066_transport_schemes_mpi.py | 78 +++ tests/test_1059_stress_transport.py | 10 +- ...st_1102_forward_nodes_rotating_gaussian.py | 4 +- 10 files changed, 504 insertions(+), 276 deletions(-) create mode 100644 tests/parallel/test_1066_transport_schemes_mpi.py diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index f5652da8..543b8a96 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -17,15 +17,13 @@ stokes.constitutive_model.Parameters.solvent_viscosity = eta_s # Oldroyd-B; o stokes.constitutive_model.Parameters.dt_elastic = dt ``` -## The five histories +## The six histories The four semi-Lagrangian names are `_`, the arguments of `uw.systems.ddt.SemiLagrangian(..., trace=, launch=)`: a backward trace follows the characteristic back from each storage point and samples the old stress at the foot; a forward trace carries the old stress from where it is known and fits the arrivals in -each cell. The fourth combination, forward from nodes, carries a field known at its -nodes (a temperature, say); a stress is formed at the integration points, so it is not -a stress history. The former names `semi_lagrangian`, `integration_point` and `forward` +each cell. The former names `semi_lagrangian`, `integration_point` and `forward` are accepted with a warning. | `stress_transport` | storage | carried by | stable at | fails by | @@ -33,6 +31,7 @@ are accepted with a warning. | `backward_nodes` (the default) | continuous P1 at the vertices | vertex trace-back, interpolation at the foot | any Courant number | excess stress in the first cells off a no-slip wall; on the confined cylinder that excess loses the conformation and the solve hangs | | `backward_integration_points` | continuous P1 store, sampled at the quadrature points | trace-back of every quadrature point | Courant near one, or below one with store smoothing | a cell-scale mode of the stress that grows below Courant one when the solvent viscosity is small | | `forward_integration_points` | discontinuous P1 per cell, fitted from the arrivals | fixed launch set of interior points (the integration points), one forward trajectory a step; the flux is read back at the launch points through a continuous P1 projection; an inflow cell's uncovered share is filled with the inflow value | the cylinder walls at dt 0.04; below Courant one with `flux_smoothing` at c = 0.023 (Waters-King 1/16, dt 0.0125: 0.9543 at t 1 and 0.5185 at t 6.5, against nodal 0.9622 and 0.5171) | the same cell-scale mode as the integration-point history without that smoothing (diverges at t 2.4 there); first order only; does not cross a periodic seam or follow a moving mesh | +| `forward_nodes` | continuous, the history's degree, at its nodes | the stress projected onto that store, launched from its nodes and from a lattice inside each element, one forward trajectory a step; a per-cell fit at the history's degree read back at the nodes | measured on transport alone (rotating diffusing Gaussian, P2: 1.78e-2 against 2.52e-2 for `backward_nodes` over half a turn) | not yet measured on a stress benchmark; first order only; does not follow a moving mesh | | `lagrangian` (particles) | a swarm the solver owns and advects, one value per particle, read through a discontinuous cells proxy | the material points themselves: the constitutive flux is evaluated at the particles each step and never projected back to the mesh; a particle that entered through an inflow takes the inflow value | any Courant number; no numerical diffusion of the history | the cost and bookkeeping of a swarm, and a proxy that needs its cells kept populated (population control refills them); the conformation check does not read a per-point tensor from it | | `eulerian` (SUPG grid) | continuous P1 | assembled transport equation with streamline upwinding | with DEVSS | without DEVSS the velocity block loses its preconditioner as the stress grows | @@ -69,6 +68,28 @@ where each flavour receives it, and anything a flavour writes from `evaluate` model with the timestep in kyr and as the same problem in plain numbers, from a moving start, at orders 1 and 2: the stores agree to solver precision. + +## In parallel + +Every history except the particle one gives the serial answer on any number of +ranks: `tests/parallel/test_1066` holds the six stress histories on a +turned-over Maxwell box and the four semi-Lagrangian value histories on a +rotating Gaussian to 1e-6 of their serial values at np 3, 4 and 6. Four things +make that so, and a new history has to respect them: + +- a field that is itself a mesh variable is recorded into a history by copying + its nodal values, never by evaluating it at its nodes; +- a trace starts from the node or point itself: a nudge toward "the nearest + cell's" centroid depends on which cells a rank holds; +- every history projection is solved to 1e-10, since its result is the carried + field and an iterative error is partition-dependent; +- a point is never given to one of two cells by a tie-break: a forward arrival + on a shared face is fitted in every cell that contains it, and a point the + parallel evaluator strands is evaluated by the rank whose cell contains it. + +The particle history differs from serial by about 1e-5 at np 4 and 6 (np 3 +matches); the cause is open. + ## The timestep is set by the wall strain rate, not the far-field Courant number The objective-rate source $L\sigma^* + \sigma^* L^T$ acts on the carried stress diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 96540e25..0c55057d 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -572,9 +572,10 @@ def global_evaluate_nd( expr, # rank whose nearest cell is globally closest, and Allreduce(SUM of # the winner-only value/flag) scatters that rank's extrapolation back. # - # A point some rank actually contains (distance ~ 0) naturally wins, so - # only genuinely-stranded points are corrected. Cost is O(boundary points) - # — no dense global tree, no exhaustive search. + # A point some rank's cell contains is then evaluated by that rank with the + # FE interpolant, not the rbf extrapolation (see the containment round + # below). Cost is O(stranded points) — no dense global tree, no + # exhaustive search. # # DEADLOCK SAFETY — read before editing. Every collective here (allgather, # Allreduce) runs unconditionally on the IDENTICAL global set on every @@ -650,6 +651,32 @@ def global_evaluate_nd( expr, best_flag = np.empty(n_ext_total, dtype=np.int32) comm.Allreduce([contrib_flag, MPI.INT], [best_flag, MPI.INT], op=MPI.SUM) + # A stranded point that a rank's cell CONTAINS is not out of the + # domain: the migration's claim (points_in_domain) is looser than + # cell containment, so a point a hair from a partition seam can be + # claimed by the neighbour, found in none of its cells and stranded. + # The rank that contains it evaluates it with the FE interpolant, + # exactly as serial evaluate() would (the rbf value above is only + # for points no rank contains). Every rank calls evaluate_nd, on + # the points it contains -- possibly none -- as the first pass does + # on the points it received, so the collectives stay in lockstep. + contains = np.asarray(mesh._robust_owning_cells(all_ext)) >= 0 + my_owner = np.where(contains, comm.rank, comm.size).astype(np.int32) + owner = np.empty(n_ext_total, dtype=np.int32) + comm.Allreduce([my_owner, MPI.INT], [owner, MPI.INT], op=MPI.MIN) + mine = owner == comm.rank + fe_vals, _fe_flag = evaluate_nd( + expr, np.ascontiguousarray(all_ext[mine]), rbf=rbf, evalf=evalf, + verbose=False, simplify=simplify, check_extrapolated=True,) + contrib_fe = np.zeros((n_ext_total,) + expr_shape, dtype=np.float64) + if mine.any(): + contrib_fe[mine] = np.asarray(fe_vals, dtype=np.float64).reshape((-1,) + expr_shape) + fe_val = np.empty_like(contrib_fe) + comm.Allreduce([contrib_fe, MPI.DOUBLE], [fe_val, MPI.DOUBLE], op=MPI.SUM) + contained = owner < comm.size + best_val[contained] = fe_val[contained] + best_flag[contained] = 0 + # Scatter this rank's segment of the global set back to its points. offset = int(counts[:comm.rank].sum()) seg = slice(offset, offset + ext_coords.shape[0]) diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 48c2d0e2..34043efc 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5923,6 +5923,10 @@ def estimate_dt(self, V_fn): import math import numpy as np + # TODO(BUG): with evalf=True evaluate() flags EVERY point as extrapolated, so + # at np > 1 this raises the 'not located on this rank' warning for particles + # that are all local (seen on the Lagrangian stress history, 2026-09-27). + # See planning file: underworld.md (Bugs section, 2026-09-27) vel = uw.function.evaluate(V_fn, self._particle_coordinates.data, evalf=True) # If vel is unit-aware (UnitAwareArray), nondimensionalise it to get diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 96d942ba..d7a034f6 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -223,6 +223,40 @@ def _history_units(psi_fn, units=None): return units if units is not None else uw.get_units(psi_fn) +# A history projection's result is the carried field itself, so it is solved +# to convergence: a mass-matrix solve left at a projection's default tolerance +# (1e-4) carries an error of that size into every step, and the error depends +# on the partition. Converging it costs a few more iterations. +_HISTORY_PROJECTION_TOLERANCE = 1.0e-10 + + +def _tracked_field(mesh, psi_fn, store): + """The mesh variable whose values ``psi_fn`` is, when it is laid out like + ``store`` (same degree, continuity and components), else ``None``. + + Recording such a field into a history is a copy of its nodal values: + exact, and the same on any partition. Evaluating it at the nodes instead + locates each node in a cell, and a node on a partition seam is located in + a different cell on each rank. ``psi_fn`` is the variable itself or its + symbol (``T.sym``, or ``T.sym[0]`` for a one-component field). + """ + if isinstance(psi_fn, uw.discretisation.MeshVariable): + field = psi_fn + else: + hit = uw.discretisation.meshVariable_lookup_by_symbol(mesh, psi_fn) + if hit is None and isinstance(psi_fn, sympy.MatrixBase) and psi_fn.shape == (1, 1): + hit = uw.discretisation.meshVariable_lookup_by_symbol(mesh, psi_fn[0, 0]) + if hit is None: + return None + field, component = hit + if component != -1 and field.num_components != 1: + return None + if (field.degree, field.continuous, field.num_components) != ( + store.degree, store.continuous, store.num_components): + return None + return field + + def _write_evaluated(var, values): """Write an evaluation of a store's quantity into the store. @@ -1128,6 +1162,33 @@ def _nondim_timestep(self, dt): reduced = _as_float(dt) return dt if reduced is None else reduced + def _project_nodally(self, expr, name="flux", smoothing=0.0, verbose=False): + """L2 projection of the stored components of ``expr`` onto a field of the + history's degree and continuity; returns that (1, ncomponents) field. + A flux holds gradients, which are discontinuous across elements, so it + is projected rather than read at nodes. Each ``name`` keeps its own + projection, so a history that projects two expressions does not + recompile one projection back and forth.""" + if not hasattr(self, "_nodal_projections"): + self._nodal_projections = {} + expr = sympy.Matrix(expr) + columns = _storage_components(self.vtype, expr.shape) + if name not in self._nodal_projections: + target = uw.discretisation.MeshVariable( + f"{name}_nodal_{self.instance_number}", self.mesh, (1, len(columns)), + vtype=uw.VarType.MATRIX, degree=self.degree, + continuous=self.continuous, + varsymbol=rf"{{{name}^{{\mathrm{{nodal}}}}_{{{self.instance_number}}}}}") + projection = uw.systems.solvers.SNES_MultiComponent_Projection( + self.mesh, u_Field=target, n_components=len(columns), verbose=self.verbose) + projection.tolerance = _HISTORY_PROJECTION_TOLERANCE + self._nodal_projections[name] = (target, projection) + target, projection = self._nodal_projections[name] + projection.uw_function = sympy.Matrix([[expr[i, j] for (i, j) in columns]]) + projection.smoothing = smoothing + projection.solve(verbose=verbose) + return target + def _write_inflow(self, var, coords, rows): """Overwrite ``rows`` of ``var`` with :attr:`inflow_value` evaluated at ``coords`` (the positions of ALL the points, so that the read, which @@ -1699,6 +1760,7 @@ def _setup_projections(self): verbose=False, ) self._psi_star_use_multicomponent = True + self._psi_star_projection_solver.tolerance = _HISTORY_PROJECTION_TOLERANCE self._psi_star_projection_solver.uw_function = self._build_projection_source( self.psi_fn) @@ -2487,15 +2549,13 @@ def _segment_exprs(self, seg): half = sympy.Rational(1, 2) return self.level_expr(k - 1), (self.level_expr(k - 1) + self.level_expr(k)) * half - def departure_points(self, key, X0, segments, evalf=False, X_eval=None, + def departure_points(self, key, X0, segments, evalf=False, clamp_final=True, subtract_v_mesh=False, v_mesh_var=None): r"""Trace ``X0`` back through ``segments`` (RK2 midpoint each): ``x_mid = x - dt/2 v_start(x)``, ``x_dep = x - dt v_mid(x_mid)``. - ``key`` names the launch node set (a variable's ``_basis_key`` plus - a tag for the nudge); ``X_eval`` are the points where the first - segment's start velocity is evaluated when they differ from ``X0`` - (the centroid-nudged nodes of the nodal history). Midpoints are + ``key`` names the launch point set (a variable's ``_basis_key`` plus + a tag). Midpoints are clamped to the domain; the last point is clamped unless ``clamp_final`` is False (old-frame reach-back). """ @@ -2513,8 +2573,7 @@ def departure_points(self, key, X0, segments, evalf=False, X_eval=None, continue kind, k, dt = segments[j] v_start, v_mid = self._segment_exprs(segments[j]) - X_start = X_eval if (j == 0 and X_eval is not None) else X - v0 = self.velocity_at(v_start, X_start, use_global=j > 0, evalf=evalf, + v0 = self.velocity_at(v_start, X, use_global=j > 0, evalf=evalf, subtract_v_mesh=subtract_v_mesh, v_mesh_var=v_mesh_var) Xm = X - v0 * (0.5 * dt) if clamp is not None: @@ -2969,6 +3028,7 @@ def __init__( verbose=False, ) self._psi_star_use_multicomponent = True + self._psi_star_projection_solver.tolerance = _HISTORY_PROJECTION_TOLERANCE # We should find a way to add natural bcs here # (self.Unknowns.u carried as a symbol from solver to solver) @@ -3283,40 +3343,15 @@ def _consume_ale_pulse(self): """ self._pending_v_mesh_disp = None - def _record_psi_star_from_field_data(self): - """Parallel-safe 'record current field into psi_star[0]'. - - The default record step evaluates ``psi_fn`` at its own node - coordinates, which under MPI mis-locates on-vertex points at a - process seam (first-pass ``get_closest_cells`` + FE extrapolation), - seeding a spurious history value. When ``psi_fn`` is a single - mesh-variable component living on this mesh with the same nodal - layout as ``psi_star[0]``, "evaluate at own nodes" is exactly that - variable's nodal data, so we copy it directly — no point location. - - Returns an array shaped like ``psi_star[0].array`` for that case, or - ``None`` (caller falls back to ``evaluate``) for non-scalar or - expression ``psi_fn`` (e.g. a flux with derivatives). - """ - try: - comps = list(self.psi_fn) # sympy Matrix, row-major - if len(comps) != 1: # scoped to scalar fields - return None - hit = uw.discretisation.meshVariable_lookup_by_symbol( - self.mesh, comps[0]) - if hit is None: - return None - var, comp = hit - vflat = np.asarray(var.array) - vflat = vflat.reshape(vflat.shape[0], -1) - out = np.array(np.asarray(self.psi_star[0].array)) - oflat = out.reshape(out.shape[0], -1) - if vflat.shape[0] != oflat.shape[0] or oflat.shape[1] != 1: - return None - oflat[:, 0] = vflat[:, comp] - return out - except Exception: - return None + def _copy_tracked_field(self): + """Record the current field into ``psi_star[0]`` by copying its nodal + values, when ``psi_fn`` is a mesh variable laid out like the store + (see :func:`_tracked_field`). Returns whether it did.""" + field = _tracked_field(self.mesh, self.psi_fn, self.psi_star[0]) + if field is None: + return False + self.psi_star[0].data[...] = field.data[...] + return True def _midtime_velocity_expr(self): r"""Velocity at :math:`t^{n+1/2}` for the mid-point stage of the @@ -3337,14 +3372,6 @@ def _record_velocity_history(self): if getattr(self, "_owns_characteristics", True): self.characteristics.finish_step() - def _centroid_shifted_var_coords(self, var): - """ND node coordinates of ``var`` nudged 0.1 % toward their cell - centroids (see :meth:`_centroid_shifted_node_coords`).""" - coords = np.asarray(var.coords_nd) - cellid = self.mesh.get_closest_cells(coords).reshape(-1) - cent = np.asarray(self.mesh._centroids)[cellid] - return 0.999 * coords + 0.001 * cent - def _velocity_nd_at( self, coords, @@ -3463,83 +3490,28 @@ def carried_tensors(self, level: int = 0): points = _to_nondim_ndarray(history.coords).reshape(-1, self.mesh.cdim) return values, points - def _centroid_shifted_node_coords(self): - r"""ND node coordinates of ``psi_star[0]``, nudged toward cell centroids. - - Point-location and FE interpolation are ambiguous exactly on - element edges/vertices (worst on quad meshes at the domain - boundary), so the sample points are moved 0.1 % of the way toward - the centroid of their owning cell: far enough to make cell - ownership unambiguous, close enough not to bias the sampled - values. Coordinates are plain non-dimensional arrays — never raw - ``.magnitude``, which would be dimensional metres (see - ``_to_nondim_ndarray`` and issue #267). - """ - psi_star_0_coords_nd = _to_nondim_ndarray(self.psi_star[0].coords) - - cellid = self.mesh.get_closest_cells( - psi_star_0_coords_nd, - ) - centroid_coords = self.mesh._centroids[cellid] - - shift = 0.001 - return (1.0 - shift) * psi_star_0_coords_nd[:, :] + shift * centroid_coords[ - :, : - ] - def _record_current_field_into_history( - self, node_coords_nd, evalf, verbose, oldframe_active + self, node_coords_nd, evalf, verbose ): r"""Record the current value of :math:`\psi` into ``psi_star[0]``. Three routes, in order of preference: - 1. direct nodal copy of the tracked field's data (parallel, or - old-frame reach-back); + 1. direct nodal copy of the tracked field's data, when ``psi_fn`` + is a mesh variable laid out like the store (:func:`_tracked_field`); 2. pointwise evaluation of ``psi_fn`` at the centroid-shifted - node coordinates (the validated serial path); + node coordinates; 3. an L2 projection for expressions that ``evaluate`` cannot handle (e.g. the NS viscous flux, which contains derivatives). """ try: - # Use shifted ND coords to avoid quad mesh boundary issues - # node_coords_nd is slightly shifted toward cell centroids - # evaluate() treats plain numpy as ND [0-1] coordinates. - # - # PARALLEL band-aid (parallel-singular-corruption, 2026-05): - # this "record current field into psi_star" step samples psi_fn - # at its OWN node coords. On-vertex sampling + first-pass - # get_closest_cells mis-locates at a process seam under MPI, - # recording a spurious history value that the implicit solve - # then propagates (the seam spike in adaptive advection- - # diffusion). When psi_fn is a single mesh-variable component on - # this mesh (the SLCN adv-diff case), "evaluate at own nodes" == - # the field's nodal data, so under MPI copy it directly (exact, - # no point location). Serial keeps the validated shifted- - # evaluate path bit-identically; non-scalar / expression psi_fn - # falls back to evaluate(). Proper fix (remap-on-adapt / ALE) - # tracked separately. - # Old-frame: record the history by a DIRECT nodal carry - # of the field rather than re-evaluating psi_fn at the - # (centroid-shifted) nodes of the DEFORMED mesh. The - # re-evaluate injects boundary-layer interpolation error - # that grows with mesh distortion and then rides the - # old-geometry sample below — the exact value we want is - # the carried nodal value (cf. the lagged-clone "store - # primitives" principle). Reuses the parallel direct-copy - # path, which returns None for non-scalar / expression - # psi_fn (those fall back to evaluate). - _direct = (self._record_psi_star_from_field_data() - if (uw.mpi.size > 1 or oldframe_active) else None) - if _direct is not None: - eval_result = _direct - else: + if not self._copy_tracked_field(): eval_result = uw.function.evaluate( self.psi_fn, node_coords_nd, evalf=evalf, ) - _write_evaluated(self.psi_star[0], eval_result) + _write_evaluated(self.psi_star[0], eval_result) except Exception: # Fallback to projection solver for expressions that can't be directly evaluated @@ -3566,7 +3538,7 @@ def _record_current_field_into_history( self.psi_star[0].array[:, j, i] = vals def _trace_departure_points( - self, i, node_coords_nd, dt_for_calc, evalf, subtract_v_mesh, oldframe_active + self, i, dt_for_calc, evalf, subtract_v_mesh, oldframe_active ): r"""RK2 midpoint trace-back: departure points for history slot ``i``. @@ -3586,16 +3558,14 @@ def _trace_departure_points( any foot that falls outside the old mesh, matching the validated prototype, which omits this clamp). """ - # One RK2 segment from the true nodes, the start velocity taken at - # the centroid-nudged coordinates (node_coords_nd). Served from the - # shared trace: a second history on the same nodes, or an older slot - # of this one, reuses the departure points computed here. + # One RK2 segment from the nodes. Served from the shared trace: a + # second history on the same nodes, or an older slot of this one, + # reuses the departure points computed here. return self.characteristics.departure_points( - (_basis_key_of(self.psi_star[i]), "nudged"), + (_basis_key_of(self.psi_star[i]), "nodes"), np.asarray(self.psi_star[i].coords_nd), (("first", 0, dt_for_calc),), evalf=evalf, - X_eval=node_coords_nd, clamp_final=not oldframe_active, subtract_v_mesh=subtract_v_mesh, v_mesh_var=getattr(self, "_v_mesh_var", None), @@ -3755,11 +3725,15 @@ def update_pre_solve( # When store_result=False (e.g. VE stress history), skip this — # psi_star[0] already contains the projected actual stress from # the previous solve and we want to advect *that*, not the flux. - node_coords_nd = self._centroid_shifted_node_coords() + # The nodes themselves: a node on a no-slip wall has zero velocity and + # departs from where it is. (A 0.1% nudge toward "the closest cell's" + # centroid gave it a velocity, and on a partition seam a different + # nudge on each rank.) + node_coords_nd = np.asarray(self.psi_star[0].coords_nd) if store_result: self._record_current_field_into_history( - node_coords_nd, evalf, verbose, _oldframe_active + node_coords_nd, evalf, verbose ) # 3. Trace the characteristics back and sample each history slot @@ -3793,7 +3767,7 @@ def update_pre_solve( _ale_this_iter = _ale_active and not (store_result and i == 0) end_pt_coords = self._trace_departure_points( - i, node_coords_nd, dt_for_calc, evalf, + i, dt_for_calc, evalf, _ale_this_iter, _oldframe_active, ) self._sample_history_at_departure( @@ -4696,6 +4670,43 @@ def update_post_solve( +def _hand_arrivals_to_owners(X, columns, locate): + """Give every arrival to the rank whose cell it landed in. + + ``locate`` returns the owning local cell of each point, -1 where the point + is in none of this rank's cells. Only the points that left their rank + travel: each rank offers its leavers to everyone and keeps the offered + points that land in its own cells (at Courant one that is the seam layer, + not the whole set). A point no rank takes has left the domain and is + dropped. The locators are local; the exchange is the only collective and + every rank makes it, a rank with no cells contributing and taking nothing + (comm.allgather rather than gather_data: the rows are vectors, and + gather_data flattens). + + Returns the kept points, their ``columns`` (arrays with a row per point), + their cells, and how many arrived from another rank. + """ + def owner(P): + return (np.asarray(locate(P), dtype=np.int64).reshape(-1) if P.shape[0] + else np.zeros(0, dtype=np.int64)) + + X = np.asarray(X) + columns = [np.asarray(c) for c in columns] + own = owner(X) + stay = own >= 0 + if uw.mpi.size == 1: + return X[stay], [c[stay] for c in columns], own[stay], 0 + comm = uw.mpi.comm + offered = np.concatenate(comm.allgather(X[~stay]), axis=0) + offered_columns = [np.concatenate(comm.allgather(c[~stay]), axis=0) for c in columns] + taken = owner(offered) + take = taken >= 0 + return (np.concatenate([X[stay], offered[take]], axis=0), + [np.concatenate([c[stay], o[take]], axis=0) for c, o in zip(columns, offered_columns)], + np.concatenate([own[stay], taken[take]]), + int(take.sum())) + + def _psi_shape_for(vtype, cdim): """The symbolic shape a history of ``vtype`` must have.""" if vtype == uw.VarType.SCALAR: @@ -4784,9 +4795,6 @@ class BackwardIntegrationPointsSemiLagrangian(_DDtBase): #: there when one is set (#745); without one they sample the edge. applies_inflow_value = True - _commit_projection = None - _commit_flat = None - def commit_flux_to_history(self, flux, verbose=False): """Project the new flux into the nodal snapshot and read it at the points, then shift both ladders. @@ -4806,22 +4814,8 @@ def commit_flux_to_history(self, flux, verbose=False): explicit and needs no snapshot substitution. """ history = self.psi_star[0] - flux = sympy.Matrix(flux) - columns = _storage_components(self.vtype, flux.shape) - - if self._commit_projection is None: - self._commit_flat = uw.discretisation.MeshVariable( - f"flux_nodal_{self.instance_number}", self.mesh, (1, len(columns)), - vtype=uw.VarType.MATRIX, degree=self.degree, - continuous=self.continuous, - varsymbol=rf"{{F^{{\mathrm{{nodal}}}}_{{{self.instance_number}}}}}") - self._commit_projection = uw.systems.solvers.SNES_MultiComponent_Projection( - self.mesh, u_Field=self._commit_flat, n_components=len(columns), - verbose=self.verbose) - self._commit_projection.uw_function = sympy.Matrix( - [[flux[i, j] for (i, j) in columns]]) - self._commit_projection.smoothing = self._store_smoothing_alpha() - self._commit_projection.solve(verbose=verbose) + columns = _storage_components(self.vtype, sympy.Matrix(flux).shape) + nodal_flux = self._project_nodally(flux, smoothing=self._store_smoothing_alpha(), verbose=verbose) # Oldest first, so each level reads the one above before it is written. for level in range(self.order - 1, 0, -1): @@ -4830,10 +4824,10 @@ def commit_flux_to_history(self, flux, verbose=False): points = np.asarray(history.integration_points).reshape(-1, self.mesh.cdim) for column in range(len(columns)): - nodal = self._commit_flat.data[:, column] + nodal = nodal_flux.data[:, column] self.psi_snap[0].data[:, column] = np.asarray(nodal).reshape(-1) history.data[:, column] = np.asarray(uw.function.evaluate( - self._commit_flat.sym[0, column], points)).reshape(-1) + nodal_flux.sym[0, column], points)).reshape(-1) self._history_committed = True @@ -5180,25 +5174,14 @@ def _object_viewer(self): display(Latex(rf"$\quad$History steps = {self.order} (at the integration points)")) # ------------------------------------------------------------------ - def _nudged_node_coords(self, var): - """ND node coordinates of ``var`` moved 0.1 % toward their cell - centroids so boundary nodes locate unambiguously (see - :meth:`BackwardNodesSemiLagrangian._centroid_shifted_node_coords`).""" - coords = np.asarray(var.coords_nd) - cellid = self.mesh.get_closest_cells(coords).reshape(-1) - cent = np.asarray(self.mesh._centroids)[cellid] - return 0.999 * coords + 0.001 * cent - def _record_current(self): """Snapshot slot 0 <- the current solution and velocity.""" ps = self.psi_snap[0] - if self._psi_meshVar is not None and ( - self._psi_meshVar.degree == ps.degree - and self._psi_meshVar.continuous == ps.continuous - ): - ps.data[...] = self._psi_meshVar.data[...] + field = _tracked_field(self.mesh, self._psi_meshVar or self.psi_fn, ps) + if field is not None: + ps.data[...] = field.data[...] else: - self._write_components(ps, self.psi_fn, self._nudged_node_coords(ps)) + self._write_components(ps, self.psi_fn, np.asarray(ps.coords_nd)) def _write_components(self, var, expr, coords, evaluate=None, **kwargs): """Evaluate ``expr`` at ``coords`` and store it component by component. @@ -5281,6 +5264,7 @@ def update_forcing_history(self, forcing_fn=None, evalf=False, verbose=False): self._forcing_projection = uw.systems.solvers.SNES_MultiComponent_Projection( self.mesh, u_Field=self._forcing_flat, n_components=len(columns), verbose=self.verbose) + self._forcing_projection.tolerance = _HISTORY_PROJECTION_TOLERANCE self._forcing_projection.uw_function = sympy.Matrix( [[forcing[i, j] for (i, j) in columns]]) self._forcing_projection.smoothing = 0.0 @@ -5595,6 +5579,7 @@ def _evaluate_at_launch(self, expr): if self._flux_projection is None: self._flux_projection = uw.systems.solvers.SNES_MultiComponent_Projection( self.mesh, u_Field=self._flux_var, n_components=self.num_components) + self._flux_projection.tolerance = _HISTORY_PROJECTION_TOLERANCE self._flux_projection.uw_function = sympy.Matrix([[expr[i, j] for (i, j) in self._components]]) self._flux_projection.smoothing = self.flux_smoothing self._flux_projection.solve() @@ -5636,38 +5621,15 @@ def _fit_arrivals(self, X, values, cell=None, inflow=None): ncell = self._cell_measure.size w = self._launch_weights if cell is None: - if uw.mpi.size > 1: - # Only the points that left this rank's partition travel: each rank - # offers its leavers to everyone and keeps the offered points that - # land in its own cells. At Courant one that is the seam layer, not - # the whole set. A rank with no cells owns nothing and keeps nothing. - # Ownership is by strict containment (face tolerance zero), not the - # evaluation locator's slab: a point a hair across a seam face would - # otherwise be kept by the rank it left and fitted into the wrong - # cell. The locators are local; the exchange is the only collective - # and every rank makes it, a rank with no cells contributing and - # taking nothing. (comm.allgather rather than gather_data: the rows - # are vectors, and gather_data flattens.) - X, values = np.asarray(X), np.asarray(values) - own = np.asarray(mesh._get_closest_local_cells_internal(X, tol=0.0), dtype=int).reshape(-1) - stay = own >= 0 - comm = uw.mpi.comm - offered = np.concatenate(comm.allgather(X[~stay]), axis=0) - offered_values = np.concatenate(comm.allgather(values[~stay]), axis=0) - offered_w = np.concatenate(comm.allgather(w[~stay]), axis=0) - taken = np.asarray(mesh._get_closest_local_cells_internal(offered, tol=0.0), dtype=int).reshape(-1) - take = taken >= 0 - self._n_relocated = int(take.sum()) - X = np.concatenate([X[stay], offered[take]], axis=0) - values = np.concatenate([values[stay], offered_values[take]], axis=0) - w = np.concatenate([w[stay], offered_w[take]], axis=0) - cell = np.concatenate([own[stay], taken[take]]) - else: - # the same strict rule as the parallel path, so a partition does - # not change which cell a point on a face is fitted into - cell = np.asarray(mesh._get_closest_local_cells_internal(X, tol=0.0), dtype=int).reshape(-1) - inside = cell >= 0 - X, values, w, cell = X[inside], values[inside], w[inside], cell[inside] + # Ownership is by strict containment (face tolerance zero), not the + # evaluation locator's slab: a point a hair across a seam face would + # otherwise be kept by the rank it left and fitted into the wrong + # cell; serially the same rule, so a partition does not change which + # cell a point on a face is fitted into. + def strict(P): + return np.asarray(mesh._get_closest_local_cells_internal(P, tol=0.0), dtype=int).reshape(-1) + X, (values, w), cell, self._n_relocated = _hand_arrivals_to_owners( + X, (values, w), strict) centroid = np.asarray(mesh._centroids)[:, :d] h = np.sqrt(self._cell_measure) if d == 2 else np.cbrt(self._cell_measure) # centred on the cell and scaled by its size, so the constant is c0 and the @@ -5785,22 +5747,32 @@ class ForwardNodesSemiLagrangian(_DDtBase): A field is known exactly at its own nodes (they are its unknowns) and, since its interpolant is a polynomial inside each element, at any point of an element's interior. Each step launches the values at the nodes and at a - lattice inside every element (the points of the discontinuous basis one - degree up), carries each one step forward along the characteristic, fits + lattice inside every element (the points of the discontinuous basis two + degrees up), carries each one step forward along the characteristic, fits the arrivals in each cell to a polynomial of the field's own degree, and - reads that fit back at the nodes. The fit never reaches across an element - boundary, where the interpolant has a kink. A cell with too few arrivals, or - arrivals on a line, falls back to a linear fit over the nearest arrivals; a - cell nothing reached keeps the field it launched (see + projects the per-cell fits (L2) onto the continuous store. The fit never + reaches across an element boundary, where the interpolant has a kink; the + projection weighs every cell that shares a node. A cell with too few arrivals + for that takes a linear fit to its own arrivals, or their mean; a cell + nothing reached keeps the field it launched (see :class:`~underworld3.utilities.cell_polynomial_projection.CellPolynomialProjector`). Arrivals that leave the domain are dropped; a node the flow reached from outside takes :attr:`inflow_value` when one is set. - The integration-point counterpart is :class:`ForwardIntegrationPointsSemiLagrangian`, which - launches from the quadrature points, where a flux (a stress) is formed, and - is the forward scheme for one: this class carries a field and refuses a flux. - First order; serial (a node near a partition seam needs arrivals from the - other rank); a fixed mesh. + A flux (a viscoelastic stress) is committed by L2 projection onto the + continuous store, since it holds gradients that are discontinuous across + elements; the next step launches from that projected field, which is a + polynomial inside each element like any other. The integration-point + counterpart, :class:`ForwardIntegrationPointsSemiLagrangian`, launches a + flux from the quadrature points where it is formed. + + An arrival on a face or vertex shared by several cells (a node on a + no-slip wall arrives where it started) is fitted in every one of them. In + parallel a node shared by two ranks is launched once, by its owner, and an + arrival that crosses a partition seam, or lands on one, is handed to every + rank with a cell that contains it, so each cell is fitted from everything + that reached it and the result does not depend on the partition. First + order; a fixed mesh. Parameters ---------- @@ -5826,15 +5798,13 @@ def __init__(self, mesh, psi_fn, V_fn, vtype=VarType.SCALAR, degree=1, varsymbol raise NotImplementedError("ForwardNodesSemiLagrangian carries one level; order must be 1") if mesh.cdim != mesh.dim: raise NotImplementedError("ForwardNodesSemiLagrangian fits in the embedding coordinates; no manifolds") - if uw.mpi.size > 1: - raise NotImplementedError("ForwardNodesSemiLagrangian is serial: the fit at a node near a " - "partition seam needs the arrivals on the other rank") if _unsupported: warnings.warn(f"ForwardNodesSemiLagrangian ignores {sorted(_unsupported)}", stacklevel=2) self.vtype = vtype self.mesh = mesh self.degree = int(degree) self.continuous = True + self.verbose = False self.order = 1 self.theta = float(theta) self.V_fn = V_fn @@ -5855,8 +5825,18 @@ def __init__(self, mesh, psi_fn, V_fn, vtype=VarType.SCALAR, degree=1, varsymbol self._fit_var = uw.discretisation.MeshVariable( f"fit_fwn_{self.instance_number}", mesh, vtype=vtype, degree=self.degree, continuous=False) self._projector = CellPolynomialProjector(self._fit_var) - self._interior = np.asarray(mesh._get_coords_for_basis(self.degree + 1, continuous=False) + # two degrees up: enough arrivals that every cell of the rotating + # Gaussian fits at the field's degree (one degree up left ~8% of the + # cells to the linear fallback; three gains nothing) + self._interior = np.asarray(mesh._get_coords_for_basis(self.degree + 2, continuous=False) ).reshape(-1, mesh.cdim) + # A node on a partition seam is on both ranks and is launched once, by its + # owner: writing the rank into the store and reading it back after the + # ghost update leaves every copy holding its owner's rank. + self.psi_star[0].data[:, 0] = uw.mpi.rank + self._owned = np.asarray(self.psi_star[0].data[:, 0]) == uw.mpi.rank + self.psi_star[0].data[:, :] = 0.0 + self._n_relocated = 0 self._n_v = 2 self._init_coefficient_expressions(1, self.theta, with_exp=True) self._register_with_default_model() @@ -5893,26 +5873,73 @@ def _values_at(self, expr, X): np.asarray(_to_nondim_ndarray(uw.function.evaluate(expr[i, j], X))).reshape(-1) for (i, j) in self._components]) - def _reconstruct(self, arrivals, values, nodes): - """Fit the arrivals in each cell at the field's degree; the fit at the nodes. - A cell nothing reached keeps the field it launched.""" - inside = np.asarray(self.mesh.points_in_domain(arrivals), dtype=bool) - launched = self._values_at(self._psi_fn, np.asarray(self._fit_var.coords_nd)) - fit = self._projector.fit(arrivals[inside], values[inside], old=launched) - at_nodes = self._projector.interpolate(fit, nodes) - if np.isnan(at_nodes).any(): - raise RuntimeError("ForwardNodesSemiLagrangian: a node lies in no cell of the fit") - return at_nodes + _FACE_TOLERANCE = 1.0e-9 + + def _arrivals_by_cell(self, X, values): + """Every arrival with every cell that contains it, on the rank that owns + the cell: one row per (arrival, cell). + + An arrival strictly inside a cell of this rank belongs to that cell + alone. One on a face or vertex (a node on a no-slip wall arrives where + it started) belongs to every cell that shares it, and one in none of + this rank's cells has left the rank; both are offered to every rank, + which takes them into each of its own cells that contains them. No + tie-break between cells is made, so the assignment, and the fit, are the + same on any partition. + """ + pj = self._projector + tol = self._FACE_TOLERANCE + point, cell, lam = pj.containing_cells(X, tol) + strict = np.zeros(X.shape[0], dtype=bool) + strict[point[lam > tol]] = True + keep = strict[point] + offered_X, offered_v = X[~strict], values[~strict] + if uw.mpi.size > 1: + comm = uw.mpi.comm + parts_X = comm.allgather(offered_X) + parts_v = comm.allgather(offered_v) + offered_X = np.concatenate(parts_X, axis=0) + offered_v = np.concatenate(parts_v, axis=0) + mine_from = sum(p.shape[0] for p in parts_X[:uw.mpi.rank]) + n_mine = parts_X[uw.mpi.rank].shape[0] + else: + mine_from, n_mine = 0, offered_X.shape[0] + q, qcell, _ = pj.containing_cells(offered_X, tol) + own = (q >= mine_from) & (q < mine_from + n_mine) + self._n_relocated = int(np.unique(q[~own]).size) + return (np.concatenate([X[point[keep]], offered_X[q]], axis=0), + np.concatenate([values[point[keep]], offered_v[q]], axis=0), + np.concatenate([cell[keep], qcell])) + + def _reconstruct(self, arrivals, values, source): + """Fit the arrivals in each cell at the field's degree, and project the + per-cell fits onto the continuous store. A node is shared by several + cells whose fits differ slightly; the projection weighs them all, where + reading one cell's fit would depend on which cell (and in parallel which + rank) the node was located in. A cell nothing reached keeps the field it + launched (``source``).""" + arrivals, values, cells = self._arrivals_by_cell(arrivals, values) + launched = self._values_at(source, np.asarray(self._fit_var.coords_nd)) + self._fit_var.data[:, :] = self._projector.fit(arrivals, values, old=launched, + cell_local=True, cells=cells) + return np.array(self._project_nodally(self._fit_var.sym, name="fit").data) # ------------------------------------------------------------------ def initialise_history(self): - """Start from the current field at the nodes.""" + """Start from the current field at the nodes. A history already placed + by :meth:`commit_flux_to_history` is the start, and is kept.""" self.characteristics.initialise_levels(self._n_v) - self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) + if not self._history_committed: + self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) self._history_initialised = True - def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): - """Carry the field forward one step from the nodes and rebuild it there.""" + def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, **_ignored): + """Carry the field forward one step from the nodes and rebuild it there. + + ``store_result=False`` says the store already holds the values to launch + (a flux placed by :meth:`commit_flux_to_history`); otherwise the tracked + field is read first. + """ self._dt = dt = self._nondim_timestep(dt) if self._projector.mesh_version != self.mesh._mesh_version: raise NotImplementedError( @@ -5926,19 +5953,27 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): if self._owns_characteristics: trace.begin_step(dt) nodes = self._nodes() - launch = np.vstack([nodes, self._interior]) - values = self._values_at(self._psi_fn, launch) + owned = nodes[self._owned] + launch = np.vstack([owned, self._interior]) + # the tracked field, or the committed store, which is a polynomial + # inside each element, so its interior values are exact + source = self._psi_fn if store_result else self.psi_star[0].sym + at_owned = (self._values_at(source, owned) if store_result + else np.asarray(self.psi_star[0].data)[self._owned]) + values = np.vstack([at_owned, self._values_at(source, self._interior)]) key = (_basis_key_of(self.psi_star[0]), "launch") arrivals = np.asarray(trace.departure_points(key, launch, (("first", 0, -float(dt)),), evalf=evalf, clamp_final=False)) - self.psi_star[0].data[:, :] = self._reconstruct(arrivals, values, nodes) - arrivals = arrivals[:nodes.shape[0]] + carried = self._reconstruct(arrivals, values, source) if self._inflow_value is not None: - # a node whose back-trace leaves the domain holds fluid that entered this step - back = 2.0 * nodes - arrivals - entered = ~np.asarray(self.mesh.points_in_domain(back), dtype=bool) - if entered.any(): - self._write_inflow(self.psi_star[0], nodes, entered) + # an owned node whose back-trace leaves the domain holds fluid that + # entered this step (a ghost copy takes its owner's value) + back = 2.0 * owned - arrivals[:owned.shape[0]] + entered = np.zeros(nodes.shape[0], dtype=bool) + entered[np.flatnonzero(self._owned)] = ~np.asarray(self.mesh.points_in_domain(back), dtype=bool) + inflow = self._values_at(self._inflow_record(), nodes) + carried[entered] = inflow[entered] + self.psi_star[0].data[:, :] = carried if self._owns_characteristics: trace.finish_step() @@ -5952,11 +5987,10 @@ def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): self._n_solves_completed += 1 def commit_flux_to_history(self, flux, verbose=False): - """Refused: a flux is formed at the integration points, not known at the nodes.""" - raise NotImplementedError( - "ForwardNodesSemiLagrangian carries a field known at its nodes; a flux (a stress) " - "is formed at the integration points, so carry it with " - "ForwardIntegrationPointsSemiLagrangian") + """Project the new flux onto the store, where the next step launches it.""" + projected = self._project_nodally(flux, verbose=verbose) + self.psi_star[0].data[:, :] = np.asarray(projected.data) + self._history_committed = True _SEMI_LAGRANGIAN_SCHEMES = { diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index d8b80618..507ae941 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -462,16 +462,11 @@ def _invalidate_solution_cache(u): from .ddt import Symbolic as Symbolic_DDt # The semi-Lagrangian schemes a solver can build its history with, named -# "_" after the arguments of ddt.SemiLagrangian. A stress is -# formed at the integration points, so it is not carried forward from nodes. +# "_" after the arguments of ddt.SemiLagrangian. _SEMI_LAGRANGIAN_TRANSPORTS = ( "backward_nodes", "backward_integration_points", "forward_integration_points", "forward_nodes", ) -_STRESS_TRANSPORTS = ( - "backward_nodes", "backward_integration_points", "forward_integration_points", - "lagrangian", "eulerian", -) _RENAMED_TRANSPORTS = { "semi_lagrangian": "backward_nodes", "integration_point": "backward_integration_points", @@ -1565,9 +1560,10 @@ def set_jacobian_F1_source(self, F1_source, linesearch="cp"): def stress_transport(self) -> str: """How a viscoelastic stress history is carried: ``"backward_nodes"`` (default), ``"backward_integration_points"``, - ``"forward_integration_points"``, ``"lagrangian"`` or ``"eulerian"``. + ``"forward_integration_points"``, ``"forward_nodes"``, ``"lagrangian"`` + or ``"eulerian"``. - The first three are semi-Lagrangian schemes of + The first four are the semi-Lagrangian schemes of :func:`~underworld3.systems.ddt.SemiLagrangian`, named by the direction of the trace and the points the history is held at. A backward trace follows the characteristic back from each storage point and samples the @@ -1583,9 +1579,9 @@ def stress_transport(self) -> str: through a continuous P1 projection; it holds the Maxwell start-up below Courant one where the backward integration-point history rings (see :class:`~underworld3.systems.ddt.ForwardIntegrationPointsSemiLagrangian`). - The fourth semi-Lagrangian scheme, forward from nodes, carries a field - known at its nodes; a stress is formed at the integration points, so it - is not offered here. + ``"forward_nodes"`` launches the stress projected onto the continuous + history space from its nodes and from a lattice inside each element (see + :class:`~underworld3.systems.ddt.ForwardNodesSemiLagrangian`). ``"eulerian"`` transports the stress on the grid with the same streamline-upwind stabilisation the Eulerian solvers use, and gives the @@ -1610,9 +1606,10 @@ def stress_transport(self, value): f"stress_transport={value!r} is now {_RENAMED_TRANSPORTS[value]!r}", FutureWarning, stacklevel=2) value = _RENAMED_TRANSPORTS[value] - if value not in _STRESS_TRANSPORTS: + if value not in _SEMI_LAGRANGIAN_TRANSPORTS + ("lagrangian", "eulerian"): raise ValueError( - f"stress_transport must be one of {_STRESS_TRANSPORTS}, not {value!r}.") + f"stress_transport must be one of {_SEMI_LAGRANGIAN_TRANSPORTS} or " + f"'lagrangian' or 'eulerian', not {value!r}.") if self.Unknowns.DFDt is not None: raise RuntimeError( "the stress history already exists: set stress_transport before the " @@ -1836,17 +1833,19 @@ def _create_stress_history_ddt(self, order=2): **ddt_kwargs, **{k: v for k, v in common.items() if k != "smoothing"}, ) - elif self.stress_transport == "forward_integration_points": + elif self.stress_transport.startswith("forward_"): if ddt_kwargs: raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " - f"{sorted(ddt_kwargs)}, which the forward flavour does not provide; " + f"{sorted(ddt_kwargs)}, which the forward flavours do not provide; " "use stress_transport='backward_nodes' for it.") - self.Unknowns.DFDt = uw.systems.ddt.ForwardIntegrationPointsSemiLagrangian( + self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( self.mesh, sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), self.u.sym, - vtype=common["vtype"], degree=common["degree"], varsymbol=common["varsymbol"], + common["vtype"], + trace="forward", launch=self.stress_transport[len("forward_"):], + degree=common["degree"], varsymbol=common["varsymbol"], order=order, units=common["units"], ) elif self.stress_transport == "lagrangian": @@ -4463,8 +4462,7 @@ class SNES_AdvectionDiffusion(SNES_Scalar): does not take (``monotone_mode`` on a forward scheme, say) is refused. Only the value history is chosen: the diffusive flux history ``DFDt`` is always backward from the nodes. ``"forward_integration_points"`` - fits a linear polynomial per cell and so needs a degree-1 field; - ``"forward_nodes"`` runs in serial only. + fits a linear polynomial per cell and so needs a degree-1 field. old_frame_traceback : bool, default=False Use the old-frame semi-Lagrangian reach-back for the advective ``DuDt`` history on a moving mesh (free surface or interior-node diff --git a/src/underworld3/utilities/cell_polynomial_projection.py b/src/underworld3/utilities/cell_polynomial_projection.py index c07b04ad..e1456004 100644 --- a/src/underworld3/utilities/cell_polynomial_projection.py +++ b/src/underworld3/utilities/cell_polynomial_projection.py @@ -106,6 +106,40 @@ def reference_coords(self, coords, cells): """Reference coordinates (PETSc's [-1, 1] frame) of points in their cells.""" return np.einsum("cij,cj->ci", self.invJ[cells], coords - self.v0[cells]) - 1.0 + def containing_cells(self, coords, tol=1.0e-9): + """Every local cell that contains each point. + + Returns ``(point, cell, lam)``: one row per (point, containing cell) + pair, with ``lam`` the point's distance inside that cell's nearest face + in reference units (the smallest barycentric coordinate of a simplex; + ``1 - max |xi|`` of a quadrilateral or hexahedron, through the same + affine cell map the fit uses). A point on a face or vertex shared by + several cells is in all of them (``lam`` within ``tol`` of zero); a + point strictly inside a cell is in that one only; a point in no local + cell has no row. + """ + coords = np.asarray(coords, dtype=np.float64).reshape(-1, self.dim) + if coords.shape[0] == 0 or self.ncells == 0: + empty = np.zeros(0, dtype=np.int64) + return empty, empty, np.zeros(0) + if getattr(self, "_centroid_tree", None) is None: + self._centroid_tree = uw.kdtree.KDTree(self.centroids) + # enough neighbours to hold every cell around a vertex + k = min(self.ncells, 16 if self.dim == 2 else 64) + _, near = self._centroid_tree.query(coords, k=k) + near = np.asarray(near, dtype=np.int64).reshape(coords.shape[0], k) + point = np.repeat(np.arange(coords.shape[0]), k) + cell = near.reshape(-1) + xi = self.reference_coords(np.repeat(coords, k, axis=0), cell) + if self.mesh.isSimplex: + # PETSc's reference simplex has vertices at -1 and +1 on each axis + lam_axes = 0.5 * (xi + 1.0) + lam = np.minimum(lam_axes.min(axis=1), 1.0 - lam_axes.sum(axis=1)) + else: + lam = 0.5 * (1.0 - np.abs(xi).max(axis=1)) + inside = lam >= -tol + return point[inside], cell[inside], lam[inside] + def locate(self, coords): """Owning local cell of each point (-1 when not on this rank) and its reference coordinates.""" coords = np.asarray(coords, dtype=np.float64) @@ -118,7 +152,8 @@ def locate(self, coords): # -- the fit ------------------------------------------------------------ - def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e6): + def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e6, + cell_local=False, cells=None): """Fit every cell; returns nodal values shaped like ``meshVar.data``. Parameters @@ -137,11 +172,23 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e advection lie on a line, and the P2 fit of a line is singular (measured: condition 1e300 at 92 particles, garbage that grew by 1e12 in ten steps through the read-back). + cell_local : a thin cell takes a lower-degree fit to its OWN points + (linear, or their mean when too few or collinear for a gradient) + instead of the linear patch fit to the nearest points, and a cell + with no points keeps ``old`` from the first fit on. No cell then + reads a point outside it, so the fit is the same on any partition. + cells : the local cell of each row, when the caller has assigned them + (see :meth:`containing_cells`); otherwise each point is located. """ coords = np.asarray(coords, dtype=np.float64).reshape(-1, self.dim) values = np.asarray(values, dtype=np.float64).reshape(coords.shape[0], -1) nc = values.shape[1] - cells, ok, xi = self.locate(coords) + if cells is None: + cells, ok, xi = self.locate(coords) + else: + cells = np.asarray(cells, dtype=np.int64) + ok = np.ones(cells.shape[0], dtype=bool) + xi = self.reference_coords(coords, cells) if cells.shape[0] else np.zeros_like(coords) c = cells[ok] B = self._scalar_basis(xi[ok]) # (Np, Nb) psi = values[ok] # (Np, nc) @@ -170,12 +217,14 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e self.n_empty = int((npc == 0).sum()) held = np.zeros(self.ncells, dtype=bool) - if old is not None and getattr(self, "_has_fit", False): + if old is not None and (cell_local or getattr(self, "_has_fit", False)): held = npc == 0 U[held] = np.asarray(old, dtype=np.float64).reshape(self.ncells, self.Nb, nc)[held] thin = np.nonzero(~dense & ~held)[0] self.n_thin = int(thin.shape[0]) - if thin.shape[0] > 0 and c.shape[0] > 0: + if cell_local and thin.shape[0] > 0: + U[thin] = self._cell_linear_fit(thin, c, xi[ok], psi, cond_max) + elif thin.shape[0] > 0 and c.shape[0] > 0: # Linear fit (monomials 1, xi_1, ..., xi_dim in the cell's frame) # to the nearest particles, evaluated at the cell's dof nodes. Xp = coords[ok] @@ -203,6 +252,29 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e self._has_fit = True return U.reshape(self.ncells * self.Nb, nc) + def _cell_linear_fit(self, cells, c, xi, psi, cond_max): + """Linear least-squares fit (1, xi_1, ..., xi_dim) of each of ``cells`` + to its own points, at the cell's dof nodes; the points' mean where they + cannot carry a gradient (fewer than dim + 2, or collinear).""" + npar = self.dim + 1 + slot = np.full(self.ncells, -1, dtype=np.int64) + slot[cells] = np.arange(cells.shape[0]) + mine = slot[c] >= 0 + s, A, p = slot[c[mine]], np.concatenate([np.ones((int(mine.sum()), 1)), xi[mine]], axis=1), psi[mine] + G = np.zeros((cells.shape[0], npar, npar)) + R = np.zeros((cells.shape[0], npar, p.shape[1])) + np.add.at(G, s, A[:, :, None] * A[:, None, :]) + np.add.at(R, s, A[:, :, None] * p[:, None, :]) + count = np.bincount(s, minlength=cells.shape[0]) + coef = np.zeros_like(R) + coef[:, 0, :] = R[:, 0, :] / np.maximum(count, 1)[:, None] # the mean + ev = np.linalg.eigvalsh(G) + linear = (count >= self.dim + 2) & (ev[:, -1] <= cond_max * np.maximum(ev[:, 0], 1e-300)) + if linear.any(): + coef[linear] = np.linalg.solve(G[linear], R[linear]) + Adof = np.concatenate([np.ones((self.Nb, 1)), self.xi_dof], axis=1) # (Nb, dim+1) + return np.einsum("ba,cak->cbk", Adof, coef) + def interpolate(self, U, coords): """The fitted polynomials evaluated at points (NaN off-rank): the FLIP read-back.""" U = np.asarray(U).reshape(self.ncells, self.Nb, -1) diff --git a/tests/parallel/test_1062_forward_stress_history_mpi.py b/tests/parallel/test_1062_forward_stress_history_mpi.py index 8327099b..8ef5dcf4 100644 --- a/tests/parallel/test_1062_forward_stress_history_mpi.py +++ b/tests/parallel/test_1062_forward_stress_history_mpi.py @@ -17,7 +17,7 @@ pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, pytest.mark.mpi(min_size=2), pytest.mark.timeout(600)] # BASELINES: the serial values (see the ledger) -XY_AT_A, XY_AT_B = -0.1190789, 0.0702257 +XY_AT_A, XY_AT_B = -0.1190800, 0.0702081 POINTS = np.array([[0.5, 0.2], [-0.3, -0.35]]) @@ -50,10 +50,8 @@ def test_the_forward_history_gives_the_serial_stress_on_every_rank(): kind, values, relocated = turned_over_maxwell_box() assert kind == "ForwardIntegrationPointsSemiLagrangian" # The serial values are the hard baseline (they pin the physics; regenerate - # them if a default changes). Parallel matches them to 1e-4, not to - # round-off: a cell's arrivals are summed into its least-squares fit in an - # order the partition sets, and ten steps of that reordering reach ~2e-5 on - # a mildly conditioned fit. A dropped or misrouted arrival would be far larger. - assert abs(values[0] - XY_AT_A) < 1.0e-4 and abs(values[1] - XY_AT_B) < 1.0e-4, values + # them if a default changes). Parallel matches them to the recorded figures + # (see test_1066); a dropped or misrouted arrival costs 1e-5 or more. + assert abs(values[0] - XY_AT_A) < 1.0e-6 and abs(values[1] - XY_AT_B) < 1.0e-6, values # the seam must actually have been crossed for this to test the parallel path assert relocated > 0 diff --git a/tests/parallel/test_1066_transport_schemes_mpi.py b/tests/parallel/test_1066_transport_schemes_mpi.py new file mode 100644 index 00000000..1e302347 --- /dev/null +++ b/tests/parallel/test_1066_transport_schemes_mpi.py @@ -0,0 +1,78 @@ +"""Every stress and advection history in parallel: np >= 3 equals serial. + +Stress: the turned-over Maxwell box of test_1062 (a shear modulus varying in x, +two counter-rotating cells between no-slip walls, an unstructured mesh), so +the stress is non-uniform and the flow crosses seams of either orientation. +Advection: a rotating Gaussian in a square box (the flow crosses every wall), +P2, a quarter turn, with each of the four semi-Lagrangian value histories. The values at two points after the +run must be the serial ones. At least three ranks: two ranks meet only along +one seam, while three or more also meet at points, where an arrival can be +handed to either of two other ranks. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from test_1062_forward_stress_history_mpi import turned_over_maxwell_box + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, pytest.mark.mpi(min_size=3), pytest.mark.timeout(1200)] + +# BASELINES: the serial values at the points of each test (2026-09-27) +STRESS_XY = { + "backward_nodes": (-0.1001236, 0.1000386), + "backward_integration_points": (-0.1166604, 0.0682694), + "forward_integration_points": (-0.1190800, 0.0702081), + "forward_nodes": (-0.1190115, 0.0701718), + "lagrangian": (-0.1214269, 0.0672075), + "eulerian": (-0.1186881, 0.0699410), +} +ADVECTED_T = { + "backward_nodes": (0.7259154, 0.0751478), + "backward_integration_points": (0.7533908, 0.0739654), + "forward_integration_points": (0.6347829, 0.0836136), + "forward_nodes": (0.7694358, 0.0721413), +} +# np 3, 4 and 6 give the serial values to the 7 figures recorded: the history +# projections are converged to 1e-10 and no stage depends on which rank, or +# which of two cells, a point is assigned to. A misrouted or dropped value +# costs 1e-5 or more. +ATOL = 1.0e-6 +T_POINTS = np.array([[0.0, 0.5], [-0.2, 0.35]]) + + +def rotating_gaussian(transport, steps=16, dt=np.pi / 32): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), + cellSize=0.1, qdegree=3, regular=False) + x, y = mesh.X + degree = 1 if transport == "forward_integration_points" else 2 + T = uw.discretisation.MeshVariable("T_rg", mesh, 1, degree=degree) + T.array[:, 0, 0] = uw.function.evaluate( + sympy.exp(-((x - 0.5) ** 2 + y ** 2) / (2 * 0.1 ** 2)), T.coords).reshape(-1) + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=sympy.Matrix([[-y, x]]), order=1, + transport=transport) + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = 1.0e-3 + for wall in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, wall) + adv.tolerance = 1.0e-10 + for _ in range(steps): + adv.solve(timestep=dt) + return np.asarray(uw.function.global_evaluate(T.sym[0], T_POINTS)).reshape(-1) + + +@pytest.mark.parametrize("transport", [ + pytest.param(t, marks=pytest.mark.xfail( + uw.mpi.size >= 4, strict=True, reason="TODO(BUG): the particle history differs from serial by 1.2e-5 " + "at np >= 4 (np 3 matches); the particle values, their ownership and the cell " + "proxy's fit are partition-independent -- cause not yet found")) + if t == "lagrangian" else t for t in STRESS_XY]) +def test_every_stress_history_gives_the_serial_stress_on_every_rank(transport): + _kind, values, _relocated = turned_over_maxwell_box(transport) + assert np.allclose(values, STRESS_XY[transport], atol=ATOL), (transport, values) + + +@pytest.mark.parametrize("transport", list(ADVECTED_T)) +def test_every_value_history_gives_the_serial_field_on_every_rank(transport): + values = rotating_gaussian(transport) + assert np.allclose(values, ADVECTED_T[transport], atol=ATOL), (transport, values) diff --git a/tests/test_1059_stress_transport.py b/tests/test_1059_stress_transport.py index 099d1409..12625e13 100644 --- a/tests/test_1059_stress_transport.py +++ b/tests/test_1059_stress_transport.py @@ -226,6 +226,7 @@ def _maxwell_shear(transport, order, steps=20, dt=0.1, integrator="bdf", solver= "backward_nodes": "BackwardNodesSemiLagrangian", "backward_integration_points": "BackwardIntegrationPointsSemiLagrangian", "forward_integration_points": "ForwardIntegrationPointsSemiLagrangian", + "forward_nodes": "ForwardNodesSemiLagrangian", "lagrangian": "Lagrangian", "eulerian": "EulerianSUPG", } @@ -292,9 +293,6 @@ def test_the_former_stress_transport_names_still_select_their_scheme(): with pytest.warns(FutureWarning, match=new): stokes.stress_transport = old assert stokes.stress_transport == new - # a stress is formed at the integration points: it is not carried from the nodes - with pytest.raises(ValueError, match="stress_transport must be"): - stokes.stress_transport = "forward_nodes" def test_semi_lagrangian_selects_its_scheme_by_trace_and_launch(): @@ -314,8 +312,6 @@ def test_semi_lagrangian_selects_its_scheme_by_trace_and_launch(): # an option the chosen scheme does not take is refused, not dropped with pytest.raises(TypeError, match="monotone_mode"): uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, trace="forward", degree=1, monotone_mode="clamp") - with pytest.raises(NotImplementedError, match="integration points"): - uw.systems.ddt.SemiLagrangian(mesh, T.sym, V, trace="forward", degree=1).commit_flux_to_history(T.sym) with pytest.warns(FutureWarning, match="ForwardIntegrationPointsSemiLagrangian"): assert uw.systems.ddt.ForwardSemiLagrangian is uw.systems.ddt.ForwardIntegrationPointsSemiLagrangian @@ -766,11 +762,11 @@ def test_the_store_smoothing_is_the_coefficient_times_the_local_cell_size_square stokes = _one_shear_step("backward_integration_points", dt=1.0) history = stokes.DFDt assert history.store_smoothing == 0.0 - assert history._commit_projection.smoothing == 0.0 + assert history._nodal_projections["flux"][1].smoothing == 0.0 history.store_smoothing = 0.05 stokes.solve(timestep=1.0, zero_init_guess=False) # the projection's smoothing is now a field: c times the cell-size field squared - alpha = history._commit_projection.smoothing + alpha = history._nodal_projections["flux"][1].smoothing x0 = np.array([[0.1, 0.1]]) h = float(np.asarray(uw.function.evaluate(stokes.mesh.cell_size(), x0)).reshape(-1)[0]) a = float(np.asarray(uw.function.evaluate(alpha, x0)).reshape(-1)[0]) diff --git a/tests/test_1102_forward_nodes_rotating_gaussian.py b/tests/test_1102_forward_nodes_rotating_gaussian.py index 9a990ba9..650d92d8 100644 --- a/tests/test_1102_forward_nodes_rotating_gaussian.py +++ b/tests/test_1102_forward_nodes_rotating_gaussian.py @@ -41,7 +41,7 @@ def _run(mesh, walls): def test_forward_from_nodes_on_the_disc(): uw.reset_default_model() err, peak = _run(_disc(24), ("Upper",)) - assert abs(err - 0.01996) < 0.002, err # BASELINE (2026-09-26) + assert abs(err - 0.01780) < 0.0018, err # BASELINE (2026-09-27) assert abs(peak - EXACT_PEAK) < 0.002, (peak, EXACT_PEAK) @@ -50,5 +50,5 @@ def test_forward_from_nodes_holds_the_box(): mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=2.0 / 24, qdegree=3, regular=False) err, peak = _run(mesh, ("Left", "Right", "Top", "Bottom")) - assert abs(err - 0.0617) < 0.006, err # BASELINE (2026-09-26) + assert abs(err - 0.0657) < 0.0065, err # BASELINE (2026-09-27) assert abs(peak - EXACT_PEAK) < 0.005, (peak, EXACT_PEAK) From 98c4193617f76fbde62e7bcf09ad85191d04527c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 27 Sep 2026 01:55:53 -0700 Subject: [PATCH 4/5] Review (b20c992e): inflow across seams, NaN fallback, reproducible history projections - forward from nodes: an inflow node is one whose back-reflected point the global restore moves (a per-rank "in domain" test called every partition face a boundary); ownership of seam nodes read from the global section (_owned_rows) instead of a write-and-read-back; the tracked field copied, not evaluated, at the nodes; a strictly-inside arrival fitted in its one cell only; containing_cells falls back to the locator for a point whose cell is not among the nearest centroids; quadrilateral centroids corrected. - global_evaluate: the containment round keeps a finite rbf value where the FE value is NaN, and runs only where the cell hint is authoritative (no DMLocatePoints collective inside it); the deadlock note says so. - every history projection is a CG/Jacobi solve (linear_solver moved to the projection mixin) to 1e-12, from zero: the warm start depended on a state a restart does not restore. - Eulerian record uses _tracked_field; _tracked_field refuses another mesh's variable; stress_transport names map through one (trace, launch) table. - tests: test_1066 resets the model per case, asserts the forward exchanges ran, gives the value histories an inflow, marks the particle case #797 (xfail, np >= 4); test_1063 covers forward_nodes; test_1102 held to 1e-4. test_1066 at np 4 is the regression test for the containment round. - TODO(DESIGN) notes: face points offered to every rank; the two forward flavours' arrival rules. np 3/4/6: 10 passed (+1 xfail at 4/6); serial: 388 passed. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- .../design/REMESH_FIELD_TRANSFER_DESIGN.md | 4 +- .../design/lagged-clone-sl-history.md | 2 +- src/underworld3/function/_function.pyx | 56 +++++--- src/underworld3/swarm.py | 2 +- src/underworld3/systems/ddt.py | 126 +++++++++++------- src/underworld3/systems/solvers.py | 120 +++++++++-------- .../utilities/cell_polynomial_projection.py | 37 +++-- .../test_1066_transport_schemes_mpi.py | 19 ++- tests/test_1063_stress_history_restart.py | 2 +- ...st_1102_forward_nodes_rotating_gaussian.py | 4 +- 10 files changed, 226 insertions(+), 146 deletions(-) diff --git a/docs/developer/design/REMESH_FIELD_TRANSFER_DESIGN.md b/docs/developer/design/REMESH_FIELD_TRANSFER_DESIGN.md index 73cf7058..0fd2cb20 100644 --- a/docs/developer/design/REMESH_FIELD_TRANSFER_DESIGN.md +++ b/docs/developer/design/REMESH_FIELD_TRANSFER_DESIGN.md @@ -204,7 +204,7 @@ REMAP. ## 7. The band-aid (already landed) and its relationship to the true fix -`ddt.py` `SemiLagrangian._record_psi_star_from_field_data()` + a guarded call: +`ddt.py` `BackwardNodesSemiLagrangian._copy_tracked_field()` (formerly the parallel-only `_record_psi_star_from_field_data`; since 2026-09-27 used in serial too) + a guarded call: under MPI, when `psi_fn` is a single mesh-variable component on this mesh, the per-step "record current field into `psi_star[0]`" copies the field's nodal data **directly** instead of evaluating it at its own (on-vertex) coords. @@ -255,6 +255,6 @@ Relationship to the true fix: |---|---|---| | `discretisation/discretisation_mesh.py` | `_deform_mesh` @2001, `nuke_coords_and_rebuild` @1757 (per-var refill @1890), `vars` registry @3082 | coords/cache rebuild; scoped-refill change (§5) | | `meshing/smoothing.py` | movers @1570/1861/2639 (per-outer `_deform_mesh`); `smooth_mesh_interior` @2683; `OT_adapt`, `follow_metric` | adapt op owns transfer; mover sweep refill scope | -| `systems/ddt.py` | `SemiLagrangian` (history @119, shift @2035, re-record @2085, band-aid `_record_psi_star_from_field_data`); flavors @98/108/119/146 | history policy = REMAP/ALE; `on_remesh` hook | +| `systems/ddt.py` | `SemiLagrangian` (history @119, shift @2035, re-record @2085, nodal copy `_copy_tracked_field`, formerly `_record_psi_star_from_field_data`); flavors @98/108/119/146 | history policy = REMAP/ALE; `on_remesh` hook | | `systems/solvers.py` | `AdvDiffusionSLCN`, NS, VE (`DuDt`/`DFDt` @997/1354) | use DDt; nothing solver-specific for ALE | | `swarm.py` | `_proxy_stale` @309, `_update` @1019/2129 | proxy REINIT + adapt-triggered staleness (§8) | diff --git a/docs/developer/design/lagged-clone-sl-history.md b/docs/developer/design/lagged-clone-sl-history.md index 04ba8bd2..26b19530 100644 --- a/docs/developer/design/lagged-clone-sl-history.md +++ b/docs/developer/design/lagged-clone-sl-history.md @@ -208,7 +208,7 @@ despite old-frame being active, because the standard `store_result` path **re-records** `psi_star[0]` by evaluating `psi_fn` on the *deformed* mesh at centroid-shifted nodes — injecting boundary-layer interpolation error that grows with `h_max` and then rides the old-geometry sample. Recording the history by a -**direct nodal carry** (reusing the parallel `_record_psi_star_from_field_data` +**direct nodal carry** (reusing the nodal copy `_copy_tracked_field`, formerly `_record_psi_star_from_field_data` path) restores the prototype's exact behaviour. This is the "store primitives, not re-derived values" principle of invariant 3, in miniature. diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 0c55057d..e35354b1 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -580,12 +580,13 @@ def global_evaluate_nd( expr, # DEADLOCK SAFETY — read before editing. Every collective here (allgather, # Allreduce) runs unconditionally on the IDENTICAL global set on every # rank, so all ranks stay in lockstep (n_ext_total is itself a reduced - # value, so the `> 0` guard is taken identically everywhere). The per-rank - # value MUST come from the LOCAL rbf path (rbf=True): the FE interpolation - # path (petsc_interpolate / DMInterpolation) is itself collective and would - # desync here, because each rank classifies the same global set against its - # own domain (different interior-point counts) → hang. Never route the - # fallback value through FE interpolation. + # value, so the `> 0` guard is taken identically everywhere). The + # best-claim value comes from the LOCAL rbf path (rbf=True), evaluated on + # the whole global set. The FE path is used only in the containment round, + # only on meshes whose cell hint is authoritative (no DMLocatePoints, so no + # collective inside it), and every rank calls it on the points it contains + # -- possibly none. Never call the FE path on the global set, and never on + # a mesh that needs DMLocatePoints. # # Serial is left untouched (the serial path above already extrapolates from # the true nearest cell). Escape hatch: GE_LOCAL_FALLBACK=0 restores the @@ -660,22 +661,33 @@ def global_evaluate_nd( expr, # for points no rank contains). Every rank calls evaluate_nd, on # the points it contains -- possibly none -- as the first pass does # on the points it received, so the collectives stay in lockstep. - contains = np.asarray(mesh._robust_owning_cells(all_ext)) >= 0 - my_owner = np.where(contains, comm.rank, comm.size).astype(np.int32) - owner = np.empty(n_ext_total, dtype=np.int32) - comm.Allreduce([my_owner, MPI.INT], [owner, MPI.INT], op=MPI.MIN) - mine = owner == comm.rank - fe_vals, _fe_flag = evaluate_nd( - expr, np.ascontiguousarray(all_ext[mine]), rbf=rbf, evalf=evalf, - verbose=False, simplify=simplify, check_extrapolated=True,) - contrib_fe = np.zeros((n_ext_total,) + expr_shape, dtype=np.float64) - if mine.any(): - contrib_fe[mine] = np.asarray(fe_vals, dtype=np.float64).reshape((-1,) + expr_shape) - fe_val = np.empty_like(contrib_fe) - comm.Allreduce([contrib_fe, MPI.DOUBLE], [fe_val, MPI.DOUBLE], op=MPI.SUM) - contained = owner < comm.size - best_val[contained] = fe_val[contained] - best_flag[contained] = 0 + # + # Only where the cell hint is authoritative: there the FE path + # locates without DMLocatePoints, so a rank holding none of the + # points takes no collective. Elsewhere (warped quads/hexes) the + # rbf value stands, as before. + all_continuous = all( + getattr(varfn.meshvar(), "continuous", True) for varfn in varfns) + if mesh._hint_is_authoritative(all_continuous): + contains = np.asarray(mesh._robust_owning_cells(all_ext)) >= 0 + my_owner = np.where(contains, comm.rank, comm.size).astype(np.int32) + owner = np.empty(n_ext_total, dtype=np.int32) + comm.Allreduce([my_owner, MPI.INT], [owner, MPI.INT], op=MPI.MIN) + mine = owner == comm.rank + fe_vals, _fe_flag = evaluate_nd( + expr, np.ascontiguousarray(all_ext[mine]), rbf=rbf, evalf=evalf, + verbose=False, simplify=simplify, check_extrapolated=True,) + contrib_fe = np.zeros((n_ext_total,) + expr_shape, dtype=np.float64) + if mine.any(): + contrib_fe[mine] = np.asarray(fe_vals, dtype=np.float64).reshape((-1,) + expr_shape) + fe_val = np.empty_like(contrib_fe) + comm.Allreduce([contrib_fe, MPI.DOUBLE], [fe_val, MPI.DOUBLE], op=MPI.SUM) + # a located point whose FE value is NaN is why the fallback + # exists (see above): its finite rbf value stands + use = (owner < comm.size) & np.isfinite( + fe_val.reshape(n_ext_total, -1)).all(axis=1) + best_val[use] = fe_val[use] + best_flag[use] = 0 # Scatter this rank's segment of the global set back to its points. offset = int(counts[:comm.rank].sum()) diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 34043efc..bf3f56c0 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5926,7 +5926,7 @@ def estimate_dt(self, V_fn): # TODO(BUG): with evalf=True evaluate() flags EVERY point as extrapolated, so # at np > 1 this raises the 'not located on this rank' warning for particles # that are all local (seen on the Lagrangian stress history, 2026-09-27). - # See planning file: underworld.md (Bugs section, 2026-09-27) + # See issue #798 vel = uw.function.evaluate(V_fn, self._particle_coordinates.data, evalf=True) # If vel is unit-aware (UnitAwareArray), nondimensionalise it to get diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index d7a034f6..4eba59a2 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -224,10 +224,28 @@ def _history_units(psi_fn, units=None): # A history projection's result is the carried field itself, so it is solved -# to convergence: a mass-matrix solve left at a projection's default tolerance -# (1e-4) carries an error of that size into every step, and the error depends -# on the partition. Converging it costs a few more iterations. -_HISTORY_PROJECTION_TOLERANCE = 1.0e-10 +# to convergence, by CG with a Jacobi preconditioner (the matrix is SPD), and +# from zero: a mass-matrix solve left at a projection's default tolerance (1e-4) +# carries an error of that size into every step, with an aggregating multigrid +# that error depends on the partition, and a solve warmed from the previous +# step's field depends on a state a restart does not restore. +_HISTORY_PROJECTION_TOLERANCE = 1.0e-12 + + +def _owned_rows(var): + """Rows of ``var.data`` this rank owns; a ghost row copies another rank's + (its point has no place in the global section).""" + _is, subdm = var.mesh.dm.createSubDM(var.field_id) + local, global_ = subdm.getLocalSection(), subdm.getGlobalSection() + nc = var.num_components + owned = np.zeros(np.asarray(var.data).shape[0], dtype=bool) + p_start, p_end = local.getChart() + for p in range(p_start, p_end): + ndof = local.getDof(p) + if ndof and global_.getOffset(p) >= 0: + offset = local.getOffset(p) + owned[offset // nc:(offset + ndof) // nc] = True + return owned def _tracked_field(mesh, psi_fn, store): @@ -242,10 +260,10 @@ def _tracked_field(mesh, psi_fn, store): """ if isinstance(psi_fn, uw.discretisation.MeshVariable): field = psi_fn + if field.mesh is not mesh: + return None else: hit = uw.discretisation.meshVariable_lookup_by_symbol(mesh, psi_fn) - if hit is None and isinstance(psi_fn, sympy.MatrixBase) and psi_fn.shape == (1, 1): - hit = uw.discretisation.meshVariable_lookup_by_symbol(mesh, psi_fn[0, 0]) if hit is None: return None field, component = hit @@ -1057,14 +1075,14 @@ def commit_flux_to_history(self, flux, verbose=False): # is a one-shot Galerkin projection and not a fixed-point iteration # (which at a yield kink admits the wrong branch). self._psi_star_projection_solver.smoothing = 0.0 - self._psi_star_projection_solver.solve(verbose=verbose) + self._psi_star_projection_solver.solve(verbose=verbose, zero_init_guess=True) for k, (i, j) in enumerate(self._psi_star_indep_indices): values = np.asarray(self._psi_star_flat_var.data[:, k]) level_0.data[:, level_0._data_layout(i, j)] = values else: self._psi_star_projection_solver.uw_function = flux self._psi_star_projection_solver.smoothing = 0.0 - self._psi_star_projection_solver.solve(verbose=verbose) + self._psi_star_projection_solver.solve(verbose=verbose, zero_init_guess=True) for level in range(self.order - 1, 0, -1): self.psi_star[level].data[...] = ( @@ -1181,12 +1199,12 @@ def _project_nodally(self, expr, name="flux", smoothing=0.0, verbose=False): varsymbol=rf"{{{name}^{{\mathrm{{nodal}}}}_{{{self.instance_number}}}}}") projection = uw.systems.solvers.SNES_MultiComponent_Projection( self.mesh, u_Field=target, n_components=len(columns), verbose=self.verbose) - projection.tolerance = _HISTORY_PROJECTION_TOLERANCE + projection.linear_solver(rtol=_HISTORY_PROJECTION_TOLERANCE) self._nodal_projections[name] = (target, projection) target, projection = self._nodal_projections[name] projection.uw_function = sympy.Matrix([[expr[i, j] for (i, j) in columns]]) projection.smoothing = smoothing - projection.solve(verbose=verbose) + projection.solve(verbose=verbose, zero_init_guess=True) return target def _write_inflow(self, var, coords, rows): @@ -1760,7 +1778,7 @@ def _setup_projections(self): verbose=False, ) self._psi_star_use_multicomponent = True - self._psi_star_projection_solver.tolerance = _HISTORY_PROJECTION_TOLERANCE + self._psi_star_projection_solver.linear_solver(rtol=_HISTORY_PROJECTION_TOLERANCE) self._psi_star_projection_solver.uw_function = self._build_projection_source( self.psi_fn) @@ -1775,16 +1793,12 @@ def update_history_fn(self): evaluation of ``psi_fn``; an L2 projection for expressions that ``evaluate`` cannot handle (e.g. containing derivatives). """ - if self._psi_meshVar is not None: - try: - self.psi_star[0].data[...] = self._psi_meshVar.data[...] - return - except ValueError: - # Sanctioned fallthrough: the tracked variable's nodal - # layout differs from psi_star's (different degree / - # continuity), so the direct copy cannot broadcast — - # evaluate psi_fn at psi_star's own nodes instead. - pass + field = _tracked_field( + self.mesh, self._psi_meshVar if self._psi_meshVar is not None else self.psi_fn, + self.psi_star[0]) + if field is not None: + self.psi_star[0].data[...] = field.data[...] + return try: self.psi_star[0].data[...] = np.asarray(_to_nondim_ndarray(uw.function.evaluate( @@ -1797,7 +1811,7 @@ def update_history_fn(self): # expressions containing derivatives (e.g. flux terms) — # project them onto psi_star[0] instead. self._setup_projections() - self._psi_star_projection_solver.solve() + self._psi_star_projection_solver.solve(zero_init_guess=True) def initialise_history(self): r"""Initialize all history slots to the current value of :math:`\psi`. @@ -3028,7 +3042,7 @@ def __init__( verbose=False, ) self._psi_star_use_multicomponent = True - self._psi_star_projection_solver.tolerance = _HISTORY_PROJECTION_TOLERANCE + self._psi_star_projection_solver.linear_solver(rtol=_HISTORY_PROJECTION_TOLERANCE) # We should find a way to add natural bcs here # (self.Unknowns.u carried as a symbol from solver to solver) @@ -3218,7 +3232,7 @@ def initialise_history(self): # semantics are consistent. self._psi_star_projection_solver.uw_function = self._build_projection_source(self.psi_fn) self._psi_star_projection_solver.smoothing = 0.0 - self._psi_star_projection_solver.solve() + self._psi_star_projection_solver.solve(zero_init_guess=True) if getattr(self, '_psi_star_use_multicomponent', False): # Fan out flat result to tensor psi_star[0] for k, (i, j) in enumerate(self._psi_star_indep_indices): @@ -3524,7 +3538,7 @@ def _record_current_field_into_history( self.psi_fn ) self._psi_star_projection_solver.smoothing = 0.0 - self._psi_star_projection_solver.solve(verbose=verbose) + self._psi_star_projection_solver.solve(verbose=verbose, zero_init_guess=True) # For tensor vtypes the projection writes into the flat (1, Nc) variable, # so we must fan it back out to psi_star[0] — otherwise subsequent @@ -4670,6 +4684,11 @@ def update_post_solve( +# TODO(DESIGN): the forward integration-point history gives an arrival on a +# shared face to one cell (nearest centroid) through this helper, where the +# forward nodal history (_arrivals_by_cell) fits it in every containing cell. +# Integration points almost never arrive on a face, so the two agree in +# practice; one rule for both would remove the difference. def _hand_arrivals_to_owners(X, columns, locate): """Give every arrival to the rank whose cell it landed in. @@ -5177,7 +5196,8 @@ def _object_viewer(self): def _record_current(self): """Snapshot slot 0 <- the current solution and velocity.""" ps = self.psi_snap[0] - field = _tracked_field(self.mesh, self._psi_meshVar or self.psi_fn, ps) + field = _tracked_field( + self.mesh, self._psi_meshVar if self._psi_meshVar is not None else self.psi_fn, ps) if field is not None: ps.data[...] = field.data[...] else: @@ -5264,11 +5284,11 @@ def update_forcing_history(self, forcing_fn=None, evalf=False, verbose=False): self._forcing_projection = uw.systems.solvers.SNES_MultiComponent_Projection( self.mesh, u_Field=self._forcing_flat, n_components=len(columns), verbose=self.verbose) - self._forcing_projection.tolerance = _HISTORY_PROJECTION_TOLERANCE + self._forcing_projection.linear_solver(rtol=_HISTORY_PROJECTION_TOLERANCE) self._forcing_projection.uw_function = sympy.Matrix( [[forcing[i, j] for (i, j) in columns]]) self._forcing_projection.smoothing = 0.0 - self._forcing_projection.solve(verbose=verbose) + self._forcing_projection.solve(verbose=verbose, zero_init_guess=True) points = np.asarray(self.forcing_star.integration_points).reshape(-1, self.mesh.cdim) for column in range(len(columns)): self.forcing_snap.data[:, column] = np.asarray( @@ -5579,10 +5599,10 @@ def _evaluate_at_launch(self, expr): if self._flux_projection is None: self._flux_projection = uw.systems.solvers.SNES_MultiComponent_Projection( self.mesh, u_Field=self._flux_var, n_components=self.num_components) - self._flux_projection.tolerance = _HISTORY_PROJECTION_TOLERANCE + self._flux_projection.linear_solver(rtol=_HISTORY_PROJECTION_TOLERANCE) self._flux_projection.uw_function = sympy.Matrix([[expr[i, j] for (i, j) in self._components]]) self._flux_projection.smoothing = self.flux_smoothing - self._flux_projection.solve() + self._flux_projection.solve(zero_init_guess=True) out = np.empty_like(self._launch_values) for k in range(self.num_components): out[:, k] = _to_nondim_ndarray(uw.function.evaluate(self._flux_var.sym[0, k], self._launch)).reshape(-1) @@ -5830,12 +5850,8 @@ def __init__(self, mesh, psi_fn, V_fn, vtype=VarType.SCALAR, degree=1, varsymbol # cells to the linear fallback; three gains nothing) self._interior = np.asarray(mesh._get_coords_for_basis(self.degree + 2, continuous=False) ).reshape(-1, mesh.cdim) - # A node on a partition seam is on both ranks and is launched once, by its - # owner: writing the rank into the store and reading it back after the - # ghost update leaves every copy holding its owner's rank. - self.psi_star[0].data[:, 0] = uw.mpi.rank - self._owned = np.asarray(self.psi_star[0].data[:, 0]) == uw.mpi.rank - self.psi_star[0].data[:, :] = 0.0 + # a node on a partition seam is on both ranks and is launched once, by its owner + self._owned = _owned_rows(self.psi_star[0]) self._n_relocated = 0 self._n_v = 2 self._init_coefficient_expressions(1, self.theta, with_exp=True) @@ -5873,8 +5889,6 @@ def _values_at(self, expr, X): np.asarray(_to_nondim_ndarray(uw.function.evaluate(expr[i, j], X))).reshape(-1) for (i, j) in self._components]) - _FACE_TOLERANCE = 1.0e-9 - def _arrivals_by_cell(self, X, values): """Every arrival with every cell that contains it, on the rank that owns the cell: one row per (arrival, cell). @@ -5888,11 +5902,17 @@ def _arrivals_by_cell(self, X, values): same on any partition. """ pj = self._projector - tol = self._FACE_TOLERANCE - point, cell, lam = pj.containing_cells(X, tol) + tol = pj.FACE_TOLERANCE + point, cell, lam = pj.containing_cells(X) + inside = lam > tol strict = np.zeros(X.shape[0], dtype=bool) - strict[point[lam > tol]] = True - keep = strict[point] + strict[point[inside]] = True + # a strictly-inside point belongs to its one cell, even where a larger + # neighbour puts it within the tolerance of a face + keep = inside + # TODO(DESIGN): every face point is offered, including those on faces + # inside this rank's partition (no-slip wall nodes, stagnant regions); + # only faces on a partition seam need to travel. offered_X, offered_v = X[~strict], values[~strict] if uw.mpi.size > 1: comm = uw.mpi.comm @@ -5904,7 +5924,7 @@ def _arrivals_by_cell(self, X, values): n_mine = parts_X[uw.mpi.rank].shape[0] else: mine_from, n_mine = 0, offered_X.shape[0] - q, qcell, _ = pj.containing_cells(offered_X, tol) + q, qcell, _ = pj.containing_cells(offered_X) own = (q >= mine_from) & (q < mine_from + n_mine) self._n_relocated = int(np.unique(q[~own]).size) return (np.concatenate([X[point[keep]], offered_X[q]], axis=0), @@ -5925,12 +5945,20 @@ def _reconstruct(self, arrivals, values, source): return np.array(self._project_nodally(self._fit_var.sym, name="fit").data) # ------------------------------------------------------------------ + def _field_at_nodes(self): + """The tracked field at the nodes: its own nodal values when it is a mesh + variable laid out like the store, else evaluated there.""" + field = _tracked_field(self.mesh, self._psi_fn, self.psi_star[0]) + if field is not None: + return np.array(field.data) + return self._values_at(self._psi_fn, self._nodes()) + def initialise_history(self): """Start from the current field at the nodes. A history already placed by :meth:`commit_flux_to_history` is the start, and is kept.""" self.characteristics.initialise_levels(self._n_v) if not self._history_committed: - self.psi_star[0].data[:, :] = self._values_at(self._psi_fn, self._nodes()) + self.psi_star[0].data[:, :] = self._field_at_nodes() self._history_initialised = True def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, **_ignored): @@ -5958,8 +5986,8 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** # the tracked field, or the committed store, which is a polynomial # inside each element, so its interior values are exact source = self._psi_fn if store_result else self.psi_star[0].sym - at_owned = (self._values_at(source, owned) if store_result - else np.asarray(self.psi_star[0].data)[self._owned]) + at_owned = (self._field_at_nodes() if store_result + else np.asarray(self.psi_star[0].data))[self._owned] values = np.vstack([at_owned, self._values_at(source, self._interior)]) key = (_basis_key_of(self.psi_star[0]), "launch") arrivals = np.asarray(trace.departure_points(key, launch, (("first", 0, -float(dt)),), @@ -5968,9 +5996,13 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** if self._inflow_value is not None: # an owned node whose back-trace leaves the domain holds fluid that # entered this step (a ghost copy takes its owner's value) + # (the global-domain test the other flavours use: restoring the + # point moves it; a per-rank "in domain" test calls a partition + # face a boundary) back = 2.0 * owned - arrivals[:owned.shape[0]] + restored = np.asarray(self.mesh.return_coords_to_bounds(back.copy())).reshape(back.shape) entered = np.zeros(nodes.shape[0], dtype=bool) - entered[np.flatnonzero(self._owned)] = ~np.asarray(self.mesh.points_in_domain(back), dtype=bool) + entered[np.flatnonzero(self._owned)] = np.any(np.abs(restored - back) > 0.0, axis=1) inflow = self._values_at(self._inflow_record(), nodes) carried[entered] = inflow[entered] self.psi_star[0].data[:, :] = carried diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 507ae941..0c6cb931 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -463,10 +463,12 @@ def _invalidate_solution_cache(u): # The semi-Lagrangian schemes a solver can build its history with, named # "_" after the arguments of ddt.SemiLagrangian. -_SEMI_LAGRANGIAN_TRANSPORTS = ( - "backward_nodes", "backward_integration_points", - "forward_integration_points", "forward_nodes", -) +_SEMI_LAGRANGIAN_TRANSPORTS = { + "backward_nodes": ("backward", "nodes"), + "backward_integration_points": ("backward", "integration_points"), + "forward_integration_points": ("forward", "integration_points"), + "forward_nodes": ("forward", "nodes"), +} _RENAMED_TRANSPORTS = { "semi_lagrangian": "backward_nodes", "integration_point": "backward_integration_points", @@ -1606,9 +1608,9 @@ def stress_transport(self, value): f"stress_transport={value!r} is now {_RENAMED_TRANSPORTS[value]!r}", FutureWarning, stacklevel=2) value = _RENAMED_TRANSPORTS[value] - if value not in _SEMI_LAGRANGIAN_TRANSPORTS + ("lagrangian", "eulerian"): + if value not in (*_SEMI_LAGRANGIAN_TRANSPORTS, "lagrangian", "eulerian"): raise ValueError( - f"stress_transport must be one of {_SEMI_LAGRANGIAN_TRANSPORTS} or " + f"stress_transport must be one of {tuple(_SEMI_LAGRANGIAN_TRANSPORTS)} or " f"'lagrangian' or 'eulerian', not {value!r}.") if self.Unknowns.DFDt is not None: raise RuntimeError( @@ -1833,7 +1835,7 @@ def _create_stress_history_ddt(self, order=2): **ddt_kwargs, **{k: v for k, v in common.items() if k != "smoothing"}, ) - elif self.stress_transport.startswith("forward_"): + elif _SEMI_LAGRANGIAN_TRANSPORTS.get(self.stress_transport, ("",))[0] == "forward": if ddt_kwargs: raise NotImplementedError( f"{type(cm).__name__} asks its stress history for " @@ -1844,7 +1846,7 @@ def _create_stress_history_ddt(self, order=2): sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), self.u.sym, common["vtype"], - trace="forward", launch=self.stress_transport[len("forward_"):], + trace="forward", launch=_SEMI_LAGRANGIAN_TRANSPORTS[self.stress_transport][1], degree=common["degree"], varsymbol=common["varsymbol"], order=order, units=common["units"], ) @@ -3493,7 +3495,9 @@ class _SmoothingLengthMixin: ``_smoothing_is_dimensional`` flag that lets a dimensional input round-trip as a Pint Quantity while plain-float input round-trips as a plain float. Subclasses keep their own property docstrings as thin - wrappers delegating here. + wrappers delegating here. It also holds :meth:`linear_solver`, the switch + to a CG solve that every projection (scalar, vector, multi-component) can + take. """ def _set_smoothing(self, value): @@ -3540,6 +3544,53 @@ def _set_smoothing_length(self, L): self._smoothing_is_dimensional = is_dim self._smoothing = sympify(L_nd) ** 2 + def linear_solver(self, pc="jacobi", rtol=1.0e-10): + """Switch this projector to a lightweight *linear* (SPD) solve. + + An L2 projection (and the screened-Poisson smoother) is a **linear, + symmetric-positive-definite** problem, so the inherited + ``newtonls / gmres / gamg`` default is unnecessarily heavy — GAMG + setup/repartition dominates cost and memory at MPI scale, which is the + bottleneck for repeated post-processing projections (UW3 issue #156). + This replaces it with ``ksponly + CG + a cheap preconditioner`` — the + right tool for the mass/Helmholtz matrix — and removes the now-unused + GAMG options. + + Opt-in: the default ``SNES_Scalar`` solver stack is unchanged for code + that relies on it. Call this on a projector used purely for output / + post-processing. + + Parameters + ---------- + pc : str, default "jacobi" + Preconditioner. ``"jacobi"`` is fine for a well-conditioned mass + matrix; use ``"bjacobi"`` or ``"icc"`` if CG iteration counts climb + on distorted or high-degree meshes. + rtol : float, default 1e-10 + KSP relative tolerance. + + Returns + ------- + self (so the call can be chained). + """ + self.petsc_options["snes_type"] = "ksponly" + self.petsc_options["ksp_type"] = "cg" + self.petsc_options["pc_type"] = pc + self.petsc_options["ksp_rtol"] = rtol + # GAMG-specific options are now unused; remove them to avoid PETSc + # "unused option" warnings (and any stale AMG configuration). + for _k in ( + "pc_gamg_type", + "pc_gamg_repartition", + "pc_gamg_agg_nsmooths", + "pc_mg_type", + ): + try: + self.petsc_options.delValue(_k) + except Exception: + pass + return self + class SNES_Projection(_SmoothingLengthMixin, SNES_Scalar): r""" @@ -3691,53 +3742,6 @@ def __init__( # Use SymbolicProperty for automatic unwrapping uw_function = SymbolicProperty(matrix_wrap=True, doc="Function to project onto mesh") - def linear_solver(self, pc="jacobi", rtol=1.0e-10): - """Switch this projector to a lightweight *linear* (SPD) solve. - - An L2 projection (and the screened-Poisson smoother) is a **linear, - symmetric-positive-definite** problem, so the inherited - ``newtonls / gmres / gamg`` default is unnecessarily heavy — GAMG - setup/repartition dominates cost and memory at MPI scale, which is the - bottleneck for repeated post-processing projections (UW3 issue #156). - This replaces it with ``ksponly + CG + a cheap preconditioner`` — the - right tool for the mass/Helmholtz matrix — and removes the now-unused - GAMG options. - - Opt-in: the default ``SNES_Scalar`` solver stack is unchanged for code - that relies on it. Call this on a projector used purely for output / - post-processing. - - Parameters - ---------- - pc : str, default "jacobi" - Preconditioner. ``"jacobi"`` is fine for a well-conditioned mass - matrix; use ``"bjacobi"`` or ``"icc"`` if CG iteration counts climb - on distorted or high-degree meshes. - rtol : float, default 1e-10 - KSP relative tolerance. - - Returns - ------- - self (so the call can be chained). - """ - self.petsc_options["snes_type"] = "ksponly" - self.petsc_options["ksp_type"] = "cg" - self.petsc_options["pc_type"] = pc - self.petsc_options["ksp_rtol"] = rtol - # GAMG-specific options are now unused; remove them to avoid PETSc - # "unused option" warnings (and any stale AMG configuration). - for _k in ( - "pc_gamg_type", - "pc_gamg_repartition", - "pc_gamg_agg_nsmooths", - "pc_mg_type", - ): - try: - self.petsc_options.delValue(_k) - except Exception: - pass - return self - @property def smoothing(self): r"""Smoothing coefficient :math:`\alpha` of the screened-Poisson form. @@ -4562,7 +4566,7 @@ def __init__( ## at the various resolutions tested. if transport not in _SEMI_LAGRANGIAN_TRANSPORTS: - raise ValueError(f"transport must be one of {_SEMI_LAGRANGIAN_TRANSPORTS}, " + raise ValueError(f"transport must be one of {tuple(_SEMI_LAGRANGIAN_TRANSPORTS)}, " f"not {transport!r}") if DuDt is not None and transport != "backward_nodes": raise ValueError("transport chooses the DuDt the solver builds; it cannot " @@ -4593,7 +4597,7 @@ def __init__( # so only those that were asked for are passed on asked = {k: v for k, v in (("monotone_mode", monotone_mode), ("old_frame_traceback", old_frame_traceback)) if v} - trace, launch = transport.split("_", 1) + trace, launch = _SEMI_LAGRANGIAN_TRANSPORTS[transport] self.Unknowns.DuDt = uw.systems.ddt.SemiLagrangian( self.mesh, u_Field.sym, diff --git a/src/underworld3/utilities/cell_polynomial_projection.py b/src/underworld3/utilities/cell_polynomial_projection.py index e1456004..28169a65 100644 --- a/src/underworld3/utilities/cell_polynomial_projection.py +++ b/src/underworld3/utilities/cell_polynomial_projection.py @@ -64,10 +64,12 @@ def __init__(self, meshVar): self.ncells = self.detJ.shape[0] probe = tabulate(self.fe, np.zeros((1, self.dim))) self.Nb = probe.shape[1] // self.num_components # scalar basis size - # Reference-cell centroid, mapped: xi_c + 1 = 2 / (dim + 1) on every axis. + # Reference-cell centroid, mapped: xi_c + 1 = 2 / (dim + 1) on every axis + # of a simplex, 1 (the origin of [-1, 1]^dim) of a quadrilateral or hexahedron. J = np.linalg.inv(self.invJ) if self.ncells else self.invJ self.centroids = self.v0 + np.einsum( - "cij,j->ci", J, np.full(self.dim, 2.0 / (self.dim + 1)) + "cij,j->ci", J, + np.full(self.dim, 2.0 / (self.dim + 1) if mesh.isSimplex else 1.0) ) self._check_layout() # Reference coordinates of the cell's dof nodes (the same in every @@ -106,7 +108,9 @@ def reference_coords(self, coords, cells): """Reference coordinates (PETSc's [-1, 1] frame) of points in their cells.""" return np.einsum("cij,cj->ci", self.invJ[cells], coords - self.v0[cells]) - 1.0 - def containing_cells(self, coords, tol=1.0e-9): + FACE_TOLERANCE = 1.0e-9 + + def containing_cells(self, coords, tol=FACE_TOLERANCE): """Every local cell that contains each point. Returns ``(point, cell, lam)``: one row per (point, containing cell) @@ -131,14 +135,29 @@ def containing_cells(self, coords, tol=1.0e-9): point = np.repeat(np.arange(coords.shape[0]), k) cell = near.reshape(-1) xi = self.reference_coords(np.repeat(coords, k, axis=0), cell) + lam = self._reference_distance(xi) + inside = lam >= -tol + point, cell, lam = point[inside], cell[inside], lam[inside] + # a point whose containing cell is not among the nearest centroids (a + # large cell next to small ones): ask the locator + missed = np.setdiff1d(np.arange(coords.shape[0]), point) + if missed.size: + found = np.asarray(self.mesh._robust_owning_cells(coords[missed]), dtype=np.int64) + ok = found >= 0 + if ok.any(): + xm = self.reference_coords(coords[missed[ok]], found[ok]) + point = np.concatenate([point, missed[ok]]) + cell = np.concatenate([cell, found[ok]]) + lam = np.concatenate([lam, self._reference_distance(xm)]) + return point, cell, lam + + def _reference_distance(self, xi): + """How far inside its cell's nearest face a point is, in reference units.""" if self.mesh.isSimplex: # PETSc's reference simplex has vertices at -1 and +1 on each axis lam_axes = 0.5 * (xi + 1.0) - lam = np.minimum(lam_axes.min(axis=1), 1.0 - lam_axes.sum(axis=1)) - else: - lam = 0.5 * (1.0 - np.abs(xi).max(axis=1)) - inside = lam >= -tol - return point[inside], cell[inside], lam[inside] + return np.minimum(lam_axes.min(axis=1), 1.0 - lam_axes.sum(axis=1)) + return 0.5 * (1.0 - np.abs(xi).max(axis=1)) def locate(self, coords): """Owning local cell of each point (-1 when not on this rank) and its reference coordinates.""" @@ -222,6 +241,8 @@ def fit(self, coords, values, nmin=None, patch_nnn=None, old=None, cond_max=1.0e U[held] = np.asarray(old, dtype=np.float64).reshape(self.ncells, self.Nb, nc)[held] thin = np.nonzero(~dense & ~held)[0] self.n_thin = int(thin.shape[0]) + if cell_local and old is None: + raise ValueError("cell_local needs old: the value a cell nothing reached keeps") if cell_local and thin.shape[0] > 0: U[thin] = self._cell_linear_fit(thin, c, xi[ok], psi, cond_max) elif thin.shape[0] > 0 and c.shape[0] > 0: diff --git a/tests/parallel/test_1066_transport_schemes_mpi.py b/tests/parallel/test_1066_transport_schemes_mpi.py index 1e302347..a02bb03f 100644 --- a/tests/parallel/test_1066_transport_schemes_mpi.py +++ b/tests/parallel/test_1066_transport_schemes_mpi.py @@ -8,6 +8,12 @@ run must be the serial ones. At least three ranks: two ranks meet only along one seam, while three or more also meet at points, where an arrival can be handed to either of two other ranks. + +The value histories that apply an inflow value are given one, so the inflow +detection runs across seams too. At np 4 the backward integration-point stress +history has three departure points that the parallel evaluator used to strand +and fill by rbf extrapolation (7.5e-5 before the 2026-09-27 fix): this file is +the regression test for global_evaluate's containment round as well. """ import numpy as np import pytest @@ -51,6 +57,8 @@ def rotating_gaussian(transport, steps=16, dt=np.pi / 32): sympy.exp(-((x - 0.5) ** 2 + y ** 2) / (2 * 0.1 ** 2)), T.coords).reshape(-1) adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=sympy.Matrix([[-y, x]]), order=1, transport=transport) + if adv.DuDt.applies_inflow_value: + adv.DuDt.inflow_value = sympy.Matrix([[0.0]]) adv.constitutive_model = uw.constitutive_models.DiffusionModel adv.constitutive_model.Parameters.diffusivity = 1.0e-3 for wall in ("Left", "Right", "Top", "Bottom"): @@ -63,16 +71,19 @@ def rotating_gaussian(transport, steps=16, dt=np.pi / 32): @pytest.mark.parametrize("transport", [ pytest.param(t, marks=pytest.mark.xfail( - uw.mpi.size >= 4, strict=True, reason="TODO(BUG): the particle history differs from serial by 1.2e-5 " - "at np >= 4 (np 3 matches); the particle values, their ownership and the cell " - "proxy's fit are partition-independent -- cause not yet found")) + uw.mpi.size >= 4, strict=False, reason="#797: the particle history differs from " + "serial by 1.2e-5 at np 4 and 6 (np 3 matches)")) if t == "lagrangian" else t for t in STRESS_XY]) def test_every_stress_history_gives_the_serial_stress_on_every_rank(transport): - _kind, values, _relocated = turned_over_maxwell_box(transport) + uw.reset_default_model() + _kind, values, relocated = turned_over_maxwell_box(transport) + if transport.startswith("forward_"): + assert relocated > 0 # the seams were crossed, so the exchange ran assert np.allclose(values, STRESS_XY[transport], atol=ATOL), (transport, values) @pytest.mark.parametrize("transport", list(ADVECTED_T)) def test_every_value_history_gives_the_serial_field_on_every_rank(transport): + uw.reset_default_model() values = rotating_gaussian(transport) assert np.allclose(values, ADVECTED_T[transport], atol=ATOL), (transport, values) diff --git a/tests/test_1063_stress_history_restart.py b/tests/test_1063_stress_history_restart.py index 2c634071..543ab231 100644 --- a/tests/test_1063_stress_history_restart.py +++ b/tests/test_1063_stress_history_restart.py @@ -42,7 +42,7 @@ def _stress_at_origin(stokes): return float(np.asarray(uw.function.evaluate(stokes.DFDt.psi_star[0].sym[0, 1], np.array([[0.0, 0.0]]))).reshape(-1)[0]) -@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points", "lagrangian"]) +@pytest.mark.parametrize("transport", ["backward_nodes", "backward_integration_points", "forward_integration_points", "forward_nodes", "lagrangian"]) def test_a_restored_history_continues_where_it_left_off(transport): orchestration_model, stokes, dt = _shear_box(transport) for _ in range(6): diff --git a/tests/test_1102_forward_nodes_rotating_gaussian.py b/tests/test_1102_forward_nodes_rotating_gaussian.py index 650d92d8..91c16a45 100644 --- a/tests/test_1102_forward_nodes_rotating_gaussian.py +++ b/tests/test_1102_forward_nodes_rotating_gaussian.py @@ -41,7 +41,7 @@ def _run(mesh, walls): def test_forward_from_nodes_on_the_disc(): uw.reset_default_model() err, peak = _run(_disc(24), ("Upper",)) - assert abs(err - 0.01780) < 0.0018, err # BASELINE (2026-09-27) + assert abs(err - 0.017804) < 1.0e-4, err # BASELINE (2026-09-27) assert abs(peak - EXACT_PEAK) < 0.002, (peak, EXACT_PEAK) @@ -50,5 +50,5 @@ def test_forward_from_nodes_holds_the_box(): mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=2.0 / 24, qdegree=3, regular=False) err, peak = _run(mesh, ("Left", "Right", "Top", "Bottom")) - assert abs(err - 0.0657) < 0.0065, err # BASELINE (2026-09-27) + assert abs(err - 0.065672) < 1.0e-4, err # BASELINE (2026-09-27) assert abs(peak - EXACT_PEAK) < 0.005, (peak, EXACT_PEAK) From 28314aa34b2b2962f1de0278cef9a71512677529 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 27 Sep 2026 09:33:19 -0700 Subject: [PATCH 5/5] The monotone clamp is applied where a point is evaluated; the Navier-Stokes velocity history offers every semi-Lagrangian scheme global_evaluate applied monotone="clamp" on the rank that asked for a point, bounding it by that rank's nearest nodes; for a departure point evaluated on another rank those are the wrong neighbourhood. The bound now runs on the evaluating rank (first pass and containment round), on non-dimensional values like the nodal data it compares with. Serial results are unchanged. #682's level-set reproduction (LeVeque swirl, 128^2 quad, P2, SLCN with clamp): volume at np 1/2/4/8 now 0.0706771442/415/438/412 (np 8 was 1.6% off). NavierStokesSLCN(velocity_transport=...) chooses the velocity history among the four semi-Lagrangian schemes, through one _value_history helper shared with AdvDiffusionSLCN(transport=...). Forward schemes need order=1; the forward integration-point fit is linear and refuses a P2 velocity. tests: test_1103 (lid-driven cavity, Re 100, three schemes against recorded values and each other); test_1066 adds the cavity (equal to serial at np 3/4/6). tests/parallel at np 4: 168 passed, 1 xfailed, 1 failed (test_1063_constrained_freeslip_parallel[ti], fails identically on #795's 05d52480). Serial set: 393 passed. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- docs/developer/subsystems/stress-transport.md | 9 ++ src/underworld3/function/_function.pyx | 7 ++ .../function/functions_unit_system.py | 27 +++-- src/underworld3/systems/solvers.py | 104 +++++++++--------- .../test_1066_transport_schemes_mpi.py | 16 ++- ...t_1103_navier_stokes_velocity_transport.py | 70 ++++++++++++ 6 files changed, 172 insertions(+), 61 deletions(-) create mode 100644 tests/test_1103_navier_stokes_velocity_transport.py diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index 543b8a96..7f505cae 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -86,6 +86,15 @@ make that so, and a new history has to respect them: - a point is never given to one of two cells by a tie-break: a forward arrival on a shared face is fitted in every cell that contains it, and a point the parallel evaluator strands is evaluated by the rank whose cell contains it. +- a monotone bound (`monotone_mode="clamp"`) is applied on the rank that + evaluates the point, among the point's own nodes: applied on the rank that + asked, it bounded a departure point on another rank by the wrong + neighbourhood (#682, 1.6% of a level set's volume at np 8). + +The same holds for the value histories of advection-diffusion +(`AdvDiffusionSLCN(transport=...)`) and for the Navier-Stokes velocity history +(`NavierStokesSLCN(velocity_transport=...)`; the forward integration-point fit +is linear, so it refuses a P2 velocity). The particle history differs from serial by about 1e-5 at np 4 and 6 (np 3 matches); the cause is open. diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index e35354b1..67137c69 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -370,6 +370,7 @@ def global_evaluate_nd( expr, force_l2=False, smoothing=1e-6, local_fallback=True, + limit=None, ): """ @@ -518,6 +519,10 @@ def global_evaluate_nd( expr, # sympy.simplify on every call for any expression holding a mesh variable # (14 of 25 s in a semi-Lagrangian step with a tanh velocity, 2026-09-08). values, extrapolated = evaluate_nd(expr, local_coords, rbf=rbf, evalf=evalf, verbose=verbose, check_extrapolated=True, simplify=simplify,) + # ``limit(coords, values)`` (a monotone bound) runs on the rank that + # evaluated the points, where their neighbourhood is local + if limit is not None and local_coords.shape[0] > 0: + values = limit(local_coords, values) if local_coords.shape[0] > 0: data_container.array[...] = values[...] @@ -677,6 +682,8 @@ def global_evaluate_nd( expr, fe_vals, _fe_flag = evaluate_nd( expr, np.ascontiguousarray(all_ext[mine]), rbf=rbf, evalf=evalf, verbose=False, simplify=simplify, check_extrapolated=True,) + if limit is not None and mine.any(): + fe_vals = limit(np.ascontiguousarray(all_ext[mine]), fe_vals) contrib_fe = np.zeros((n_ext_total,) + expr_shape, dtype=np.float64) if mine.any(): contrib_fe[mine] = np.asarray(fe_vals, dtype=np.float64).reshape((-1,) + expr_shape) diff --git a/src/underworld3/function/functions_unit_system.py b/src/underworld3/function/functions_unit_system.py index a92926bf..7f5930a7 100644 --- a/src/underworld3/function/functions_unit_system.py +++ b/src/underworld3/function/functions_unit_system.py @@ -361,6 +361,7 @@ def _global_evaluate_impl( rbf=None, force_l2=None, local_fallback=True, + limit=None, ): """ Global evaluate with automatic unit-aware results. @@ -542,6 +543,7 @@ def _global_evaluate_impl( force_l2=force_l2_flag, smoothing=smoothing, local_fallback=local_fallback, + limit=limit, ) # Step 2: Re-dimensionalize and wrap with units (GATEWAY PRINCIPLE) @@ -692,12 +694,11 @@ def _apply_monotone_limit( psi_coords_nd = np.asarray(psi_coords_nd.magnitude) # --- kNN neighbour stats from the source nodal data ------------------ - # TODO(parallel): the KDTree is built from rank-local `var.coords_nd`, - # so near a partition seam the neighbour stats bound against a - # truncated neighbourhood. This matches the validated SL behaviour. For - # full parallel correctness the bound should include halo / global DOF - # neighbours -- see the nav-only overlap-clone machinery - # (project_parallel_point_eval_decision) as the hook if hardened. + # The KDTree is built from rank-local `var.coords_nd`. global_evaluate + # calls this on the rank that evaluated each point, so the point's own + # cell and its nodes are local; a point's nearest nodes can still lie in a + # cell the rank does not hold only when it sits against a partition seam. + # TODO(parallel): a halo of the neighbouring cells' nodes would close that. # TODO(units): nbr bounds come from `var.data` (always non-dimensional) # while `value` is dimensional in a units-active run -- a pre-existing # latent mismatch (scaling is inactive in the validated baseline so it @@ -1009,6 +1010,17 @@ def global_evaluate( "check_extrapolated in global_evaluate." ) + # The clamp bounds a value by the nodal data around the point, so it is + # applied where the point is evaluated -- on the rank whose cells hold it, + # among that point's own nodes -- not on the rank that asked, whose nodes + # near a departure point on another rank are the wrong neighbourhood + # (#682: 1.4% of a level set's volume at np 8). The values there are + # non-dimensional, as the nodal data are. + limit = None + if monotone_mode == "clamp": + def limit(local_coords, values): + return _apply_monotone_limit(expr, local_coords, values, "clamp") + result = _global_evaluate_impl( expr, coords=coords, @@ -1024,9 +1036,10 @@ def global_evaluate( rbf=rbf, force_l2=force_l2, local_fallback=local_fallback, + limit=limit, ) - if monotone_mode is None: + if monotone_mode is None or limit is not None: return result return _apply_monotone_limit( diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 0c6cb931..a16ba146 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -469,6 +469,35 @@ def _invalidate_solution_cache(u): "forward_integration_points": ("forward", "integration_points"), "forward_nodes": ("forward", "nodes"), } +def _value_history(transport, mesh, field, V_fn, vtype, order, nodal_options, requested, + **common): + """The semi-Lagrangian history a solver builds for the value of ``field``. + + ``transport`` names the scheme (a key of ``_SEMI_LAGRANGIAN_TRANSPORTS``). + ``nodal_options`` go to the backward nodal scheme only (its boundary + conditions, smoothing, ...); ``requested`` are options the user set, which + go to whichever scheme is chosen and are refused by one that does not + take them; ``common`` go to every scheme. + """ + if transport not in _SEMI_LAGRANGIAN_TRANSPORTS: + raise ValueError(f"transport must be one of {tuple(_SEMI_LAGRANGIAN_TRANSPORTS)}, " + f"not {transport!r}") + if transport == "backward_nodes": + return BackwardNodesSemiLagrangian( + mesh, field.sym, V_fn, vtype=vtype, degree=field.degree, + continuous=field.continuous, varsymbol=field.symbol, order=order, + **nodal_options, **requested, **common) + if not field.continuous: + raise NotImplementedError( + f"transport={transport!r} holds a continuous history; " + "use transport='backward_nodes' for a discontinuous field") + trace, launch = _SEMI_LAGRANGIAN_TRANSPORTS[transport] + return uw.systems.ddt.SemiLagrangian( + mesh, field.sym, V_fn, vtype, trace=trace, launch=launch, + degree=field.degree, varsymbol=field.symbol, order=order, + **requested, **common) + + _RENAMED_TRANSPORTS = { "semi_lagrangian": "backward_nodes", "integration_point": "backward_integration_points", @@ -4565,51 +4594,16 @@ def __init__( ## NB - Smoothing is generally required for stability. 0.0001 is effective ## at the various resolutions tested. - if transport not in _SEMI_LAGRANGIAN_TRANSPORTS: - raise ValueError(f"transport must be one of {tuple(_SEMI_LAGRANGIAN_TRANSPORTS)}, " - f"not {transport!r}") if DuDt is not None and transport != "backward_nodes": raise ValueError("transport chooses the DuDt the solver builds; it cannot " "apply to a DuDt that is supplied") - if DuDt is None and transport == "backward_nodes": - self.Unknowns.DuDt = BackwardNodesSemiLagrangian( - self.mesh, - u_Field.sym, # Symbolic expression - SemiLagrangian evaluates this at each update - self._V_fn, - vtype=uw.VarType.SCALAR, - degree=u_Field.degree, - continuous=u_Field.continuous, - varsymbol=u_Field.symbol, - verbose=verbose, - bcs=self.essential_bcs, - order=1, - smoothing=0.0, - monotone_mode=monotone_mode, - theta=theta, - old_frame_traceback=old_frame_traceback, - ) - elif DuDt is None: - if not u_Field.continuous: - raise NotImplementedError( - f"transport={transport!r} holds a continuous history; " - "use transport='backward_nodes' for a discontinuous field") - # options a scheme does not take are refused by ddt.SemiLagrangian, - # so only those that were asked for are passed on - asked = {k: v for k, v in (("monotone_mode", monotone_mode), - ("old_frame_traceback", old_frame_traceback)) if v} - trace, launch = _SEMI_LAGRANGIAN_TRANSPORTS[transport] - self.Unknowns.DuDt = uw.systems.ddt.SemiLagrangian( - self.mesh, - u_Field.sym, - self._V_fn, - uw.VarType.SCALAR, - trace=trace, - launch=launch, - degree=u_Field.degree, - varsymbol=u_Field.symbol, - order=1, + if DuDt is None: + self.Unknowns.DuDt = _value_history( + transport, self.mesh, u_Field, self._V_fn, uw.VarType.SCALAR, order=1, + nodal_options=dict(verbose=verbose, bcs=self.essential_bcs, smoothing=0.0), + requested={k: v for k, v in (("monotone_mode", monotone_mode), + ("old_frame_traceback", old_frame_traceback)) if v}, theta=theta, - **asked, ) else: @@ -5452,6 +5446,14 @@ class SNES_NavierStokes(SNES_Stokes_SaddlePt): Time derivative operator for velocity. DFDt : BackwardNodesSemiLagrangian or Lagrangian_DDt, optional Time derivative operator for stress. + velocity_transport : str, default="backward_nodes" + The semi-Lagrangian scheme of the internally-constructed velocity + history: ``"backward_nodes"``, ``"backward_integration_points"``, + ``"forward_integration_points"`` or ``"forward_nodes"``, named by the + ``trace`` and ``launch`` arguments of + :func:`~underworld3.systems.ddt.SemiLagrangian`. The forward schemes + carry one level, so they need ``order=1``; the viscoelastic stress + history is chosen separately, by :attr:`stress_transport`. Notes ----- @@ -5498,6 +5500,7 @@ def __init__( verbose: Optional[bool] = False, DuDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, DFDt: Union[BackwardNodesSemiLagrangian, Lagrangian_DDt] = None, + velocity_transport: str = "backward_nodes", ): ## Parent class will set up default values and load u_Field into the solver super().__init__( @@ -5534,20 +5537,15 @@ def __init__( ### sets up DuDt and DFDt ## ._setup_history_terms() - # If DuDt is not provided, then we can build a SLCN version + if DuDt is not None and velocity_transport != "backward_nodes": + raise ValueError("velocity_transport chooses the DuDt the solver builds; it " + "cannot apply to a DuDt that is supplied") if self.Unknowns.DuDt is None: - self.Unknowns.DuDt = uw.systems.ddt.SemiLagrangian( - self.mesh, - self.u.sym, # Symbolic expression - SemiLagrangian evaluates this at each update - self.u.sym, - vtype=uw.VarType.VECTOR, - degree=self.u.degree, - continuous=self.u.continuous, - varsymbol=self.u.symbol, - verbose=self.verbose, - bcs=self.essential_bcs, + self.Unknowns.DuDt = _value_history( + velocity_transport, self.mesh, self.u, self.u.sym, uw.VarType.VECTOR, order=self._order, - smoothing=0.0001, + nodal_options=dict(verbose=self.verbose, bcs=self.essential_bcs, smoothing=0.0001), + requested={}, ) # F (at least for N-S) is a nodal point variable so there is no benefit diff --git a/tests/parallel/test_1066_transport_schemes_mpi.py b/tests/parallel/test_1066_transport_schemes_mpi.py index a02bb03f..cbd2c58e 100644 --- a/tests/parallel/test_1066_transport_schemes_mpi.py +++ b/tests/parallel/test_1066_transport_schemes_mpi.py @@ -4,7 +4,8 @@ two counter-rotating cells between no-slip walls, an unstructured mesh), so the stress is non-uniform and the flow crosses seams of either orientation. Advection: a rotating Gaussian in a square box (the flow crosses every wall), -P2, a quarter turn, with each of the four semi-Lagrangian value histories. The values at two points after the +P2, a quarter turn, with each of the four semi-Lagrangian value histories. +Momentum: the lid-driven cavity of test_1103, with each velocity history. The values at two points after the run must be the serial ones. At least three ranks: two ranks meet only along one seam, while three or more also meet at points, where an arrival can be handed to either of two other ranks. @@ -19,9 +20,15 @@ import pytest import sympy +import sys +from pathlib import Path + import underworld3 as uw from test_1062_forward_stress_history_mpi import turned_over_maxwell_box +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from test_1103_navier_stokes_velocity_transport import CAVITY_U, lid_driven_cavity # noqa: E402 + pytestmark = [pytest.mark.level_2, pytest.mark.tier_b, pytest.mark.mpi(min_size=3), pytest.mark.timeout(1200)] # BASELINES: the serial values at the points of each test (2026-09-27) @@ -87,3 +94,10 @@ def test_every_value_history_gives_the_serial_field_on_every_rank(transport): uw.reset_default_model() values = rotating_gaussian(transport) assert np.allclose(values, ADVECTED_T[transport], atol=ATOL), (transport, values) + + +@pytest.mark.parametrize("transport", list(CAVITY_U)) +def test_every_velocity_history_gives_the_serial_flow_on_every_rank(transport): + uw.reset_default_model() + _kind, values = lid_driven_cavity(transport) + assert np.allclose(values, CAVITY_U[transport], atol=ATOL), (transport, values) diff --git a/tests/test_1103_navier_stokes_velocity_transport.py b/tests/test_1103_navier_stokes_velocity_transport.py new file mode 100644 index 00000000..e3884b5b --- /dev/null +++ b/tests/test_1103_navier_stokes_velocity_transport.py @@ -0,0 +1,70 @@ +"""The Navier-Stokes velocity history carried by each semi-Lagrangian scheme. + +A lid-driven cavity at Reynolds number 100 (lid speed 1, viscosity 0.01, +density 1), first order, ten steps of 0.05 from rest: the momentum is advected, +so the velocity history is doing work. There is no closed form; each scheme is +held to its own recorded value. The forward integration-point history fits a +linear polynomial per cell and refuses the P2 velocity. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +POINTS = np.array([[0.5, 0.75], [0.3, 0.5]]) +# BASELINES: horizontal velocity at POINTS after ten steps (2026-09-27) +CAVITY_U = { + "backward_nodes": (-0.1302979, -0.0553993), + "backward_integration_points": (-0.1300480, -0.0554152), + "forward_nodes": (-0.1300408, -0.0554839), +} + + +def lid_driven_cavity(transport, steps=10, dt=0.05, cell_size=1.0 / 12): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=cell_size, qdegree=3, regular=False) + v = uw.discretisation.MeshVariable("U_cav", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_cav", mesh, 1, degree=1) + ns = uw.systems.NavierStokesSLCN(mesh, velocityField=v, pressureField=p, rho=1.0, order=1, + velocity_transport=transport) + ns.constitutive_model = uw.constitutive_models.ViscousFlowModel + ns.constitutive_model.Parameters.shear_viscosity_0 = 0.01 + ns.add_dirichlet_bc((1.0, 0.0), "Top") + for wall in ("Bottom", "Left", "Right"): + ns.add_dirichlet_bc((0.0, 0.0), wall) + ns.bodyforce = sympy.Matrix([[0.0, 0.0]]) + ns.tolerance = 1.0e-8 + for _ in range(steps): + ns.solve(timestep=dt) + values = np.asarray(uw.function.global_evaluate(v.sym[0], POINTS)).reshape(-1) + return type(ns.DuDt).__name__, values + + +@pytest.mark.parametrize("transport", list(CAVITY_U)) +def test_each_velocity_history_gives_its_recorded_cavity_flow(transport): + uw.reset_default_model() + kind, values = lid_driven_cavity(transport) + assert kind == "".join(w.capitalize() for w in transport.split("_")) + "SemiLagrangian" + assert np.allclose(values, CAVITY_U[transport], atol=1.0e-6), (transport, values) + # the schemes differ by their transport error, a few parts in a thousand here + assert np.allclose(values, CAVITY_U["backward_nodes"], rtol=5.0e-3), (transport, values) + + +def test_the_forward_integration_point_fit_refuses_a_p2_velocity(): + uw.reset_default_model() + with pytest.raises(NotImplementedError, match="degree must be 1"): + lid_driven_cavity("forward_integration_points", steps=0) + + +def test_the_forward_schemes_refuse_order_two(): + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25) + v = uw.discretisation.MeshVariable("U_o2", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_o2", mesh, 1, degree=1) + with pytest.raises(NotImplementedError, match="order must be 1"): + uw.systems.NavierStokesSLCN(mesh, velocityField=v, pressureField=p, order=2, + velocity_transport="forward_nodes")