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
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()