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/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/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..7f505cae 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 @@ -17,13 +17,21 @@ 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 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 | +| `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 | @@ -60,6 +68,37 @@ 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. +- 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. + ## 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/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/function/_function.pyx b/src/underworld3/function/_function.pyx index 96540e25..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[...] @@ -572,19 +577,21 @@ 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 # 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 @@ -650,6 +657,45 @@ 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. + # + # 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,) + 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) + 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()) seg = slice(offset, offset + ext_coords.shape[0]) 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/swarm.py b/src/underworld3/swarm.py index 48c2d0e2..bf3f56c0 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 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/__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 b21f6a02..4eba59a2 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 = "" @@ -222,6 +223,58 @@ 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, 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): + """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 + if field.mesh is not mesh: + return None + else: + hit = uw.discretisation.meshVariable_lookup_by_symbol(mesh, psi_fn) + 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. @@ -569,7 +622,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 +1042,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 @@ -1022,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[...] = ( @@ -1073,8 +1126,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 +1155,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 @@ -1127,6 +1180,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.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, zero_init_guess=True) + 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 @@ -1698,6 +1778,7 @@ def _setup_projections(self): verbose=False, ) self._psi_star_use_multicomponent = True + 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) @@ -1712,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( @@ -1734,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`. @@ -1917,7 +1994,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. @@ -2486,15 +2563,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). """ @@ -2512,8 +2587,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: @@ -2559,12 +2633,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 +2678,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 +2785,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=[], @@ -2968,6 +3042,7 @@ def __init__( verbose=False, ) self._psi_star_use_multicomponent = True + 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) @@ -3157,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): @@ -3282,40 +3357,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 @@ -3336,14 +3386,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, @@ -3462,83 +3504,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 @@ -3551,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 @@ -3565,7 +3552,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``. @@ -3585,16 +3572,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), @@ -3754,11 +3739,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 @@ -3792,7 +3781,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( @@ -4695,6 +4684,48 @@ 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. + + ``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: @@ -4734,7 +4765,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 +4774,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 +4797,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 @@ -4783,9 +4814,6 @@ class IntegrationPointSemiLagrangian(_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. @@ -4805,22 +4833,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): @@ -4829,10 +4843,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 @@ -4862,7 +4876,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 +4937,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 +5084,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 +5093,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 +5135,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 +5151,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 +5160,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 +5178,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.", @@ -5179,25 +5193,15 @@ 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:`SemiLagrangian._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 if self._psi_meshVar is not None else 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. @@ -5280,10 +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.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( @@ -5351,7 +5356,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 +5408,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 +5417,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 +5436,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 +5461,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 +5539,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 +5582,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): @@ -5590,9 +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.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) @@ -5631,38 +5641,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 @@ -5729,7 +5716,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: @@ -5766,3 +5753,355 @@ 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 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 + 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. + + 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 + ---------- + 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 _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 + 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) + # 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 + 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) + 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(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.""" + 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 _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 = pj.FACE_TOLERANCE + point, cell, lam = pj.containing_cells(X) + inside = lam > tol + strict = np.zeros(X.shape[0], dtype=bool) + 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 + 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) + 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 _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._field_at_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 + (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( + "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) + _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() + 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._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)),), + evalf=evalf, clamp_final=False)) + carried = self._reconstruct(arrivals, values, source) + 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.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 + 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): + """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 = { + ("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..a16ba146 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,55 @@ 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. +_SEMI_LAGRANGIAN_TRANSPORTS = { + "backward_nodes": ("backward", "nodes"), + "backward_integration_points": ("backward", "integration_points"), + "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", + "forward": "forward_integration_points", +} + class _ConstitutiveModelStateMixin: """Single definition of the constitutive-model readiness flag. @@ -504,9 +549,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 +572,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 +1439,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 +1495,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 +1589,58 @@ 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"``, ``"forward_nodes"``, ``"lagrangian"`` + or ``"eulerian"``. + + 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 + 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`). + ``"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 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 (*_SEMI_LAGRANGIAN_TRANSPORTS, "lagrangian", "eulerian"): 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 {tuple(_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 " @@ -1790,51 +1850,53 @@ 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 _SEMI_LAGRANGIAN_TRANSPORTS.get(self.stress_transport, ("",))[0] == "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; " - "use stress_transport='semi_lagrangian' for it.") - self.Unknowns.DFDt = uw.systems.ddt.ForwardSemiLagrangian( + f"{sorted(ddt_kwargs)}, which the forward flavours do not provide; " + "use stress_transport='backward_nodes' for it.") + self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( 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"], + common["vtype"], + trace="forward", launch=_SEMI_LAGRANGIAN_TRANSPORTS[self.stress_transport][1], + 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 +1918,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 +2869,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 +2997,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, @@ -3462,7 +3524,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): @@ -3509,6 +3573,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""" @@ -3660,53 +3771,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. @@ -4382,13 +4446,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 +4475,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 +4486,16 @@ 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. 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 +4556,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,22 +4594,16 @@ def __init__( ## NB - Smoothing is generally required for stability. 0.0001 is effective ## at the various resolutions tested. + 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: - self.Unknowns.DuDt = SemiLagrangian_DDt( - 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, + 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, - old_frame_traceback=old_frame_traceback, ) else: @@ -4554,7 +4623,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 +4654,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 +5116,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 +5150,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,10 +5442,18 @@ 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. + 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 ----- @@ -5421,8 +5498,9 @@ 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, + velocity_transport: str = "backward_nodes", ): ## Parent class will set up default values and load u_Field into the solver super().__init__( @@ -5459,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 @@ -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/src/underworld3/utilities/cell_polynomial_projection.py b/src/underworld3/utilities/cell_polynomial_projection.py index c07b04ad..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,6 +108,57 @@ 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 + 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) + 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) + 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) + 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.""" coords = np.asarray(coords, dtype=np.float64) @@ -118,7 +171,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 +191,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 +236,16 @@ 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 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: # 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 +273,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 db421412..8ef5dcf4 100644 --- a/tests/parallel/test_1062_forward_stress_history_mpi.py +++ b/tests/parallel/test_1062_forward_stress_history_mpi.py @@ -17,11 +17,11 @@ 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]]) -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,12 +48,10 @@ 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 - # 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..cbd2c58e --- /dev/null +++ b/tests/parallel/test_1066_transport_schemes_mpi.py @@ -0,0 +1,103 @@ +"""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. +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. + +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 +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) +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) + 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"): + 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=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): + 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) + + +@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_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..12625e13 100644 --- a/tests/test_1059_stress_transport.py +++ b/tests/test_1059_stress_transport.py @@ -223,9 +223,10 @@ 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", + "forward_nodes": "ForwardNodesSemiLagrangian", "lagrangian": "Lagrangian", "eulerian": "EulerianSUPG", } @@ -247,8 +248,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 +279,41 @@ 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 + + +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.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 +359,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 +403,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 +464,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 +477,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 +491,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 +547,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 +558,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 +583,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 +608,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 +661,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 +706,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 +728,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 +747,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,14 +759,14 @@ 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 + 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_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..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", ["semi_lagrangian", "integration_point", "forward", "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_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 new file mode 100644 index 00000000..91c16a45 --- /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. 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). +""" + +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]]) + 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: + 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.017804) < 1.0e-4, err # BASELINE (2026-09-27) + 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.065672) < 1.0e-4, err # BASELINE (2026-09-27) + assert abs(peak - EXACT_PEAK) < 0.005, (peak, EXACT_PEAK) 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")