diff --git a/docs/advanced/troubleshooting.md b/docs/advanced/troubleshooting.md index ddd7eaf96..588cbd99e 100644 --- a/docs/advanced/troubleshooting.md +++ b/docs/advanced/troubleshooting.md @@ -1,13 +1,52 @@ --- -title: "trouuleshooting" +title: "Troubleshooting" --- -# trouuleshooting +# Troubleshooting -```{note} Documentation In Progress -Content being developed for advanced users. -``` +Notes on behaviour changes and common pitfalls that show up as "my old script +runs but gives different answers / different performance". -## Overview +## Stokes penalty is now viscosity-scaled (June 2026) -Advanced topic documentation coming soon. +The Stokes augmented-Lagrangian (grad-div) incompressibility penalty changed +from a bare constant to a **viscosity-scaled** term. The operator changed +from + +$$ +\sigma + \lambda \, (\nabla \cdot \mathbf{u}) \, I +$$ + +to + +$$ +\sigma + \lambda \, \mu \, (\nabla \cdot \mathbf{u}) \, I +$$ + +where $\mu$ is the local viscosity (`constitutive_model.K`) and $\lambda$ is +the value you set with `stokes.penalty`. + +Why: with spatially variable viscosity, a constant penalty is enormous +relative to the stress in low-viscosity regions (over-stiffening them into +velocity locking) and negligible in high-viscosity regions. Scaling by the +local viscosity keeps the penalty proportional to the local stress scale +everywhere. + +**What to change in older scripts:** + +1. `stokes.penalty` is now a *dimensionless* number of order 1. Scripts that + tuned a large constant against the viscosity magnitude (say `penalty=1e6` + for a model with $\mu \sim 1$, or a value chosen to sit above the largest + viscosity in a contrast model) should drop that tuning — use `1.0` (or a + modest multiple) instead. Keeping the old large value multiplies it by + the viscosity again and can lock the velocity field or destroy the + conditioning of the velocity block. +2. A penalty is still **off by default** (`penalty = 0`) and usually + unnecessary: the pressure Schur complement is already preconditioned with + the local $1/\mu$ scaling. +3. The penalty term is part of the operator, so converged pressures include + the grad-div contribution; diagnostics that compared pressure against a + constant-penalty run will differ at the level of the divergence residual. + +See the `Stokes.penalty` property docstring for the full description, and the +`CONSTRAINED_FREESLIP_MULTIPLIER` design note for the derivation. diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index e39058ce2..86d694ef2 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -199,6 +199,11 @@ in `Stokes_Constrained` (#224), then made parallel-correct. - `selfp` Schur preconditioner default, viscosity-scaled penalty, and nullspace re-setup fix (#229); over-conservative serial guard removed (#240); gauge, convergence, knockout, and rotation-gauge fixes (#265). +- The main `Stokes.penalty` (augmented-Lagrangian grad-div) is likewise + viscosity-scaled since June 2026: the parameter is now a dimensionless + O(1) number, not a large constant tuned against the viscosity magnitude. + Migration note for older scripts in `docs/advanced/troubleshooting.md` + (#292). ### Boundary Conditions: Local-h Nitsche and Boundary-Slip Surfaces (June 2026) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 70592a4a0..b5149a0b4 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2406,9 +2406,11 @@ class SNES_Scalar(SolverBaseClass): ## Todo: some validity checking on the size / type of u_Field supplied if u_Field is None: - # TODO(BUG): num_components=mesh.dim for a SCALAR unknown looks wrong - # (a scalar has one component); kept as-is pending review (READ-22). - self.Unknowns.u = uw.discretisation.MeshVariable( mesh=mesh, num_components=mesh.dim, + # A scalar unknown has one component. (MeshVariable already forced + # this: vtype=SCALAR overrides num_components, so the old + # num_components=mesh.dim here was ignored, never over-allocated — + # issue #367.) + self.Unknowns.u = uw.discretisation.MeshVariable( mesh=mesh, num_components=1, varname="Us{}".format(SNES_Scalar._obj_count), vtype=uw.VarType.SCALAR, degree=degree, ) else: diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 8a2515246..32840c9e9 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1896,11 +1896,21 @@ def set_pressure_gauge(self, boundary, reference=0.0): boundary : str Mesh boundary label over which the mean pressure is fixed. reference : float, default 0.0 - Target mean pressure on ``boundary``. + Target mean pressure on ``boundary``, in the model (non-dimensional) + frame. Returns ------- the registered callback (so it can be identified/removed later). + + Notes + ----- + For **pressure-dependent plasticity** (e.g. Drucker-Prager, + :math:`\sigma_y = C + p\sin\varphi`) a pressure gauge is + counter-productive: shifting the pressure moves the yield surface at + every iteration and the nonlinear solve oscillates instead of + converging. The pressure level there is physical, not a free gauge — + set the level through the boundary conditions instead. """ p = self.Unknowns.p area = uw.maths.BdIntegral(self.mesh, 1.0, boundary).evaluate() @@ -1908,6 +1918,13 @@ def set_pressure_gauge(self, boundary, reference=0.0): def _pressure_gauge(solver, iteration): mean = p_surface_integral.evaluate() / area + # This callback fires inside the SNES non-dimensionalisation + # cordon, where p.data holds non-dimensional values. Under an + # active units model the boundary integral is a dimensional + # quantity, so take its non-dimensional value before applying + # the shift (issue #271). + if isinstance(mean, uw.function.quantities.UWQuantity): + mean = mean.data p.data[...] -= (mean - reference) return self.add_update_callback(_pressure_gauge) diff --git a/src/underworld3/units.py b/src/underworld3/units.py index 3714442d0..089fab46b 100644 --- a/src/underworld3/units.py +++ b/src/underworld3/units.py @@ -876,15 +876,16 @@ def non_dimensionalise(expression, model=None) -> Any: try: scale = model.get_scale_for_dimensionality(dimensionality) - result_qty = expression / scale - # TODO(BUG): UWQuantity.__init__ has no `dimensionality` - # kwarg, so this raises TypeError whenever a raw - # pint.Quantity is non-dimensionalised with active scaling - # and a resolvable scale. Found while verifying BF-12 - # (2026-07 audit). The uw.quantity (UWQuantity) path works - # and is the workaround. - # Return UWQuantity to preserve dimensionality - return UWQuantity(float(result_qty.magnitude), units="dimensionless", dimensionality=dimensionality) + # Convert to SI base units BEFORE dividing, matching the + # UnitAwareArray path above: the scale is computed in SI, so + # a non-SI input (km, cm/year, ...) must cancel through SI + # or the magnitude keeps the unit prefactor. + result_qty = expression.to_base_units() / scale + # A non-dimensionalised scalar is a plain float, matching the + # .magnitude returns above. (A dimensionless UWQuantity cannot + # carry the source dimensionality; use dimensionalise() with + # an explicit target_dimensionality to round-trip.) + return float(result_qty.magnitude) except ValueError as e: raise ValueError(f"Cannot non-dimensionalise Pint quantity: {e}") except ImportError: diff --git a/tests/parallel/ptest_004_nodal_swarm_advection.py b/tests/parallel/ptest_004_nodal_swarm_advection.py deleted file mode 100644 index 60a2c54db..000000000 --- a/tests/parallel/ptest_004_nodal_swarm_advection.py +++ /dev/null @@ -1,97 +0,0 @@ -import underworld3 as uw -from mpi4py import MPI - -# A simple parallel test of the SemiLagrangian's nodal swarm advection. -# This determines if nodal swarm particles are lost during -# large (going to different processors) advection steps. -# NOTE: Make sure to run this with 4 processors to ensure -# that the test scenario (particles going beyond neighboring processors) happens. - -dt = 1 -maxsteps = 1 -res = 16 -velocity = 1.25 # seems to take a long time if a higher velocity is set - -xmin, xmax = 0.0, 4.0 -ymin, ymax = 0.0, 1.0 - -mesh = uw.meshing.UnstructuredSimplexBox( - minCoords=(xmin, ymin), - maxCoords=(xmax, ymax), - cellSize=1 / res, - regular=False, - qdegree=3, - refinement=0, -) - -v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) - -# vector being advected -vec_tst = uw.discretisation.MeshVariable("Vn", mesh, mesh.dim, degree=2) - -DuDt = uw.systems.ddt.SemiLagrangian( - mesh, - vec_tst.sym, - v.sym, - vtype=uw.VarType.VECTOR, - degree=2, - continuous=vec_tst.continuous, - varsymbol=vec_tst.symbol, - verbose=False, - bcs=None, - order=1, -) - -# initialize variables -# TODO: DELETE remove swarm.access / data, replace with direct array assignment -# with mesh.access(v): -# v.data[:, 0] = velocity - -v.array[:, 0, 0] = velocity - -# we are only interested in monitoring the number of nodal swarm particles -# before and after advection -# TODO: DELETE remove swarm.access / data, replace with direct array assignment -# with mesh.access(vec_tst): -# vec_tst.data[:, 0] = 1. - -vec_tst.array[:, 0, 0] = 1.0 - -# get the number of nodal swarm particles BEFORE advection -# TODO: DELETE remove swarm.access / data, replace with direct array assignment -# with DuDt._nswarm_psi.access(DuDt._nswarm_psi): -# before_adv_swarm_num = len(DuDt._nswarm_psi.data) - -before_adv_swarm_num = len(DuDt._nswarm_psi.data) - -comm = uw.mpi.comm - -before_total = comm.allreduce(before_adv_swarm_num, op=MPI.SUM) - -# print(f"Before advection; rank {uw.mpi.rank} particles: {before_adv_swarm_num}", flush = True) - -# do one huge timestep -DuDt.update_pre_solve(dt, verbose=False, evalf=False) -# TODO: DELETE remove swarm.access / data, replace with direct array assignment -# with mesh.access(vec_tst): -# vec_tst.data[...] = DuDt.psi_star[0].data[...] - -vec_tst.array[...] = DuDt.psi_star[0].array[...] - -# get the number of nodal swarm particles AFTER advection -# TODO: DELETE remove swarm.access / data, replace with direct array assignment -# with DuDt._nswarm_psi.access(DuDt._nswarm_psi): -# after_adv_swarm_num = len(DuDt._nswarm_psi.data) - -after_adv_swarm_num = len(DuDt._nswarm_psi.data) - -after_total = comm.allreduce(after_adv_swarm_num, op=MPI.SUM) - -# print(f"After advection; rank {uw.mpi.rank} particles: {after_adv_swarm_num}", flush = True) - -if uw.mpi.rank == 0: - print(f"Before advection; Total particles in all ranks: {before_total}", flush=True) - print(f"After advection; Total particles in all ranks: {after_total}", flush=True) - -if uw.mpi.rank == 0: - assert after_total == before_total, "Error: Nodal swarm particles lost during advection." diff --git a/tests/parallel/test_0700_basic_parallel_operations.py b/tests/parallel/test_0700_basic_parallel_operations.py index 9236037ee..5bae24889 100644 --- a/tests/parallel/test_0700_basic_parallel_operations.py +++ b/tests/parallel/test_0700_basic_parallel_operations.py @@ -242,8 +242,10 @@ def test_nodal_swarm_advection_basic(): """ Test basic nodal swarm advection in parallel. - This is a simplified version of ptest_004 that verifies the advection - mechanism works without running a full simulation. + This carries the coverage of the retired ptest_004 script (issue #340): + it verifies the advection mechanism works without running a full + simulation. (The particle-count assertion ptest_004 made is obsolete — + the SL trace-back no longer allocates a nodal swarm.) Requires 4 ranks to test cross-processor particle movement. """ @@ -261,12 +263,14 @@ def test_nodal_swarm_advection_basic(): # Vector being advected vec_tst = uw.discretisation.MeshVariable("Vn", mesh, mesh.dim, degree=2) - # Create semi-Lagrangian advection + # Create semi-Lagrangian advection (degree/continuous match vec_tst) DuDt = uw.systems.ddt.SemiLagrangian( mesh, vec_tst.sym, v.sym, vtype=uw.VarType.VECTOR, + degree=2, + continuous=vec_tst.continuous, order=2, smoothing=0.0, ) diff --git a/tests/parallel/test_0765_internal_boundary_integral_mpi.py b/tests/parallel/test_0765_internal_boundary_integral_mpi.py index 7849fe30c..907159b8a 100644 --- a/tests/parallel/test_0765_internal_boundary_integral_mpi.py +++ b/tests/parallel/test_0765_internal_boundary_integral_mpi.py @@ -87,11 +87,10 @@ def test_deformed_spherical_shell_boundary_area_parallel(): a = math.log(2.0) mapped = (np.exp(a * t) - 1.0) / (math.exp(a) - 1.0) new_radii = 0.5 + thickness * mapped - # TODO(BUG): fails since the _deform_mesh live-mesh guard landed (PR #326): - # a direct _deform_mesh call on a mesh that already carries a variable is - # rejected. Needs migration to mesh.deform() or a sanctioned mutation - # scope. Pre-existing on development; observed while gating the #324 fix. - mesh._deform_mesh(coords * (new_radii / radii)[:, None]) + # Public deform() (pure geometric move, no dt): the raw _deform_mesh + # primitive is rejected by the live-mesh coordinate-mutation guard when + # a variable exists (issue #331). + mesh.deform(coords * (new_radii / radii)[:, None]) lower = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Lower").evaluate()) upper = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Upper").evaluate()) diff --git a/tests/test_0819_units_nondim_regressions.py b/tests/test_0819_units_nondim_regressions.py new file mode 100644 index 000000000..0103f8595 --- /dev/null +++ b/tests/test_0819_units_nondim_regressions.py @@ -0,0 +1,79 @@ +"""Regressions at the units <-> non-dimensional boundary. + +Covers two crashes found with an active scaling model: + +- issue #328: ``uw.non_dimensionalise(pint.Quantity)`` raised ``TypeError`` + (the protocol-5 branch constructed UWQuantity with an invalid + ``dimensionality=`` keyword whenever a scale was resolvable); +- issue #271: ``set_pressure_gauge`` computed the surface-mean pressure as a + dimensional quantity but applied it to the non-dimensional ``p.data`` + inside the SNES update callback, crashing the solve + (``_UFuncOutputCastingError``). +""" + +import numpy as np +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +@pytest.fixture() +def scaling_model(): + """Active reference scaling covering length, time and mass.""" + uw.reset_default_model() + orchestration_model = uw.get_default_model() + orchestration_model.set_reference_quantities( + length=uw.quantity(500, "m"), + velocity=uw.quantity(1, "cm/year"), + viscosity=uw.quantity(1e21, "Pa*s"), + ) + yield orchestration_model + uw.reset_default_model() + + +def test_non_dimensionalise_plain_pint_quantity(scaling_model): + # issue #328: a raw pint quantity with a resolvable scale must + # non-dimensionalise to a plain float, not raise TypeError. + result = uw.non_dimensionalise(250 * uw.units("m")) + assert isinstance(result, float) + assert result == pytest.approx(0.5) # length scale is 500 m + + +def test_non_dimensionalise_pint_quantity_non_si_units(scaling_model): + # The scale is computed in SI, so a non-SI input must reduce through SI: + # 0.25 km == 250 m == 0.5 length-scales. + assert uw.non_dimensionalise(0.25 * uw.units("km")) == pytest.approx(0.5) + + +def test_pressure_gauge_solves_under_units_model(scaling_model): + # issue #271: the gauge callback fires inside the SNES non-dimensional + # cordon; the dimensional surface-mean must be non-dimensionalised + # before it is subtracted from p.data. Before the fix this solve + # crashed with an object-dtype ufunc casting error. + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + units="metre", qdegree=3) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2, units="cm/year") + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, units="Pa") + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = uw.quantity(1e21, "Pa*s") + stokes.add_dirichlet_bc((uw.quantity(1.0, "cm/year"), 0.0), "Top") + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((0.0, 0.0), "Left") + stokes.add_dirichlet_bc((0.0, 0.0), "Right") + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1.0e-8 + + stokes.set_pressure_gauge("Top") + stokes.solve() + assert stokes.snes.getConvergedReason() > 0 + + # The gauge held: the non-dimensional mean pressure on Top is zero. + area = uw.maths.BdIntegral(mesh, 1.0, "Top").evaluate() + mean_top = uw.maths.BdIntegral(mesh, p.sym[0, 0], "Top").evaluate() / area + if isinstance(mean_top, uw.function.quantities.UWQuantity): + mean_top = mean_top.data + assert np.isclose(float(mean_top), 0.0, atol=1.0e-6), f"top-mean pressure = {mean_top}"