Stress histories in a units-aware model: reduce on the way in, read back in Pa (#788)

docs/developer/subsystems/stress-transport.md

@@ -48,6 +48,18 @@ 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 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
 
 The objective-rate source $L\sigma^* + \sigma^* L^T$ acts on the carried stress

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
 
@@ -2003,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):
@@ -3757,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):

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."""

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)