diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 39bf5692..25b7fbd7 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1172,6 +1172,7 @@ def __init__( smoothing=0.0, preserve_moments=False, with_forcing_history: bool = False, + monotone_mode: Optional[str] = None, ): super().__init__() @@ -1185,6 +1186,20 @@ def __init__( self.order = order self.preserve_moments = preserve_moments self.with_forcing_history = with_forcing_history + # Monotonicity limiter for the SL trace-back result. Bound + # the FE-interpolated upstream sample to the local data range + # of psi_star at each trace-back point. Cures the FE Lagrange + # overshoot pattern in cells with sharp gradients while + # preserving FE accuracy elsewhere. + # None → pure FE (legacy; can overshoot at non-nodal + # points in cells with sharp gradients) + # "clamp" → B.2: clip FE result to [nbr_min, nbr_max] of + # the k=dim+1 nearest psi_star DOFs + # "pick" → B.1: keep FE if in nbr bounds, else re-evaluate + # via RBF (Shepard) at out-of-bounds DOFs + # Settable after construction: + # ``adv_diff.DuDt.monotone_mode = "clamp"`` + self.monotone_mode = monotone_mode # Forcing-history storage. Allocated only if requested. Populated # each step via update_forcing_history(forcing_fn) — used by ETD-2 @@ -1636,6 +1651,7 @@ def update_pre_solve( verbose: Optional[bool] = False, dt_physical: Optional[float] = None, store_result: Optional[bool] = True, + monotone_mode: Optional[str] = "__instance__", ): """Sample upstream values along characteristics before solve. @@ -1651,10 +1667,19 @@ def update_pre_solve( existing psi_star levels upstream. Used by VE_Stokes where psi_star[0] already contains the projected actual stress from the previous solve. + monotone_mode : str or None, optional + Override the instance ``self.monotone_mode`` for this call. + Default ``"__instance__"`` (sentinel) means "use whatever is + on the instance". Pass ``None``, ``"clamp"``, or ``"pick"`` + to force a particular mode for one call. """ self._dt = dt + # Resolve monotone_mode: explicit kwarg overrides instance attr. + if monotone_mode == "__instance__": + monotone_mode = getattr(self, "monotone_mode", None) + if not self._history_initialised: self.initialise_history() @@ -1903,6 +1928,7 @@ def update_pre_solve( v_mid_result = uw.function.global_evaluate( self.V_fn, mid_pt_coords, + evalf=evalf, ) # CRITICAL: Preserve UnitAwareArray through slicing @@ -1967,9 +1993,16 @@ def update_pre_solve( # Evaluate psi_star at upstream coordinates # global_evaluate now returns dimensional results (gateway fix 2025-11-28) + # When evalf=True, route through RBF (Shepard, bounded by + # neighbour values) instead of FE shape functions. FE + # Lagrange P3 can overshoot at non-nodal upstream points + # in cells with sharp gradients — observed as the 'pepper' + # DOF scatter that ignites catastrophic ringing on free- + # surface convection at high Ra. value_at_end_points = uw.function.global_evaluate( expr_to_evaluate, end_pt_coords, + evalf=evalf, ) # CRITICAL FIX (2025-11-27): If psi_star has units, ensure the assigned @@ -1978,6 +2011,87 @@ def update_pre_solve( 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) + # Monotonicity limiter (B.1 / B.2). Bound the FE/RBF + # trace-back result to the local data range of psi_star + # near each upstream coord. + # monotone_mode = "clamp" → B.2: clip to [nbr_min, nbr_max] + # monotone_mode = "pick" → B.1: keep FE in bounds, else RBF + if monotone_mode in ("clamp", "pick"): + # Convert end_pt_coords to ND space matching the + # psi_star DOF coords so the kdtree query compares + # like-with-like. Same pattern as the + # psi_star_0_coords → psi_star_0_coords_nd conversion + # above. In unit-aware runs `end_pt_coords` is + # dimensional (e.g. metres) while + # `psi_star[i].coords_nd` is [0,1] non-dimensional, + # so without this conversion the kdtree picks wrong + # neighbours. + if hasattr(end_pt_coords, "magnitude"): + epc_raw = uw.non_dimensionalise(end_pt_coords) + if isinstance(epc_raw, UnitAwareArray): + epc_nd = np.array(epc_raw) + elif hasattr(epc_raw, "magnitude"): + epc_nd = epc_raw.magnitude + else: + epc_nd = np.asarray(epc_raw) + else: + epc_nd = np.asarray(end_pt_coords) + psi_coords_nd = np.asarray(self.psi_star[i].coords_nd) + if hasattr(psi_coords_nd, "magnitude"): + psi_coords_nd = np.asarray(psi_coords_nd.magnitude) + nnn = self.mesh.dim + 1 + kdt = uw.kdtree.KDTree( + np.ascontiguousarray(psi_coords_nd)) + _, idxs = kdt.query( + np.ascontiguousarray(epc_nd), k=nnn, + sqr_dists=False) + psi_data = np.asarray(self.psi_star[i].data) + # Flatten trailing dims to compute per-coord nbr stats, + # then restore the original shape afterwards. + psi_data_flat = psi_data.reshape(psi_data.shape[0], -1) + nbr_vals = psi_data_flat[idxs] + nbr_min = nbr_vals.min(axis=1) + nbr_max = nbr_vals.max(axis=1) + veep_np = np.asarray(value_at_end_points) + orig_shape = veep_np.shape + veep_flat = veep_np.reshape(nbr_min.shape) + if monotone_mode == "clamp": + veep_lim = np.clip(veep_flat, nbr_min, nbr_max) + else: + # B.1 "pick": re-evaluate via RBF only at the + # subset of upstream coords whose FE result is + # out-of-bounds. Keeps pick-mode cost dominated + # by the cheap FE pass when most DOFs are + # in-bounds; the RBF eval is the expensive part + # so we don't want it for points we're going to + # keep anyway. + out_of_bounds = ((veep_flat < nbr_min) + | (veep_flat > nbr_max)) + oob_mask = out_of_bounds.any(axis=tuple( + range(1, out_of_bounds.ndim))) + veep_lim = veep_flat.copy() + if oob_mask.any(): + oob_coords = end_pt_coords[oob_mask] + value_rbf_oob = uw.function.global_evaluate( + expr_to_evaluate, oob_coords, evalf=True) + vrbf_flat = np.asarray(value_rbf_oob).reshape( + (-1,) + veep_flat.shape[1:]) + # Within the oob subset, only overwrite + # entries that are individually out of bounds + # (multi-component case). + sub_oob = out_of_bounds[oob_mask] + veep_sub = veep_lim[oob_mask] + veep_sub = np.where( + sub_oob, vrbf_flat, veep_sub) + veep_lim[oob_mask] = veep_sub + value_at_end_points = veep_lim.reshape(orig_shape) + # Re-wrap units after numpy ops, if needed. + 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 # disable this for now - Compute moments before update