From 4462b7d045663d7e07ef96f8ebb25e7491bb43c8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 25 Sep 2026 18:10:36 -0700 Subject: [PATCH 1/2] Stress histories in a units-aware model: stores are non-dimensional work arrays, one dimensional read (#788) Every history store is a non-dimensional work array behind the units boundary and carries no units of its own. In a model with reference quantities the histories failed in four ways: the nodal commit and the SUPG carry copied through unit-aware .array views (np.copy raised); the particle flavours wrote the dimensional values evaluate returns straight into storage and ran away (1e148 after 20 steps); the integration-point trace cache keyed on an unhashable unit-carrying timestep; and nothing read back in pascals. - The timestep is reduced where every flavour receives it. - Anything written from evaluate is reduced on the way in. - Copies between stores go through .data with the component layout (EnhancedMeshVariable gains the layout delegate). - DFDt.carried(level) is the one way to read a history back: the constitutive model supplies the map from the stored value to the stress, the value times the model's unit of stress, so the read is in pascals when reference scales are set and unchanged when they are not. Test: the Maxwell shear box with reference quantities, all five stress_transport flavours, read through carried and held to the loading curve in Pa at the order-1 tolerance. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- docs/developer/subsystems/stress-transport.md | 13 +++ src/underworld3/constitutive_models.py | 22 +++++ .../discretisation/enhanced_variables.py | 14 +++ src/underworld3/systems/ddt.py | 94 ++++++++++++------- src/underworld3/systems/solvers.py | 6 +- tests/test_1064_stress_history_units.py | 90 ++++++++++++++++++ 6 files changed, 205 insertions(+), 34 deletions(-) create mode 100644 tests/test_1064_stress_history_units.py diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index 0459a49b..8dbd5f0d 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -48,6 +48,19 @@ 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. +### Units + +Every history store is a non-dimensional work array behind the units boundary, +whichever flavour holds it, and carries no units of its own. What enters is +reduced on the way in: the timestep where each flavour receives it, and anything +a flavour writes from `evaluate` (which returns dimensional values), such as the +flux at the particles or an inflow datum. Copies between stores go through +`.data`. The one way to read a history back is `DFDt.carried(level)`: the +constitutive model supplies the map from the stored value to the stress, so the +read is in pascals when reference scales are set. `DFDt.psi_star[level]` is the +raw store. The Maxwell shear box with reference quantities (`test_1064`) holds +all five flavours to the same loading curve in Pa. + ## The timestep is set by the wall strain rate, not the far-field Courant number The objective-rate source $L\sigma^* + \sigma^* L^T$ acts on the carried stress diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 0abf2239..75212587 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1626,6 +1626,17 @@ def _object_viewer(self): return +def _reference_stress(): + """The model's unit of stress as an expression (non-dimensional value one), + or ``None`` when the model has no reference scales.""" + orchestration_model = uw.get_default_model() + if not orchestration_model.has_units(): + return None + scale = orchestration_model.get_scale_for_dimensionality((1 * uw.units.Pa).dimensionality) + return expression(r"\sigma_{\mathrm{ref}}", uw.quantity(float(scale.magnitude), str(scale.units)), + "the model's unit of stress") + + class ViscoElasticPlasticFlowModel(ViscousFlowModel): r""" Viscoelastic-plastic flow constitutive model. @@ -2085,6 +2096,17 @@ def stress_history_ddt_kwargs(self): return {"with_forcing_history": True} return {} + def decode_history(self, stored): + r"""The stress a stored history value stands for, in the model's units. + + The history stores the stress non-dimensionally; this multiplies it by + the model's unit of stress (an expression whose non-dimensional value is + one), so a read through ``DFDt.carried`` comes back in pascals when + reference scales are set, and unchanged when they are not. + """ + unit = _reference_stress() + return stored if unit is None else stored * unit + # The following should have no setters @property def stress_star(self): diff --git a/src/underworld3/discretisation/enhanced_variables.py b/src/underworld3/discretisation/enhanced_variables.py index 9a1388fa..420809ee 100644 --- a/src/underworld3/discretisation/enhanced_variables.py +++ b/src/underworld3/discretisation/enhanced_variables.py @@ -325,6 +325,20 @@ def units(self): """Units for this variable.""" return self._base_var.units + def _data_layout(self, i, j=None): + """Column of a component in the flat ``.data`` storage. + + Parameters + ---------- + i, j : int + Component indices; ``j`` is omitted for a vector. + + Returns + ------- + int + """ + return self._base_var._data_layout(i, j) + @property def has_units(self) -> bool: """Check if this variable has units.""" diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 4c1aa4e7..66742672 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -987,7 +987,10 @@ def commit_flux_to_history(self, flux, verbose=False): if not hasattr(self, "_psi_star_projection_solver"): self._setup_projections() - transported = np.copy(self.psi_star[0].array[...]) + # The stores are non-dimensional work arrays: copies between them go + # through .data, never the unit-aware .array (#788). + level_0 = self.psi_star[0] + transported = np.array(level_0.data) if getattr(self, "_psi_star_use_multicomponent", False): # The snapshot machinery has frozen the projection's input, so this @@ -996,18 +999,18 @@ def commit_flux_to_history(self, flux, verbose=False): 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 + values = np.asarray(self._psi_star_flat_var.data[:, k]) + level_0.data[:, level_0._data_layout(i, j)] = values if i != j: - self.psi_star[0].array[:, j, i] = values + level_0.data[:, level_0._data_layout(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.psi_star[level].data[...] = ( + transported if level == 1 else np.asarray(self.psi_star[level - 1].data)) self._history_committed = True @@ -1084,6 +1087,30 @@ def inflow_value(self, value): #: Whether this flavour compiles :attr:`inflow_value` into its transport. applies_inflow_value = False + def carried(self, level: int = 0): + """The carried quantity of history ``level``, as an expression in its + own units. + + The history's stores are non-dimensional work arrays; this is the one + way to read them back. The owner of the quantity (a solver's + constitutive model) supplies the map from what is stored to what it + stands for (``_record_unmap``): the stored value times the quantity's + unit, or the decoding of a stored encoding. Without an owner the + stored value is returned as it is. + + Parameters + ---------- + level : int, default 0 + 0 is the newest level. + + Returns + ------- + sympy.Matrix + """ + stored = self.psi_star[level].sym + unmap = self.__dict__.get("_record_unmap") + return stored if unmap is None else unmap(stored) + def _nondim_timestep(self, dt): r"""Reduce ``dt`` to a plain non-dimensional model-time value. @@ -1399,7 +1426,7 @@ def update_pre_solve( verbose: Optional[bool] = False, ): """Pre-solve update hook. Auto-initialises history on first call.""" - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if not self._history_initialised: self.initialise_history() @@ -1417,7 +1444,7 @@ def update_post_solve( verbose: Optional[bool] = False, ): r"""Shift history chain after solve: :math:`\psi^{*n} \leftarrow \psi^{*(n-1)}`.""" - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if verbose: print(f"Updating history for ψ = {self.psi_fn}", flush=True) @@ -1804,7 +1831,7 @@ def update_pre_solve( grid-based advection correction so that bdf() approximates the material derivative Dφ/Dt rather than ∂φ/∂t. """ - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if not self._history_initialised: self.initialise_history() @@ -1859,7 +1886,7 @@ def update_post_solve( level it should hold. Invisible at order 1, a 150-fold error at order 2 on the analytic Maxwell shear box. """ - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if verbose and uw.mpi.rank == 0: print(f"Update {self.psi_fn}", flush=True) @@ -2270,18 +2297,19 @@ def _transport_history(self, dt, verbose=False): 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[...] + flat.data[:, k] = np.asarray(history.data[:, history._data_layout(i, j)]) + previous.data[...] = np.asarray(flat.data) self._transport_solver.solve(verbose=verbose) for k, (i, j) in enumerate(indices): - values = flat.array[:, 0, k] - history.array[:, i, j] = values + values = np.asarray(flat.data[:, k]) + history.data[:, history._data_layout(i, j)] = values if i != j: - history.array[:, j, i] = values + history.data[:, history._data_layout(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.""" + dt = self._nondim_timestep(dt) 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: @@ -3395,7 +3423,7 @@ def update_post_solve( dt_physical: Optional[float] = None, ): """Post-solve: record timestep and increment solve counter.""" - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) # Record timestep history for variable-dt BDF for i in range(self.order - 1, 0, -1): @@ -3699,7 +3727,7 @@ def update_pre_solve( to force a particular mode for one call. """ - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) # Resolve monotone_mode: explicit kwarg overrides instance attr. if monotone_mode == "__instance__": @@ -4137,9 +4165,9 @@ def initialise_history(self): ij = psi_star_0._data_layout(i, j) if ij in updated: continue - updated[ij] = np.asarray( + updated[ij] = np.asarray(_to_nondim_ndarray( uw.function.evaluate(self.psi_fn[i, j], coords) - ).reshape(-1) + )).reshape(-1) for ij, vals in updated.items(): psi_star_0.data[:, ij] = vals @@ -4180,7 +4208,7 @@ def update_pre_solve( ``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 + self._dt = dt = self._nondim_timestep(dt) if not self._history_initialised: self.initialise_history() @@ -4263,7 +4291,7 @@ def update_post_solve( **_ignored, ): """Shift history chain and advect swarm after solve.""" - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) # Record timestep history for variable-dt BDF for i in range(self.order - 1, 0, -1): @@ -4297,9 +4325,9 @@ def update_post_solve( ij = psi_star_0._data_layout(i, j) if ij in updated: continue - updated[ij] = np.asarray( + updated[ij] = np.asarray(_to_nondim_ndarray( uw.function.evaluate(self.psi_fn[i, j], coords, evalf=evalf) - ).reshape(-1) + )).reshape(-1) for ij, vals in updated.items(): psi_star_0.data[:, ij] = vals @@ -4599,7 +4627,7 @@ def update_pre_solve( arguments (the nodal manager's ``store_result``, ``dt_physical``, ``monotone_mode``) are accepted and ignored so a solver written for the nodal history can drive this one.""" - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if not self._history_initialised: self.initialise_history() @@ -4630,9 +4658,9 @@ def _proxy_values_at_particles(self, slot, coords, evalf): for i in range(slot.shape[0]): for j in range(slot.shape[1]): ij = slot._data_layout(i, j) - out[:, ij] = np.asarray( + out[:, ij] = np.asarray(_to_nondim_ndarray( uw.function.evaluate(mv.sym[i, j], coords, evalf=evalf) - ).reshape(-1) + )).reshape(-1) return out def update_post_solve( @@ -4653,7 +4681,7 @@ def update_post_solve( new (audit SWARM-06); a shift before the evaluation would hand the stress expression the wrong levels. """ - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) phi = 1 / self.step_averaging psi_star_0 = self.psi_star[0] @@ -4668,9 +4696,9 @@ def update_post_solve( ij = psi_star_0._data_layout(i, j) if ij in updated: continue # symmetric storage: one evaluation per slot - updated[ij] = np.asarray( + updated[ij] = np.asarray(_to_nondim_ndarray( uw.function.evaluate(self.psi_fn[i, j], coords, evalf=evalf) - ).reshape(-1) + )).reshape(-1) # Record timestep history for variable-dt BDF for i in range(self.order - 1, 0, -1): @@ -5330,7 +5358,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, 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 + self._dt = dt = self._nondim_timestep(dt) if not self._history_initialised: self.initialise_history() _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) @@ -5351,7 +5379,7 @@ 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 = dt = self._nondim_timestep(dt) for i in range(self.order - 1, 0, -1): self._dt_history[i] = self._dt_history[i - 1] self._dt_history[0] = dt @@ -5736,7 +5764,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, store_result=True, ** :meth:`commit_flux_to_history` (the viscoelastic case); otherwise the tracked field is read at the launch points first. """ - self._dt = dt + self._dt = dt = self._nondim_timestep(dt) if self._geometry_stamp() != self._launch_geometry: raise NotImplementedError( "ForwardSemiLagrangian: the launch set, cell measures and boundary faces were " @@ -5765,7 +5793,7 @@ 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 = dt = self._nondim_timestep(dt) self._dt_history[0] = dt if self._n_solves_completed < self.order: self._n_solves_completed += 1 diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index fd4b4ba0..3859cb5f 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -88,7 +88,11 @@ def expression(*args, **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.""" + whole flux. Transposed to the history's row layout. The model also hands + the history its map back to a dimensional stress (``DFDt.carried``).""" + DDt = constitutive_model.Unknowns.DFDt + if DDt is not None: + DDt._record_unmap = getattr(constitutive_model, "decode_history", None) if hasattr(constitutive_model, "history_flux"): return constitutive_model.history_flux.T return constitutive_model.flux.T diff --git a/tests/test_1064_stress_history_units.py b/tests/test_1064_stress_history_units.py new file mode 100644 index 00000000..a5b9eb33 --- /dev/null +++ b/tests/test_1064_stress_history_units.py @@ -0,0 +1,90 @@ +"""Every stress history carries a viscoelastic stress through a units-aware model. + +The Maxwell shear box of test_1059, with reference quantities set: the +velocity is given in km/Myr, the viscosity in Pa s, the modulus in Pa, the +timestep in Myr. The history stores are non-dimensional work arrays behind +the units boundary; what enters them must be reduced on the way in, and +what a user reads back through ``.array`` or ``evaluate`` must come out in +Pa. The analytic loading curve +:math:`\\sigma_{xy} = \\eta\\dot{\\gamma}\\,(1 - e^{-t\\mu/\\eta})` +is the baseline, in Pa (#788). +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +VISCOSITY_PA_S = 1.0e21 +RELAXATION_MYR = 1.0 +SPEED_KM_MYR, HEIGHT_KM, WIDTH_KM = 0.5, 1.0, 2.0 +MYR_S = 3.15576e13 + +TRANSPORTS = ["semi_lagrangian", "integration_point", "forward", "lagrangian", "eulerian"] + + +def _units_model(): + uw.reset_default_model() + orchestration_model = uw.get_default_model() + orchestration_model.set_reference_quantities( + length=uw.quantity(1.0, "km"), + viscosity=uw.quantity(VISCOSITY_PA_S, "Pa*s"), + time=uw.quantity(1.0, "Myr"), + ) + return orchestration_model + + +def _maxwell_shear_with_units(transport, steps=20, dt_myr=0.1): + _units_model() + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-WIDTH_KM / 2, -HEIGHT_KM / 2), + maxCoords=(WIDTH_KM / 2, HEIGHT_KM / 2)) + v = uw.discretisation.MeshVariable(f"U_{transport}", mesh, mesh.dim, degree=2, units="km/Myr") + p = uw.discretisation.MeshVariable(f"P_{transport}", mesh, 1, degree=1, units="Pa") + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.stress_transport = transport + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1) + parameters = stokes.constitutive_model.Parameters + parameters.shear_viscosity_0 = uw.quantity(VISCOSITY_PA_S, "Pa*s") + parameters.shear_modulus = uw.quantity(VISCOSITY_PA_S / (RELAXATION_MYR * MYR_S), "Pa") + dt = uw.quantity(dt_myr, "Myr") + parameters.dt_elastic = dt + stokes.add_dirichlet_bc((uw.quantity(SPEED_KM_MYR, "km/Myr"), 0.0), "Top") + stokes.add_dirichlet_bc((uw.quantity(-SPEED_KM_MYR, "km/Myr"), 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) + return stokes, dt + + +def _exact_pa(t_myr): + shear_rate_per_s = 2.0 * SPEED_KM_MYR / HEIGHT_KM / MYR_S + return VISCOSITY_PA_S * shear_rate_per_s * (1.0 - np.exp(-t_myr / RELAXATION_MYR)) + + +def _carried_pa(stokes): + """The history's first level at the origin, read through ``carried``, in Pa.""" + carried = uw.function.evaluate(stokes.DFDt.carried(0)[0, 1], np.array([[0.0, 0.0]])) + quantity = uw.units.Quantity(float(np.asarray(carried).reshape(-1)[0]), carried.units) + return float(quantity.to("Pa").magnitude) + + +@pytest.mark.parametrize("transport", TRANSPORTS) +def test_stress_history_loads_in_pascals(transport): + """The carried stress, read back through ``DFDt.carried``, is the Maxwell + curve in Pa: the same 2% order-1 tolerance as the non-dimensional box + (test_1059), on 2.74e7 Pa at t = 2 Myr. The store itself is a + non-dimensional work array (sigma/G of order one here).""" + steps, dt_myr = 20, 0.1 + stokes, dt = _maxwell_shear_with_units(transport, steps, dt_myr) + exact_pa = _exact_pa(steps * dt_myr) + assert abs(_carried_pa(stokes) - exact_pa) / exact_pa < 0.02, transport + assert np.abs(np.asarray(stokes.DFDt.psi_star[0].data)).max() < 10.0 From fd01e6167cca1de1fe82adf6c58b9cd0dac990ae Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 25 Sep 2026 21:12:35 -0700 Subject: [PATCH 2/2] Stress histories: the stores carry their quantity's units; drop the separate dimensional read (#788) Adversarial review of the previous commit found the separate read (carried, decode_history, a reference-stress expression, a private map handed from the solver) reimplements what MeshVariable units already do, and found defects the test could not see. Now: - A solver builds its stress history with the stress's units (the flux it hands over at construction is a zero placeholder and cannot say so). .data stays non-dimensional; .array and evaluate read in Pa. One helper, _history_units, replaces the derivation pasted into four flavours, and Eulerian/EulerianSUPG stores take units too. - One rule for writing an evaluation into a store (_write_evaluated): a result that carries units is reduced, one that does not is already non-dimensional; the re-wraps that assumed an unlabelled result was dimensional are gone, and so are the units attached before reducing in the inflow, integration-point and forward writes. The Eulerian and Lagrangian_Swarm initialisations were writing dimensional values unreduced. - The BDF-2 step-ratio guard compared a dimensional dt_elastic with the non-dimensional history and silently fell back to BDF-1 whenever the time unit differed from the reference; both are compared non-dimensionally. - The mirror write in the component copies (wrong for a full tensor, a no-op for a symmetric one) is removed; one timestep reducer, built on _as_float. - Swarm advection wrote the dimensional velocity global_evaluate returns into its model-unit arithmetic, so in a units model the particles did not move (every swarm-carried quantity, not only stress); the reads are reduced. Test: every flavour, in a units model with the timestep in kyr (not the reference unit) and as the same problem in plain numbers, from a moving start, orders 1 and 2: the stores agree to 1e-6 and the stress reads back in Pa as the non-dimensional value times the stress scale. The previous test was blind to the timestep, the order-2 guard, the initial writes and the swarm motion. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_017kSkAq7oJ5J3XisuvLBovo --- docs/developer/subsystems/stress-transport.md | 19 +- src/underworld3/constitutive_models.py | 32 +--- src/underworld3/swarm.py | 14 +- src/underworld3/systems/ddt.py | 169 ++++++------------ src/underworld3/systems/solvers.py | 10 +- tests/test_1064_stress_history_units.py | 121 ++++++------- 6 files changed, 146 insertions(+), 219 deletions(-) diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index 8dbd5f0d..cd197add 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/docs/developer/subsystems/stress-transport.md @@ -50,16 +50,15 @@ is wrong: sub-cell layers and no-slip walls. ### Units -Every history store is a non-dimensional work array behind the units boundary, -whichever flavour holds it, and carries no units of its own. What enters is -reduced on the way in: the timestep where each flavour receives it, and anything -a flavour writes from `evaluate` (which returns dimensional values), such as the -flux at the particles or an inflow datum. Copies between stores go through -`.data`. The one way to read a history back is `DFDt.carried(level)`: the -constitutive model supplies the map from the stored value to the stress, so the -read is in pascals when reference scales are set. `DFDt.psi_star[level]` is the -raw store. The Maxwell shear box with reference quantities (`test_1064`) holds -all five flavours to the same loading curve in Pa. +Every history store holds non-dimensional values in `.data`, whichever flavour +holds it; that is what the solver reads and what copies between stores go +through. A solver builds its stress history with the stress's units, so the +store's `.array`, and `evaluate` of its symbol, read back in pascals when +reference scales are set. What enters is reduced on the way in: the timestep +where each flavour receives it, and anything a flavour writes from `evaluate` +(which returns dimensional values). `test_1064` runs every flavour in a units +model with the timestep in kyr and as the same problem in plain numbers, from a +moving start, at orders 1 and 2: the stores agree to solver precision. ## The timestep is set by the wall strain rate, not the far-field Courant number diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 75212587..0aa15984 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -48,7 +48,7 @@ from underworld3.swarm import IndexSwarmVariable from underworld3.discretisation import MeshVariable from underworld3.systems.ddt import SemiLagrangian as SemiLagrangian_DDt -from underworld3.systems.ddt import _bdf_coefficients +from underworld3.systems.ddt import _bdf_coefficients, _as_float from underworld3.function.quantities import UWQuantity from underworld3.systems.ddt import Lagrangian as Lagrangian_DDt @@ -1626,17 +1626,6 @@ def _object_viewer(self): return -def _reference_stress(): - """The model's unit of stress as an expression (non-dimensional value one), - or ``None`` when the model has no reference scales.""" - orchestration_model = uw.get_default_model() - if not orchestration_model.has_units(): - return None - scale = orchestration_model.get_scale_for_dimensionality((1 * uw.units.Pa).dimensionality) - return expression(r"\sigma_{\mathrm{ref}}", uw.quantity(float(scale.magnitude), str(scale.units)), - "the model's unit of stress") - - class ViscoElasticPlasticFlowModel(ViscousFlowModel): r""" Viscoelastic-plastic flow constitutive model. @@ -2014,7 +2003,9 @@ def _update_bdf_coefficients(self): dt_history = self.Unknowns.DFDt._dt_history if order >= 2 and len(dt_history) > 0 and dt_history[0] is not None: try: - ratio = float(dt_current) / float(dt_history[0]) + # both as non-dimensional model time: the history keeps its + # steps reduced, while dt_elastic may be a quantity + ratio = _as_float(dt_current) / _as_float(dt_history[0]) if ratio > self._max_dt_ratio_for_higher_order: order = 1 except (TypeError, ZeroDivisionError): @@ -2096,17 +2087,6 @@ def stress_history_ddt_kwargs(self): return {"with_forcing_history": True} return {} - def decode_history(self, stored): - r"""The stress a stored history value stands for, in the model's units. - - The history stores the stress non-dimensionally; this multiplies it by - the model's unit of stress (an expression whose non-dimensional value is - one), so a read through ``DFDt.carried`` comes back in pascals when - reference scales are set, and unchanged when they are not. - """ - unit = _reference_stress() - return stored if unit is None else stored * unit - # The following should have no setters @property def stress_star(self): @@ -3779,7 +3759,9 @@ def _update_bdf_coefficients(self): dt_history = self.Unknowns.DFDt._dt_history if order >= 2 and len(dt_history) > 0 and dt_history[0] is not None: try: - ratio = float(dt_current) / float(dt_history[0]) + # both as non-dimensional model time: the history keeps its + # steps reduced, while dt_elastic may be a quantity + ratio = _as_float(dt_current) / _as_float(dt_history[0]) if ratio > self._max_dt_ratio_for_higher_order: order = 1 except (TypeError, ZeroDivisionError): diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index b24792ce..48c2d0e2 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5717,6 +5717,12 @@ def advection( import underworld3 as uw delta_t_model = uw.scaling.non_dimensionalise(delta_t) + # The particle arithmetic is in model units: a velocity read by + # global_evaluate comes back dimensional in a units model and is reduced. + from underworld3.systems.ddt import _to_nondim_ndarray + + def _nondim_velocity(value): + return np.asarray(_to_nondim_ndarray(value))[:, 0, :] dt_limit = self.estimate_dt(V_fn) @@ -5821,9 +5827,9 @@ def advection( # rank-local evaluation silently extrapolates wrong values # for it (SWARM-16 / BF-16). - v_at_Vpts[...] = uw.function.global_evaluate( + v_at_Vpts[...] = _nondim_velocity(uw.function.global_evaluate( V_fn_matrix, self._particle_coordinates.data - )[:, 0, :] + )) mid_pt_coords = ( self._particle_coordinates.data[...] @@ -5838,7 +5844,7 @@ def advection( # (since the mid-points might have moved off-proc) # - v_at_Vpts[...] = uw.function.global_evaluate(v_mid_matrix, mid_pt_coords)[:, 0, :] + v_at_Vpts[...] = _nondim_velocity(uw.function.global_evaluate(v_mid_matrix, mid_pt_coords)) new_coords = X0.array[:, 0, :] + delta_t_model * v_at_Vpts / substeps @@ -5858,7 +5864,7 @@ def advection( print(f"1. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", flush=True) v_at_Vpts = np.zeros_like(coords) - v_at_Vpts[...] = uw.function.global_evaluate(V_fn_matrix, coords[...])[:, 0, :] + v_at_Vpts[...] = _nondim_velocity(uw.function.global_evaluate(V_fn_matrix, coords[...])) if self.verbose: print(f"2. Advection (1st): {coords.shape} v {self.local_size} - swarm point shape", flush=True) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 66742672..7739c992 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -211,6 +211,31 @@ def _as_float(value): return None +def _history_units(psi_fn, units=None): + """The units a history's stores are built with: ``units`` when given (a + solver knows what its history carries), else those of ``psi_fn``; ``None`` + outside a model with reference quantities. The stores hold non-dimensional + values in ``.data`` whatever the units; the units only say what ``.array`` + and ``evaluate`` read back as.""" + if not uw.get_default_model().has_units(): + return None + return units if units is not None else uw.get_units(psi_fn) + + +def _write_evaluated(var, values): + """Write an evaluation of a store's quantity into the store. + + ``values`` is what ``evaluate`` or ``global_evaluate`` returned: dimensional + when it carries units (it is reduced), non-dimensional when it does not. + It is written component by component into ``.data``, the store's + non-dimensional storage, whatever units the store is built with. + """ + shape = tuple(var.sym.shape) + values = np.asarray(_to_nondim_ndarray(values)).reshape(-1, *shape) + for (i, j) in _storage_components(var.vtype, shape): + var.data[:, var._data_layout(i, j)] = values[:, i, j] + + def _to_nondim_ndarray(value, units=None): """Reduce a possibly unit-carrying array to a plain non-dimensional ndarray. @@ -1001,8 +1026,6 @@ def commit_flux_to_history(self, flux, verbose=False): for k, (i, j) in enumerate(self._psi_star_indep_indices): values = np.asarray(self._psi_star_flat_var.data[:, k]) level_0.data[:, level_0._data_layout(i, j)] = values - if i != j: - level_0.data[:, level_0._data_layout(j, i)] = values else: self._psi_star_projection_solver.uw_function = flux self._psi_star_projection_solver.smoothing = 0.0 @@ -1087,69 +1110,24 @@ def inflow_value(self, value): #: Whether this flavour compiles :attr:`inflow_value` into its transport. applies_inflow_value = False - def carried(self, level: int = 0): - """The carried quantity of history ``level``, as an expression in its - own units. - - The history's stores are non-dimensional work arrays; this is the one - way to read them back. The owner of the quantity (a solver's - constitutive model) supplies the map from what is stored to what it - stands for (``_record_unmap``): the stored value times the quantity's - unit, or the decoding of a stored encoding. Without an owner the - stored value is returned as it is. - - Parameters - ---------- - level : int, default 0 - 0 is the newest level. - - Returns - ------- - sympy.Matrix - """ - stored = self.psi_star[level].sym - unmap = self.__dict__.get("_record_unmap") - return stored if unmap is None else unmap(stored) - def _nondim_timestep(self, dt): - r"""Reduce ``dt`` to a plain non-dimensional model-time value. - - The semi-Lagrangian trace-back is performed ENTIRELY in the mesh's - NON-DIMENSIONAL (DM) coordinate space: evaluate()/global_evaluate - treat plain arrays as DM coords and the DM point-location uses DM - values (0..L_model, NOT dimensional metres). So coords, velocity - AND dt are all reduced to non-dimensional values, whether or not - the model carries units. (Previously the has_units branch kept - dimensional coords/velocity and left dt unitless -> a 'meter' vs - 'meter/second' subtraction crash and mislocation against the ND - DM; UW3 issue #267.) - """ - if hasattr(dt, "magnitude") or hasattr(dt, "value"): - # dt carries units -> non-dimensionalise it - dt_nondim = uw.non_dimensionalise(dt, uw.get_default_model()) - if hasattr(dt_nondim, "magnitude"): - return float(dt_nondim.magnitude) - elif hasattr(dt_nondim, "value"): - return float(dt_nondim.value) - else: - return float(dt_nondim) - else: - # already non-dimensional model-time - return dt - + """The timestep as a non-dimensional model time (:func:`_as_float`); a + symbolic timestep passes through unchanged.""" + reduced = _as_float(dt) + return dt if reduced is None else reduced def _write_inflow(self, var, coords, rows): """Overwrite ``rows`` of ``var`` with :attr:`inflow_value` evaluated at ``coords`` (the positions of ALL the points, so that the read, which is collective when the value holds a field, is made on every rank; - only ``rows`` are written). Storage is non-dimensional, so the value - is reduced through the history's units. A flavour that calls this - sets ``_components`` (its stored columns) and ``_psi_units``.""" + only ``rows`` are written). Storage is non-dimensional, so a value + that evaluates with units is reduced. A flavour that calls this sets + ``_components`` (its stored columns).""" expr = 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] + _to_nondim_ndarray(vals)).reshape(-1)[rows] def _unknown_shape(self): """Shape of the unknown as a matrix (``Symbolic`` stores ``_shape`` as data).""" @@ -1552,6 +1530,7 @@ def __init__( order=1, smoothing=0.0, num_components=None, + units=None, ): super().__init__() @@ -1605,6 +1584,7 @@ def __init__( degree=degree, continuous=continuous, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", + units=_history_units(self._psi_fn, units), ) ) @@ -1733,11 +1713,11 @@ def update_history_fn(self): pass try: - self.psi_star[0].data[...] = uw.function.evaluate( + self.psi_star[0].data[...] = np.asarray(_to_nondim_ndarray(uw.function.evaluate( self.psi_fn, self.psi_star[0].coords, evalf=self.evalf, - ).reshape(-1, max(self.psi_fn.shape)) + ))).reshape(-1, max(self.psi_fn.shape)) except Exception: # Sanctioned fallback: evaluate() cannot interpolate # expressions containing derivatives (e.g. flux terms) — @@ -1992,6 +1972,7 @@ def __init__( peclet_weight: float = 4.0, num_components=None, transport_on_update: bool = False, + units=None, ): order = int(order) if order not in (1, 2, 3): @@ -2012,7 +1993,7 @@ def __init__( super().__init__( mesh, psi_fn, vtype, degree, continuous, V_fn=None, theta=theta, varsymbol=varsymbol, verbose=verbose, bcs=[] if bcs is None else bcs, - order=order, smoothing=smoothing, num_components=num_components, + order=order, smoothing=smoothing, num_components=num_components, units=units, ) self._advection_mode = "assembled" self._integrator = "am" if order == 1 else "bdf" @@ -2303,8 +2284,6 @@ def _transport_history(self, dt, verbose=False): for k, (i, j) in enumerate(indices): values = np.asarray(flat.data[:, k]) history.data[:, history._data_layout(i, j)] = values - if i != j: - history.data[:, history._data_layout(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 @@ -2730,6 +2709,7 @@ def __init__( theta: float = 0.5, old_frame_traceback: bool = False, midtime_velocity: bool = True, + units=None, ): super().__init__() @@ -2845,16 +2825,7 @@ def __init__( psi_star = [] self.psi_star = psi_star - # Propagate units from psi_fn to psi_star if the model supports units. - # Internal psi_star variables should match the user's variable units when possible, - # but if no reference quantities are set, use unitless variables to avoid strict mode errors. - psi_units = uw.get_units(psi_fn) - - # Check if the model can handle units (has reference quantities set) - model = uw.get_default_model() - if psi_units is not None and not model.has_units(): - # Model doesn't have reference quantities - don't propagate units to internal vars - psi_units = None + psi_units = _history_units(psi_fn, units) for i in range(order): self.psi_star.append( @@ -2865,7 +2836,7 @@ def __init__( degree=self.degree, continuous=self.continuous, varsymbol=rf"{{ {varsymbol}^{{ {'*'*(i+1)} }} }}", - units=psi_units, # Inherit units from psi_fn (or None if model has no units) + units=psi_units, ) ) @@ -3165,11 +3136,7 @@ def initialise_history(self): coords_nd = _to_nondim_ndarray(self.psi_star[0].coords) try: - eval_result = uw.function.evaluate(self.psi_fn, coords_nd) - psi_units = self.psi_star[0].units - if psi_units is not None and not isinstance(eval_result, UnitAwareArray): - eval_result = UnitAwareArray(eval_result, units=psi_units) - self.psi_star[0].array[...] = eval_result + _write_evaluated(self.psi_star[0], uw.function.evaluate(self.psi_fn, coords_nd)) except Exception: # Fallback: project psi_fn onto psi_star[0] via the SNES projector. # Route through the shared builder so snapshot substitution @@ -3554,12 +3521,7 @@ def _record_current_field_into_history( node_coords_nd, evalf=evalf, ) - # Wrap result with units if psi_star has units but eval didn't return UnitAwareArray - psi_star_units = self.psi_star[0].units - if psi_star_units is not None and not isinstance(eval_result, UnitAwareArray): - eval_result = UnitAwareArray(eval_result, units=psi_star_units) - - self.psi_star[0].array[...] = eval_result + _write_evaluated(self.psi_star[0], eval_result) except Exception: # Fallback to projection solver for expressions that can't be directly evaluated @@ -3677,13 +3639,7 @@ def _sample_history_at_departure( monotone=monotone_mode, ) - # CRITICAL FIX (2025-11-27): If psi_star has units, ensure the assigned - # value also has units. global_evaluate may return plain arrays. - psi_star_units = self.psi_star[i].units - if psi_star_units is not None and not isinstance(value_at_end_points, UnitAwareArray): - value_at_end_points = UnitAwareArray(value_at_end_points, units=psi_star_units) - - self.psi_star[i].array[...] = value_at_end_points + _write_evaluated(self.psi_star[i], value_at_end_points) # TODO(DESIGN): a moment-preserving correction (restore mean and L2 # moment of psi_star after the semi-Lagrangian update) was removed @@ -3791,7 +3747,7 @@ def update_pre_solve( # 3. Trace the characteristics back and sample each history slot # at its departure points. Work from the oldest slot backwards # so we don't overwrite history terms we still need to sample. - dt_for_calc = self._nondim_timestep(dt) + dt_for_calc = dt # Phase-2 ALE: if an adapt stashed Δx, build v_mesh = Δx / dt as # a per-DDt MeshVariable now so the trace-back below can use @@ -4027,6 +3983,7 @@ def __init__( fill_param=3, proxy_location="cells", proxy_sampling="reconstruct", + units=None, ): super().__init__() @@ -4039,12 +3996,7 @@ def __init__( self.V_fn = V_fn self.verbose = verbose self.order = order - # Particle storage is non-dimensional; an inflow datum with units is - # reduced through these before it is written (#783). - psi_units = uw.get_units(psi_fn) - if psi_units is not None and not uw.get_default_model().has_units(): - psi_units = None - self._psi_units = psi_units + psi_units = _history_units(psi_fn, units) self._components = _storage_components(vtype, tuple(sympy.Matrix(psi_fn).shape)) self._init_history_tracking(order) @@ -4063,6 +4015,7 @@ def __init__( proxy_location=proxy_location, proxy_sampling=proxy_sampling, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", + units=psi_units, ) ) @@ -4316,10 +4269,6 @@ def update_post_solve( psi_star_0 = self.psi_star[0] coords = np.asarray(self.swarm._particle_coordinates.data) updated = {} - # TODO(BUG): the evaluated psi_fn is written into non-dimensional - # storage without reduction through _psi_units (see _write_inflow); a - # psi_fn carrying units lands at its physical magnitude. Same in - # initialise_history and in Lagrangian_Swarm. for i in range(psi_star_0.shape[0]): for j in range(psi_star_0.shape[1]): ij = psi_star_0._data_layout(i, j) @@ -4589,7 +4538,7 @@ def initialise_history(self): coords, ) psi_star_0.data[:, psi_star_0._data_layout(i, j)] = np.asarray( - updated_psi + _to_nondim_ndarray(updated_psi) ).reshape(-1) # Copy to all other history slots @@ -4886,6 +4835,7 @@ def __init__( monotone_mode: Optional[str] = None, with_forcing_history: bool = False, store_smoothing: float = 0.0, + units=None, **_unsupported, # TODO(BUG): swallowed without a stated failure mode (Charter 5) ): super().__init__() @@ -4921,10 +4871,7 @@ def __init__( 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 + psi_units = _history_units(self._psi_fn, units) # A vector or tensor history is one dof per INDEPENDENT component per # point. The trace-back and the weighted sums are shape-agnostic, so @@ -5247,7 +5194,7 @@ def _write_components(self, var, expr, coords, evaluate=None, **kwargs): for column, (i, j) in enumerate(self._components): vals = evaluate(expr[i, j], coords, **kwargs) var.data[:, column] = np.asarray( - _to_nondim_ndarray(vals, units=self._psi_units) + _to_nondim_ndarray(vals) ).reshape(-1) def _segment_dt(self, j, dt): @@ -5442,6 +5389,7 @@ def __init__( varsymbol: Optional[str] = None, order: int = 1, theta: float = 0.5, + units=None, **_unsupported, ): super().__init__() @@ -5468,10 +5416,7 @@ def __init__( 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 + psi_units = _history_units(self._psi_fn, units) self.psi_star = [ uw.discretisation.MeshVariable( f"psi_star_fwd_{inst}", mesh, vtype=vtype, degree=1, continuous=False, @@ -5723,8 +5668,8 @@ def _fit_arrivals(self, X, values, cell=None, inflow=None): 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) + dofs[bcells].reshape(-1, mesh.cdim)) + ).reshape(-1) for (i, j) in self._components]).reshape(bcells.size, ndof, self.num_components) short = fed[bcells] if short.any(): diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 3859cb5f..8021763d 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -88,11 +88,7 @@ def expression(*args, **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. The model also hands - the history its map back to a dimensional stress (``DFDt.carried``).""" - DDt = constitutive_model.Unknowns.DFDt - if DDt is not None: - DDt._record_unmap = getattr(constitutive_model, "decode_history", None) + 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 @@ -1784,6 +1780,9 @@ def _create_stress_history_ddt(self, order=2): bcs=None, order=order, smoothing=0.0001, + # the history carries a stress; the flux handed over at construction + # is a zero placeholder, so it cannot say so itself + units=uw.units.Pa, ) if self.stress_transport == "integration_point": unsupported = set(ddt_kwargs) - {"with_forcing_history"} @@ -1810,6 +1809,7 @@ def _create_stress_history_ddt(self, order=2): sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), self.u.sym, vtype=common["vtype"], varsymbol=common["varsymbol"], order=order, + units=common["units"], ) elif self.stress_transport == "lagrangian": if ddt_kwargs: diff --git a/tests/test_1064_stress_history_units.py b/tests/test_1064_stress_history_units.py index a5b9eb33..3a8f2db1 100644 --- a/tests/test_1064_stress_history_units.py +++ b/tests/test_1064_stress_history_units.py @@ -1,13 +1,15 @@ -"""Every stress history carries a viscoelastic stress through a units-aware model. - -The Maxwell shear box of test_1059, with reference quantities set: the -velocity is given in km/Myr, the viscosity in Pa s, the modulus in Pa, the -timestep in Myr. The history stores are non-dimensional work arrays behind -the units boundary; what enters them must be reduced on the way in, and -what a user reads back through ``.array`` or ``evaluate`` must come out in -Pa. The analytic loading curve -:math:`\\sigma_{xy} = \\eta\\dot{\\gamma}\\,(1 - e^{-t\\mu/\\eta})` -is the baseline, in Pa (#788). +"""A stress history in a units-aware model is the non-dimensional history, read in Pa. + +The Maxwell shear box twice: once with reference quantities (1 km, 1e21 Pa s, +1 Myr) and every input given with units -- the timestep in kyr, deliberately not +the reference time unit -- and once as the same problem in plain numbers +(eta = G = 1, dt = 0.1, speed 0.5). Both start from the steady shear profile +already in place, so each history's first level is written from the flux of a +moving flow. The history's stores are non-dimensional work arrays: they must +agree between the two runs to solver precision, at order 1 and at order 2 (where +the order-2 weights depend on the ratio of the current step to the previous +one). Read back through the variable, the stress is in Pa: its non-dimensional +value times the stress scale, 1e21 Pa s / 1 Myr (#788). """ import numpy as np @@ -18,73 +20,66 @@ pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] -VISCOSITY_PA_S = 1.0e21 -RELAXATION_MYR = 1.0 -SPEED_KM_MYR, HEIGHT_KM, WIDTH_KM = 0.5, 1.0, 2.0 MYR_S = 3.15576e13 - -TRANSPORTS = ["semi_lagrangian", "integration_point", "forward", "lagrangian", "eulerian"] +STRESS_SCALE_PA = 1.0e21 / MYR_S +STEPS = 12 +ORIGIN = np.array([[0.0, 0.0]]) +CASES = [(t, 1) for t in ("semi_lagrangian", "integration_point", "forward", "lagrangian", "eulerian")] \ + + [(t, 2) for t in ("semi_lagrangian", "integration_point", "eulerian")] -def _units_model(): +def _shear_box(transport, order, with_units): uw.reset_default_model() - orchestration_model = uw.get_default_model() - orchestration_model.set_reference_quantities( - length=uw.quantity(1.0, "km"), - viscosity=uw.quantity(VISCOSITY_PA_S, "Pa*s"), - time=uw.quantity(1.0, "Myr"), - ) - return orchestration_model - - -def _maxwell_shear_with_units(transport, steps=20, dt_myr=0.1): - _units_model() - mesh = uw.meshing.StructuredQuadBox( - elementRes=(16, 8), minCoords=(-WIDTH_KM / 2, -HEIGHT_KM / 2), - maxCoords=(WIDTH_KM / 2, HEIGHT_KM / 2)) - v = uw.discretisation.MeshVariable(f"U_{transport}", mesh, mesh.dim, degree=2, units="km/Myr") - p = uw.discretisation.MeshVariable(f"P_{transport}", mesh, 1, degree=1, units="Pa") + if with_units: + uw.get_default_model().set_reference_quantities( + length=uw.quantity(1.0, "km"), viscosity=uw.quantity(1.0e21, "Pa*s"), + time=uw.quantity(1.0, "Myr")) + mesh = uw.meshing.StructuredQuadBox(elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5)) + tag = f"{transport[:3]}{order}{'u' if with_units else 'n'}" + v = uw.discretisation.MeshVariable(f"U_{tag}", mesh, 2, degree=2, units="km/Myr" if with_units else None) + p = uw.discretisation.MeshVariable(f"P_{tag}", mesh, 1, degree=1, units="Pa" if with_units else None) stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) stokes.stress_transport = transport - stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( - stokes.Unknowns, order=1) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=order) parameters = stokes.constitutive_model.Parameters - parameters.shear_viscosity_0 = uw.quantity(VISCOSITY_PA_S, "Pa*s") - parameters.shear_modulus = uw.quantity(VISCOSITY_PA_S / (RELAXATION_MYR * MYR_S), "Pa") - dt = uw.quantity(dt_myr, "Myr") + if with_units: + parameters.shear_viscosity_0 = uw.quantity(1.0e21, "Pa*s") + parameters.shear_modulus = uw.quantity(STRESS_SCALE_PA, "Pa") + dt = uw.quantity(100.0, "kyr") + speed = uw.quantity(0.5, "km/Myr") + else: + parameters.shear_viscosity_0 = 1.0 + parameters.shear_modulus = 1.0 + dt, speed = 0.1, 0.5 parameters.dt_elastic = dt - stokes.add_dirichlet_bc((uw.quantity(SPEED_KM_MYR, "km/Myr"), 0.0), "Top") - stokes.add_dirichlet_bc((uw.quantity(-SPEED_KM_MYR, "km/Myr"), 0.0), "Bottom") + 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.tolerance = 1.0e-8 stokes.petsc_options["snes_type"] = "newtonls" stokes.petsc_options["ksp_type"] = "fgmres" - for _ in range(steps): + # The steady shear profile already in place, written non-dimensionally (the + # coordinates come back in metres in a units model; 1 km/Myr is the velocity + # scale, so the non-dimensional values are the same in both runs). + X = np.asarray(uw.non_dimensionalise(v.coords) if with_units else v.coords) + v.data[:, 0], v.data[:, 1] = X[:, 1], 0.0 + for _ in range(STEPS): stokes.solve(timestep=dt, zero_init_guess=False) - return stokes, dt - - -def _exact_pa(t_myr): - shear_rate_per_s = 2.0 * SPEED_KM_MYR / HEIGHT_KM / MYR_S - return VISCOSITY_PA_S * shear_rate_per_s * (1.0 - np.exp(-t_myr / RELAXATION_MYR)) + return stokes -def _carried_pa(stokes): - """The history's first level at the origin, read through ``carried``, in Pa.""" - carried = uw.function.evaluate(stokes.DFDt.carried(0)[0, 1], np.array([[0.0, 0.0]])) - quantity = uw.units.Quantity(float(np.asarray(carried).reshape(-1)[0]), carried.units) - return float(quantity.to("Pa").magnitude) +@pytest.mark.parametrize("transport, order", CASES) +def test_units_model_history_is_the_nondimensional_history(transport, order): + plain = _shear_box(transport, order, with_units=False) + plain_store = np.array(plain.DFDt.psi_star[0].data) + plain_xy = float(np.asarray(uw.function.evaluate(plain.DFDt.psi_star[0].sym[0, 1], ORIGIN)).reshape(-1)[0]) + stokes = _shear_box(transport, order, with_units=True) + store = np.asarray(stokes.DFDt.psi_star[0].data) + assert np.abs(store - plain_store).max() < 1.0e-6 * np.abs(plain_store).max(), (transport, order) -@pytest.mark.parametrize("transport", TRANSPORTS) -def test_stress_history_loads_in_pascals(transport): - """The carried stress, read back through ``DFDt.carried``, is the Maxwell - curve in Pa: the same 2% order-1 tolerance as the non-dimensional box - (test_1059), on 2.74e7 Pa at t = 2 Myr. The store itself is a - non-dimensional work array (sigma/G of order one here).""" - steps, dt_myr = 20, 0.1 - stokes, dt = _maxwell_shear_with_units(transport, steps, dt_myr) - exact_pa = _exact_pa(steps * dt_myr) - assert abs(_carried_pa(stokes) - exact_pa) / exact_pa < 0.02, transport - assert np.abs(np.asarray(stokes.DFDt.psi_star[0].data)).max() < 10.0 + carried = uw.function.evaluate(stokes.DFDt.psi_star[0].sym[0, 1], ORIGIN) + value_pa = float(uw.units.Quantity(float(np.asarray(carried).reshape(-1)[0]), carried.units).to("Pa").magnitude) + expected_pa = plain_xy * STRESS_SCALE_PA + assert abs(value_pa - expected_pa) < 1.0e-6 * abs(expected_pa), (transport, order, value_pa, expected_pa)