diff --git a/docs/developer/index.md b/docs/developer/index.md index 08b180d12..a9bf7073f 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -230,6 +230,7 @@ subsystems/petsc-jacobian-layout subsystems/constitutive-models subsystems/constitutive-models-theory subsystems/constitutive-models-anisotropy +subsystems/stress-transport subsystems/swarm-system subsystems/data-access subsystems/interpolation diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md new file mode 100644 index 000000000..422422927 --- /dev/null +++ b/docs/developer/subsystems/stress-transport.md @@ -0,0 +1,160 @@ +# Viscoelastic stress transport + +A viscoelastic constitutive model carries the stress from one step to the next. The +solver owns the unknowns (velocity, pressure); the stress history is a transported +field managed by one of the `DDt` flavours, chosen with `solver.stress_transport` +before the constitutive model is assigned. This page says what each flavour does, +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.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 +stokes.constitutive_model.Parameters.shear_modulus = G +stokes.constitutive_model.Parameters.solvent_viscosity = eta_s # Oldroyd-B; omit for Maxwell +stokes.constitutive_model.Parameters.dt_elastic = dt +``` + +## The five histories + +| `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 | +| `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 | 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 | + +The nodal and integration-point flavours store the stress after every solve by +the same global L2 projection onto the continuous space. The integration-point +flavour differs only in where it samples that field: the quadrature points along +their own characteristics, rather than the vertices. Nothing is carried at the +points from one step to the next. The forward flavour is the exception among the +mesh-based ones: its launch values persist, and a cell that receives no fit keeps +its previous one. The Lagrangian flavour is a swarm history, not a mesh one: the +solver creates and advects a swarm, evaluates the constitutive flux at the +particles after each solve, and reads it back through a discontinuous cells proxy; +every tensor component is evaluated before any is written, so the proxy is never +half-updated. The same scheme on a user-supplied swarm (a material swarm that +already exists) is :class:`~underworld3.systems.ddt.Lagrangian_Swarm`, passed to +the solver as ``DFDt=``. + +The momentum equation sees the stress only through $\int \sigma : \nabla v$, so +only its per-cell P1 projection matters, and the viscous part is rebuilt from +$\nabla u$ every step while the memory part decays on the relaxation time. That +is why an interpolating history that would be too diffusive for temperature or +velocity is acceptable for stress. Where it fails is where the projection itself +is wrong: sub-cell layers and no-slip walls. + +## 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 +with the current velocity gradient, so over one step it stretches the conformation +by about $(1 + \Delta t\,\dot\gamma)^2$ before relaxation acts. When +$\Delta t\,\dot\gamma$ is of order one that update loses the positive-definiteness +of the conformation $c = \sigma^*/G + I$ in the first step, and no arrangement of +the split recovers it. On the confined cylinder at Courant one on the far-field +mesh the wall shear rate is ten times the far-field one: every history lost the +conformation at the cylinder top in step one, the nodal history then ran away and +the solve hung, and the integration-point history gave out at Wi 0.6. At a step +ten times smaller the integration-point history is admissible everywhere to Wi 0.6, +and at Wi 0.8 the conformation is mildly indefinite on two percent of its points. + +```python +dt = min(courant_dt, stokes.constitutive_model.max_elastic_timestep(safety=0.3)) +``` + +`max_elastic_timestep` is the safety factor over the largest strain-rate magnitude +$\sqrt{2\,\mathbf{D}:\mathbf{D}}$ (the shear rate in simple shear), read from the +continuous projection of the strain rate at the points of the carried stress. That +projection sits a little below the per-cell gradient at a wall, which the safety +factor covers: 0.3 kept the conformation positive on the cylinder. The value is a +physical time when reference scales are active, as `estimate_dt` is. + +## Watching the conformation + +```python +health = stokes.constitutive_model.conformation_min_eigenvalue() +# {'min': 0.29, 'max': 79.3, 'fraction_negative': 0.0, 'where': (0.001, 1.07)} +``` + +An Oldroyd-B or Maxwell stress is $G(c - I)$ with $c$ positive-definite, so the +most compressive eigenvalue of the polymer stress is bounded by $-G$, and the +symmetric part of the momentum tangent stays positive exactly as long as that +holds. A negative `min` is a discretisation defect and it tells the two failure +modes apart: a solve that has lost its preconditioner (the non-symmetric, +co-rotational part of the tangent grows with $|W||\sigma^*|/G$ and is a solver +setting) from a solve that has lost its problem (nothing recovers it). Print this +line every step on a new problem. + +## The recommended configuration + +Integration-point history, the step set by the wall strain rate +(`max_elastic_timestep`), store smoothing at c = 0.07 when the solvent viscosity +is a small fraction of the total, DEVSS off. That combination is characterised on +Waters and King (regular and irregular meshes) and on the confined cylinder, +admissible to Wi 0.6 and mildly indefinite at Wi 0.8. The coefficient 0.023 is +the least that holds the mode on a regular mesh; 0.07 holds it on an irregular +one as well and costs half a percent, so it is the recommended value. The forward flavour is the same scheme with a per-cell fit and interior +launch points, measured at 17 s a step against 28 on the cylinder, parallel by +handing the arrivals that cross a seam to the rank that owns them; it needs its read-back smoothing (`DFDt.flux_smoothing = 0.023 * mesh.cell_size()**2`). Neither transports its memory without a cell-scale mode below Courant +one on a Maxwell element: a version that did not ring turned out not to be +transporting the memory at all. + +## Store smoothing for the integration-point history below Courant one + +The store cycle of the integration-point flavour, sample at the points then +project, is a consistent-mass Galerkin transport of the carried stress. It has no +dissipation at the cell scale, so below Courant one a cell-scale mode grows from +round-off at a rate $\gamma$ set by the elastic feedback: about 2.4 per unit time +on the Maxwell Waters-King start-up, 1.8 with a solvent fraction of 0.2, and not +observed at 0.59 over the run lengths used. The vertex interpolation of the nodal +flavour damps that mode; a discontinuous store (private nodes per cell) frees it +and diverges; DEVSS does not touch it. + +A Laplacian term in the store projection does. It multiplies wavenumber $k$ by +$1/(1 + \alpha k^2)$ once per step, so the mode is held when +$\alpha \approx \gamma\,\Delta t\,(h/\pi)^2$. The coefficient is a field, so the +dose follows the local cell on a graded mesh: + +```python +stokes.DFDt.store_smoothing = 0.07 # alpha = 0.07 * mesh.cell_size()**2 +``` + +Measured on the Waters-King start-up at 1/32 and $\Delta t = 0.01$: $c = 0.023$ +($\alpha = 10^{-5}$) holds the mode for eight time units at 0.1% on the peak, +$c = 0.07$ holds it unconditionally at 0.5%, and $c = 0.23$ (ten times 0.023) +costs 3%. The irregular mesh needs 0.07. Zero, the default, is the plain projection. When the +solvent viscosity is a fair fraction of the total, as on the cylinder benchmark, +the smoothing is unnecessary and costs nothing if left on. + +## DEVSS + +`solver.devss_viscosity = eta_a` adds $2\eta_a(\mathbf{D} - \bar{\mathbf{D}})$ +to the momentum flux, with $\bar{\mathbf{D}}$ the continuous projection of the +strain rate, iterated within the solve until the pair cancels. It is a viscosity +on cell-scale velocity modes only. It rescues the grid (SUPG) history on the +cylinder at Wi 0.4 (80,000 inner solves at the iteration cap without it, none +with), is harmless on the integration-point history, does nothing for the nodal +collapse, and does nothing for the store-cycle mode above. Use it with the grid +history; treat it as optional elsewhere. + +## Waters and King as the discriminating test + +The start-up of plane Poiseuille flow of a Maxwell fluid (Waters and King 1970) +has a reference solution and separates the histories: nodal converges as the +step falls, the integration-point history rings then diverges as the step falls +without smoothing, and the grid history diverges at every step. Note that the +exact solution is independent of $x$, so $u\cdot\nabla\sigma \equiv 0$ there: +the test exercises the store cycle, not the trace-back. A pure Maxwell element is +also the most demanding case for anything explicit in the stress. + +Related: [constitutive models](constitutive-models.md), +[integration-point variables](integration-point-variables.md), +[solvers](solvers.md). Tests: `tests/test_1059_stress_transport.py` (the +histories, the conformation check, the elastic timestep), `test_1060` (store +smoothing below Courant one), `test_1061` (the forward flavour on Waters and +King), `tests/parallel/test_1062` (the forward flavour at np 2), `test_1063` +(save and restore). diff --git a/scripts/test.sh b/scripts/test.sh index 424881dc1..b6caaedda 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -200,8 +200,8 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # CI). The band runs as a whole now, so the named line is gone rather than # running it twice. - # Diffusion / Advection tests - "${PYTEST[@]}" tests/test_1100*py || status=1 + # Diffusion / Advection tests (test_110*: SUPG/SLCN/swarm adv-diff benchmarks) + "${PYTEST[@]}" tests/test_110*py || status=1 "${PYTEST[@]}" tests/test_1110*py tests/test_1120*py || status=1 # Annulus + vector SL "${PYTEST[@]}" tests/test_1450*py || status=1 diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 44424922a..7da160f34 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1732,7 +1732,7 @@ class ViscoElasticPlasticFlowModel(ViscousFlowModel): """ def __init__(self, unknowns, order=1, integrator: str = "bdf", - material_name: str = None): + material_name: str = None, objective_rate: str = "none"): """Construct a viscoelastic-plastic flow model. Parameters @@ -1779,6 +1779,17 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", f"VEP where it produces global runaway). Got order={order}." ) self._integrator = integrator + if objective_rate not in ("none", "upper_convected", "jaumann"): + raise ValueError( + f"objective_rate must be 'none', 'upper_convected' or 'jaumann', got {objective_rate!r}") + # Which objective stress rate the carried stress obeys. "none" transports + # the components as a passive tensor (a linear Maxwell fluid carried by + # the flow); "upper_convected" adds L.sigma + sigma.L^T, which makes + # the element the UCM / Oldroyd-B fluid; "jaumann" adds W.sigma - + # sigma.W (co-rotational). The term is taken on the carried stress and + # the CURRENT velocity gradient, so it is linear in the unknown and + # first order in time. + self._objective_rate = objective_rate # Store material_name before creating expressions (needed by create_unique_symbol) self._material_name = material_name @@ -1875,6 +1886,12 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias): "Shear modulus", units="Pa", ) + solvent_viscosity = api_tools.Parameter( + R"{\eta_s}", + lambda inner_self: 0, + "Solvent viscosity in parallel with the Maxwell element (Oldroyd-B when non-zero)", + units="Pa*s", + ) @property def dt_elastic(inner_self): @@ -2178,6 +2195,9 @@ def stress_2star(self): def E_eff(self): r"""Effective strain rate including elastic-history coupling. + With an objective rate the newest level also contributes the explicit + source :math:`S(\sigma^*)/(2\mu)` (see :meth:`_objective_term`). + For BDF integration: .. math:: @@ -2221,6 +2241,8 @@ def E_eff(self): (1 - phi) * E + (alpha / (2 * eta_raw)) * sigma_star + (phi - alpha) * edot_star + # objective rate: the flux gains alpha * dt * S(sigma*) + + alpha * self.Parameters.dt_elastic * self._objective_term(sigma_star) / (2 * eta_raw) ) return self._E_eff @@ -2229,9 +2251,134 @@ def E_eff(self): bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] for i in range(DDt.order): E += -bdf_cs[i] * DDt.psi_star[i].sym / (2 * mu_dt) + # objective rate on the newest level, explicit and first order: the + # constitutive law is sigma + (eta/mu)(Dsigma/Dt - S) = 2 eta E, so the + # flux gains (eta/mu) S(sigma*) and E_eff gains S/(2 mu) with weight + # ONE whatever the BDF order (the BDF weights belong to the time + # derivative, not to the source; -c1 = 2 at order 2 would double it). + E += self._objective_term(DDt.psi_star[0].sym) / (2 * self.Parameters.shear_modulus) self._E_eff.sym = E return self._E_eff + def _objective_term(self, sigma): + r"""The stretching / rotation source of the chosen objective rate, + :math:`S(\sigma)`, on the carried stress with the current velocity + gradient :math:`L_{ij} = \partial u_i / \partial x_j`: zero for + ``"none"``, :math:`L\sigma + \sigma L^T` for ``"upper_convected"``, + :math:`W\sigma - \sigma W` with :math:`W = (L - L^T)/2` for + ``"jaumann"``. Symmetric whenever ``sigma`` is.""" + dim = self.Unknowns.u.mesh.dim + if self._objective_rate == "none": + return sympy.zeros(dim, dim) + sigma = sympy.Matrix(sigma) + L = sympy.Matrix(self.Unknowns.u.sym).jacobian(self.Unknowns.u.mesh.X) + if self._objective_rate == "upper_convected": + return L * sigma + sigma * L.T + W = (L - L.T) / 2 + return W * sigma - sigma * W + + # ----- Health of the carried stress (#768) ----- + + def _carried_stress(self): + """The carried polymer stress as tensors at its own points, non-dimensional. + Only the trace-back histories carry one tensor per point; a particle + history does not, and is refused rather than misread.""" + DDt = self.Unknowns.DFDt + if DDt is None: + raise RuntimeError("the model has no stress history yet: assign it to a solver and solve once") + if not hasattr(DDt, "carried_tensors"): + raise NotImplementedError(f"{type(DDt).__name__} does not expose the carried stress as " + "one tensor per point; the conformation check reads a " + "per-point tensor history (the trace-back flavours)") + return DDt.carried_tensors() + + def max_elastic_timestep(self, safety: float = 0.3) -> float: + r"""The largest step the explicit stretching term tolerates: + ``safety`` divided by the largest strain rate anywhere. + + The objective-rate source :math:`L\sigma^* + \sigma^* L^T` is taken on + the carried stress with the current gradient, so over one step it + stretches the conformation by about :math:`(1 + \Delta t\,\dot\gamma)^2` + before relaxation acts. With :math:`\Delta t\,\dot\gamma` of order one + that update loses the positive-definiteness of the conformation in the + first step (measured on the confined cylinder at Courant one on the + far-field mesh, where the wall shear rate is ten times the far-field + one), and no operator split recovers it. A run should take + ``dt = min(courant_dt, model.max_elastic_timestep())``; 0.3 keeps the + conformation positive on the cylinder to Wi 0.6. Like ``estimate_dt``, + the value is a physical time when reference scales are active and a + plain number otherwise, so the two can be compared directly. + + The strain-rate measure is :math:`\dot\gamma = \sqrt{2\,\mathbf{D}:\mathbf{D}}`, + the shear rate in simple shear, read from the continuous projection of + the strain rate (what ``evaluate`` returns for a gradient) at the points + of the carried stress. That projection sits a little below the per-cell + gradient at a wall, which the safety factor covers. Reduced over + ranks. ``inf`` when there is no objective rate (nothing stretches) or + no flow; a velocity that is not finite raises rather than returning + a silent ``nan``. + """ + if self._objective_rate == "none": + return float("inf") + E = sympy.Matrix(self.Unknowns.E) + rate = sympy.sqrt(2 * (E.T * E).trace()) + _, points = self._carried_stress() + # evaluate is collective: every rank calls it, with its own (possibly empty) points + from underworld3.systems.ddt import _to_nondim_ndarray + values = np.asarray(_to_nondim_ndarray(uw.function.evaluate(rate, points))).reshape(-1) + local = float(np.abs(values).max()) if values.size else 0.0 + finite = bool(np.all(np.isfinite(values))) if values.size else True + peak = float(uw.mpi.comm.allreduce(local, op=uw.MPI.MAX)) + finite = bool(uw.mpi.comm.allreduce(finite, op=uw.MPI.LAND)) + if not finite: + raise RuntimeError("max_elastic_timestep: the strain rate is not finite") + if peak == 0.0: + return float("inf") + dt = float(safety) / peak + try: + return uw.dimensionalise(dt, {"[time]": 1}) + except Exception: + return dt + + def conformation_min_eigenvalue(self): + r"""Smallest eigenvalue of the conformation :math:`c = \sigma^*/G + I` on + the carried polymer stress, with where it is and how much of the field + is below zero. + + An Oldroyd-B or Maxwell stress is :math:`G(c - I)` with :math:`c` + positive-definite, so the most compressive eigenvalue of the polymer + stress is bounded by :math:`-G`. The symmetric part of the momentum + tangent stays positive exactly as long as that holds. A negative value + here is a discretisation defect (the step against the local strain rate, + or a wall-layer excess of the nodal history), and it separates a solve + that has lost its preconditioner from one that has lost its problem: + the first is a solver setting, the second is not. + + Returns a dict: ``min`` and ``max`` (reduced over ranks), + ``fraction_negative``, and ``where``, the coordinates of the minimum on + any rank whose local minimum is the global one (``None`` on the others). + Non-dimensional throughout. Needs a per-point stress history. + """ + tau, points = self._carried_stress() + dim = tau.shape[-1] + from underworld3.systems.ddt import _to_nondim_ndarray + G = np.asarray(_to_nondim_ndarray( + uw.function.evaluate(self.Parameters.shear_modulus.sym, points))).reshape(-1) + c = tau / G[:, None, None] + np.eye(dim)[None, :, :] + ev = np.linalg.eigvalsh(c) # ascending per point + lo = ev[:, 0] + n_local = lo.size + local_min = float(lo.min()) if n_local else float("inf") + local_max = float(ev[:, -1].max()) if n_local else float("-inf") + n_neg = int((lo < 0.0).sum()) + comm = uw.mpi.comm + gmin = float(comm.allreduce(local_min, op=uw.MPI.MIN)) + gmax = float(comm.allreduce(local_max, op=uw.MPI.MAX)) + n_all = int(comm.allreduce(n_local, op=uw.MPI.SUM)) + n_neg_all = int(comm.allreduce(n_neg, op=uw.MPI.SUM)) + where = tuple(float(x) for x in points[int(lo.argmin())]) if (n_local and local_min == gmin) else None + return {"min": gmin, "max": gmax, "fraction_negative": (n_neg_all / n_all) if n_all else 0.0, "where": where} + @property def E_eff_inv_II(self): r"""Second invariant of effective strain rate: :math:`\dot{\varepsilon}_{II} = \sqrt{\frac{1}{2}\dot{\varepsilon}_{ij}\dot{\varepsilon}_{ij}}`.""" @@ -2454,8 +2601,9 @@ def stress(self): wrapped effective viscosity (``ve_effective_viscosity`` for BDF, raw ``η`` for ETD-2 since the time factor is in ``E_eff``). """ + solvent = 2 * self.Parameters.solvent_viscosity * self.Unknowns.E if not self.is_elastic or self.Unknowns.DFDt is None: - return 2 * self.viscosity * self.grad_u + return 2 * self.viscosity * self.grad_u + solvent # ETD-1 (order=1) uses the same E_eff machinery but with φ=α, so # forcing_star is not required (the (φ-α)·ε̇* term zeros out). @@ -2466,12 +2614,23 @@ def stress(self): and self.Unknowns.DFDt.forcing_star is None ): raise RuntimeError( - "integrator='etd' requires a SemiLagrangian DDt with " - "with_forcing_history=True. The auto-DDt creation path " + "integrator='etd' at order 2 needs a stress history with a " + "forcing slot (with_forcing_history=True: the nodal or the " + "integration-point flavour). The auto-DDt creation path " "reads stress_history_ddt_kwargs — re-create the solver/" "model so the kwargs propagate." ) + return 2 * self.viscosity * self.E_eff.sym + solvent + + @property + def history_flux(self): + """The stress the history carries: the Maxwell element's own stress, + without the solvent's Newtonian part. The momentum flux is + :attr:`flux` = this + 2 eta_s E; committing the total would put the + solvent stress into the memory and feed it back a step later.""" + if not self.is_elastic or self.Unknowns.DFDt is None: + return 2 * self.viscosity * self.grad_u return 2 * self.viscosity * self.E_eff.sym # def eff_edot(self): diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index d9050c1e2..e51cd1d89 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3523,7 +3523,8 @@ def return_coords_to_bounds(self): * **analytic** (default) — the closure installed by the mesh constructor (radial on an annulus/sphere, face clamps on a box). Cheap and exact **while the boundary keeps the shape it was written for**. - * **general facet restore** (once :meth:`deform` has moved the geometry) — the + * **general facet restore** (once :meth:`deform` has moved the geometry, and for + any mesh with no analytic closure at all — every mesh read from a file) — the nearest point on the mesh's CURRENT boundary facets, with an outward-normal side test (:meth:`_facet_return_coords_to_bounds`). @@ -3538,7 +3539,14 @@ def return_coords_to_bounds(self): Assigning to this attribute overrides both (the setter replaces the analytic closure and is honoured until the geometry deforms). """ - if getattr(self, "_geometry_deformed", False): + if (getattr(self, "_geometry_deformed", False) + or self._analytic_return_coords_to_bounds is None): + # No analytic closure means a mesh built from a file (every gmsh mesh), which + # is where the benchmark geometries live. Without this they returned None and + # a trace-back foot leaving through an inlet was never restored: it fell + # through to the evaluator's distance-weighted fallback, which is both slow + # and wrong. The general restore already handles that case, and returns + # interior points untouched, so it is the right fallback rather than nothing. return self._facet_return_coords_to_bounds return self._analytic_return_coords_to_bounds @@ -3558,7 +3566,13 @@ def _boundary_facet_geometry(self): return cache[1], cache[2] cdim = self.cdim - facets, opp = _boundary_facets(self, cdim) + # The DM's boundary label separates domain faces from partition faces; + # the topological search (a facet in exactly one local cell) cannot, and + # on a distributed mesh would restore every foot crossing a rank seam to + # that seam. Only a mesh without the label falls back to the search. + facets, opp = self._labelled_boundary_facets() + if facets is None: + facets, opp = _boundary_facets(self, cdim) if facets is None: # non-simplicial: no general restore self._bnd_restore_cache = (stamp, None, None) return None, None @@ -3576,6 +3590,31 @@ def _boundary_facet_geometry(self): self._bnd_restore_cache = (stamp, fpts, n) return fpts, n + def _labelled_boundary_facets(self): + """Boundary facets and the opposite cell vertex from the DM's own + ``All_Boundaries`` label: ``(facets, opp)`` as :func:`_boundary_facets` + returns them (vertex indices in point order), or ``(None, None)`` when + the label is absent or the mesh is not simplicial.""" + dm = self.dm + if not dm.hasLabel("All_Boundaries"): + return None, None + label = dm.getLabel("All_Boundaries") + d = self.dim + f0, f1 = dm.getHeightStratum(1) + v0, v1 = dm.getDepthStratum(0) + facets, opp = [], [] + for f in range(f0, f1): + if dm.getSupportSize(f) != 1 or label.getValue(f) == -1: + continue + fverts = [p for p in dm.getTransitiveClosure(f)[0] if v0 <= p < v1] + cverts = [p for p in dm.getTransitiveClosure(dm.getSupport(f)[0])[0] if v0 <= p < v1] + if len(fverts) != d or len(cverts) != d + 1: + return None, None + facets.append([p - v0 for p in fverts]) + opp.append([p for p in cverts if p not in fverts][0] - v0) + return (numpy.array(facets, dtype=int).reshape(-1, d), + numpy.array(opp, dtype=int).reshape(-1)) + def _facet_return_coords_to_bounds(self, coords): """General restore: snap points lying OUTSIDE the current boundary to just inside the nearest boundary facet. Interior points are returned untouched (the @@ -3599,8 +3638,10 @@ def _facet_return_coords_to_bounds(self, coords): owner = numpy.asarray(tree.query(numpy.ascontiguousarray(closest), 1)[1]).flatten() nvec = nrm[owner] outside = numpy.einsum("ij,ij->i", pts - closest, nvec) > 0.0 + # get_min_radius is collective: read it on every rank, not only on the + # ranks that have a point to restore + eps = 1.0e-3 * float(self.get_min_radius()) if numpy.any(outside): - eps = 1.0e-3 * float(self.get_min_radius()) pts[outside] = closest[outside] - eps * nvec[outside] return pts diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index c47a33832..5ba37b93a 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -69,6 +69,7 @@ # These are now implemented the same way using the ddt module from .solvers import SNES_AdvectionDiffusion as AdvDiffusionSLCN +from .solvers import SNES_AdvectionDiffusion_Swarm as AdvDiffusionSwarm # The generic names are the composing solvers: the transport (assembled SUPG # advection, or a semi-Lagrangian history) is the DDt manager they hold. from .advection_diffusion_eulerian import SNES_AdvectionDiffusion_Composed as AdvDiffusion @@ -94,6 +95,7 @@ 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 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 cc796eecd..047bb3b83 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -144,6 +144,25 @@ class DDtSemiLagrangianState(_DDtCoreState): with_forcing_history: bool = False +@dataclass +class DDtIntegrationPointState(_DDtCoreState): + """Snapshot of an :class:`IntegrationPointSemiLagrangian` 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) + psi_snap_var_names: list[str] = field(default_factory=list) + with_forcing_history: bool = False + history_committed: bool = False + + +@dataclass +class DDtForwardState(_DDtCoreState): + """Snapshot of a :class:`ForwardSemiLagrangian` 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 = "" + + @dataclass class DDtLagrangianState(_DDtCoreState): """Snapshot of a :class:`Lagrangian` DDt instance. @@ -552,9 +571,8 @@ class _DDtBase(uw_object): - **History symbols**: Symbolic stores raw sympy matrices in ``psi_star``; the storage-backed flavors store variables and contribute ``.sym`` (see :meth:`_history_syms`). - - **ETD-2 exp coefficients** exist only on the flavors used by the - Maxwell / viscoelastic relaxation path (``with_exp=True``: - Symbolic, Eulerian, SemiLagrangian). + - **ETD-2 exp coefficients** exist on every flavour that can carry a + viscoelastic stress (``with_exp=True``: all but the particle flavours). """ @classmethod @@ -620,6 +638,18 @@ def _init_history_tracking(self, order): # History tracking: deferred initialization and effective order self._history_initialised = False self._n_solves_completed = 0 + # Snapshot substitution in the projection's source (see + # enable_source_snapshot): every flavour that projects a flux into its + # own history needs it, not only the semi-Lagrangian one. + self._psi_snapshot_enabled = False + self._psi_snapshot = None + # Set by commit_flux_to_history: the levels are already placed for this + # step, so a post-solve must not shift or re-record them again. + self._history_committed = False + # What the transported quantity is worth where the flow enters. Held on + # the base so a driver can set it whatever flavour it holds; see the + # inflow_value property for which flavours act on it. + self._inflow_value = None self._dt = None # current timestep (set by solver or update_pre_solve) self._dt_history = [None] * order # previous timesteps for variable-dt BDF @@ -668,7 +698,7 @@ def _init_coefficient_expressions(self, order, theta, with_exp): Also create the ETD-2 ``[α, φ]`` coefficients used by Maxwell-relaxation integration; values are pushed via PetscDSSetConstants every step in ``update_exp_coefficients`` - (Symbolic, Eulerian, SemiLagrangian only). + (every flavour but the particle ones, #739). """ self._bdf_coeffs = _create_coefficients(order, r"c^{\mathrm{BDF}}", self.instance_number) self._am_coeffs = _create_coefficients(order, r"a^{\mathrm{AM}}", self.instance_number) @@ -680,6 +710,28 @@ def _init_coefficient_expressions(self, order, theta, with_exp): if with_exp: _update_exp_values(self._exp_coeffs, None, None) + def update_exp_coefficients(self, dt, tau_eff): + r"""Set the exponential (ETD) coefficients for this step. + + ``self._exp_coeffs[0].sym = α = exp(-Δt/τ_eff)`` and + ``self._exp_coeffs[1].sym = φ = (1-α)/(Δt/τ_eff)``, with τ_eff the + Maxwell relaxation time :math:`\eta_\mathrm{eff}/\mu`. Called by the + constitutive model, which owns τ_eff, before each solve -- peer to the + BDF/AM coefficient updates that ``update_pre_solve`` makes itself. + Every flavour that can carry a Maxwell stress allocates these + coefficients (#739); the particle flavours do not. + """ + _update_exp_values(self._exp_coeffs, dt, tau_eff) + + @property + def _exp_alpha(self): + """The ETD ``α`` coefficient UWexpression.""" + return self._exp_coeffs[0] + + @property + def _exp_phi(self): + """The ETD ``φ`` coefficient UWexpression.""" + return self._exp_coeffs[1] def _note_history_shift(self, dt, **detail): """Tell the model's open step that this history advanced. @@ -967,6 +1019,151 @@ def initiate_history_fn(self): """Deprecated: use ``initialise_history`` instead.""" self.initialise_history() + def _build_projection_source(self, source_fn): + """Construct the row matrix used as the projection's ``uw_function``. + + Applies snapshot substitution (psi_star[0] → snap) when enabled. + Used by both ``psi_fn.setter`` and the ``initialise_history`` + fallback path so substitution semantics are consistent. + """ + if getattr(self, '_psi_star_use_multicomponent', False): + indep = self._psi_star_indep_indices + row = sympy.Matrix([[source_fn[i, j] for (i, j) in indep]]) + if self._psi_snapshot_enabled and self._psi_snapshot is not None: + ps0 = self.psi_star[0] + psi_snapshot = self._psi_snapshot + substitutions = { + ps0.sym[i, j]: psi_snapshot.sym[i, j] + for i in range(self.mesh.dim) + for j in range(self.mesh.dim) + } + row = row.subs(substitutions) + return row + else: + # Scalar / vector path: psi_star[0] is a scalar/vector field. If + # snapshot is needed for these vtypes, extend here similarly. + return source_fn + + def enable_source_snapshot(self): + """Enable snapshot substitution in the projection's source field. + + Call this once when the source expression (``psi_fn``) references + ``psi_star[0]`` itself — without it the projection's residual + ``(target − flux(psi_star[0]))·weight`` is implicit in the target + because target and source share the same data field. With Min-mode + plasticity at the yield kink, the implicit projection admits two + fixed points (elastic and yield branches); under timestep change the + iteration drifts to the elastic-branch fixed point and σ violates + the yield surface. + + The snapshot is a separate mesh variable matching ``psi_star[0]``'s + shape/vtype/degree. Each call to ``update_pre_solve`` copies + ``psi_star[0].array → psi_snapshot.array``, freezing the source's + input for the upcoming projection. Substitution makes the + projection's compiled C code read from ``psi_snapshot.array`` + instead of ``psi_star[0].array`` — there's no recompile per step, + just a memcpy. + + Idempotent: safe to call more than once. + """ + if not getattr(self, '_psi_star_use_multicomponent', False): + # Currently only wired for tensor projections (the case that + # exposed the bug). Scalar/vector extension is straightforward + # if needed later. + return + + if self._psi_snapshot is None: + ps0 = self.psi_star[0] + # NOTE: this currently registers a persistent MeshVariable in the + # mesh DM, which is overkill for a transient buffer that's only + # read by this DDt's projection. A future improvement would be + # a transient/scratch-variable mechanism (likely backed by + # PETSc's auxiliary Vec machinery — already used elsewhere in + # the codebase via DMSetAuxiliaryVec_UW) so the snapshot doesn't + # accumulate in the DM across DDt creations. See: + # docs/developer/ai-notes/historical-notes.md for the + # variable-deletion limitation context. + self._psi_snapshot = uw.discretisation.MeshVariable( + f"psi_snapshot_{self.instance_number}", + self.mesh, + ps0.shape, + vtype=ps0.vtype, + degree=ps0.degree, + continuous=ps0.continuous, + ) + # Initialise psi_snapshot's data to current psi_star[0]'s data + # so the source evaluates consistently before the first refresh. + self._psi_snapshot.data[...] = ps0.data[...] + + self._psi_snapshot_enabled = True + + # Re-run the psi_fn setter so the substitution is applied to the + # currently-installed projection source. + self.psi_fn = self._psi_fn + + def _refresh_source_snapshot(self): + """Freeze the projection's input for this step (a memcpy, no recompile). + + Routes through ``.data`` rather than ``.array`` to skip unit conversion + (both variables are non-dimensional) while keeping the callback sync + that pushes values into the underlying PETSc local vector. + """ + if self._psi_snapshot_enabled and self._psi_snapshot is not None: + self._psi_snapshot.data[...] = self.psi_star[0].data[...] + + #: A history that places the new flux on its own storage during + #: ``update_post_solve`` (the particle flavours evaluate it at their + #: particles) rather than through :meth:`commit_flux_to_history`. + commits_flux_in_post_solve = False + + #: The forcing (strain-rate) history the second-order exponential + #: integrator reads; only :class:`SemiLagrangian` allocates one, on + #: request. ``None`` means the integrator runs at first order (#739). + forcing_star = None + + def commit_flux_to_history(self, flux, verbose=False): + r"""Project ``flux`` into ``psi_star[0]`` and shift the history levels. + + What a solver does after solving with a flux history (a viscoelastic + stress, say): the flux the constitutive model has just formed becomes + the new level 0, and the level that was 0 -- already carried to the new + configuration by ``update_pre_solve``, whether by a trace-back or by an + assembled transport -- becomes level 1. The shift is the history + manager's business, so every flavour does it the same way and a solver + does not need to know which one it holds. + + ``flux`` is the expression to project (the constitutive model's flux). + It is ignored on the multi-component path, where the projection's source + was compiled once and is refreshed through the snapshot machinery (see + :meth:`enable_source_snapshot`). + """ + if not hasattr(self, "_psi_star_projection_solver"): + self._setup_projections() + + transported = np.copy(self.psi_star[0].array[...]) + + if getattr(self, "_psi_star_use_multicomponent", False): + # The snapshot machinery has frozen the projection's input, so this + # 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) + for k, (i, j) in enumerate(self._psi_star_indep_indices): + values = self._psi_star_flat_var.array[:, 0, k] + self.psi_star[0].array[:, i, j] = values + if i != j: + self.psi_star[0].array[:, j, i] = 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) + + for level in range(self.order - 1, 0, -1): + self.psi_star[level].array[...] = ( + transported if level == 1 else self.psi_star[level - 1].array[...]) + + self._history_committed = True + # ----- The transport contract ----- # # A solver that owns an unknown composes its residual from these terms @@ -985,6 +1182,57 @@ def integrator(self) -> str: """``"am"`` (the theta rule on the spatial terms) at order 1, ``"bdf"`` above.""" return "am" if self.order == 1 else "bdf" + @property + def inflow_value(self): + r"""What enters the domain where the flow comes in, or ``None``. + + A transported quantity needs data wherever the flow enters, and nowhere + else. Which parts of the boundary those are is not fixed: on a shedding + wake the outflow boundary carries reversed flow that migrates along it, + so the condition is applied by the sign of :math:`\mathbf{u}\cdot\mathbf{n}` + rather than by naming a boundary. Left ``None`` the transport is + unconstrained at an inflow, and whatever the solve produces there is + carried into the domain: measured on the viscoelastic cylinder, that is + what destroys the run (the stress maximum leaves the cylinder for the + outlet as soon as the wake reverses through it). + + Set it to an expression of the unknown's shape -- for a stress history, + 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` + fills the uncovered share of an inflow cell with it. The nodal + trace-back and the particle flavours do not use it: a departure point + or a particle that lands outside the domain is restored to the + boundary and takes the transported field's value THERE, which + constrains the inflow but is not the value set. Setting a value on + such a flavour says so once rather than dropping it in silence (#733). + """ + return self._inflow_value + + @inflow_value.setter + def inflow_value(self, value): + if value is not None: + value = sympy.Matrix(value) + if value.shape != self._unknown_shape(): + raise ValueError( + f"inflow_value has shape {value.shape}, but the transported " + f"quantity is {self._unknown_shape()}.") + if not self.applies_inflow_value: + warnings.warn( + f"{type(self).__name__} does not apply inflow_value: it " + "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 and ForwardSemiLagrangian " + "apply it (#733).", + stacklevel=2) + self._inflow_value = value + + #: Whether this flavour compiles :attr:`inflow_value` into its transport. + applies_inflow_value = False + def _unknown_shape(self): """Shape of the unknown as a matrix (``Symbolic`` stores ``_shape`` as data).""" psi = self.psi_fn @@ -1296,18 +1544,6 @@ def _history_syms(self): """Symbolic stores raw sympy matrices in ``psi_star`` — return them as-is.""" return list(self.psi_star) - def update_exp_coefficients(self, dt, tau_eff): - r"""Update the ETD-2 (exponential) coefficient values for this step. - - Sets ``self._exp_coeffs[0].sym = α`` and ``self._exp_coeffs[1].sym = φ`` - from current ``dt`` and ``tau_eff`` (Maxwell relaxation time - :math:`\tau = \eta_\mathrm{eff}/\mu`). Called by the constitutive - model (which owns τ_eff) before each solve, peer to the BDF/AM - coefficient updates that happen automatically in - ``update_pre_solve``. - """ - _update_exp_values(self._exp_coeffs, dt, tau_eff) - class Eulerian(_DDtBase): r""" @@ -1478,14 +1714,16 @@ def psi_fn(self): @psi_fn.setter def psi_fn(self, new_fn): - """Set the tracked expression.""" + """Set the tracked expression, and the source of any live projection.""" self._psi_fn = new_fn - # self._psi_star_projection_solver.uw_function = self.psi_fn - return + if getattr(self, "_psi_star_projection_solver", None) is not None: + self._psi_star_projection_solver.uw_function = self._build_projection_source(new_fn) def _setup_projections(self): - """Initialize projection solvers for history updates.""" + """Initialize projection solvers for history updates (once).""" + if getattr(self, "_psi_star_projection_solver", None) is not None: + return ### using this to store terms that can't be evaluated (e.g. derivatives) # The projection operator for mapping derivative values to the mesh - needs to be different for each variable type, unfortunately ... if self.vtype == uw.VarType.SCALAR: @@ -1539,13 +1777,8 @@ def _setup_projections(self): ) self._psi_star_use_multicomponent = True - if getattr(self, '_psi_star_use_multicomponent', False): - # Flatten tensor to (1, Nc) row for multicomponent solver - indep = self._psi_star_indep_indices - row = sympy.Matrix([[self.psi_fn[i, j] for (i, j) in indep]]) - self._psi_star_projection_solver.uw_function = row - else: - self._psi_star_projection_solver.uw_function = self.psi_fn + self._psi_star_projection_solver.uw_function = self._build_projection_source( + self.psi_fn) self._psi_star_projection_solver.bcs = self.bcs self._psi_star_projection_solver.smoothing = self.smoothing @@ -1654,9 +1887,14 @@ def update_pre_solve( dt, evalf: Optional[bool] = False, verbose: Optional[bool] = False, + store_result: Optional[bool] = True, ): """Pre-solve: auto-initialise history and apply advection correction. + ``store_result`` is accepted for interface parity with the + semi-Lagrangian flavour (which can sample without storing) and is not + used here: this flavour has nothing to sample. + On the first call, automatically initialises history from the current field values. If V_fn is set, also applies an explicit grid-based advection correction so that bdf() approximates the @@ -1668,6 +1906,8 @@ def update_pre_solve( self.initialise_history() # Update coefficient values for current effective_order and dt + self._refresh_source_snapshot() + _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) @@ -1707,7 +1947,14 @@ def update_post_solve( evalf: Optional[bool] = False, verbose: Optional[bool] = False, ): - r"""Shift history chain after solve: :math:`\psi^{*n} \leftarrow \psi^{*(n-1)}`.""" + r"""Shift history chain after solve: :math:`\psi^{*n} \leftarrow \psi^{*(n-1)}`. + + A history committed this step (a flux projected into level 0 and the + levels shifted by :meth:`commit_flux_to_history`) is already placed: + shifting again pushes the new value straight into level 1 and loses the + level it should hold. Invisible at order 1, a 150-fold error at order 2 + on the analytic Maxwell shear box. + """ self._dt = dt if verbose and uw.mpi.rank == 0: @@ -1719,22 +1966,21 @@ def update_post_solve( self._dt_history[0] = dt self._note_history_shift(dt) - ### copy values down the chain - for i in range(self.order - 1, 0, -1): - self.psi_star[i].data[...] = self.psi_star[i - 1].data[...] + if self._history_committed: + self._history_committed = False + else: + ### copy values down the chain + for i in range(self.order - 1, 0, -1): + self.psi_star[i].data[...] = self.psi_star[i - 1].data[...] - ### update the history fn - self.update_history_fn() + ### update the history fn + self.update_history_fn() if self._n_solves_completed < self.order: self._n_solves_completed += 1 return - def update_exp_coefficients(self, dt, tau_eff): - r"""Update the ETD-2 (exponential) coefficient values for this step.""" - _update_exp_values(self._exp_coeffs, dt, tau_eff) - class EulerianSUPG(Eulerian): r"""Eulerian history manager that assembles its transport: implicit advection with SUPG. @@ -1815,6 +2061,7 @@ def __init__( tau_shape: str = "inverse_sum", peclet_weight: float = 4.0, num_components=None, + transport_on_update: bool = False, ): order = int(order) if order not in (1, 2, 3): @@ -1839,6 +2086,15 @@ def __init__( ) self._advection_mode = "assembled" self._integrator = "am" if order == 1 else "bdf" + # When the unknown is not the solver's own (a stress carried by a Stokes + # solve, say) the manager transports its history itself, on the grid, + # in place of a semi-Lagrangian trace-back. See _transport_history. + self.transport_on_update = bool(transport_on_update) + self._transport_theta = 0.5 + self._transport_flat = None + self._transport_old = None + self._transport_solver = None + self._inflow_value = None self.V_fn = V_fn self.V_fn_history = None self.diffusivity = diffusivity @@ -1976,6 +2232,168 @@ def stabilisation_flux(self, R): column = R.reshape(len(R), 1) return self.tau() * (column * self.advecting_velocity(0)) + applies_inflow_value = True + + @_DDtBase.inflow_value.setter + def inflow_value(self, value): + """As the base class, and then drop the transport solver: this flavour + compiles the condition into the weak form, so the solver is stale.""" + _DDtBase.inflow_value.fset(self, value) + self._transport_solver = None + + # ----- transporting the history on the grid ----- + + def _transport_components(self): + """Independent components of the history variable as ``(i, j)`` pairs. + + A symmetric tensor contributes its upper triangle; every other shape + contributes every entry. + """ + rows, cols = self.psi_star[0].shape + if self.vtype == uw.VarType.SYM_TENSOR: + return [(i, j) for i in range(rows) for j in range(i, cols)] + return [(i, j) for i in range(rows) for j in range(cols)] + + @property + def transport_theta(self) -> float: + r"""Blend of the transport step: 0.5 Crank-Nicolson (default), 1 backward Euler. + + Distinct from :attr:`theta`, which weights the SCHEME's spatial terms at + each stored level. The transport of a history level is a time + discretisation of the same physical step as the scheme around it, so it + must be of the same order: backward Euler here is first order and costs + a factor of eight on uniform translation (0.0139 against 0.109 at + Courant 0.6). There is no reason to lower it; the setter exists to make + that measurable rather than to recommend it. + """ + return self._transport_theta + + @transport_theta.setter + def transport_theta(self, value): + value = float(value) + if not 0.0 < value <= 1.0: + raise ValueError(f"transport_theta must be in (0, 1], not {value}.") + self._transport_theta = value + + def _transport_residual(self, solver): + r"""Strong residual of one transport step, one entry per component. + + The theta rule at :attr:`transport_theta`, Crank-Nicolson by default. + """ + dim = self.mesh.dim + S, S_old = solver.u.sym, self._transport_old.sym + gradient = self.mesh.vector.gradient + a = self.advecting_velocity(0) + theta = self._transport_theta + entries = [] + for k in range(S.shape[1]): + new = sum(a[0, i] * gradient(S[0, k])[0, i] for i in range(dim)) + old = sum(a[0, i] * gradient(S_old[0, k])[0, i] for i in range(dim)) + entries.append((S[0, k] - S_old[0, k]) / self._delta_t + + theta * new + (1 - theta) * old) + return sympy.Matrix([entries]) + + def _transport_flux(self, solver): + """The SUPG flux of the transport step, one row per component.""" + return self.stabilisation_flux(self._transport_residual(solver)) + + def _build_transport_solver(self): + """The implicit SUPG solve that carries one history level over a step. + + The history's independent components are flattened onto one matrix + variable and marched together: backward Euler in time (the step is a + transport within a step, not the scheme's own time discretisation), + the manager's advecting velocity and its stabilisation parameter. No + boundary condition is imposed: the transported quantity is not the + solver's unknown and what enters at an inflow is the caller's to say. + """ + from underworld3.utilities._api_tools import Template + + components = len(self._transport_components()) + tag = self.instance_number + self._transport_flat = uw.discretisation.MeshVariable( + f"psi_transport_{tag}", self.mesh, (1, components), + vtype=uw.VarType.MATRIX, degree=self.degree, + continuous=self.continuous, varsymbol=rf"{{\psi^{{T}}_{{{tag}}}}}") + self._transport_old = uw.discretisation.MeshVariable( + f"psi_transport_old_{tag}", self.mesh, (1, components), + vtype=uw.VarType.MATRIX, degree=self.degree, + continuous=self.continuous, varsymbol=rf"{{\psi^{{T-}}_{{{tag}}}}}") + + class _HistoryTransport(uw.systems.SNES_MultiComponent): + F0 = Template(r"f_0", lambda solver: solver._manager._transport_residual(solver), + "Transport of a history level: time derivative and advection.") + F1 = Template(r"\mathbf{F}_1", lambda solver: solver._manager._transport_flux(solver), + "The SUPG flux of the transported history level.") + + solver = _HistoryTransport(self.mesh, u_Field=self._transport_flat, verbose=self.verbose) + solver._manager = self + solver.constitutive_model = uw.constitutive_models.Constitutive_Model + if self._inflow_value is not None: + # The inflow condition, weakly: on every boundary, the term is the + # NEGATIVE part of u.n, so it is active exactly where the flow enters + # and vanishes where it leaves. Walls (u.n = 0) contribute nothing, + # so no boundary needs naming and a migrating inflow patch is covered. + a = self.advecting_velocity(0) + normal_flow = sum(a[0, i] * self.mesh.Gamma[i] for i in range(self.mesh.dim)) + entering = sympy.Min(normal_flow, 0) + indices = self._transport_components() + incoming = self._inflow_value + # Sign: `entering` is non-positive, so -entering is |u.n| on the + # inflow and zero elsewhere, and the term is dissipative in + # (sigma - sigma_in). With the sign the other way it amplifies: + # measured on the cylinder, the stress reached 16 within ten steps. + condition = sympy.Matrix([[ + -entering * (self._transport_flat.sym[0, k] - incoming[i, j]) + for k, (i, j) in enumerate(indices)]]) + for boundary in self.mesh.boundaries: + solver.add_natural_bc(condition, boundary.name) + solver.petsc_options["snes_rtol"] = 1.0e-8 + solver.petsc_options["ksp_rtol"] = 1.0e-9 + solver.petsc_options["ksp_type"] = "gmres" + solver.petsc_options["pc_type"] = "asm" + solver.petsc_options["sub_pc_type"] = "ilu" + self._transport_solver = solver + return solver + + def _transport_history(self, dt, verbose=False): + """Carry every stored level forward by one step of the flow. + + Each level is transported once per step, so the level that is two + steps old has been carried twice: the grid counterpart of sampling + the semi-Lagrangian trace-back at two departure points. + """ + if self._transport_solver is None: + self._build_transport_solver() + indices = self._transport_components() + flat, previous = self._transport_flat, self._transport_old + + for level in range(self.order): + history = self.psi_star[level] + for k, (i, j) in enumerate(indices): + flat.array[:, 0, k] = history.array[:, i, j] + previous.array[...] = flat.array[...] + self._transport_solver.solve(verbose=verbose) + for k, (i, j) in enumerate(indices): + values = flat.array[:, 0, k] + history.array[:, i, j] = values + if i != j: + history.array[:, j, i] = values + + def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True): + """Refresh the scheme's coefficients and, when this manager owns the + transport of its history, carry every level forward by one step.""" + super().update_pre_solve(dt, evalf=evalf, verbose=verbose, + store_result=store_result) + if self.transport_on_update and self.V_fn is not None and dt: + self._transport_history(dt, verbose=verbose) + + def _object_viewer(self): + from IPython.display import Latex, display + + super()._object_viewer() + display(Latex(r"$\quad\mathbf{a} = $ " + self.V_fn._repr_latex_())) + display(Latex(rf"$\quad$ integrator: {self.integrator}, tau shape: {self.tau_shape}")) class CharacteristicTrace: @@ -2483,8 +2901,6 @@ def __init__( # substituted with a frozen snapshot variable that's refreshed each # step from psi_star[0]'s data array. The projection becomes a true # one-shot Galerkin projection. - self._psi_snapshot_enabled = False - self._psi_snapshot = None if varsymbol is None: @@ -2803,87 +3219,6 @@ def psi_fn(self, new_fn): self._psi_star_projection_solver.uw_function = self._build_projection_source(new_fn) return - def _build_projection_source(self, source_fn): - """Construct the row matrix used as the projection's ``uw_function``. - - Applies snapshot substitution (psi_star[0] → snap) when enabled. - Used by both ``psi_fn.setter`` and the ``initialise_history`` - fallback path so substitution semantics are consistent. - """ - if getattr(self, '_psi_star_use_multicomponent', False): - indep = self._psi_star_indep_indices - row = sympy.Matrix([[source_fn[i, j] for (i, j) in indep]]) - if self._psi_snapshot_enabled and self._psi_snapshot is not None: - ps0 = self.psi_star[0] - psi_snapshot = self._psi_snapshot - substitutions = { - ps0.sym[i, j]: psi_snapshot.sym[i, j] - for i in range(self.mesh.dim) - for j in range(self.mesh.dim) - } - row = row.subs(substitutions) - return row - else: - # Scalar / vector path: psi_star[0] is a scalar/vector field. If - # snapshot is needed for these vtypes, extend here similarly. - return source_fn - - def enable_source_snapshot(self): - """Enable snapshot substitution in the projection's source field. - - Call this once when the source expression (``psi_fn``) references - ``psi_star[0]`` itself — without it the projection's residual - ``(target − flux(psi_star[0]))·weight`` is implicit in the target - because target and source share the same data field. With Min-mode - plasticity at the yield kink, the implicit projection admits two - fixed points (elastic and yield branches); under timestep change the - iteration drifts to the elastic-branch fixed point and σ violates - the yield surface. - - The snapshot is a separate mesh variable matching ``psi_star[0]``'s - shape/vtype/degree. Each call to ``update_pre_solve`` copies - ``psi_star[0].array → psi_snapshot.array``, freezing the source's - input for the upcoming projection. Substitution makes the - projection's compiled C code read from ``psi_snapshot.array`` - instead of ``psi_star[0].array`` — there's no recompile per step, - just a memcpy. - - Idempotent: safe to call more than once. - """ - if not getattr(self, '_psi_star_use_multicomponent', False): - # Currently only wired for tensor projections (the case that - # exposed the bug). Scalar/vector extension is straightforward - # if needed later. - return - - if self._psi_snapshot is None: - ps0 = self.psi_star[0] - # NOTE: this currently registers a persistent MeshVariable in the - # mesh DM, which is overkill for a transient buffer that's only - # read by this DDt's projection. A future improvement would be - # a transient/scratch-variable mechanism (likely backed by - # PETSc's auxiliary Vec machinery — already used elsewhere in - # the codebase via DMSetAuxiliaryVec_UW) so the snapshot doesn't - # accumulate in the DM across DDt creations. See: - # docs/developer/ai-notes/historical-notes.md for the - # variable-deletion limitation context. - self._psi_snapshot = uw.discretisation.MeshVariable( - f"psi_snapshot_{self.instance_number}", - self.mesh, - ps0.shape, - vtype=ps0.vtype, - degree=ps0.degree, - continuous=ps0.continuous, - ) - # Initialise psi_snapshot's data to current psi_star[0]'s data - # so the source evaluates consistently before the first refresh. - self._psi_snapshot.data[...] = ps0.data[...] - - self._psi_snapshot_enabled = True - - # Re-run the psi_fn setter so the substitution is applied to the - # currently-installed projection source. - self.psi_fn = self._psi_fn def initialise_history(self): @@ -3212,6 +3547,15 @@ def _shift_history_with_blend(self, dt, dt_physical=None): phi * self.psi_star[i - 1].array[...] + (1 - phi) * self.psi_star[i].array[...] ) + def carried_tensors(self, level: int = 0): + """The carried history as one tensor per vertex, non-dimensional, with the + vertices: ``(values[n, d, d], coords[n, cdim])``.""" + history = self.psi_star[level] + dim = self.mesh.dim + values = np.asarray(_to_nondim_ndarray(np.asarray(history.array), units=history.units)).reshape(-1, dim, dim) + 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. @@ -3524,8 +3868,7 @@ def update_pre_solve( # variables already live in non-dimensional space) while keeping # the callback sync that pushes values into the underlying PETSc # local Vec. - if self._psi_snapshot_enabled and self._psi_snapshot is not None: - self._psi_snapshot.data[...] = self.psi_star[0].data[...] + self._refresh_source_snapshot() # Update coefficient values for current effective_order and dt _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) @@ -3609,26 +3952,6 @@ def update_pre_solve( return - def update_exp_coefficients(self, dt, tau_eff): - r"""Update the scalar ETD-2 (exponential) coefficient UWexpressions. - - Sets ``self._exp_coeffs[0].sym = α = exp(-Δt/τ_eff)`` and - ``self._exp_coeffs[1].sym = φ = (1-α)/(Δt/τ_eff)`` so the next solve - uses the correct exponential coefficients via PetscDSSetConstants - on the next ``_update_constants`` call. - """ - _update_exp_values(self._exp_coeffs, dt, tau_eff) - - @property - def _exp_alpha(self): - """Convenience accessor for the ETD-2 ``α`` coefficient UWexpression.""" - return self._exp_coeffs[0] - - @property - def _exp_phi(self): - """Convenience accessor for the ETD-2 ``φ`` coefficient UWexpression.""" - return self._exp_coeffs[1] - def update_forcing_history(self, forcing_fn=None, evalf=False, verbose=False): r"""Refresh ``forcing_star`` from ``forcing_fn`` via direct nodal evaluation. @@ -3781,6 +4104,7 @@ class Lagrangian(_DDtBase): Lagrangian_Swarm : For user-provided swarms. """ + commits_flux_in_post_solve = True instances = ( 0 # count how many of these there are in order to create unique private mesh variable ids @@ -3801,6 +4125,8 @@ def __init__( order=1, smoothing=0.0, fill_param=3, + proxy_location="cells", + proxy_sampling="reconstruct", ): super().__init__() @@ -3827,6 +4153,8 @@ def __init__( vtype=vtype, proxy_degree=degree, proxy_continuous=continuous, + proxy_location=proxy_location, + proxy_sampling=proxy_sampling, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", ) ) @@ -3836,6 +4164,21 @@ def __init__( self._init_coefficient_expressions(order, 0.5, with_exp=False) dudt_swarm.populate(fill_param) + # The class owns this swarm and carries the stress on it, so it also + # keeps it populated: without refilling starved cells, a flow that + # carries particles out through an open boundary empties the downstream + # cells and the proxy has nothing to interpolate. The control both fills + # starved cells and caps over-full ones (a bare refill-only dict grows + # the swarm without bound where particles pile up against a wall): the + # bounds are set from the initial occupancy, with a floor at the linear + # fit minimum. + _npart = int(uw.mpi.comm.allreduce(np.asarray(dudt_swarm._particle_coordinates.data).shape[0], op=uw.MPI.SUM)) + _ncell = int(uw.mpi.comm.allreduce(mesh._centroids.shape[0], op=uw.MPI.SUM)) + _mean = _npart / max(_ncell, 1) + dudt_swarm.population_control = { + "min_per_cell": max(mesh.dim + 1, int(0.5 * _mean)), + "max_per_cell": max(2 * (mesh.dim + 1), int(3.0 * _mean)), + } # Register with the active default model as a Snapshottable # state-bearer. Safe if no model is active. @@ -3891,20 +4234,21 @@ def initialise_history(self): be called manually after setting initial conditions. """ psi_star_0 = self.psi_star[0] - # Component-wise write through the canonical (N, components) storage. - # Indexing the SwarmVariable itself (``psi_star_0[i, j]``) returns a - # *symbolic* component with no ``.data`` — the modern component - # address is ``.data[:, var._data_layout(i, j)]`` (audit SWARM-06). + # Every component evaluated before any is written (audit SWARM-06): a + # partial write marks the proxy stale and a later evaluation of a psi_fn + # that reads psi_star would see a half-updated history. coords = np.asarray(self.swarm._particle_coordinates.data) + updated = {} for i in range(psi_star_0.shape[0]): for j in range(psi_star_0.shape[1]): - updated_psi = uw.function.evaluate( - self.psi_fn[i, j], - coords, - ) - psi_star_0.data[:, psi_star_0._data_layout(i, j)] = np.asarray( - updated_psi + ij = psi_star_0._data_layout(i, j) + if ij in updated: + continue + updated[ij] = np.asarray( + uw.function.evaluate(self.psi_fn[i, j], coords) ).reshape(-1) + for ij, vals in updated.items(): + psi_star_0.data[:, ij] = vals # Copy to all other history slots for k in range(1, self.order): @@ -3935,8 +4279,14 @@ def update_pre_solve( dt: float, evalf: Optional[bool] = False, verbose: Optional[bool] = False, + store_result: bool = True, + **_ignored, ): - """Pre-solve: auto-initialise history on first call.""" + """Pre-solve: auto-initialise history on first call. + + ``store_result`` is accepted for a uniform hook signature and ignored: + this flavour records the stress at its particles in the post-solve. + """ self._dt = dt if not self._history_initialised: @@ -3953,6 +4303,7 @@ def update_post_solve( dt: float, evalf: Optional[bool] = False, verbose: Optional[bool] = False, + **_ignored, ): """Shift history chain and advect swarm after solve.""" self._dt = dt @@ -3973,22 +4324,24 @@ def update_post_solve( self.psi_star[i].array[...] = self.psi_star[i - 1].array[...] - # Now update the swarm variable - + # Now update the swarm variable. psi_fn is the constitutive flux and + # reads psi_star[0] itself, so every component is evaluated BEFORE any is + # written: writing one marks the proxy stale, and a later evaluation + # would then read a history that is half new (audit SWARM-06, the reason + # Lagrangian_Swarm computes all its components first). psi_star_0 = self.psi_star[0] - # Grab the current psi values at the (pre-advection) particle - # positions via the canonical component storage (audit SWARM-06). coords = np.asarray(self.swarm._particle_coordinates.data) + updated = {} for i in range(psi_star_0.shape[0]): for j in range(psi_star_0.shape[1]): - updated_psi = uw.function.evaluate( - self.psi_fn[i, j], - coords, - evalf=evalf, - ) - psi_star_0.data[:, psi_star_0._data_layout(i, j)] = np.asarray( - updated_psi + ij = psi_star_0._data_layout(i, j) + if ij in updated: + continue + updated[ij] = np.asarray( + uw.function.evaluate(self.psi_fn[i, j], coords, evalf=evalf) ).reshape(-1) + for ij, vals in updated.items(): + psi_star_0.data[:, ij] = vals # Now update the swarm locations @@ -4097,6 +4450,7 @@ class Lagrangian_Swarm(_DDtBase): Eulerian : Pure mesh-based history (no particle tracking). """ + commits_flux_in_post_solve = True instances = ( 0 # count how many of these there are in order to create unique private mesh variable ids @@ -4445,7 +4799,7 @@ class IntegrationPointSemiLagrangian(_DDtBase): history needs. What is not here (yet): units-aware velocity reduction, ALE / old-frame - trace-back, forcing history, checkpoint state. Use + trace-back. Use :class:`SemiLagrangian` for those, or :class:`Lagrangian_Swarm` when the history should ride on particles rather than on the rule. @@ -4459,6 +4813,62 @@ class IntegrationPointSemiLagrangian(_DDtBase): velocity history caches it by evaluation at each time level. """ + #: Departure points restored to the boundary take :attr:`inflow_value` + #: 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. + + Two ladders have to move together here. ``psi_star`` holds the history + at the integration points, where the assembler reads it; ``psi_snap`` + holds it nodally, which is what the next trace-back samples at the + departure points. A flux committed only to ``psi_star`` does not + survive a step -- the fill overwrites it from ``psi_snap`` -- so the + committed stress goes to both. + + The projection, rather than a pointwise evaluation, is what the nodal + trace-back does and is what the snapshot needs: the trace-back samples + the snapshot between its nodes, so the snapshot has to be the L2 fit of + the flux on the history space, not the flux read at the nodes. Its + target is a separate field from the flux's inputs, so the projection is + 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) + + # Oldest first, so each level reads the one above before it is written. + for level in range(self.order - 1, 0, -1): + self.psi_star[level].data[...] = self.psi_star[level - 1].data[...] + self.psi_snap[level].data[...] = self.psi_snap[level - 1].data[...] + + points = np.asarray(history.integration_points).reshape(-1, self.mesh.cdim) + for column in range(len(columns)): + nodal = self._commit_flat.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) + + self._history_committed = True def __init__( self, @@ -4474,13 +4884,20 @@ def __init__( order: int = 1, theta: float = 0.5, monotone_mode: Optional[str] = None, - **_unsupported, + with_forcing_history: bool = False, + store_smoothing: float = 0.0, + **_unsupported, # TODO(BUG): swallowed without a stated failure mode (Charter 5) ): super().__init__() self.vtype = vtype self.monotone_mode = monotone_mode + self.with_forcing_history = bool(with_forcing_history) + self.store_smoothing = store_smoothing self.mesh = mesh - self.bcs = list(bcs) if bcs is not None else [] # per instance, never a shared default + if bcs: + raise ValueError("IntegrationPointSemiLagrangian applies no boundary conditions to its " + "store; an inflow is set through inflow_value") + self.bcs = [] self.verbose = verbose self.degree = degree self.continuous = continuous @@ -4551,7 +4968,104 @@ def __init__( # V_fn evaluated at the nodes at that time, so V_fn may be any # expression (variables, ramping constants, swarm proxies). self._n_v = max(order, 2) # velocity levels the segments read - self._init_coefficient_expressions(order, self.theta, with_exp=False) + self._init_coefficient_expressions(order, self.theta, with_exp=True) + self._register_with_default_model() + # The forcing (strain-rate) history the second-order exponential + # integrator reads. Unlike the nodal flavour, which re-evaluates the + # strain rate at its nodes, this one is carried along the same + # characteristic as the stress: ``forcing_star`` at the points is what + # the parcel saw a step ago, sampled from the continuous snapshot at + # the departure point. Committed from the solved strain rate in + # :meth:`update_forcing_history` (the constitutive model's post-solve + # hook), traced in :meth:`_fill_slots`. + self.forcing_star = None + self._forcing_projection = None + self._forcing_flat = None + if self.with_forcing_history: + self.forcing_star = uw.discretisation.IntegrationPointVariable( + f"forcing_star_ip_{inst}", mesh, vtype=vtype, + varsymbol=rf"{{ \dot\varepsilon^{{ * }}_{{ [{inst}] }} }}", units=None) + self.forcing_snap = uw.discretisation.MeshVariable( + f"forcing_snap_ip_{inst}", mesh, vtype=vtype, + degree=degree, continuous=continuous, + varsymbol=rf"{{ \dot\varepsilon^{{ (n) }}_{{ [{inst}] }} }}", units=None) + + @property + def state(self) -> "DDtIntegrationPointState": + return DDtIntegrationPointState( + **self._core_state_kwargs(), + psi_star_var_names=[ps.clean_name for ps in self.psi_star], + psi_snap_var_names=[ps.clean_name for ps in self.psi_snap], + with_forcing_history=bool(self.with_forcing_history), + history_committed=bool(getattr(self, "_history_committed", False)), + ) + + @state.setter + def state(self, s: "DDtIntegrationPointState") -> None: + self._validate_state_schema(s, DDtIntegrationPointState) + self._validate_psi_star_names(s.psi_star_var_names) + if s.psi_snap_var_names != [ps.clean_name for ps in self.psi_snap]: + raise ValueError("psi_snap variable names changed since snapshot") + if s.with_forcing_history != bool(self.with_forcing_history): + raise ValueError("with_forcing_history differs between snapshot and instance") + self._restore_core_state(s, am_theta=self.theta) + self._history_committed = bool(s.history_committed) + + @property + def store_smoothing(self) -> float: + r"""Coefficient :math:`c` of the Laplacian term in the store projection, + :math:`\alpha = c\,h(\mathbf{x})^2` with :math:`h` the local cell size. + + Every step the new flux is L2-projected onto the continuous snapshot and + read back at the points. That cycle is a consistent-mass Galerkin + transport of the carried stress and has no dissipation at the cell + scale, so below Courant one a cell-scale mode of the stress grows from + round-off at a rate :math:`\gamma` set by the elastic feedback (about + 2.4 per unit time on the Maxwell Waters-King start-up, 1.8 with a + solvent fraction of 0.2, and negligible at 0.59). The term + :math:`\alpha\nabla^2` in the projection multiplies wavenumber + :math:`k` by :math:`1/(1+\alpha k^2)` once per step, so the mode is held + when :math:`\alpha (\pi/h)^2 \gtrsim \gamma\,\Delta t`, i.e. + :math:`\alpha \approx \gamma\,\Delta t\,(h/\pi)^2`. Measured on the + 1/32 mesh at :math:`\Delta t = 0.01`: 1e-5 holds it for eight time + units at 0.1% on the peak, 3e-5 holds it unconditionally at 0.5%, + 1e-4 costs 3%. In units of the mesh cell-size field (RMS vertex-to- + centroid distance, about 2h/3 on triangles) that is :math:`c` between + 0.03 and 0.07; the irregular mesh needs 0.07. Zero (the default) is + the plain projection. Only the stress store is smoothed; the forcing + history is not. + """ + return self._store_smoothing + + @store_smoothing.setter + def store_smoothing(self, value): + value = float(value) + if value < 0.0: + raise ValueError(f"store_smoothing must be >= 0, got {value}") + self._store_smoothing = value + + def carried_tensors(self, level: int = 0): + """The carried history as one tensor per point, non-dimensional, with the + points: ``(values[n, d, d], coords[n, cdim])``. The integration-point + storage keeps the independent components in columns; this is the one + place that unpacks them.""" + history = self.psi_star[level] + dim = self.mesh.dim + cols = _storage_components(self.vtype, (dim, dim)) + data = np.asarray(history.data) + values = np.zeros((data.shape[0], dim, dim)) + for k, (i, j) in enumerate(cols): + values[:, i, j] = data[:, k] + values[:, j, i] = data[:, k] + points = np.asarray(history.integration_points).reshape(-1, self.mesh.cdim) + return values, points + + def _store_smoothing_alpha(self): + """The smoothing the store projection uses this step: a field, so the + dose follows the local cell on a graded mesh.""" + if self._store_smoothing <= 0.0: + return 0.0 + return self._store_smoothing * self.mesh.cell_size() ** 2 def spatial_weights(self): """As the base class, except that at ``theta = 1`` the old-level @@ -4730,6 +5244,17 @@ def _write_components(self, var, expr, coords, evaluate=None, **kwargs): _to_nondim_ndarray(vals, units=self._psi_units) ).reshape(-1) + def _write_inflow(self, var, coords, rows): + """Overwrite ``rows`` of ``var`` with :attr:`inflow_value` evaluated at + ``coords`` (the restored positions of ALL the points, so that the + read, which is collective when the value holds a field, is made on + every rank; only ``rows`` are written).""" + expr = sympy.Matrix(self._inflow_value) + for column, (i, j) in enumerate(self._components): + vals = uw.function.evaluate(expr[i, j], coords) + var.data[rows, column] = np.asarray( + _to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1)[rows] + def _segment_dt(self, j, dt): """Length of segment ``j`` (0 = the current step).""" if j == 0: @@ -4756,10 +5281,65 @@ def _fill_slots(self, dt, evalf): evaluate=uw.function.global_evaluate, evalf=evalf, monotone=self.monotone_mode, ) + if self._inflow_value is not None: + # A departure point that left the domain was restored to the + # boundary by the trace. What it should carry is the stress of + # the fluid ENTERING there, not a sample of the boundary edge: + # on a box channel that edge sample fed a growing mode in the + # inlet cell column (#745). The unclamped end point says which + # points left. + X_raw = trace.departure_points(key, X0, tuple(segments), evalf=evalf, + clamp_final=False) + left = np.any(np.abs(np.asarray(X_raw) - np.asarray(X)) > 0.0, axis=1) + # every rank decides together whether the (collective) read happens + any_left = uw.mpi.comm.allreduce(int(left.any()), op=uw.MPI.SUM) if uw.mpi.size > 1 else int(left.any()) + if any_left: + self._write_inflow(self.psi_star[k], X, left) + if k == 0 and self.forcing_star is not None: + # the strain rate the parcel saw a step ago, at the same + # departure point as its stress + self._write_components( + self.forcing_star, self.forcing_snap.sym, X, + evaluate=uw.function.global_evaluate, evalf=evalf, + ) + + def update_forcing_history(self, forcing_fn=None, evalf=False, verbose=False): + """Commit the solved strain rate as the forcing history: an L2 fit onto + the continuous snapshot, read back at the points. Same construction + as :meth:`commit_flux_to_history`, for the same reason -- the next + trace-back samples the snapshot between its nodes. A no-op unless + ``with_forcing_history`` was asked for.""" + if self.forcing_star is None or forcing_fn is None: + return + forcing = sympy.Matrix(forcing_fn) + columns = _storage_components(self.vtype, forcing.shape) + if self._forcing_projection is None: + self._forcing_flat = uw.discretisation.MeshVariable( + f"forcing_nodal_{self.instance_number}", self.mesh, (1, len(columns)), + vtype=uw.VarType.MATRIX, degree=self.degree, continuous=self.continuous, + varsymbol=rf"{{E^{{\mathrm{{nodal}}}}_{{{self.instance_number}}}}}") + 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.uw_function = sympy.Matrix( + [[forcing[i, j] for (i, j) in columns]]) + self._forcing_projection.smoothing = 0.0 + self._forcing_projection.solve(verbose=verbose) + 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( + self._forcing_flat.data[:, column]).reshape(-1) + self.forcing_star.data[:, column] = np.asarray(uw.function.evaluate( + self._forcing_flat.sym[0, column], points)).reshape(-1) def initialise_history(self): """Start every snapshot and slot from the current field, so - ``bdf()`` is zero on the first step.""" + ``bdf()`` is zero on the first step. A history already placed by + :meth:`commit_flux_to_history` is the start, and is kept.""" + if self._history_committed: + self.characteristics.initialise_levels(self._n_v) + self._history_initialised = True + return self._record_current() for k in range(1, self.order): self.psi_snap[k].data[...] = self.psi_snap[0].data[...] @@ -4770,18 +5350,32 @@ def initialise_history(self): self.psi_star[k].data[...] = self.psi_star[0].data[...] self._history_initialised = True - def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): + def update_pre_solve(self, dt, evalf=False, verbose=False, + store_result=True, **_ignored): + """Carry the history to the departure points of this step. + + ``store_result=False`` says the snapshots already hold what is to be + transported and must not be rebuilt from ``psi_fn`` -- the viscoelastic + case, where ``psi_fn`` is the constitutive flux and the flux is a + function of the history. Recording it there applies the constitutive + update a second time on a field that is already the stress: on the + analytic Maxwell shear box that overshoots the relaxed stress by 12.8% + where the nodal and grid flavours sit at 1.5% (#732). The snapshots are + placed by ``commit_flux_to_history`` instead, and the shift with them. + """ self._dt = dt if not self._history_initialised: self.initialise_history() _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) _update_am_values(self._am_coeffs, self.effective_order, self.theta) - for k in range(self.order - 1, 0, -1): - self.psi_snap[k].data[...] = self.psi_snap[k - 1].data[...] + if store_result: + for k in range(self.order - 1, 0, -1): + self.psi_snap[k].data[...] = self.psi_snap[k - 1].data[...] trace = self.characteristics if self._owns_characteristics: trace.begin_step(dt) - self._record_current() + if store_result: + self._record_current() self._fill_slots(dt, evalf) if self._owns_characteristics: trace.finish_step() @@ -4797,3 +5391,422 @@ def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): self._note_history_shift(dt) if self._n_solves_completed < self.order: self._n_solves_completed += 1 + + +class ForwardSemiLagrangian(_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. + + The carried field is known at the launch points, the mesh's integration + points with their quadrature weights scaled by the cell measure. Each step + every point is moved forward one step along the velocity, the arrivals in + each cell are fitted by weighted least squares to a linear polynomial, and + that discontinuous P1 field is ``psi_star[0]``, what the weak form reads. + After the solve the new flux is projected onto the continuous P1 space and + read at the launch points. Nothing persists at the arrivals: one fixed + point set, one velocity-dependent map, one fit, no particle state and no + repopulation. + + Why this form. The launch points are inside the cells, so no point sits + on a no-slip wall with zero velocity (the nodal history's wall-layer + defect); the arrivals are fitted per cell (measured 17 s a step on the + confined cylinder against 28 for the integration-point history, which + samples its store at every foot); and the wall stress it carries is the + more self-consistent (the drag by stress integral and by reaction agree + to 1.4% where the integration-point history has them 12% apart). + It is, like the integration-point history, a consistently transported + scheme with no dissipation of its own at the cell scale: below Courant one + on a Maxwell element a cell-scale mode grows from round-off, and the + read-back projection needs :attr:`flux_smoothing` for the same reason and + at the same dose as the integration-point store. See + :doc:`/developer/subsystems/stress-transport`. + + First order only. In parallel the arrivals that left their rank travel + with their values and weights to every rank, and each rank keeps the ones + that landed in its own cells, so a point that crosses a seam is fitted by + the rank that owns its arrival cell. A point that leaves the domain is + dropped, so a periodic seam is not crossed. The launch set and the cell + geometry are taken once from the mesh, so the flavour does not follow a + mesh that moves. An inflow cell (a boundary cell whose boundary face has + fluid entering) that receives less than it launched has the missing share + filled with :attr:`inflow_value` when one is set; a cell whose arrivals + cannot determine a linear fit keeps its previous fit, the one piece of + state carried between steps. + + TODO(DESIGN): periodic seams (wrap the end point), moving meshes (refresh + on the topology version). + """ + + applies_inflow_value = True + + def __init__( + self, + mesh, + psi_fn, + V_fn, + vtype=VarType.SCALAR, + varsymbol: Optional[str] = None, + order: int = 1, + theta: float = 0.5, + **_unsupported, + ): + super().__init__() + if order != 1: + raise NotImplementedError("ForwardSemiLagrangian carries one level; order must be 1") + if mesh.cdim != mesh.dim: + raise NotImplementedError("ForwardSemiLagrangian fits in the embedding coordinates; no manifolds") + if _unsupported: + warnings.warn(f"ForwardSemiLagrangian 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 + self.degree = 1 + self.continuous = 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]]) + 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)} " + f"but vtype={vtype} on a cdim={mesh.cdim} mesh needs {expected}") + self._init_history_tracking(1) + if varsymbol is None: + varsymbol = rf"u_{{ [{self.instance_number}] }}" + inst = self.instance_number + psi_units = uw.get_units(self._psi_fn) + if psi_units is not None and not uw.get_default_model().has_units(): + psi_units = None + self._psi_units = psi_units + self.psi_star = [ + uw.discretisation.MeshVariable( + f"psi_star_fwd_{inst}", mesh, vtype=vtype, degree=1, continuous=False, + varsymbol=rf"{{ {varsymbol}^{{ * }} }}", units=psi_units) + ] + self._components = _storage_components(vtype, tuple(self.psi_star[0].sym.shape)) + self.num_components = len(self._components) + # The launch set: the integration points, cell-major in the rule's order, + # with the rule's weights scaled by the cell measure (a share of area). + # The launch values live in an integration-point variable so a snapshot + # captures them with every other variable; the class reads them as columns. + self._launch_var = uw.discretisation.IntegrationPointVariable( + f"launch_fwd_{inst}", mesh, vtype=vtype, + varsymbol=rf"{{ {varsymbol}^{{ \ell }} }}", units=psi_units) + launch_var = self._launch_var + 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 " + 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() + self._launch_cell = np.repeat(np.arange(self._cell_measure.size), self._nq) + self._launch_weights = np.tile(w_ref, self._cell_measure.size) * self._cell_measure[self._launch_cell] / w_ref.sum() + self._n_relocated = 0 # arrivals this rank took from other ranks at the last carry + self._launch_geometry = self._geometry_stamp() + self._bface_cell, self._bface_centroid, self._bface_normal = self._boundary_faces() + # The flux read at the launch points, through a continuous P1 projection. + # `flux_smoothing` is the Laplacian coefficient of that projection, a + # number or a field (length^2): c * mesh.cell_size()**2 with c between + # 0.03 and 0.07 is the dose the integration-point store needs below + # Courant one on a Maxwell element, and this cycle needs the same. + self.flux_smoothing = 0.0 + self._flux_var = uw.discretisation.MeshVariable( + f"flux_fwd_{inst}", mesh, (1, self.num_components), vtype=VarType.MATRIX, + degree=1, continuous=True, varsymbol=rf"{{ F^{{\mathrm{{nodal}}}}_{{ [{inst}] }} }}") + self._flux_projection = None + self._n_v = 2 + self._init_coefficient_expressions(1, self.theta, with_exp=True) + self._register_with_default_model() + + @property + def _launch_values(self): + """The carried values at the launch points, one column per component.""" + return np.asarray(self._launch_var.data) + + @_launch_values.setter + def _launch_values(self, values): + # one write, one PETSc flush (a per-column write is a collective + # round trip per component) + self._launch_var.data[:, :] = np.asarray(values).reshape(self._launch.shape[0], self.num_components) + + @property + def state(self) -> "DDtForwardState": + return DDtForwardState( + **self._core_state_kwargs(), + psi_star_var_names=[ps.clean_name for ps in self.psi_star], + launch_var_name=self._launch_var.clean_name, + ) + + @state.setter + def state(self, s: "DDtForwardState") -> None: + self._validate_state_schema(s, DDtForwardState) + self._validate_psi_star_names(s.psi_star_var_names) + if s.launch_var_name != self._launch_var.clean_name: + raise ValueError("launch variable name changed since snapshot") + self._restore_core_state(s, am_theta=self.theta) + + # ------------------------------------------------------------------ + def _geometry_stamp(self): + """The mesh geometry the launch set was built for: the vertex count and + the coordinate sum (a moved or re-meshed mesh changes one of them; adding + a variable, which rebuilds the DM, changes neither).""" + coords = np.asarray(self.mesh.X.coords) + return (coords.shape, float(coords.sum())) + + def _cell_measures(self): + """Area (2-D) or volume (3-D) of every cell, in cell order.""" + dm = self.mesh.dm + c0, c1 = dm.getHeightStratum(0) + return np.array([dm.computeCellGeometryFVM(c)[0] for c in range(c0, c1)]) + + def _boundary_faces(self): + """For every face on the domain boundary: its owning cell, its centroid + and its outward normal (centroid of the face away from the centroid of + the cell: outward on a convex cell). The inflow test reads the velocity + at these centroids.""" + dm = self.mesh.dm + d = self.mesh.dim + c0, c1 = dm.getHeightStratum(0) + f0, f1 = dm.getHeightStratum(1) + cells, centroids, normals = [], [], [] + # A face with one local cell is a domain boundary face OR a partition face; + # 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 " + "partition faces cannot be told from domain faces in parallel") + for f in range(f0, f1): + if dm.getSupportSize(f) != 1: + continue + if label is not None and label.getValue(f) == -1: + continue + c = dm.getSupport(f)[0] - c0 + _, fc, fn = dm.computeCellGeometryFVM(f) + _, cc, _ = dm.computeCellGeometryFVM(c + c0) + # the face's own normal (the centroid difference is the cell's median, + # normal to the face only on a right cell), pointing out of the cell + n = np.asarray(fn[:d], dtype=float) + if np.dot(n, np.asarray(fc[:d]) - np.asarray(cc[:d])) < 0.0: + n = -n + cells.append(c); centroids.append(np.asarray(fc[:d])); normals.append(n / np.linalg.norm(n)) + return (np.array(cells, dtype=int), + np.array(centroids).reshape(-1, d), np.array(normals).reshape(-1, d)) + + def _inflow_cells(self, trace, evalf): + """Boundary cells whose boundary face has fluid entering now. The + velocity read is collective, so a rank with no boundary face still + takes part, with no points.""" + mask = np.zeros(self._cell_measure.size, dtype=bool) + v = trace.velocity_at(trace.V_matrix(), self._bface_centroid, use_global=True, evalf=evalf) + if self._bface_cell.size == 0: + return mask + v = np.asarray(v).reshape(-1, self.mesh.dim) + # a wall with u.n = 0 to round-off (free slip) must not flip in and out + entering = np.einsum("fi,fi->f", v, self._bface_normal) < -1.0e-10 * (np.abs(v).max() if v.size else 0.0) + mask[self._bface_cell[entering]] = True + return mask + + @property + def psi_fn(self): + r"""Current symbolic expression :math:`\psi` being tracked.""" + return self._psi_fn + + @psi_fn.setter + 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}") + self._psi_fn = new_fn + + def _object_viewer(self): + from IPython.display import Latex, display + super()._object_viewer() + display(Latex(r"$\quad\psi = $ " + self.psi_fn._repr_latex_())) + display(Latex(r"$\quad\mathbf{v} = $ " + sympy.Matrix(self.V_fn)._repr_latex_())) + display(Latex(r"$\quad$Carried forward from the integration points, fitted per cell")) + + def _evaluate_at_launch(self, expr): + """Every stored component of ``expr`` at the launch points, as columns, + through the continuous P1 projection.""" + expr = sympy.Matrix(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.uw_function = sympy.Matrix([[expr[i, j] for (i, j) in self._components]]) + self._flux_projection.smoothing = self.flux_smoothing + self._flux_projection.solve() + 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) + return out + + def carried_tensors(self, level: int = 0): + """The carried values as one tensor per launch point, non-dimensional, + with the points: ``(values[n, d, d], coords[n, cdim])``.""" + dim = self.mesh.dim + values = np.zeros((self._launch.shape[0], dim, dim)) + for k, (i, j) in enumerate(self._components): + values[:, i, j] = self._launch_values[:, k] + if self.vtype == VarType.SYM_TENSOR: + values[:, j, i] = self._launch_values[:, k] + return values, self._launch + + # ------------------------------------------------------------------ + def _fit_arrivals(self, X, values, cell=None, inflow=None): + """Weighted least-squares linear fit of the carried values at their + arrival points, cell by cell, written into ``psi_star[0]``. + + ``cell`` is the owning cell of each arrival when it is known exactly + (the launch points themselves). Otherwise the arrivals of every rank + are gathered with their values and weights, each rank keeps the ones + in its partition and locates them in its cells; a point outside the + domain (it left through an outflow) is dropped. + A boundary cell that lost more points than it received has the + missing share filled with the inflow value at its own dofs, weighted + by that share: the state of the part of the cell nothing has reached + is the incoming fluid. A cell whose arrivals cannot determine a linear + fit keeps its previous one. + """ + mesh = self.mesh + d = mesh.dim + npar = d + 1 + 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] + 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 + # moment matrix is well conditioned exactly when the arrivals span the cell + A = np.concatenate([np.ones((X.shape[0], 1)), (X[:, :d] - centroid[cell]) / h[cell, None]], axis=1) + M = np.zeros((ncell, npar, npar)) + np.add.at(M, cell, w[:, None, None] * A[:, :, None] * A[:, None, :]) + R = np.zeros((ncell, npar, self.num_components)) + np.add.at(R, cell, w[:, None, None] * A[:, :, None] * values[:, None, :]) + received = np.bincount(cell, weights=w, minlength=ncell) + dofs = np.asarray(self.psi_star[0].coords_nd).reshape(-1, mesh.cdim) + ndof = dofs.shape[0] // ncell if ncell else 0 + dofs = dofs.reshape(ncell, ndof, mesh.cdim) + Adof = np.concatenate([np.ones((ncell, ndof, 1)), (dofs[:, :, :d] - centroid[:, None, :]) / h[:, None, None]], axis=2) + if self._inflow_value is not None and inflow is not None: + deficit = np.clip(self._cell_measure - received, 0.0, None) * inflow + fed = deficit > 0.0 + # evaluate is collective: every rank reads the inflow value at every + # boundary-cell dof, whether or not any of its cells is short + bcells = np.flatnonzero(np.isin(np.arange(ncell), self._bface_cell)) + filled = np.column_stack([ + _to_nondim_ndarray(uw.function.evaluate(self._inflow_value[i, j], + dofs[bcells].reshape(-1, mesh.cdim)), + units=self._psi_units).reshape(-1) + for (i, j) in self._components]).reshape(bcells.size, ndof, self.num_components) + short = fed[bcells] + if short.any(): + sel = bcells[short] + wi = np.broadcast_to((deficit[sel] / ndof)[:, None], (sel.size, ndof)) # per dof + M[sel] += np.einsum("cq,cqi,cqj->cij", wi, Adof[sel], Adof[sel]) + R[sel] += np.einsum("cq,cqi,cqk->cik", wi, Adof[sel], filled[short]) + # with the columns scaled, the eigenvalue ratio of the moment matrix is a + # conditioning number: 1e-6 rejects a cell whose arrivals sit on a line + ev = np.linalg.eigvalsh(M) + fit_ok = ev[:, 0] > 1.0e-6 * np.maximum(ev[:, -1], 1.0e-300) + beta = np.zeros_like(R) + beta[fit_ok] = np.linalg.solve(M[fit_ok], R[fit_ok]) + fitted = np.einsum("cqi,cik->cqk", Adof, beta).reshape(ncell * ndof, self.num_components) + rows_ok = np.repeat(fit_ok, ndof) + for k in range(self.num_components): + column = np.array(self.psi_star[0].data[:, k]) + column[rows_ok] = fitted[rows_ok, k] + self.psi_star[0].data[:, k] = column + + def initialise_history(self): + """Start from the current field: its values at the launch points, and + their fit, so ``bdf()`` is zero on the first step. A history already + placed by :meth:`commit_flux_to_history` is the start, and is kept.""" + self.characteristics.initialise_levels(self._n_v) + if self._history_committed: + self._history_initialised = True + return + self._launch_values = self._evaluate_at_launch(self._psi_fn) + self._fit_arrivals(self._launch, self._launch_values, cell=self._launch_cell) + self._history_initialised = True + + def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, **_ignored): + """Carry the launch values forward one step and fit the arrivals. + + ``store_result=False`` says the launch values were placed by + :meth:`commit_flux_to_history` (the viscoelastic case); otherwise the + tracked field is read at the launch points first. + """ + self._dt = dt + if self._geometry_stamp() != self._launch_geometry: + raise NotImplementedError( + "ForwardSemiLagrangian: 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: + 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) + if store_result: + self._launch_values = self._evaluate_at_launch(self._psi_fn) + # forward along the same characteristic the backward flavours trace, with + # the step reversed: start velocity v^n, mid-time velocity at n+1/2. The + # end point is not restored to the domain: a point that leaves has left. + key = (_basis_key_of(self.psi_star[0]), "launch") + X = trace.departure_points(key, self._launch, (("first", 0, -float(dt)),), + evalf=evalf, clamp_final=False) + self._fit_arrivals(np.asarray(X), self._launch_values, inflow=self._inflow_cells(trace, evalf)) + 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._dt_history[0] = dt + if self._n_solves_completed < self.order: + self._n_solves_completed += 1 + + def commit_flux_to_history(self, flux, verbose=False): + """Read the new flux at the launch points, and leave its fit in the + slot until the next carry, as the other flavours do.""" + self._launch_values = self._evaluate_at_launch(flux) + self._fit_arrivals(self._launch, self._launch_values, cell=self._launch_cell) + self._history_committed = True diff --git a/src/underworld3/systems/navier_stokes_eulerian.py b/src/underworld3/systems/navier_stokes_eulerian.py index 98763dd0c..bcce239f3 100644 --- a/src/underworld3/systems/navier_stokes_eulerian.py +++ b/src/underworld3/systems/navier_stokes_eulerian.py @@ -175,9 +175,9 @@ def __init__( ): if DFDt is not None: raise ValueError( - "AdvDiffusion-style Navier-Stokes carries no stress history: " - "the viscous stress at earlier levels is rebuilt from the stored " - "velocity. Do not pass DFDt." + "Do not pass DFDt: a stress history is created by assigning a " + "constitutive model that asks for one (a viscoelastic model), and " + "`stress_transport` chooses the flavour that carries it." ) if restore_points_func is not None: warnings.warn( @@ -409,13 +409,38 @@ def _viscous_stress(self, u_row): return 2 * eta * sympy.Matrix(self.mesh.vector.strain_tensor(u_row)) def _viscous_flux(self): + r"""The flux of the time scheme: the new stress, blended with the stored levels. + + The theta rule weights the flux across time levels, so it needs the stress + at the stored levels as well as the new one. For a viscous fluid that is + rebuilt from the stored velocity, :math:`2\eta\dot\varepsilon(\mathbf{u}^{n})`. + For a viscoelastic one that rebuild is wrong -- the stress there is not a + viscous stress -- and the right object is already held by the stress + history, which is the same quantity the theta rule is asking for. So the + two histories meet here: the flux history of the time scheme IS the + elastic stress history. BDF puts every spatial term at the new level and + the question does not arise. + """ states = self.DuDt.states() weights = self.DuDt.spatial_weights() total = weights[0] * self.stress_deviator - for w, u_k in zip(weights[1:], states[1:]): + stress_history = self.Unknowns.DFDt + for level, (w, u_k) in enumerate(zip(weights[1:], states[1:])): if w == 0: continue - total = total + w * self._viscous_stress(u_k) + if stress_history is None: + total = total + w * self._viscous_stress(u_k) + elif level < len(stress_history.psi_star): + # the history carries the memory part only; a solvent viscosity + # is rebuilt from the stored velocity, as the new level has it + eta_s = getattr(self.constitutive_model.Parameters, "solvent_viscosity", 0) + solvent = 2 * eta_s * sympy.Matrix(self.mesh.vector.strain_tensor(u_k)) + total = total + w * (sympy.Matrix(stress_history.psi_star[level].sym) + solvent) + else: + raise ValueError( + f"the time scheme weights the flux at level {level + 1}, but the " + f"stress history holds {len(stress_history.psi_star)} level(s): " + "give the constitutive model a higher order, or the solver a lower one.") return total def _stabilisation_flux(self): @@ -442,7 +467,7 @@ def F1(self): self.p.sym[0] - self.penalty * self.constitutive_model.K * self.div_u) F1 = public_expression( r"\mathbf{F}_1\left( \mathbf{u} \right)", - self._viscous_flux() - sympy.eye(dim) * mechanical_pressure + self._viscous_flux() + self._devss_flux() - sympy.eye(dim) * mechanical_pressure + self._stabilisation_flux(), "Navier-Stokes SUPG: viscous flux of the time scheme, pressure, tau R (x) a", ) @@ -523,8 +548,18 @@ def solve( # The base _build resolves the preconditioner choice against the mesh # before the SNES reads its options; the setup stages must not be run # directly here (they mark the solver set up first, #683). + carries_stress = self.Unknowns.DFDt is not None + if carries_stress: + # BEFORE the build: the order ramp decides whether the compiled + # functions must be rewired, and the build reads that flag (#727). + self._stress_history_prepare(dt) + self._build(verbose) + if carries_stress: + # AFTER the build, once per step -- not once per pass. + self._stress_history_advance(dt, verbose=verbose, evalf=False) + self._prime_history() u_n = np.array(self.u.array[...]) if self._advection_mode == "extrapolated": @@ -546,6 +581,7 @@ def solve( self, zero_init_guess if k == 0 else False, _force_setup=_force_setup if k == 0 else False, verbose=verbose, picard=0, divergence_retries=divergence_retries, + _skip_stress_history=True, ) # The reductions run on every pass, outside any branch: a rank must # never skip a collective its peers take (tests/test_0052). @@ -564,6 +600,9 @@ def solve( local = float(change.max()) if change.size else 0.0 self._last_change_rate = comm.allreduce(local, op=MPI.MAX) / dt + if carries_stress: + self._stress_history_post_solve(dt, verbose=verbose, evalf=False) + # Shift the extrapolation level, then the history. self._u_prev.array[...] = self.DuDt.psi_star[0].array[...] self.DuDt.update_post_solve(dt, verbose=verbose) diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index c12b206cf..977115898 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -85,6 +85,15 @@ def expression(*args, **kwargs): return public_expression(*args, _unique_name_generation=True, **kwargs) +def _history_psi_fn(constitutive_model): + """The flux a stress history carries: the model's memory part when it + separates one out (a solvent viscosity is rebuilt each step), else the + whole flux. Transposed to the history's row layout.""" + if hasattr(constitutive_model, "history_flux"): + return constitutive_model.history_flux.T + return constitutive_model.flux.T + + def _as_scalar(value): """Collapse a zero-dimensional array to a plain scalar, leave the rest. @@ -449,6 +458,7 @@ def _invalidate_solution_cache(u): from .ddt import SemiLagrangian as SemiLagrangian_DDt 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 @@ -1122,7 +1132,7 @@ def solve( if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - self.DFDt.psi_fn = self.constitutive_model.flux.T + self.DFDt.psi_fn = _history_psi_fn(self.constitutive_model) if not self.is_setup: self._setup_pointwise_functions(verbose) @@ -1439,6 +1449,12 @@ class SNES_Stokes(_ConstitutiveModelStateMixin, SNES_Stokes_SaddlePt): >>> stokes.bodyforce = [0, -1] # gravity >>> stokes.solve() """ + #: DEVSS lag iterations on a plain Stokes solve: D is lagged data, so it has to + #: catch up with the strain rate before the added and subtracted terms cancel + #: (#754). Two passes suffice on a linear problem; the cap bounds a nonlinear one. + _DEVSS_MAX_LAG_ITERATIONS = 4 + _DEVSS_LAG_TOLERANCE = 1.0e-8 + instances = 0 @@ -1470,6 +1486,15 @@ def __init__( self._Estar = None + # DEVSS (Guenette & Fortin 1995): an artificial viscosity added to the + # momentum flux on the velocity and subtracted through a projected + # strain rate, so the velocity always sees an elliptic operator when the + # stress is elastic and lives in a coarser, discontinuous space. Off + # unless devss_viscosity is set; see _devss_flux. + self._devss_viscosity = None + self._devss_D = None + self._devss_flat = None + self._devss_projection = None self._penalty = expression(R"\uplambda", 0, "Numerical Penalty") # Whether `penalty` still holds its automatic value. An explicit # assignment latches this False and the auto default stands down -- @@ -1534,12 +1559,226 @@ def set_jacobian_F1_source(self, F1_source, linesearch="cp"): if F1_source is not None and linesearch is not None: self.petsc_options["snes_linesearch_type"] = linesearch + @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 + 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. + """ + return getattr(self, "_stress_transport", "semi_lagrangian") + + @stress_transport.setter + def stress_transport(self, value): + value = str(value) + if value not in ("semi_lagrangian", "integration_point", "forward", + "lagrangian", "eulerian"): + raise ValueError( + "stress_transport must be 'semi_lagrangian', 'integration_point', " + f"'forward', '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 " + "constitutive model that asks for one.") + self._stress_transport = value + + # ----- DEVSS: stabilising a discontinuous elastic stress ----- + + @property + def devss_viscosity(self): + r"""Artificial viscosity :math:`\eta_a` of the DEVSS split, or ``None``. + + With a transported elastic stress the momentum flux is + :math:`2\eta_{\mathrm{eff}}\dot\varepsilon(\mathbf u) + c\,\sigma^*`, and as the + Weissenberg number rises the elastic term dominates: the velocity is + driven by the divergence of a field that is coarser than + :math:`\dot\varepsilon(\mathbf u)`, discontinuous across elements, and not + derived from the velocity. The discrete operator loses its viscous + character and the stress's inter-element jumps force the velocity at + mesh scale with nothing to damp them -- measured on the viscoelastic + cylinder with an integration-point stress history, 14% of the stress + by rms, every step. + + DEVSS (Guenette & Fortin 1995) adds and subtracts a viscous term, + + .. math:: + \mathbf F_1 \mathrel{+}= 2\eta_a\,(\dot\varepsilon(\mathbf u) - \mathbf D), + + with :math:`\mathbf D` the projection of the strain rate onto the + stress history's space, lagged one step. At convergence the two cancel + to projection error, so nothing physical is added; but the first term + is implicit in the velocity and the second is data, so a mesh-scale + velocity response -- which the projection does not carry -- sees the + full viscosity :math:`\eta_a` while smooth modes see none. It acts on + the RESPONSE to the discontinuous stress, not on the stress, so the + history keeps everything it carries. A natural value for our + discretised Maxwell element is the viscosity the time discretisation + removed, :math:`\eta - \eta_{\mathrm{eff}}`. + """ + return self._devss_viscosity + + @devss_viscosity.setter + def devss_viscosity(self, value): + self._devss_viscosity = value + # The flux expression changes: the compiled functions must be rewired. + self._needs_function_rewire = True + cm = getattr(self, "_constitutive_model", None) + if cm is not None: + cm._solver_is_setup = False + + def _devss_space(self): + """Build the projected strain rate and its projection on first use, in + the stress history's space (degree u-1, continuous).""" + if self._devss_D is not None: + return + dim = self.mesh.dim + columns = [(i, j) for i in range(dim) for j in range(i, dim)] + self._devss_columns = columns + self._devss_D = uw.discretisation.MeshVariable( + f"devss_D_{self.instance_number}", self.mesh, (dim, dim), + vtype=uw.VarType.SYM_TENSOR, degree=self.u.degree - 1, continuous=True, + varsymbol=rf"{{\mathbf{{D}}_{{{self.instance_number}}}}}") + self._devss_flat = uw.discretisation.MeshVariable( + f"devss_flat_{self.instance_number}", self.mesh, (1, len(columns)), + vtype=uw.VarType.MATRIX, degree=self.u.degree - 1, continuous=True) + self._devss_projection = SNES_MultiComponent_Projection( + self.mesh, u_Field=self._devss_flat, n_components=len(columns), + verbose=self.verbose) + self._devss_projection.smoothing = 0.0 + + def _devss_flux(self): + """The DEVSS term of the momentum flux, or a zero matrix when off.""" + dim = self.mesh.dim + if self._devss_viscosity is None: + return sympy.zeros(dim, dim) + self._devss_space() + return 2 * self._devss_viscosity * ( + sympy.Matrix(self.strainrate) - sympy.Matrix(self._devss_D.sym)) + + def _devss_refresh(self, verbose=False): + """Lag the projected strain rate: D <- projection of strain(u) now.""" + if self._devss_viscosity is None: + return + self._devss_space() + E = sympy.Matrix(self.strainrate) + self._devss_projection.uw_function = sympy.Matrix( + [[E[i, j] for (i, j) in self._devss_columns]]) + self._devss_projection.solve(verbose=verbose) + for k, (i, j) in enumerate(self._devss_columns): + values = self._devss_flat.array[:, 0, k] + self._devss_D.array[:, i, j] = values + if i != j: + self._devss_D.array[:, j, i] = values + + def _stress_history_prepare(self, timestep, _force_setup=False): + """Set the elastic timestep and the flags a rebuild depends on. + + Runs BEFORE the solver is built. The effective order of the stress + history ramps over the opening steps, and when it changes the compiled + functions must be rewired; setting that flag after the build has already + decided whether to set up leaves the managed multigrid block asking a + preconditioner for sub-solvers it has not created (#727). + """ + # dt_elastic must always equal the solve timestep. The constitutive + # model's VE formulas (eta_eff, stress history terms) all reference + # Parameters.dt_elastic. If it differs from the actual timestep, + # the stress computation is inconsistent with the time integration. + self.constitutive_model.Parameters.dt_elastic = timestep + # The integrator coefficients must be current BEFORE the history is + # first carried: a trace-back history initialises its first level from + # the constitutive flux of the velocity it finds, and with the + # exponential integrator that flux read the viscous-limit coefficients + # (alpha = phi = 0) until the update that used to follow the carry + # (#740). The BDF coefficients are refreshed again after the carry, as + # before, since they read the step history the carry updates. + self.constitutive_model._update_history_coefficients() + + if _force_setup: + self._needs_function_rewire = True + + # Re-setup when effective_order changes (DDt history ramp-up) + _current_eff_order = self.constitutive_model.effective_order + if _current_eff_order != self._prev_effective_order: + self._needs_function_rewire = True + self.constitutive_model._solver_is_setup = False + self._prev_effective_order = _current_eff_order + + if not self.constitutive_model._solver_is_setup: + self._needs_function_rewire = True + self.DFDt.psi_fn = _history_psi_fn(self.constitutive_model) + # D starts from the velocity as it is now, so the DEVSS pair + # cancels on the first step as it does on every later one. + self._devss_refresh() + + def _stress_history_advance(self, timestep, verbose=False, evalf=False): + """Carry the stress history to where the momentum solve will read it. + + Runs AFTER the solver is built, and once per step: a solver that takes + several passes over the momentum equation (the Navier-Stokes one, with + its Picard corrections) must not advance the history once per pass. + """ + if uw.mpi.rank == 0 and verbose: + print("Stokes solver - carry the stress history", flush=True) + + self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, + store_result=False) + # Uniform pre-solve coefficient hook: VEP delegates to + # _update_bdf_coefficients(); MaxwellExponentialFlowModel updates + # α, φ on the DDt via _update_exp_coefficients(). No isinstance + # checks at the solver layer. + self.constitutive_model._update_history_coefficients() + + def _stress_history_post_solve(self, timestep, verbose=False, evalf=False): + """Commit the stress the solve produced and shift the history levels.""" + # The history manager places the new stress in level 0 and shifts the + # levels. A particle-carried history does that itself in its post-solve, + # by evaluating the new stress at its own particles. + if not self.DFDt.commits_flux_in_post_solve: + self.DFDt.commit_flux_to_history( + getattr(self.constitutive_model, "history_flux", self.constitutive_model.flux), + verbose=verbose) + + self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf) + self._devss_refresh(verbose=verbose) + + # Uniform post-solve hook for any extra integrator-state storage. + # VEP: no-op. ETD-2 / MaxwellExponentialFlowModel: refresh + # forcing_star with current ε̇^{n+1} so the next step's history + # term has access to ε̇ⁿ. + self.constitutive_model._update_history_post_solve() + + self.is_setup = True + self.constitutive_model._solver_is_setup = True + def _create_stress_history_ddt(self, order=2): """Create DFDt for stress history tracking (VE/VEP models). Called automatically when a constitutive model with ``requires_stress_history = True`` is assigned. Can also be called explicitly to pre-create the DFDt with a specific order. + :attr:`stress_transport` chooses which flavour carries it. Constitutive models can inject extra SemiLagrangian kwargs via the ``stress_history_ddt_kwargs`` property — used e.g. by @@ -1556,10 +1795,7 @@ def _create_stress_history_ddt(self, order=2): if cm is not None: ddt_kwargs = dict(getattr(cm, "stress_history_ddt_kwargs", {})) - self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( - self.mesh, - sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), - self.u.sym, + common = dict( vtype=uw.VarType.SYM_TENSOR, degree=self.u.degree - 1, continuous=True, @@ -1568,8 +1804,88 @@ def _create_stress_history_ddt(self, order=2): bcs=None, order=order, smoothing=0.0001, - **ddt_kwargs, ) + if self.stress_transport == "integration_point": + 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( + 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": + 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( + self.mesh, + sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), + self.u.sym, + vtype=common["vtype"], varsymbol=common["varsymbol"], order=order, + ) + 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.") + # 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 " + "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.") + # 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. + # A cells proxy (discontinuous, reconstructed from the particles in + # each cell) is what the particle history is validated on + # (test_0070); the continuous nodal store the mesh flavours use is + # not the right target for a swarm-carried field. + lag_common = {k: v for k, v in common.items() if k != "continuous"} + self.Unknowns.DFDt = uw.systems.ddt.Lagrangian( + self.mesh, + sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), + self.u.sym, + continuous=False, + proxy_location="cells", + **lag_common, + ) + elif self.stress_transport == "eulerian": + if ddt_kwargs: + 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.") + self.Unknowns.DFDt = uw.systems.ddt.EulerianSUPG( + self.mesh, + sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), + self.u.sym, + transport_on_update=True, + **common, + ) + else: + self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( + self.mesh, + sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), + self.u.sym, + **ddt_kwargs, + **common, + ) # Stress flux = 2·viscosity·E_eff references psi_star[0] in E_eff's # history term — without snapshot substitution the projection of # flux→psi_star[0] becomes implicit in psi_star[0] and Min-mode at @@ -1591,6 +1907,7 @@ def solve( order=None, picard: int = 0, divergence_retries: int = 0, + _skip_stress_history: bool = False, homotopy: bool = False, homotopy_options: dict = None, ): @@ -1684,7 +2001,9 @@ def solve( # bet is confirmed after setup, on either path. self._apply_automatic_penalty() - has_stress_history = self.Unknowns.DFDt is not None + # A solver that takes several passes over the momentum equation drives + # the stress history itself, once per step, and asks to be left alone. + has_stress_history = self.Unknowns.DFDt is not None and not _skip_stress_history if has_stress_history: if timestep is None: @@ -1693,46 +2012,18 @@ def solve( "Call stokes.solve(timestep=dt)" ) - # dt_elastic must always equal the solve timestep. The constitutive - # model's VE formulas (eta_eff, stress history terms) all reference - # Parameters.dt_elastic. If it differs from the actual timestep, - # the stress computation is inconsistent with the time integration. - self.constitutive_model.Parameters.dt_elastic = timestep + self._stress_history_prepare(timestep, _force_setup=_force_setup) if order is None or order > self._order: order = self._order - if _force_setup: - self._needs_function_rewire = True - - # Re-setup when effective_order changes (DDt history ramp-up) - _current_eff_order = self.constitutive_model.effective_order - if _current_eff_order != self._prev_effective_order: - self._needs_function_rewire = True - self.constitutive_model._solver_is_setup = False - self._prev_effective_order = _current_eff_order - - if not self.constitutive_model._solver_is_setup: - self._needs_function_rewire = True - self.DFDt.psi_fn = self.constitutive_model.flux.T - if not self.is_setup: self._setup_pointwise_functions(verbose) self._setup_discretisation(verbose) self._setup_solver(verbose) self._check_velocity_preconditioner() - # 1. ADVECT stress history along characteristics - if uw.mpi.rank == 0 and verbose: - print(f"Stokes solver - advect stress history", flush=True) - - self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, - store_result=False) - # Uniform pre-solve coefficient hook: VEP delegates to - # _update_bdf_coefficients(); MaxwellExponentialFlowModel updates - # α, φ on the DDt via _update_exp_coefficients(). No isinstance - # checks at the solver layer. - self.constitutive_model._update_history_coefficients() + self._stress_history_advance(timestep, verbose=verbose, evalf=evalf) # 2. SOLVE if uw.mpi.rank == 0 and verbose: @@ -1751,53 +2042,7 @@ def solve( if uw.mpi.rank == 0 and verbose: print(f"Stokes solver - store stress and shift history", flush=True) - # A particle-carried history (Lagrangian_Swarm) evaluates the new - # stress at its particles and shifts its own chain in - # update_post_solve; the projection and shift below are the - # nodal semi-Lagrangian history's. - if isinstance(self.DFDt, SemiLagrangian_DDt): - _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) - - if getattr(self.DFDt, '_psi_star_use_multicomponent', False): - # Multi-component projection of flux → psi_star[0]. - # - # The DFDt's source-snapshot machinery (enabled once in - # _create_stress_history_ddt) intercepts psi_fn assignment - # to substitute psi_star[0] symbols with a frozen - # psi_snapshot variable, refreshed each step in - # update_pre_solve. So the projection's compiled source - # reads from psi_snapshot (not psi_star[0] itself) and is a - # true one-shot Galerkin projection — no implicit - # fixed-point iteration. - self.DFDt._psi_star_projection_solver.smoothing = 0.0 - self.DFDt._psi_star_projection_solver.solve(verbose=verbose) - # Fan flat result back to psi_star[0] tensor variable - for k, (i, j) in enumerate(self.DFDt._psi_star_indep_indices): - vals = self.DFDt._psi_star_flat_var.array[:, 0, k] - self.DFDt.psi_star[0].array[:, i, j] = vals - if i != j: - self.DFDt.psi_star[0].array[:, j, i] = vals - else: - self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux - self.DFDt._psi_star_projection_solver.smoothing = 0.0 - self.DFDt._psi_star_projection_solver.solve(verbose=verbose) - - for i in range(self.DFDt.order - 1, 0, -1): - if i == 1: - self.DFDt.psi_star[i].array[...] = _advected_sigma_star - else: - self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] - - self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf) - - # Uniform post-solve hook for any extra integrator-state storage. - # VEP: no-op. ETD-2 / MaxwellExponentialFlowModel: refresh - # forcing_star with current ε̇^{n+1} so the next step's history - # term has access to ε̇ⁿ. - self.constitutive_model._update_history_post_solve() - - self.is_setup = True - self.constitutive_model._solver_is_setup = True + self._stress_history_post_solve(timestep, verbose=verbose, evalf=evalf) else: # Plain Stokes — no stress history @@ -1809,6 +2054,42 @@ def solve( time=time, divergence_retries=divergence_retries, ) + # DEVSS adds 2 eta_a (E - D) and relies on D tracking E for the pair to + # cancel, leaving only the unrepresentable part of the strain rate. The + # catch-up lived ONLY in the stress-history post-solve, so a plain Stokes + # solve never refreshed D: the term stayed a bare 2 eta_a E and the run + # silently used eta + eta_a, forever (#754). + # + # Refreshing after the solve is not enough on its own: D starts at zero, so + # the FIRST solve is still wrong and only a second call would be right. A + # steady problem is solved once. So lag-iterate here — refresh D from the + # velocity just found and solve again — until the pair has settled. On a + # time-stepping viscoelastic run the same catch-up happens across steps and + # this loop exits after its check solve. + # A composed solver (Navier-Stokes) makes several passes through here + # per step and refreshes D itself in its post-solve; the lag loop + # here is for the single-pass solve. + if self._devss_viscosity is not None and not _skip_stress_history: + for _ in range(self._DEVSS_MAX_LAG_ITERATIONS): + before = np.array(self._devss_D.array, copy=True) + self._devss_refresh(verbose=verbose) + # the exit test is collective: every rank must take the same + # number of solves, and a rank may hold no dofs at all + D = np.asarray(self._devss_D.array) + moved = float(np.abs(D - before).max()) if D.size else 0.0 + scale = float(np.abs(D).max()) if D.size else 0.0 + moved = uw.mpi.comm.allreduce(moved, op=uw.MPI.MAX) + scale = max(uw.mpi.comm.allreduce(scale, op=uw.MPI.MAX), 1.0e-30) + super().solve( + zero_init_guess=False, + _force_setup=False, + verbose=verbose, + picard=picard, + time=time, + divergence_retries=divergence_retries, + ) + if moved <= self._DEVSS_LAG_TOLERANCE * scale: + break # Confirm the preconditioner the automatic penalty was chosen for. self._check_velocity_preconditioner() @@ -1851,7 +2132,7 @@ def tau(self): F1 = Template( r"\mathbf{F}_1\left( \mathbf{u} \right)", - lambda self: self.stress, + lambda self: self.stress + self._devss_flux(), r"""Velocity equation flux/stress term (pointwise). The $\mathbf{F}_1$ tensor represents the stress response of the fluid, @@ -4598,7 +4879,7 @@ def solve( if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - self.DFDt.psi_fn = self.constitutive_model.flux.T + self.DFDt.psi_fn = _history_psi_fn(self.constitutive_model) if not self.is_setup: self._setup_pointwise_functions(verbose) @@ -4631,6 +4912,154 @@ def solve( return +class SNES_AdvectionDiffusion_Swarm(SNES_AdvectionDiffusion): + r"""Advection-diffusion with the value history carried on a swarm. + + The particle counterpart of :class:`SNES_AdvectionDiffusion` (semi-Lagrangian, + ``AdvDiffusionSLCN``) and :class:`SNES_AdvectionDiffusion_Composed` + (streamline-upwind, ``AdvDiffusion``). The advected quantity's history rides + on a swarm of material points through a :class:`~underworld3.systems.ddt.Lagrangian_Swarm` + manager: after each solve the scalar is read at the particles and, on the next + advection, carried with them. + + The swarm is **supplied by the caller and advected by its owner** — it is the + material swarm of a coupled model, not a private one created here (creating a + private swarm would defeat the purpose of a particle solver, which is to share + the one swarm every field rides on). This solver does not advect the swarm: in + a coupled model the flow solver advects it once a step; in a standalone run, + advect it yourself before each ``solve()``. + + ``particle_update`` and ``step_averaging`` set how the mesh solution returns to + the particles. The default, PIC with ``step_averaging=1``, gives every particle + the full mesh solution each step, so the diffusion the mesh applied is captured + (a half-blend keeps the particle's old, sharper value and under-diffuses). Use + ``particle_update="flip"`` with ``residual_retention`` near + ``exp(-kappa dt pi^2 / h^2)`` to keep the sub-cell sharpness of a + weakly-diffusing field while still diffusing it correctly. + + Geometry note: where the flow crosses the domain boundary, particles + advecting out are clamped against the wall and pile up in a thin layer, and a + high-degree per-cell projection of that layer overshoots (see ``proxy_degree``). + With the default degree-1 projection the solver holds a rotating square (which + the flow crosses on all four sides) as well as a disc, though a domain the flow + keeps well filled — a disc or annulus under rotation — is still the accurate + choice: there the rotating diffusing Gaussian matches the SLCN and SUPG solvers + over a full revolution. + + Parameters + ---------- + mesh, u_Field, V_fn, order, theta, restore_points_func, verbose + As for :class:`SNES_AdvectionDiffusion`. + swarm : underworld3.swarm.Swarm + The material swarm the history rides on. Required. Pass it EMPTY: the + solver declares its history variable on it, and swarm variables must be + declared before ``swarm.populate()``. Construct the solver, then populate + the swarm, then advect it each step (the flow solver does this in a + coupled model). The solver does NOT advect the swarm and does NOT diffuse + anything until a :class:`~underworld3.constitutive_models.DiffusionModel` + is assigned to ``constitutive_model`` (as for :class:`SNES_AdvectionDiffusion`). + particle_update : {"pic", "flip"}, default "pic" + How the mesh solution updates the particles after a solve. + step_averaging : int, default 1 + PIC blend length; 1 takes the whole mesh solution each step. + residual_retention : float, default 1.0 + FLIP residual scale (1 = full FLIP, 0 = PIC). + proxy_location : {"cells", "nodes", "integration_points"}, default "cells" + Where the swarm history's proxy mesh variable lives. + proxy_degree : int, default 1 + Polynomial degree of the per-cell fit that projects the particle history + onto the mesh. Kept LOW on purpose: a high-degree per-cell least-squares + fit needs its particles to span the cell in every direction, and where + the flow clamps particles into a thin layer against a wall (an outflow + boundary) they lie on a line, so a degree-2 fit is only mildly + ill-conditioned — under the projector's guard — and overshoots the nodal + value, which then feeds the solve and diverges. Degree 1 needs three + spanning points, is well conditioned on the clamped layer, and holds a + rotating square (where the flow crosses the boundary) that degree 2 blows + up. Raise it for extra accuracy only where the flow keeps every cell's + particles well spread (a disc under rotation). + """ + + @timing.routine_timer_decorator + def __init__( + self, + mesh: uw.discretisation.Mesh, + u_Field: uw.discretisation.MeshVariable, + V_fn, + swarm: uw.swarm.Swarm, + order: int = 1, + particle_update: str = "pic", + step_averaging: int = 1, + residual_retention: float = 1.0, + proxy_location: str = "cells", + proxy_degree: int = 1, + restore_points_func: Callable = None, + verbose=False, + theta: float = 0.5, + ): + if swarm is None: + raise ValueError( + "SNES_AdvectionDiffusion_Swarm needs a swarm to carry the history on; " + "pass the material swarm (this solver does not create a private one). " + "Use AdvDiffusionSLCN for the mesh-based semi-Lagrangian scheme.") + if int(step_averaging) < 1: + raise ValueError(f"step_averaging must be >= 1, not {step_averaging!r}") + if abs(float(theta) - 0.5) > 1e-12: + warnings.warn( + "theta only sets the diffusive-flux time integration here; the swarm " + "value history is fixed Crank-Nicolson (theta=0.5). The two will be " + "inconsistent for theta != 0.5.", stacklevel=2) + DuDt = Lagrangian_Swarm_DDt( + swarm=swarm, + psi_fn=u_Field.sym, + vtype=uw.VarType.SCALAR, + degree=proxy_degree, + continuous=u_Field.continuous, + varsymbol=u_Field.symbol, + verbose=verbose, + order=order, + step_averaging=step_averaging, + proxy_location=proxy_location, + proxy_sampling="reconstruct", + particle_update=particle_update, + residual_retention=residual_retention, + ) + super().__init__( + mesh, + u_Field, + V_fn, + order=order, + restore_points_func=restore_points_func, + verbose=verbose, + DuDt=DuDt, + theta=theta, + ) + self.swarm = swarm + self._last_swarm_key = None + self._warned_static_swarm = False + + def _swarm_position_key(self): + c = np.asarray(self.swarm._particle_coordinates.data) + return (c.shape, float(c.sum()), float((c * c).sum())) + + @timing.routine_timer_decorator + def solve(self, *args, **kwargs): + # The solver does not advect the swarm (its owner does). If the swarm has + # not moved since the last solve, the material history is not being + # transported: warn once rather than return a plausible, un-advected field. + key = self._swarm_position_key() + if (self._last_swarm_key is not None and key == self._last_swarm_key + and not self._warned_static_swarm): + warnings.warn( + "AdvDiffusionSwarm: the swarm has not moved since the last solve, so the " + "material history is not being transported. Advect the swarm before each " + "solve() (swarm.advection(V, dt) in a standalone run; the flow solver in a " + "coupled model). This warning is issued once.", stacklevel=2) + self._warned_static_swarm = True + self._last_swarm_key = key + return super().solve(*args, **kwargs) + + class SNES_Diffusion(SNES_Scalar): r""" Diffusion equation solver using mesh-based finite elements. @@ -4912,7 +5341,7 @@ def solve( if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - self.DFDt.psi_fn = self.constitutive_model.flux.T + self.DFDt.psi_fn = _history_psi_fn(self.constitutive_model) # self._flux = self.constitutive_model.flux.T # self._flux_star = self._flux.copy() @@ -5129,7 +5558,7 @@ def F1(self): if DFDt is not None: # We can flag to only do this if the constitutive model has been updated - DFDt.psi_fn = self._constitutive_model.flux.T + DFDt.psi_fn = getattr(self._constitutive_model, 'history_flux', self._constitutive_model.flux).T F1 = expression( r"\mathbf{F}_1\left( \mathbf{u} \right)", @@ -5331,7 +5760,18 @@ def solve( if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - self.DFDt.psi_fn = self.constitutive_model.flux.T + self.DFDt.psi_fn = _history_psi_fn(self.constitutive_model) + + # A viscoelastic constitutive model integrates its stress over the + # solve step: it has to be told the step, and its integrator + # coefficients refreshed, exactly as the Stokes family does in + # _stress_history_prepare / _stress_history_advance. Without this the + # memory term was silently absent here (dt_elastic never set) and the + # exponential integrator ran in its viscous limit (#741). + _cm = self.constitutive_model + if getattr(_cm, "requires_stress_history", False) and hasattr(_cm.Parameters, "dt_elastic"): + _cm.Parameters.dt_elastic = timestep + _cm._update_history_coefficients() if not self.is_setup: self._setup_pointwise_functions(verbose) @@ -5353,6 +5793,8 @@ def solve( self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=_evalf) if trace is not None: trace.finish_step() + if getattr(_cm, "requires_stress_history", False): + _cm._update_history_coefficients() # BDF reads the carried step history # Override AM coefficients if flux_order is explicitly set if self._flux_order is not None: diff --git a/tests/parallel/test_1062_forward_stress_history_mpi.py b/tests/parallel/test_1062_forward_stress_history_mpi.py new file mode 100644 index 000000000..db421412f --- /dev/null +++ b/tests/parallel/test_1062_forward_stress_history_mpi.py @@ -0,0 +1,59 @@ +"""The forward stress history in parallel: np >= 2 equals serial. + +A Maxwell fluid whose shear modulus varies in x, driven by a body force that +turns over two counter-rotating cells between no-slip walls, on a gmsh +(file-read) mesh. The stress is then non-uniform in both directions, so an +arrival handed across a seam with the wrong value, fitted into the wrong +cell, or dropped shows in the answer; and the flow crosses a seam of either +orientation. The stress at two points after ten steps must be the serial +value, and arrivals must actually have changed rank. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +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 +POINTS = np.array([[0.5, 0.2], [-0.3, -0.35]]) + + +def turned_over_maxwell_box(transport="forward", 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) + x, _y = mesh.X + v = uw.discretisation.MeshVariable("U_to", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("P_to", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=1) + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.constitutive_model.Parameters.shear_modulus = 1.0 + 0.5 * sympy.sin(sympy.pi * x / Lx) + stokes.constitutive_model.Parameters.dt_elastic = dt + for wall in ("Top", "Bottom", "Left", "Right"): + stokes.add_dirichlet_bc((0.0, 0.0), wall) + stokes.bodyforce = sympy.Matrix([[0.0, 4.0 * sympy.sin(sympy.pi * x / Lx)]]) + stokes.tolerance = 1.0e-8 + relocated = 0 + for _ in range(steps): + stokes.solve(timestep=dt, zero_init_guess=False) + relocated += getattr(stokes.DFDt, "_n_relocated", 0) + values = np.asarray(uw.function.global_evaluate(stokes.DFDt.psi_star[0].sym[0, 1], POINTS)).reshape(-1) + return type(stokes.DFDt).__name__, values, uw.mpi.comm.allreduce(relocated, op=uw.MPI.SUM) + + +def test_the_forward_history_gives_the_serial_stress_on_every_rank(): + kind, values, relocated = turned_over_maxwell_box() + assert kind == "ForwardSemiLagrangian" + # 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 + # the seam must actually have been crossed for this to test the parallel path + assert relocated > 0 diff --git a/tests/test_0074_return_to_bounds_on_a_file_mesh.py b/tests/test_0074_return_to_bounds_on_a_file_mesh.py new file mode 100644 index 000000000..c7a0d0cb4 --- /dev/null +++ b/tests/test_0074_return_to_bounds_on_a_file_mesh.py @@ -0,0 +1,33 @@ +"""A mesh read from a file has no analytic closure for return_coords_to_bounds. + +Every gmsh geometry (the benchmark meshes) is read from a file. Before the fix the +property returned None for them, so a trace-back foot leaving through an inlet was +never restored to the boundary and fell through to the evaluator's distance-weighted +fallback. The general facet restore already handles the case and leaves interior +points untouched, so it is the fallback. +""" +import os +import numpy as np +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def test_a_file_mesh_restores_outside_points_to_its_boundary(tmp_path): + path = os.path.join(str(tmp_path), "box.msh") + uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=0.25, filename=path) + mesh = uw.discretisation.Mesh(path) + assert mesh._analytic_return_coords_to_bounds is None + restore = mesh.return_coords_to_bounds + assert restore is not None + + pts = np.array([[1.3, 0.5], [0.5, 0.5], [-0.2, 1.4]]) + out = np.asarray(restore(pts)) + # outside points land on the boundary (the facet restore steps 1e-4 inside) + assert abs(out[0, 0] - 1.0) < 1.0e-3 and abs(out[0, 1] - 0.5) < 1.0e-9 + assert abs(out[2, 0] - 0.0) < 1.0e-3 and abs(out[2, 1] - 1.0) < 1.0e-3 + # an interior point is untouched + assert np.array_equal(out[1], [0.5, 0.5]) diff --git a/tests/test_1059_stress_transport.py b/tests/test_1059_stress_transport.py new file mode 100644 index 000000000..f9d14cc55 --- /dev/null +++ b/tests/test_1059_stress_transport.py @@ -0,0 +1,740 @@ +"""Transporting a stress history on the grid instead of along characteristics. + +A viscoelastic solve carries a stress that is not its own unknown. The history +manager owns that transport: the semi-Lagrangian flavour traces it back along +characteristics, and ``EulerianSUPG`` with ``transport_on_update`` assembles it +implicitly with the same SUPG stabilisation the solvers use. These are the +transport tests the existing viscoelastic benchmarks cannot give us, because +those are all spatially uniform and so transport nothing. + +Run: pixi run python -m pytest tests/test_1059_stress_transport.py -v +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +import itertools + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] # several hundred solves: minutes + +_uid = itertools.count() # every helper names its variables uniquely: a name reused on a new mesh is a trap + +# BASELINES for the cases without a closed form, filled from the baseline run (see the ledger) +VARYING_SL2, VARYING_EU2, VARYING_EU2_REFINED = 0.87325, 0.87206, 0.87563 +NS_ORDER2_XY, CN_ORDER1_XY = 0.70726, 0.64668 + +COMPONENTS = ((0, 0), (0, 1), (1, 1)) +AMPLITUDES = (1.0, 0.5, -1.0) + + +def _blob(x, y, x0, y0=0.0, width=0.04): + return sympy.exp(-((x - x0) ** 2 + (y - y0) ** 2) / width) + + +def _plant(var, expression): + """Set every independent component of a symmetric tensor from one shape.""" + values = uw.function.evaluate(expression, var.coords).reshape(-1) + for (i, j), amplitude in zip(COMPONENTS, AMPLITUDES): + var.array[:, i, j] = amplitude * values + if i != j: + var.array[:, j, i] = amplitude * values + return values + + +def _error(var, expression): + """Relative L2 error of the tensor against one shape times the amplitudes.""" + exact = uw.function.evaluate(expression, var.coords).reshape(-1) + numerator = denominator = 0.0 + for (i, j), amplitude in zip(COMPONENTS, AMPLITUDES): + target = amplitude * exact + numerator += float(np.sum((np.asarray(var.array[:, i, j]) - target) ** 2)) + denominator += float(np.sum(target ** 2)) + return np.sqrt(numerator / denominator) + + +def _grid_manager(mesh, tag, velocity, order=1): + """A stress variable and the Eulerian-SUPG history that transports it.""" + stress = uw.discretisation.MeshVariable( + f"S_{tag}", mesh, vtype=uw.VarType.SYM_TENSOR, degree=2) + return stress, uw.systems.ddt.EulerianSUPG( + mesh, stress, velocity, vtype=uw.VarType.SYM_TENSOR, degree=2, + continuous=True, order=order, transport_on_update=True) + + +def _traced_manager(mesh, tag, velocity, order=1): + """The same, with the semi-Lagrangian trace-back carrying the stress. + + Each manager needs its OWN stress variable: they read the field back on + every step, so sharing one couples the two transports. + """ + stress = uw.discretisation.MeshVariable( + f"S_{tag}", mesh, vtype=uw.VarType.SYM_TENSOR, degree=2) + return stress, uw.systems.ddt.SemiLagrangian( + mesh, stress.sym, velocity, vtype=uw.VarType.SYM_TENSOR, degree=2, + continuous=True, order=order, varsymbol=rf"S^{{{tag}}}") + + +def _march(stress, manager, dt, steps): + """The step a solver takes: carry the history, then commit the new value. + + With no constitutive update the commit is the identity, which is what + makes this a pure transport comparison. + """ + for _ in range(steps): + manager.update_pre_solve(dt) + stress.array[...] = manager.psi_star[0].array[...] + + +def test_a_stress_blob_is_carried_by_a_uniform_flow(): + """Uniform translation has an exact answer: the blob arrives where the flow + puts it, with its components unchanged (no rotation in a uniform flow).""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 24, qdegree=3) + x, y = mesh.X + speed, dt, steps, start = 0.5, 0.05, 8, -0.5 + velocity = sympy.Matrix([[speed, 0.0]]) + grid_stress, eulerian = _grid_manager(mesh, "uniform_g", velocity) + traced_stress, lagrangian = _traced_manager(mesh, "uniform_s", velocity) + assert eulerian.transport_on_update and eulerian._advection_mode == "assembled" + + for stress, manager in ((grid_stress, eulerian), (traced_stress, lagrangian)): + _plant(stress, _blob(x, y, start)) + manager.initialise_history() + _march(stress, manager, dt, steps) + + exact = _blob(x, y, start + speed * dt * steps) + grid = _error(eulerian.psi_star[0], exact) + traced = _error(lagrangian.psi_star[0], exact) + # Uniform translation is the trace-back's best case (the departure point is + # exact), so it sets the bar; the grid transport must be of the same order. + assert traced < 0.02, traced + assert grid < 0.05, grid + # negative control: the field really moved + assert _error(eulerian.psi_star[0], _blob(x, y, start)) > 0.5 + # symmetry is preserved component by component + carried = np.asarray(eulerian.psi_star[0].array) + assert np.allclose(carried[:, 0, 1], carried[:, 1, 0]) + + +def test_the_components_ride_round_a_rigid_rotation_unchanged(): + """With no rotation terms in the transport the tensor components are + advected as scalars: after a quarter turn each component sits where the + flow carried it, with the same amplitude.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 24, qdegree=3) + x, y = mesh.X + stress, eulerian = _grid_manager(mesh, "rot", sympy.Matrix([[-y, x]])) + + radius, dt, steps = 0.5, np.pi / 2 / 40, 40 # a quarter turn + _plant(stress, _blob(x, y, radius)) + eulerian.initialise_history() + _march(stress, eulerian, dt, steps) + + turned = _blob(x, y, 0.0, radius) # a quarter turn from (r, 0) + assert _error(eulerian.psi_star[0], turned) < 0.2 + peak = float(np.abs(np.asarray(eulerian.psi_star[0].array[:, 0, 0])).max()) + assert 0.6 < peak < 1.05, peak # amplitude carried, not created + + +def test_transport_is_off_unless_asked_for(): + """A manager built for a solver's own unknown assembles its advection in + that solver's residual and must not also move its history.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 8, qdegree=3) + x, y = mesh.X + T = uw.discretisation.MeshVariable("T_off", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(_blob(x, y, 0.3), T.coords).reshape(-1) + manager = uw.systems.ddt.EulerianSUPG( + mesh, T, sympy.Matrix([[1.0, 0.0]]), vtype=uw.VarType.SCALAR, + degree=2, continuous=True) + manager.initialise_history() + before = np.array(manager.psi_star[0].array) + manager.update_pre_solve(0.05) + assert manager.transport_on_update is False + assert np.array_equal(np.asarray(manager.psi_star[0].array), before) + assert manager._transport_solver is None + + +def _maxwell_shear(transport, order, steps=20, dt=0.1, integrator="bdf", solver="stokes", initial_velocity=False, + objective_rate="none", solvent=0.0): + """The analytic Maxwell shear box, with the stress history of one's choosing. + + Simple shear of a Maxwell material: sigma_xy = eta gammadot (1 - exp(-t/t_r)). + The stress is spatially uniform, so its transport is a no-op and the two + flavours must agree exactly. That is what makes this the correctness check + on the plumbing rather than on the transport. + """ + eta = shear_modulus = 1.0 + speed, height, width = 0.5, 1.0, 2.0 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-width / 2, -height / 2), + maxCoords=(width / 2, height / 2)) + tag = f"{transport[0]}{order}_{next(_uid)}" + v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1) + if solver == "stokes": + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, verbose=False) + else: + # The trace-back Navier-Stokes solver at negligible inertia: the same box. + stokes = uw.systems.NavierStokesSLCN(mesh, v, p, rho=1.0e-6, order=1) + stokes.bodyforce = sympy.Matrix([[0.0, 0.0]]) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=order, integrator=integrator, objective_rate=objective_rate) + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.shear_modulus = shear_modulus + stokes.constitutive_model.Parameters.solvent_viscosity = solvent + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((speed, 0.0), "Top") + stokes.add_dirichlet_bc((-speed, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + + if initial_velocity: + # The steady shear profile, already in place before the first solve. + v.array[:, 0, :] = uw.function.evaluate( + sympy.Matrix([[2.0 * speed * mesh.X[1] / height, 0.0]]), v.coords).reshape(-1, 2) + for _ in range(steps): + stokes.solve(timestep=dt, zero_init_guess=False) + + # After a solve the Stokes family has committed the new stress into the + # history's first level; the trace-back Navier-Stokes solver records it at + # the NEXT carry, so there its first level still holds the previous step + # and the stress just solved for is the constitutive flux (#742). + latest = stokes.DFDt.psi_star[0].sym if solver == "stokes" else stokes.constitutive_model.flux + origin = np.array([[0.0, 0.0]]) + stress = float(np.asarray(uw.function.evaluate(latest[0, 1], origin)).reshape(-1)[0]) + rate = 2.0 * speed / height + exact = eta * rate * (1.0 - np.exp(-steps * dt * shear_modulus / eta)) + if objective_rate != "none" or solvent: + n1 = float(np.asarray(uw.function.evaluate(latest[0, 0] - latest[1, 1], origin)).reshape(-1)[0]) + # the momentum flux and the polymer part of it, both formed on the + # just-committed level (one step ahead of `latest`, see #742): their + # difference is the solvent's stress alone + total_xy = float(np.asarray(uw.function.evaluate(stokes.constitutive_model.flux[0, 1], origin)).reshape(-1)[0]) + polymer_flux_xy = float(np.asarray(uw.function.evaluate(stokes.constitutive_model.history_flux[0, 1], origin)).reshape(-1)[0]) + return type(stokes.DFDt).__name__, stress, exact, n1, total_xy - polymer_flux_xy + return type(stokes.DFDt).__name__, stress, exact + + +KINDS = { + "semi_lagrangian": "SemiLagrangian", + "integration_point": "IntegrationPointSemiLagrangian", + "forward": "ForwardSemiLagrangian", + "lagrangian": "Lagrangian", + "eulerian": "EulerianSUPG", +} + + +@pytest.mark.parametrize("order, integrator, tolerance", + [(1, "bdf", 0.02), (2, "bdf", 0.002), (1, "etd", 1e-4), (2, "etd", 0.01)]) +def test_every_stress_history_solves_the_maxwell_shear_box(order, integrator, tolerance): + """A Stokes solve carries its viscoelastic stress with the history its + `stress_transport` names, and on a uniform stress every flavour agrees exactly. + + The analytic tolerance is what catches a history that is rebuilt from its + own flux instead of carried: that applies the constitutive update twice a + step and lands at 12.8% on this box at order 1 (#732). The exponential + integrator is exact for a constant strain rate, and it must run on every + flavour, not only the nodal one (#739). + """ + results = {} + 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 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) + kind, stress, exact = _maxwell_shear(transport, order, integrator=integrator) + assert kind == expected_kind + assert abs(stress - exact) / exact < tolerance, (transport, stress, exact) + results[transport] = stress + # Transport is a no-op on a uniform field, so the flavours differ only in + # how they carry it: the plumbing must not add anything of its own. The + # particle flavour reconstructs its history from a swarm and still agrees to + # 2e-8 on this field, so it is held to the same floor as the mesh flavours. + spread = max(results.values()) - min(results.values()) + assert spread < 1e-6 * abs(exact), results + + +def test_stress_transport_is_validated_and_fixed_once_the_history_exists(): + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("U_val", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P_val", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + with pytest.raises(ValueError, match="stress_transport must be"): + stokes.stress_transport = "particles" + stokes.stress_transport = "eulerian" + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1) + stokes.constitutive_model.Parameters.shear_modulus = 1.0 + 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" + + +def _sheared_varying_modulus(transport, order, steps=10, dt=0.1, res=6): + """Simple shear of a Maxwell material whose shear modulus varies in x. + + The stress is then non-uniform and the shear carries it, so the transport + term is genuinely active: the case the uniform benchmarks cannot provide. + There is no closed form, so the schemes are judged against each other and + against their own behaviour under refinement. + """ + eta, modulus, speed, height, width = 1.0, 1.0, 0.5, 1.0, 2.0 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(2 * res, res), minCoords=(-width / 2, -height / 2), + maxCoords=(width / 2, height / 2)) + x, _y = mesh.X + tag = f"{transport[0]}{order}{steps}_{next(_uid)}" + v = uw.discretisation.MeshVariable(f"Uv_{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"Pv_{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, verbose=False) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=order) + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.shear_modulus = ( + modulus * (1 + 0.5 * sympy.sin(2 * sympy.pi * x / width))) + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((speed, 0.0), "Top") + stokes.add_dirichlet_bc((-speed, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + for _ in range(steps): + stokes.solve(timestep=dt, zero_init_guess=False, evalf=False) + + shear = stokes.DFDt.psi_star[0].sym[0, 1] + norm = float(np.sqrt(uw.maths.Integral(mesh, shear ** 2).evaluate())) + slope = float(np.sqrt(uw.maths.Integral(mesh, shear.diff(x) ** 2).evaluate())) + return norm, slope + + +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) + 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 + # is held to its own number, not only to the other + assert abs(traced - VARYING_SL2) < 1.0e-3 * VARYING_SL2, traced + assert abs(grid - VARYING_EU2) < 1.0e-3 * VARYING_EU2, grid + assert abs(grid - traced) / traced < 5.0e-3, (grid, traced) + + # halving the step moves each scheme by no more than they differ from + # each other: the gap between them is discretisation, not a defect + refined, _ = _sheared_varying_modulus("eulerian", 2, steps=20, dt=0.05) + assert abs(refined - VARYING_EU2_REFINED) < 1.0e-3 * VARYING_EU2_REFINED, refined + assert abs(refined - grid) / grid < 5.0e-3, (refined, grid) + + +def test_navier_stokes_carries_a_viscoelastic_stress_either_way(): + """With inertia the momentum solve takes several passes over a step, so the + stress history must advance once per step, not once per pass. Both flavours + must give the same answer; the shear box with inertia has no closed form + (the velocity is far from steady simple shear at this effective viscosity), + so they are judged against each other.""" + def run(transport): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)) + v = uw.discretisation.MeshVariable(f"Un_{transport[0]}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"Pn_{transport[0]}", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=2) + ns.stress_transport = transport + ns.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + ns.Unknowns, order=2) + ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + ns.constitutive_model.Parameters.shear_modulus = 1.0 + ns.constitutive_model.Parameters.dt_elastic = 0.1 + ns.add_dirichlet_bc((0.5, 0.0), "Top") + ns.add_dirichlet_bc((-0.5, 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([[0.0, 0.0]]) + for _ in range(10): + ns.solve(timestep=0.1) + 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") + grid_kind, grid = run("eulerian") + assert traced_kind == "SemiLagrangian" 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) + + +def test_the_theta_rule_takes_its_stored_flux_from_the_stress_history(): + """Crank-Nicolson weights the momentum flux across time levels. For a + viscous fluid the stored level is rebuilt from the stored velocity; for a + viscoelastic one that rebuild is wrong, and the stress history already holds + exactly what the theta rule is asking for. Check the flux reads it.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Uf", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pf", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=1) + ns.constitutive_model = uw.constitutive_models.ViscousFlowModel + ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + viscous_flux = ns._viscous_flux() + + ve = uw.systems.NavierStokes( + mesh, uw.discretisation.MeshVariable("Ug", mesh, mesh.dim, degree=2), + uw.discretisation.MeshVariable("Pg", mesh, 1, degree=1), rho=1.0, order=1) + ve.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + ve.Unknowns, order=1) + ve.constitutive_model.Parameters.shear_modulus = 1.0 + ve.constitutive_model.Parameters.dt_elastic = 0.1 + assert ve.integrator == "am", "order 1 is the theta rule" + stored = set(sympy.Matrix(ve.DFDt.psi_star[0].sym).atoms(sympy.Function)) + assert stored & set(ve._viscous_flux().atoms(sympy.Function)), \ + "the theta rule must read the stress history at the stored level" + assert not (stored & set(viscous_flux.atoms(sympy.Function))), \ + "a viscous fluid has no stress history to read" + + +def test_crank_nicolson_carries_a_viscoelastic_stress_either_way(): + """The scheme the Navier-Stokes benchmarks use, order 1, with a stress + history: both transports must agree.""" + def run(transport): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)) + v = uw.discretisation.MeshVariable(f"Uc_{transport[0]}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"Pc_{transport[0]}", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=1) + ns.stress_transport = transport + ns.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + ns.Unknowns, order=1) + ns.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + ns.constitutive_model.Parameters.shear_modulus = 1.0 + ns.constitutive_model.Parameters.dt_elastic = 0.1 + ns.add_dirichlet_bc((0.5, 0.0), "Top") + ns.add_dirichlet_bc((-0.5, 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([[0.0, 0.0]]) + for _ in range(10): + ns.solve(timestep=0.1) + 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") + # 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) + + +def test_devss_is_off_by_default_and_vanishes_on_a_uniform_strain_rate(monkeypatch): + """DEVSS adds 2 eta_a (edot - D) to the momentum flux with D the projected + strain rate. On the uniform Maxwell shear box D equals edot exactly, so the + term must vanish and the answer must not move for any eta_a; and it must be + 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) + assert abs(off - exact) / exact < 0.02 + + def with_devss(builder, *args, **kw): + # the helpers build their own solver; re-run them with the term on by + # patching the class default for the duration of the call + original = uw.systems.Stokes.__init__ + def patched(self, *a, **k): + original(self, *a, **k) + self.devss_viscosity = 1.0 + with monkeypatch.context() as m: + m.setattr(uw.systems.Stokes, "__init__", patched) + return builder(*args, **kw) + + _kind, on, _ = with_devss(_maxwell_shear, "integration_point", 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) + moved = abs(norm_on - norm_off) / norm_off + assert 1e-6 < moved < 5e-2, moved # live, and only projection-sized + + +def test_devss_cancels_on_a_plain_stokes_solve_too(): + """DEVSS must not change the rheology when there is NO stress history. + + The pair 2 eta_a (edot - D) cancels only if D is refreshed to track edot. That + refresh lived solely in the stress-history post-solve, so a plain Stokes solve + never performed it: D stayed at its initial value and the term was a bare + 2 eta_a edot -- the run silently used eta + eta_a, forever, and repeated solves + did not heal it (#754). + + Pinned against the EXACT fully developed value rather than a comparison, so it + fails the moment the cancellation stops happening: plane Poiseuille between + plates at y = 0, h has dp/dx = -3 eta U / (h/2)^2 with U the mean speed.""" + eta, eta_a, h, u_max = 1.0, 0.25, 1.0, 1.5 + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(4.0, h), cellSize=0.2, qdegree=3) + x, y = mesh.X + v = uw.discretisation.MeshVariable("U_devss", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("P_devss", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.devss_viscosity = eta_a + stokes.add_dirichlet_bc((u_max * (1.0 - ((y - h / 2) / (h / 2)) ** 2), 0.0), "Left") + stokes.add_dirichlet_bc((0.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + + xs = np.linspace(1.5, 3.5, 60) + pts = np.column_stack([xs, np.full_like(xs, h / 2)]) + exact = -3.0 * eta * (2.0 / 3.0 * u_max) / (h / 2) ** 2 + for _ in range(2): # the defect also survived repeated solves + stokes.solve(zero_init_guess=False) + pr = np.asarray(uw.function.evaluate(p.sym[0], pts)).reshape(-1) + grad = float(np.polyfit(xs, pr, 1)[0]) + assert abs(grad - exact) / abs(exact) < 2.0e-3, (grad, exact) + + +def test_the_exponential_integrator_runs_on_the_trace_back_navier_stokes(): + """The trace-back Navier-Stokes solver has its own history path. It must + tell a viscoelastic model the step and refresh the integrator coefficients + 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") + assert abs(stress - exact) / exact < tolerance, (integrator, stress, exact) + + +def test_a_preset_velocity_gives_both_integrators_the_same_first_stress(): + """A trace-back history initialises its first level from the constitutive + flux of the velocity it finds. With a velocity already in place that flux + must be formed with the integrator's coefficients for this step -- the + 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) + # 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 + # ETD-1: eta gdot (1 - e^{-dt/lambda}) (1 + e^{-dt/lambda}) = 1 - e^{-0.2} + # The profile is linear, so P2 holds it exactly; only the solve tolerance is left. + assert abs(bdf - 21.0 / 121.0) < 1.0e-4, bdf + assert abs(etd - (1.0 - np.exp(-0.2))) < 1.0e-4, etd + + +def test_the_integration_point_history_takes_the_inflow_value_at_an_inlet(): + """A departure point that leaves through the inlet is restored to the + boundary by the trace. Left to sample the boundary edge, the + integration-point history on a box channel grew a mode in the inlet cell + column at Wi 1 and 5 (#745); given an inflow value it takes that instead. + Uniform flow into a box carrying zero stress: after k steps the inflow + value has entered a distance k * speed * dt, and no further.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=1 / 24, qdegree=3) + speed, dt, steps = 0.5, 0.05, 8 + velocity = sympy.Matrix([[speed, 0.0]]) + 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( + 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 + if inflow is not None: + manager.inflow_value = inflow + manager.initialise_history() + for _ in range(steps): + manager.update_pre_solve(dt, store_result=False) + manager.commit_flux_to_history(manager.psi_star[0].sym) # identity commit + points = np.asarray(manager.psi_star[0].integration_points).reshape(-1, 2) + xy = np.asarray(manager.psi_star[0].data)[:, manager._components.index((0, 1))] + entered = points[:, 0] < -1.0 + 0.6 * speed * dt * steps + untouched = points[:, 0] > -1.0 + 1.5 * speed * dt * steps + results[tag] = (xy[entered], xy[untouched]) + with_in, with_out = results["with"] + assert np.allclose(with_in, 0.25, atol=0.02), (with_in.min(), with_in.max()) + assert np.abs(with_out).max() < 0.02 + # negative control: sampling the edge of a zero field brings nothing in + without_in, _ = results["without"] + assert np.abs(without_in).max() < 1e-6 + + + + +@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point"]) +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 + N1 = sigma_xx - sigma_yy = 2 eta lambda gdot^2 (1 - e^{-t/lambda}(1 + t/lambda)). + The passive Maxwell element makes no normal stress at all; the objective + rate is what makes it. First order in time on the stretching term, so the + tolerance is loose; the shear stress is unchanged by the term.""" + eta = lam = 1.0; gdot = 1.0; steps, dt = 20, 0.1 + t = steps * dt + _, xy, exact_xy, n1, _ = _maxwell_shear(transport, 1, steps=steps, dt=dt, objective_rate="upper_convected") + n1_exact = 2 * eta * lam * gdot ** 2 * (1 - np.exp(-t / lam) * (1 + t / lam)) + assert abs(xy - exact_xy) / exact_xy < 0.02, (xy, exact_xy) + assert abs(n1 - n1_exact) / n1_exact < 0.10, (n1, n1_exact) + _, _, _, n1_passive, _ = _maxwell_shear(transport, 1, steps=steps, dt=dt, objective_rate="none", solvent=1e-12) + assert abs(n1_passive) < 1e-6 * n1_exact + + +def _oldroyd_channel_first_step(eta_s, dt=0.1): + """Plane Poiseuille flow of an Oldroyd-B fluid from rest, one BDF-1 step: + the stress history is zero, so the momentum equation sees the viscosity + eta_s + eta_p / (1 + lambda / dt) and the centreline speed is f H^2 / 8 + over it. A P2 velocity holds the parabola exactly.""" + eta_p, G, f, H = 1.0, 1.0, 1.0, 1.0 + mesh = uw.meshing.StructuredQuadBox(elementRes=(8, 8), minCoords=(-1.0, -H / 2), maxCoords=(1.0, H / 2)) + tag = next(_uid) + v = uw.discretisation.MeshVariable(f"U_ch_{tag}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P_ch_{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=1) + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_p + stokes.constitutive_model.Parameters.shear_modulus = G + stokes.constitutive_model.Parameters.solvent_viscosity = eta_s + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((0.0, 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.bodyforce = sympy.Matrix([[f, 0.0]]) + stokes.tolerance = 1.0e-10 + stokes.solve(timestep=dt, zero_init_guess=False) + u_c = float(np.asarray(uw.function.evaluate(v.sym[0], np.array([[0.0, 0.0]]))).reshape(-1)[0]) + eta_eff = eta_p / (1.0 + (eta_p / G) / dt) + return u_c, f * H * H / (8.0 * (eta_s + eta_eff)) + + +def test_a_solvent_viscosity_adds_its_newtonian_stress(): + """Oldroyd-B in shear: the total shear stress is the solvent's eta_s gdot at + once plus the polymer's eta_p gdot (1 - e^{-t/lambda}) building up; and the + solvent's stress has to be in the ASSEMBLED momentum flux, which a uniform + 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) + 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 + with_solvent, exact_with = _oldroyd_channel_first_step(eta_s) + without, exact_without = _oldroyd_channel_first_step(0.0) + assert abs(with_solvent - exact_with) < 1.0e-6 * exact_with, (with_solvent, exact_with) + assert abs(without - exact_without) < 1.0e-6 * exact_without, (without, exact_without) + assert without > 5.0 * with_solvent # the solvent term is live in the assembly + + +# ----- Health of the carried stress and the store smoothing (#737, #768) ----- + +def _one_shear_step(transport, dt=1.0, objective_rate="none"): + """One BDF-1 step of the Maxwell shear box from rest, returning the solver. + + From rest there is no history and no stretching, so the first stress is + 2 eta_eff D exactly, with eta_eff = eta G dt / (eta + G dt). With eta = G = 1 + and dt = 1 that is D: the shear stress is gammadot / 2 = 0.5 and the + conformation tau/G + I has eigenvalues 1 +- 0.5. + """ + eta = shear_modulus = 1.0 + speed, height, width = 0.5, 1.0, 2.0 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-width / 2, -height / 2), + maxCoords=(width / 2, height / 2)) + tag = f"{transport[0]}_{next(_uid)}" + v = uw.discretisation.MeshVariable(f"U_h_{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"P_h_{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, verbose=False) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1, integrator="bdf", objective_rate=objective_rate) + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.shear_modulus = shear_modulus + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((speed, 0.0), "Top") + stokes.add_dirichlet_bc((-speed, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-8 + stokes.solve(timestep=dt, zero_init_guess=False) + return stokes + + +@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward"]) +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() + # eta_eff = 1/2 at dt = 1: tau_xy = 0.5, eigenvalues of tau/G + I are 0.5 and 1.5 + assert abs(health["min"] - 0.5) < 1.0e-6 + assert abs(health["max"] - 1.5) < 1.0e-6 + assert health["fraction_negative"] == 0.0 + + +def test_the_conformation_check_sees_a_lost_conformation(): + # dt = 3: eta_eff = 3/4, tau_xy = 0.75 ... still fine; the linear first step + # loses positivity when eta_eff gammadot / G exceeds 1, which needs gammadot > 1 + # at any dt. Drive it with a faster wall: gammadot = 4 -> tau_xy = 2 eta_eff. + eta = shear_modulus = 1.0 + speed, height, width = 2.0, 1.0, 2.0 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-width / 2, -height / 2), + maxCoords=(width / 2, height / 2)) + 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.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1, integrator="bdf") + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.shear_modulus = shear_modulus + stokes.constitutive_model.Parameters.dt_elastic = 1.0 + stokes.add_dirichlet_bc((speed, 0.0), "Top") + stokes.add_dirichlet_bc((-speed, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-8 + stokes.solve(timestep=1.0, zero_init_guess=False) + health = stokes.constitutive_model.conformation_min_eigenvalue() + # gammadot = 4, eta_eff = 1/2: tau_xy = 2, eigenvalues 1 -+ 2 -> -1 and 3, everywhere + assert abs(health["min"] + 1.0) < 1.0e-6 + assert health["fraction_negative"] == 1.0 + assert health["where"] is not None + + +@pytest.mark.parametrize("transport", ["semi_lagrangian", "integration_point", "forward"]) +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 + # read from a projection converged to 1e-6, hence the tolerance. + assert abs(stokes.constitutive_model.max_elastic_timestep(0.3) - 0.3) < 1.0e-5 + assert abs(stokes.constitutive_model.max_elastic_timestep(1.0) - 1.0) < 1.0e-5 + # nothing stretches without an objective rate, so there is no cap + assert _one_shear_step(transport, dt=1.0).constitutive_model.max_elastic_timestep() == float("inf") + + +def test_the_store_smoothing_is_the_coefficient_times_the_local_cell_size_squared(): + stokes = _one_shear_step("integration_point", dt=1.0) + history = stokes.DFDt + assert history.store_smoothing == 0.0 + assert history._commit_projection.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 + 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]) + assert abs(a - 0.05 * h * h) < 1.0e-12 * max(1.0, h * h) + with pytest.raises(ValueError): + history.store_smoothing = -1.0 diff --git a/tests/test_1060_stress_store_smoothing.py b/tests/test_1060_stress_store_smoothing.py new file mode 100644 index 000000000..90b93db0c --- /dev/null +++ b/tests/test_1060_stress_store_smoothing.py @@ -0,0 +1,104 @@ +"""The store smoothing of the integration-point stress history (#737). + +Waters and King start-up below Courant one on a pure Maxwell element: the case +that rings. Minutes, so level 2. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_c] # 15 min: a slow characterisation with a hard baseline, reported not gated (as test_1064 is) + + + +def _cell_scale_content(history, mesh): + """RMS of the part of the carried point values no per-cell P1 function + can represent, relative to the field: the mode the store cycle grows.""" + values, pts = history.carried_tensors() + nq = int(history.psi_star[0].num_points_per_cell) + ncell = pts.shape[0] // nq + w = np.asarray(mesh.integration_rule.getData()[1]).reshape(-1); w = w / w.sum() + A = np.concatenate([np.ones((ncell, nq, 1)), pts.reshape(ncell, nq, mesh.cdim)], axis=2) + M = np.einsum("cqi,q,cqj->cij", A, w, A) + raw = values.reshape(values.shape[0], -1) + S = raw.reshape(ncell, nq, -1) + beta = np.linalg.solve(M, np.einsum("cqi,q,cqk->cik", A, w, S)) + fit = np.einsum("cqi,cik->cqk", A, beta).reshape(raw.shape) + rms = lambda a: float(np.sqrt((a ** 2).mean())) + 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", + 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 + of the carried stress at t 1 and at t_end, and u at t_end.""" + h, Lx, eta, lam, G = 1.0, 1.0, 1.0, 1.0, 1.0 + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(-Lx, -h), maxCoords=(Lx, h), + cellSize=h / res, qdegree=3, regular=True) + v = uw.discretisation.MeshVariable(f"U_wk{store_smoothing}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P_wk{store_smoothing}", mesh, 1, degree=1) + ns = uw.systems.NavierStokes(mesh, v, p, rho=1.0, order=1) + ns.stress_transport = transport + ns.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + ns.Unknowns, order=1, integrator="bdf") + ns.constitutive_model.Parameters.shear_viscosity_0 = eta + ns.constitutive_model.Parameters.shear_modulus = eta / lam + ns.constitutive_model.Parameters.dt_elastic = dt + 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": + ns.DFDt.store_smoothing = store_smoothing + elif transport == "forward": + 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": + carry = ns.DFDt.update_pre_solve + def carry_and_measure(*args, **kwargs): + out = carry(*args, **kwargs) + latest["content"] = _cell_scale_content(ns.DFDt, mesh) + return out + ns.DFDt.update_pre_solve = carry_and_measure + centre = np.array([[0.0, 0.0]]) + content = {} + u1 = None + for step in range(int(round(t_end / dt))): + ns.solve(timestep=dt, zero_init_guess=False) + t = (step + 1) * dt + if abs(t - 1.0) < dt / 2: + u1 = float(np.asarray(uw.function.evaluate(v.sym[0], centre)).reshape(-1)[0]) + content[1.0] = latest["content"] + content[t_end] = latest["content"] + u_end = float(np.asarray(uw.function.evaluate(v.sym[0], centre)).reshape(-1)[0]) + if return_kind: + return u1, content, u_end, type(ns.DFDt).__name__ + return u1, content, u_end + + +def test_the_store_smoothing_holds_the_cell_scale_mode_of_the_integration_point_history(): + """Waters and King, 1/16, dt 0.0125 (Courant 0.2), pure Maxwell. + + The mode grows from round-off, so its AMPLITUDE is the platform's (a + checkerboard seeded into the store is projected away within a few steps + and does not set it); its GROWTH RATE is the scheme's, about 2.2 per unit + time here: measured without smoothing, the cell-scale content of the + carried stress goes from 5.8e-6 at t 1 to 5.1e-5 at t 2 (a factor 8.8) + and rings by t 5. With c = 0.07 in the store it decays instead + (1.3e-6 at t 1, 2.3e-7 at t 2). The cost on the centre velocity at t 1 + (0.9617 plain) is two percent here and scales with h^2. + """ + 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" + 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 + assert smooth[2.0] < 1.0e-6, smooth # and at the round-off floor + assert abs(u_plain - 0.9617) < 0.002, u_plain + assert abs(u_smooth - 0.9429) < 0.002, u_smooth diff --git a/tests/test_1061_stress_forward_history.py b/tests/test_1061_stress_forward_history.py new file mode 100644 index 000000000..195cacc26 --- /dev/null +++ b/tests/test_1061_stress_forward_history.py @@ -0,0 +1,24 @@ +"""The forward semi-Lagrangian stress history on the case that rings. + +Waters and King start-up below Courant one on a pure Maxwell element, to t 6.5, +where the unsmoothed integration-point history has rung (0.626 at 1/16) and the +nodal one sits at 0.5171. The forward history transports its memory +consistently and has the same cell-scale mode: unsmoothed it diverges at t 2.4 +here. With its read-back smoothing at c = 0.023 (the least dose that holds the +integration-point store on a regular mesh; 0.07 is the recommended one) it runs +clean to t 8. Twenty minutes: level 3. +""" +import pytest + +from test_1060_stress_store_smoothing import waters_king_start_up + +pytestmark = [pytest.mark.level_3, pytest.mark.tier_c] # 20 min on Waters-King: reported not gated, as test_1064 is + + +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" + # 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 + assert abs(u65 - 0.5185) < 0.003 diff --git a/tests/test_1063_stress_history_restart.py b/tests/test_1063_stress_history_restart.py new file mode 100644 index 000000000..24db68a26 --- /dev/null +++ b/tests/test_1063_stress_history_restart.py @@ -0,0 +1,66 @@ +"""A stress history survives a snapshot and restore. + +Run the Maxwell shear box for six steps, snapshot the model, run six more and +record the stress; restore, run the same six again: the stress must be the +same to round-off. For every history that carries the stress. A history +whose bookkeeping (step history, initialisation flag, launch values) were not +captured would either re-initialise from zero or carry the wrong step. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_a] # solves: not level 1 + + +def _shear_box(transport): + uw.reset_default_model() + orchestration_model = uw.get_default_model() + eta = G = 1.0 + speed, height, width, dt = 0.5, 1.0, 2.0, 0.1 + mesh = uw.meshing.StructuredQuadBox(elementRes=(16, 8), minCoords=(-width / 2, -height / 2), + maxCoords=(width / 2, height / 2)) + v = uw.discretisation.MeshVariable(f"U_r_{transport}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P_r_{transport}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=1) + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.shear_modulus = G + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.add_dirichlet_bc((speed, 0.0), "Top") + stokes.add_dirichlet_bc((-speed, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-10 + return orchestration_model, stokes, dt + + +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"]) +def test_a_restored_history_continues_where_it_left_off(transport): + orchestration_model, stokes, dt = _shear_box(transport) + for _ in range(6): + stokes.solve(timestep=dt, zero_init_guess=False) + snap = orchestration_model.save_state() + for _ in range(6): + stokes.solve(timestep=dt, zero_init_guess=False) + straight = _stress_at_origin(stokes) + # A swarm history carries its stress on particles whose count and order need + # not be reproduced identically across a restore (population control refills + # cells), so the raw per-particle array is compared only for the mesh + # flavours; the carried stress read at a point is compared for all. + straight_field = None if transport == "lagrangian" else np.array(stokes.DFDt.psi_star[0].data) + orchestration_model.load_state(snap) + for _ in range(6): + stokes.solve(timestep=dt, zero_init_guess=False) + assert abs(_stress_at_origin(stokes) - straight) < 1.0e-12 + if straight_field is not None: + assert np.allclose(np.asarray(stokes.DFDt.psi_star[0].data), straight_field, rtol=0.0, atol=1.0e-12) + # and it is a real twelve-step state, not a re-initialised six-step one + assert abs(straight - (1.0 - 1.1 ** -12)) < 1.0e-6 diff --git a/tests/test_1101_advdiff_swarm_rotating_gaussian.py b/tests/test_1101_advdiff_swarm_rotating_gaussian.py new file mode 100644 index 000000000..d84172fae --- /dev/null +++ b/tests/test_1101_advdiff_swarm_rotating_gaussian.py @@ -0,0 +1,133 @@ +"""The swarm (Lagrangian) advection-diffusion solver against the rotating Gaussian. + +AdvDiffusionSwarm carries the scalar's history on a user swarm and reads the +mesh solution back onto the particles each step. On the diffusing rotating +Gaussian (exact at every time) on a disc it is as accurate as the +semi-Lagrangian and streamline-upwind solvers over a half revolution, and it +diffuses correctly: the trap of a particle scheme is to keep the particle's old, +sharp value and under-diffuse, which shows as a peak above the exact one. + +The disc is the accurate geometry (its rotation keeps every cell well filled), +but the solver must also hold a domain the flow crosses. On a square the flow +crosses all four walls and clamps out-flowing particles into a thin boundary +layer; a high-degree per-cell projection of that layer overshoots and diverges, +which is why the history proxy defaults to degree 1 (well conditioned on the +clamped layer). Both geometries are exercised here. + +The swarm is supplied by the caller and is NOT private to the solver; the solver +declares its history variable, so it is built before the swarm is populated, and +the caller advects the swarm each step. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +SIGMA = 0.12 +KAPPA = 0.01 +T_END = float(np.pi) # half revolution +EXACT_PEAK = SIGMA ** 2 / (SIGMA ** 2 + 2 * KAPPA * T_END) + + +def _disc(res=24): + return uw.meshing.Annulus(radiusInner=0.0, radiusOuter=1.2, cellSize=2.4 / res, qdegree=3) + + +def _run(scheme, particle_update="pic", step_averaging=1, res=24): + uw.reset_default_model() + mesh = _disc(res) + 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]]) + swarm = None + if scheme == "supg": + adv = uw.systems.AdvDiffusion(mesh, T, V, order=1) + elif scheme == "slcn": + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, order=1) + elif scheme == "swarm": + swarm = uw.swarm.Swarm(mesh) + adv = uw.systems.AdvDiffusionSwarm(mesh, T, V, swarm=swarm, order=1, + particle_update=particle_update, step_averaging=step_averaging) + swarm.populate(fill_param=4) + swarm.population_control = dict() + if scheme != "supg": + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = KAPPA + adv.add_dirichlet_bc(0.0, "Upper") + dt = 0.02 + nsteps = int(round(T_END / dt)); dt = T_END / nsteps + for _ in range(nsteps): + if swarm is not None: + swarm.advection(V, dt, order=2) + adv.solve(timestep=dt) + err = float(sol.error(sol.at(T_END), T, norm="integral")) + peak = float(np.asarray(T.data).max()) + return adv, err, peak + + +def test_the_swarm_solver_is_as_accurate_as_the_mesh_solvers_over_a_half_revolution(): + adv, err, peak = _run("swarm") + assert type(adv.DuDt).__name__ == "Lagrangian_Swarm" + assert adv.swarm is not None + # BASELINE: half revolution, kappa 0.01, dt 0.02, disc res 24, degree-1 + # history proxy (the default; see the ledger). SUPG and SLCN give 1.7e-2 and + # 2.5e-2 here; the swarm solver sits between them, held to a hard number. + assert abs(err - 0.0204) < 0.004, err + assert abs(peak - EXACT_PEAK) < 0.01, (peak, EXACT_PEAK) # it diffuses to the exact peak + + +def test_the_half_blend_under_diffuses_a_particle_scheme_trap(): + """step_averaging=2 keeps half the particle's old sharp value each step, so + the field does not diffuse to the exact peak: the trap the default avoids.""" + _adv, err, peak = _run("swarm", step_averaging=2) + assert peak > EXACT_PEAK + 0.1, peak + assert err > 0.2, err + + +def test_it_will_not_run_without_a_swarm(): + uw.reset_default_model() + mesh = _disc(8) + x, y = mesh.X + T = uw.discretisation.MeshVariable("Tn", mesh, 1, degree=2) + with pytest.raises(ValueError, match="needs a swarm"): + uw.systems.AdvDiffusionSwarm(mesh, T, sympy.Matrix([[-y, x]]), swarm=None) + + +def test_the_solver_holds_a_domain_the_flow_crosses_a_square_box(): + """The degree-1 history proxy holds a rotating SQUARE, where the flow crosses + all four walls and clamps out-flowing particles into a thin boundary layer. + A degree-2 proxy overshoots the projection of that layer and diverges by a + half turn; degree 1 stays bounded and diffuses to the exact peak. + """ + 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) + x, y = mesh.X + sol = uw.analytic.RotatingGaussian(mesh, sigma=SIGMA, centre_radius=0.5, omega=1.0, diffusivity=KAPPA) + T = uw.discretisation.MeshVariable("Tb", mesh, 1, degree=2) + T.array[:, 0, 0] = uw.function.evaluate(sol.at(0.0), T.coords).reshape(-1) + V = sympy.Matrix([[-y, x]]) + swarm = uw.swarm.Swarm(mesh) + adv = uw.systems.AdvDiffusionSwarm(mesh, T, V, swarm=swarm, order=1) # proxy_degree=1 default + swarm.populate(fill_param=4) + swarm.population_control = dict() + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = KAPPA + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + dt = 0.02 + nsteps = int(round(T_END / dt)); dt = T_END / nsteps + for _ in range(nsteps): + swarm.advection(V, dt, order=2) + adv.solve(timestep=dt) + peak = float(np.asarray(T.data).max()) + err = float(sol.error(sol.at(T_END), T, norm="integral")) + assert peak < 0.5, peak # bounded: degree 2 reaches ~13 here + assert abs(peak - EXACT_PEAK) < 0.02, (peak, EXACT_PEAK) # and diffuses to the exact peak + # BASELINE: 6.7e-2 (the box is harder than the disc; the point is it holds) + assert abs(err - 0.067) < 0.02, err