diff --git a/docs/developer/subsystems/stress-transport.md b/docs/developer/subsystems/stress-transport.md index 0459a49b..cd197add 100644 --- a/docs/developer/subsystems/stress-transport.md +++ b/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 diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 7da160f3..53334563 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 @@ -2086,7 +2086,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): @@ -3755,7 +3757,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/discretisation/enhanced_variables.py b/src/underworld3/discretisation/enhanced_variables.py index 19f58d22..311ef6e6 100644 --- a/src/underworld3/discretisation/enhanced_variables.py +++ b/src/underworld3/discretisation/enhanced_variables.py @@ -344,6 +344,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/swarm.py b/src/underworld3/swarm.py index 2e9cfdc6..4fc831e2 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5736,6 +5736,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) @@ -5845,9 +5851,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[...] @@ -5862,7 +5868,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 @@ -5882,7 +5888,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 375569bf..cd395b66 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. @@ -1140,7 +1165,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 @@ -1149,18 +1177,16 @@ 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 - if i != j: - self.psi_star[0].array[:, j, i] = values + values = np.asarray(self._psi_star_flat_var.data[:, k]) + level_0.data[:, level_0._data_layout(i, j)] = values else: self._psi_star_projection_solver.uw_function = flux self._psi_star_projection_solver.smoothing = 0.0 self._psi_star_projection_solver.solve(verbose=verbose) 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 @@ -1238,44 +1264,23 @@ def inflow_value(self, value): applies_inflow_value = False 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).""" @@ -1545,7 +1550,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() @@ -1563,7 +1568,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) @@ -1673,6 +1678,7 @@ def __init__( order=1, smoothing=0.0, num_components=None, + units=None, ): super().__init__() @@ -1726,6 +1732,7 @@ def __init__( degree=degree, continuous=continuous, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", + units=_history_units(self._psi_fn, units), ) ) @@ -1846,11 +1853,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) — @@ -1944,7 +1951,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() @@ -1999,7 +2006,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) @@ -2106,6 +2113,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): @@ -2126,7 +2134,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" @@ -2415,18 +2423,17 @@ 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 - if i != j: - history.array[:, j, i] = values + values = np.asarray(flat.data[:, k]) + history.data[:, history._data_layout(i, j)] = 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: @@ -2848,6 +2855,7 @@ def __init__( theta: float = 0.5, old_frame_traceback: bool = False, midtime_velocity: bool = True, + units=None, ): super().__init__() @@ -2963,16 +2971,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( @@ -2983,7 +2982,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, ) ) @@ -3276,11 +3275,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 @@ -3534,7 +3529,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): @@ -3676,12 +3671,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 @@ -3799,13 +3789,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 @@ -3849,7 +3833,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__": @@ -3913,7 +3897,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 @@ -4149,6 +4133,7 @@ def __init__( fill_param=3, proxy_location="cells", proxy_sampling="reconstruct", + units=None, ): super().__init__() @@ -4161,12 +4146,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) @@ -4185,6 +4165,7 @@ def __init__( proxy_location=proxy_location, proxy_sampling=proxy_sampling, varsymbol=rf"{varsymbol}^{{ {'*'*(i+1)} }}", + units=psi_units, ) ) @@ -4276,9 +4257,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 @@ -4319,7 +4300,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() @@ -4402,7 +4383,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): @@ -4428,18 +4409,14 @@ 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) 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 @@ -4690,7 +4667,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 @@ -4728,7 +4705,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() @@ -4759,9 +4736,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( @@ -4782,7 +4759,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] @@ -4797,9 +4774,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): @@ -4988,6 +4965,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__() @@ -5023,10 +5001,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 @@ -5343,7 +5318,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): @@ -5454,7 +5429,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) @@ -5475,7 +5450,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 @@ -5539,6 +5514,7 @@ def __init__( varsymbol: Optional[str] = None, order: int = 1, theta: float = 0.5, + units=None, **_unsupported, ): super().__init__() @@ -5565,10 +5541,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, @@ -5820,8 +5793,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(): @@ -5861,7 +5834,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 " @@ -5890,7 +5863,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 97711589..731ac32d 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1804,6 +1804,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"} @@ -1830,6 +1833,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 new file mode 100644 index 00000000..3a8f2db1 --- /dev/null +++ b/tests/test_1064_stress_history_units.py @@ -0,0 +1,85 @@ +"""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 +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +MYR_S = 3.15576e13 +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 _shear_box(transport, order, with_units): + uw.reset_default_model() + 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=order) + parameters = stokes.constitutive_model.Parameters + 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((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.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + # 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 + + +@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) + + 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)