From ba4f349df0462e411db2b9462abe36f82b2bc36f Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Thu, 26 Mar 2026 13:33:36 +1100 Subject: [PATCH 1/2] Add solver barrier for VEP and fix is_viscoplastic comparison - Add requires_stress_history property to Constitutive_Model base (False) and override in ViscoElasticPlasticFlowModel (True). - Stokes solver setter raises TypeError when assigned a constitutive model that requires stress history but DFDt is not available. Prevents silent failure where VEP on plain Stokes drops all history terms. - Fix is_viscoplastic: was comparing UWexpression == sympy.oo (always False due to type mismatch); now uses .sym is sympy.oo. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 17 ++++++++++++++++- .../cython/petsc_generic_snes_solvers.pyx | 10 ++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 47293a680..76ea0d7be 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -579,6 +579,16 @@ def _reset(self): return + @property + def requires_stress_history(self): + """Whether this model needs DFDt stress history tracking. + + Models that return True require a solver with stress history + management (e.g. VE_Stokes). Assigning such a model to a plain + Stokes solver will raise an error. + """ + return False + def _build_c_tensor(self): """Return the identity tensor of appropriate rank (e.g. for projections)""" @@ -1588,6 +1598,11 @@ def _object_viewer(self): ## Todo: add all the other properties in here ) + @property + def requires_stress_history(self): + """VEP models always require stress history tracking.""" + return True + @property def is_elastic(self): """True if elastic behavior is active (finite dt_elastic and shear_modulus).""" @@ -1604,7 +1619,7 @@ def is_elastic(self): @property def is_viscoplastic(self): """True if plastic yielding is active (finite yield_stress).""" - if self.Parameters.yield_stress == sympy.oo: + if self.Parameters.yield_stress.sym is sympy.oo: return False return True diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c4b3377e0..da048e46f 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -950,6 +950,16 @@ class SolverBaseClass(uw_object): "constitutive_model must be a valid class or instance of a valid class" ) + # Check that the solver can support this constitutive model's requirements. + # Models with stress history (VEP) need a solver that manages DFDt — e.g. VE_Stokes. + # Using them on a plain Stokes solver silently drops the history terms. + if self._constitutive_model.requires_stress_history and self.Unknowns.DFDt is None: + raise TypeError( + f"{type(self._constitutive_model).__name__} requires stress history tracking " + f"(DFDt). Use uw.systems.VE_Stokes instead of uw.systems.Stokes, or provide " + f"a DFDt object when constructing the solver." + ) + # May not work due to flux being incomplete if self.Unknowns.DFDt is not None: self.Unknowns.DFDt.psi_fn = self._constitutive_model.flux.T From b5d57a23acbedbf968e9f71fbd01f4f0e2dd4dde Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Mon, 30 Mar 2026 15:28:37 +1100 Subject: [PATCH 2/2] Update test_1050 to use VE_Stokes for viscoelastic models The test used plain Stokes with ViscoElasticPlasticFlowModel, which silently dropped stress history terms. Now uses VE_Stokes with timestep=0.1, matching the dt_elastic=1/10 already set in the test. Required by PR #95 (VEP solver barrier) and PR #97 (solver unification) which correctly reject VEP models on plain Stokes. Underworld development team with AI support from Claude Code --- tests/test_1050_VEstokesCart.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_1050_VEstokesCart.py b/tests/test_1050_VEstokesCart.py index 26f7e63ec..4fb5edf50 100644 --- a/tests/test_1050_VEstokesCart.py +++ b/tests/test_1050_VEstokesCart.py @@ -52,7 +52,7 @@ def test_stokes_boxmesh(mesh): ) p = uw.discretisation.MeshVariable(r"mathbf{p}", mesh, 1, vtype=uw.VarType.SCALAR, degree=1) - stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p) + stokes = uw.systems.VE_Stokes(mesh, velocityField=u, pressureField=p, order=1) stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 stokes.constitutive_model.Parameters.shear_modulus = 1 @@ -101,7 +101,7 @@ def test_stokes_boxmesh(mesh): stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Front") stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Back") - stokes.solve() + stokes.solve(timestep=0.1) print(f"Mesh dimensions {mesh.dim}", flush=True) stokes.dm.ds.view() @@ -229,7 +229,7 @@ def test_stokes_boxmesh_bc_failure(mesh): ) p = uw.discretisation.MeshVariable(r"mathbf{p}", mesh, 1, vtype=uw.VarType.SCALAR, degree=1) - stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p) + stokes = uw.systems.VE_Stokes(mesh, velocityField=u, pressureField=p, order=1) stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 stokes.constitutive_model.Parameters.shear_modulus = 1 @@ -275,7 +275,7 @@ def test_stokes_boxmesh_bc_failure(mesh): stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Front") stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Back") - stokes.solve() + stokes.solve(timestep=0.1) print(f"Mesh dimensions {mesh.dim}", flush=True) stokes.dm.ds.view()