Forward-from-nodes history; semi-Lagrangian schemes named by trace and launch; every history equals serial in parallel

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`

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) |

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.
 

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

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 `<trace>_<launch>`, 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

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,15 +459,15 @@ 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
 
     @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

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,
                 ):