diff --git a/docs/advanced/benchmarks/index.md b/docs/advanced/benchmarks/index.md
new file mode 100644
index 000000000..56ee9ae04
--- /dev/null
+++ b/docs/advanced/benchmarks/index.md
@@ -0,0 +1,14 @@
+---
+title: "Benchmarks"
+---
+
+# Solver Benchmarks
+
+Validation benchmarks comparing Underworld3 solvers against analytical solutions.
+
+```{toctree}
+:maxdepth: 1
+
+ve-oscillatory-shear
+ve-square-wave-shear
+```
diff --git a/docs/advanced/benchmarks/run_ve_square_wave.py b/docs/advanced/benchmarks/run_ve_square_wave.py
new file mode 100644
index 000000000..56a988b82
--- /dev/null
+++ b/docs/advanced/benchmarks/run_ve_square_wave.py
@@ -0,0 +1,244 @@
+"""Variable-timestep VE benchmark: square-wave forcing.
+
+Maxwell material under square-wave shear rate (Fourier series, N harmonics):
+
+    γ̇(t) = (4γ̇₀/π) Σ_{k=1..N} sin((2k-1)ωt) / (2k-1)
+
+Since the Maxwell equation is linear, the analytical stress is the
+superposition of single-frequency solutions:
+
+    σ(t) = Σ_{k=1..N} maxwell_oscillatory(t, η, μ, aₖ, ωₖ)
+
+where aₖ = 4γ̇₀/(π(2k-1)) and ωₖ = (2k-1)ω.
+
+The sharp transitions demand small dt at the edges, large dt in the
+flat regions — testing variable-dt BDF correctness.
+
+Usage:
+    python tests/run_ve_square_wave.py
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+from underworld3.function import expression
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for single-frequency oscillatory Maxwell shear."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def square_wave_analytical(t, eta, mu, gamma_dot_0, omega, n_harmonics=20):
+    """Analytical stress for square-wave forcing via Fourier superposition."""
+    sigma = np.zeros_like(t)
+    for k in range(1, n_harmonics + 1):
+        n = 2 * k - 1  # odd harmonics: 1, 3, 5, ...
+        a_k = 4.0 * gamma_dot_0 / (np.pi * n)
+        omega_k = n * omega
+        sigma += maxwell_oscillatory(t, eta, mu, a_k, omega_k)
+    return sigma
+
+
+def square_wave_shear_rate(t, gamma_dot_0, omega, n_harmonics=20):
+    """Square-wave shear rate via truncated Fourier series."""
+    rate = np.zeros_like(t)
+    for k in range(1, n_harmonics + 1):
+        n = 2 * k - 1
+        rate += 4.0 * gamma_dot_0 / (np.pi * n) * np.sin(n * omega * t)
+    return rate
+
+
+def adaptive_dt(t_current, omega, dt_min, dt_max):
+    """Adaptive timestep: small near square-wave transitions, large on plateaux.
+
+    Transitions occur at t = (2m+1)·π/(2ω) for integer m, i.e. at odd
+    multiples of quarter-period. We use distance to nearest transition
+    to interpolate between dt_min and dt_max.
+    """
+    half_period = np.pi / omega
+    # Phase within half-period [0, half_period)
+    phase = t_current % half_period
+    # Distance to nearest transition (0 or half_period boundary)
+    dist = min(phase, half_period - phase)
+    # Normalise to [0, 1] where 0 = at transition, 1 = mid-plateau
+    frac = dist / (half_period / 2.0)
+    # Smooth interpolation
+    return dt_min + (dt_max - dt_min) * frac**2
+
+
+def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr,
+                    n_harmonics=20, uniform=False):
+    """Run VE square-wave shear box with adaptive or uniform timestep."""
+
+    ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0
+    t_r = ETA / MU
+    omega = De / t_r
+    V0 = 0.5
+    gamma_dot_0 = 2.0 * V0 / H
+    dt_min = dt_min_over_tr * t_r
+    dt_max = dt_max_over_tr * t_r
+    T = 2.0 * np.pi / omega
+    t_end = n_periods * T
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W / 2, -H / 2), maxCoords=(W / 2, H / 2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+
+    # Boundary conditions: simple shear driven by top/bottom velocity
+    # V_top updated numerically each timestep to produce square-wave γ̇
+    V_top = expression(R"V_{\mathrm{top}}", sympy.Float(0.0), "Top boundary velocity")
+
+    stokes.add_dirichlet_bc((V_top, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_top, 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-6
+
+    # Time loop
+    times = []
+    numerical_stress = []
+    timesteps_used = []
+    t_current = 0.0
+    step = 0
+
+    while t_current < t_end:
+        # Determine timestep
+        if uniform:
+            dt = dt_min
+        else:
+            dt = adaptive_dt(t_current, omega, dt_min, dt_max)
+
+        t_next = t_current + dt
+
+        # Update boundary velocity: V_top(t) such that γ̇ = 2·V_top/H = square wave
+        # square_wave_shear_rate returns the actual γ̇, so V_top = γ̇ · H/2
+        gamma_dot_t = square_wave_shear_rate(
+            np.array([t_next]), gamma_dot_0, omega, n_harmonics
+        )[0]
+        V_top.sym = sympy.Float(gamma_dot_t * H / 2.0)
+
+        stokes.constitutive_model.Parameters.dt_elastic = dt
+
+        stokes.solve(zero_init_guess=False, timestep=dt)
+
+        t_current += dt
+        step += 1
+
+        # Extract stress at mesh centre
+        centre = np.array([[0.0, 0.0]])
+        tau_xy = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(tau_xy.flatten()[0])
+
+        times.append(t_current)
+        numerical_stress.append(sigma_xy)
+        timesteps_used.append(dt)
+
+        if step % 50 == 0:
+            ana = square_wave_analytical(
+                np.array([t_current]), ETA, MU, gamma_dot_0, omega, n_harmonics
+            )[0]
+            print(f"  Step {step:4d}: t/t_r = {t_current / t_r:.3f}, "
+                  f"dt/t_r = {dt / t_r:.4f}, "
+                  f"σ_xy = {sigma_xy:.6f}, ana = {ana:.6f}")
+
+    times = np.array(times)
+    numerical_stress = np.array(numerical_stress)
+    timesteps_used = np.array(timesteps_used)
+
+    # Analytical solution
+    analytical_stress = square_wave_analytical(times, ETA, MU, gamma_dot_0, omega, n_harmonics)
+
+    # Error metrics (skip first period for startup transient)
+    mask = times > T
+    if mask.sum() > 0:
+        l2_err = np.sqrt(np.mean((numerical_stress[mask] - analytical_stress[mask]) ** 2))
+        linf_err = np.max(np.abs(numerical_stress[mask] - analytical_stress[mask]))
+    else:
+        l2_err = linf_err = np.nan
+
+    return {
+        "times": times,
+        "numerical": numerical_stress,
+        "analytical": analytical_stress,
+        "timesteps": timesteps_used,
+        "l2_error": l2_err,
+        "linf_error": linf_err,
+        "n_steps": step,
+    }
+
+
+if __name__ == "__main__":
+    De = 1.5
+    order = 2
+    n_periods = 3
+    n_harmonics = 10
+
+    print("=" * 60)
+    print(f"Square-wave VE benchmark: De={De}, order={order}")
+    print(f"  {n_harmonics} Fourier harmonics, {n_periods} periods")
+    print("=" * 60)
+
+    # Run with adaptive dt
+    print("\n--- Adaptive timestep ---")
+    t0 = timer.time()
+    result_adaptive = run_square_wave(
+        order=order, De=De, n_periods=n_periods,
+        dt_min_over_tr=0.02, dt_max_over_tr=0.15,
+        n_harmonics=n_harmonics, uniform=False,
+    )
+    t_adaptive = timer.time() - t0
+    print(f"  {result_adaptive['n_steps']} steps in {t_adaptive:.1f}s")
+    print(f"  L2 error:   {result_adaptive['l2_error']:.6e}")
+    print(f"  Linf error: {result_adaptive['linf_error']:.6e}")
+    print(f"  dt range:   [{result_adaptive['timesteps'].min():.4f}, "
+          f"{result_adaptive['timesteps'].max():.4f}]")
+
+    # Run with uniform dt (reference)
+    print("\n--- Uniform timestep (dt_min) ---")
+    t0 = timer.time()
+    result_uniform = run_square_wave(
+        order=order, De=De, n_periods=n_periods,
+        dt_min_over_tr=0.02, dt_max_over_tr=0.02,
+        n_harmonics=n_harmonics, uniform=True,
+    )
+    t_uniform = timer.time() - t0
+    print(f"  {result_uniform['n_steps']} steps in {t_uniform:.1f}s")
+    print(f"  L2 error:   {result_uniform['l2_error']:.6e}")
+    print(f"  Linf error: {result_uniform['linf_error']:.6e}")
+
+    # Summary
+    print("\n" + "=" * 60)
+    print("Summary:")
+    print(f"  Adaptive: {result_adaptive['n_steps']} steps, "
+          f"L2 = {result_adaptive['l2_error']:.2e}")
+    print(f"  Uniform:  {result_uniform['n_steps']} steps, "
+          f"L2 = {result_uniform['l2_error']:.2e}")
+    ratio = result_adaptive['l2_error'] / result_uniform['l2_error']
+    print(f"  Error ratio (adaptive/uniform): {ratio:.2f}")
+    print(f"  Step savings: {1 - result_adaptive['n_steps'] / result_uniform['n_steps']:.0%}")
+
+    # Save results
+    np.savez(
+        f"tests/ve_square_wave_{n_harmonics}h.npz",
+        adaptive_times=result_adaptive["times"],
+        adaptive_numerical=result_adaptive["numerical"],
+        adaptive_analytical=result_adaptive["analytical"],
+        adaptive_timesteps=result_adaptive["timesteps"],
+        uniform_times=result_uniform["times"],
+        uniform_numerical=result_uniform["numerical"],
+        uniform_analytical=result_uniform["analytical"],
+        uniform_timesteps=result_uniform["timesteps"],
+    )
+    print(f"\nResults saved to tests/ve_square_wave_{n_harmonics}h.npz")
diff --git a/docs/advanced/benchmarks/ve-oscillatory-shear.md b/docs/advanced/benchmarks/ve-oscillatory-shear.md
new file mode 100644
index 000000000..03e95c437
--- /dev/null
+++ b/docs/advanced/benchmarks/ve-oscillatory-shear.md
@@ -0,0 +1,108 @@
+---
+title: "Viscoelastic Oscillatory Shear Benchmark"
+---
+
+# Maxwell Oscillatory Shear
+
+This benchmark validates the viscoelastic Stokes solver against the analytical
+solution for a Maxwell material under oscillatory simple shear.
+
+## Problem Setup
+
+A box with height $H$ and width $2H$ is sheared by imposing time-dependent
+velocities on the top and bottom boundaries:
+
+$$v_x(y=\pm H/2, t) = \pm V_0 \sin(\omega t)$$
+
+The left and right boundaries are free-slip (no vertical velocity). The shear
+rate is $\dot\gamma(t) = \dot\gamma_0 \sin(\omega t)$ where $\dot\gamma_0 = 2V_0/H$.
+
+## Analytical Solution
+
+The Maxwell constitutive law gives:
+
+$$\dot\sigma_{xy} + \frac{\sigma_{xy}}{t_r} = \mu \dot\gamma_0 \sin(\omega t)$$
+
+where $t_r = \eta/\mu$ is the relaxation time. With $\sigma(0) = 0$, the full
+solution (including the startup transient) is:
+
+$$\sigma_{xy}(t) = \frac{\eta \dot\gamma_0}{1 + \text{De}^2}
+\left[\sin(\omega t) - \text{De}\cos(\omega t) + \text{De}\,e^{-t/t_r}\right]$$
+
+where $\text{De} = \omega t_r$ is the Deborah number.
+
+**Steady-state properties** (after transient decays):
+
+- Amplitude: $A = \eta \dot\gamma_0 / \sqrt{1 + \text{De}^2}$
+- Phase lag: $\delta = \arctan(\text{De})$
+
+At $\text{De} = 0$ (viscous limit): $A = \eta\dot\gamma_0$, $\delta = 0$.
+At $\text{De} \to \infty$ (elastic limit): $A \to 0$, $\delta \to 90°$.
+
+## Convergence with BDF Order
+
+The VE solver uses BDF-$k$ time integration ($k = 1, 2, 3$). The convergence
+study (constant shear, $\text{De} = 1$) shows:
+
+| $\Delta t / t_r$ | BDF-1 error | BDF-2 error | BDF-3 error |
+|-------------------|-------------|-------------|-------------|
+| 0.200 | 3.0e-02 | 4.0e-03 | 6.4e-03 |
+| 0.100 | 1.5e-02 | 9.3e-04 | 1.7e-03 |
+| 0.050 | 7.8e-03 | 2.3e-04 | 4.3e-04 |
+| 0.020 | 3.1e-03 | 3.7e-05 | — |
+
+BDF-2 achieves second-order convergence (~4x error reduction per halving) and
+is the recommended default. BDF-1 is first-order. BDF-3 converges at nearly
+second order but with a larger error constant.
+
+## Resolution Study (Oscillatory, De = 5)
+
+At high Deborah number, the oscillation period is short relative to the
+relaxation time, requiring fine time resolution. The plot below shows the
+effect of timestep size at $\text{De} = 5$ ($\omega t_r = 5$, phase lag = 79°):
+
+- **63 pts/period** ($\Delta t/t_r = 0.02$): both orders match analytical
+- **31 pts/period** ($\Delta t/t_r = 0.04$): O1 shows slight amplitude reduction, O2 still accurate
+- **16 pts/period** ($\Delta t/t_r = 0.08$): O1 amplitude visibly damped, O2 remains good
+
+```{note}
+The amplitude reduction at coarse timesteps is numerical dissipation from
+the BDF-1 discrete transfer function, not a cumulative error. The discrete
+steady-state amplitude is a fixed fraction of the analytical amplitude,
+determined by $\omega \Delta t$.
+```
+
+## Running the Benchmarks
+
+```bash
+# Oscillatory validation (De=1.5, order 1 and 2)
+python tests/plot_ve_oscillatory_validation.py
+
+# Resolution study (De=5, three timestep sizes)
+# Saves .npz data files for re-analysis
+python tests/plot_ve_oscillatory_validation.py
+
+# Replot from saved data (no re-running)
+python tests/plot_ve_oscillatory_validation.py --replot
+```
+
+## Notes on `dt_elastic`
+
+The parameter `dt_elastic` on the constitutive model is the elastic relaxation
+timescale used in the BDF discretisation. It controls the effective viscosity
+$\eta_{\text{eff}}$ and the stress history weighting:
+
+- BDF-1: $\eta_{\text{eff}} = \eta\mu\Delta t_e / (\eta + \mu\Delta t_e)$
+- BDF-2: $\eta_{\text{eff}} = 2\eta\mu\Delta t_e / (3\eta + 2\mu\Delta t_e)$
+- BDF-3: $\eta_{\text{eff}} = 6\eta\mu\Delta t_e / (11\eta + 6\mu\Delta t_e)$
+
+When `timestep` is passed to `VE_Stokes.solve()`, it controls the advection
+step for semi-Lagrangian history transport. It does **not** overwrite
+`dt_elastic` — these are independent parameters.
+
+A running-average approach for accumulating history when $\Delta t \ll \Delta t_e$
+was investigated but found to be extremely diffusive for semi-Lagrangian
+transport and is not implemented. To prevent runaway or unstable behaviour
+when timesteps become small (e.g. due to CFL constraints or failure events),
+we advise limiting the minimum effective viscosity, in line with the physics
+of the problem.
diff --git a/docs/advanced/benchmarks/ve-square-wave-shear.md b/docs/advanced/benchmarks/ve-square-wave-shear.md
new file mode 100644
index 000000000..6d80998a4
--- /dev/null
+++ b/docs/advanced/benchmarks/ve-square-wave-shear.md
@@ -0,0 +1,79 @@
+---
+title: "Viscoelastic Square-Wave Shear Benchmark"
+---
+
+# Maxwell Square-Wave Shear
+
+This benchmark validates the viscoelastic Stokes solver with **variable timesteps**
+against an analytical solution for a Maxwell material under square-wave shear forcing.
+
+It tests both the variable-dt BDF-2 coefficients and the PetscDS constants mechanism
+that routes these coefficients to the compiled pointwise functions at runtime.
+
+## Problem Setup
+
+Same geometry as the oscillatory shear benchmark: a box with height $H$ and width $2H$,
+sheared by top/bottom boundary velocities. The shear rate is a truncated Fourier series
+approximation of a square wave:
+
+$$\dot\gamma(t) = \frac{4\dot\gamma_0}{\pi}
+\sum_{k=1}^{N} \frac{\sin\bigl((2k-1)\omega t\bigr)}{2k-1}$$
+
+The sharp transitions between positive and negative shear demand small timesteps
+near the transition points, while the plateaux can use much larger steps. This
+makes it a natural test for adaptive (variable) timestepping.
+
+## Analytical Solution
+
+Since the Maxwell equation is linear, the stress is the superposition of
+single-frequency Maxwell solutions at each Fourier harmonic:
+
+$$\sigma_{xy}(t) = \sum_{k=1}^{N} \sigma_k(t)$$
+
+where each $\sigma_k$ is the oscillatory Maxwell solution with amplitude
+$a_k = 4\dot\gamma_0 / (\pi(2k-1))$ and frequency $\omega_k = (2k-1)\omega$:
+
+$$\sigma_k(t) = \frac{\eta\, a_k}{1 + \text{De}_k^2}
+\left[\sin(\omega_k t) - \text{De}_k\cos(\omega_k t)
++ \text{De}_k\,e^{-t/t_r}\right]$$
+
+with $\text{De}_k = \omega_k t_r$.
+
+## Adaptive Timestep Strategy
+
+The timestep varies between `dt_min` near transitions and `dt_max` on plateaux,
+based on distance to the nearest square-wave transition point:
+
+$$\Delta t = \Delta t_{\min} + (\Delta t_{\max} - \Delta t_{\min})\, f^2$$
+
+where $f \in [0,1]$ measures the normalised distance from the nearest transition.
+
+## Variable-dt BDF Coefficients
+
+With uniform timesteps, BDF-2 uses constant coefficients $[3/2, -2, 1/2]$.
+With variable timesteps (ratio $r = \Delta t_n / \Delta t_{n-1}$), the
+coefficients become:
+
+$$c_0 = \frac{1+2r}{1+r}, \quad c_1 = -(1+r), \quad c_2 = \frac{r^2}{1+r}$$
+
+These coefficients are stored as UWexpressions and updated each step via
+`_update_bdf_coefficients()`, flowing through PetscDS `constants[]` to the
+compiled pointwise functions without JIT recompilation.
+
+## Results (De=1.5, BDF-2, 10 harmonics)
+
+| Run | Steps | L2 Error | Ratio |
+|-----|-------|----------|-------|
+| Adaptive dt | 295 | 9.55e-04 | 1.78x |
+| Uniform dt | 629 | 5.35e-04 | 1.0x |
+
+The adaptive run uses 53% fewer steps at only 1.78x the error.
+
+## Running the Benchmark
+
+```bash
+pixi run -e default python docs/advanced/benchmarks/run_ve_square_wave.py
+```
+
+The script runs both adaptive and uniform timestep cases, prints convergence
+data, and saves results to `.npz` files.
diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py
index 8d02acde9..47293a680 100644
--- a/src/underworld3/constitutive_models.py
+++ b/src/underworld3/constitutive_models.py
@@ -48,6 +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.function.quantities import UWQuantity
 from underworld3.systems.ddt import Lagrangian as Lagrangian_DDt
 
@@ -971,12 +972,12 @@ def viscosity(self):
         # Keep this as an sub-expression for clarity
 
         if inner_self.shear_viscosity_min.sym != -sympy.oo:
-            self._plastic_eff_viscosity._sym = sympy.simplify(
-                sympy.Max(effective_viscosity, inner_self.shear_viscosity_min)
+            self._plastic_eff_viscosity._sym = sympy.Max(
+                effective_viscosity, inner_self.shear_viscosity_min
             )
 
         else:
-            self._plastic_eff_viscosity._sym = sympy.simplify(effective_viscosity)
+            self._plastic_eff_viscosity._sym = effective_viscosity
 
         # Returns an expression that has a different description
         return self._plastic_eff_viscosity
@@ -1086,6 +1087,14 @@ def __init__(self, unknowns, order=1, material_name: str = None):
 
         self._order = order
 
+        # BDF coefficients as UWexpressions — route through PetscDS constants[].
+        # Updated each step by _update_bdf_coefficients() before solve.
+        # Initialised to BDF-1 values: [1, -1, 0, 0].
+        self._bdf_c0 = expression(r"{c_0^{\mathrm{BDF}}}", sympy.Integer(1), "BDF leading coefficient")
+        self._bdf_c1 = expression(r"{c_1^{\mathrm{BDF}}}", sympy.Integer(-1), "BDF history coefficient 1")
+        self._bdf_c2 = expression(r"{c_2^{\mathrm{BDF}}}", sympy.Integer(0), "BDF history coefficient 2")
+        self._bdf_c3 = expression(r"{c_3^{\mathrm{BDF}}}", sympy.Integer(0), "BDF history coefficient 3")
+
         self._reset()
 
         super().__init__(unknowns, material_name=material_name)
@@ -1198,32 +1207,15 @@ def ve_effective_viscosity(inner_self):
             if inner_self.shear_modulus == sympy.oo:
                 return inner_self.shear_viscosity_0
 
-            # Note, 1st order only here but we should add higher order versions of this
-
-            # 1st Order version (default)
-            if inner_self._owning_model.order != 2:
-                el_eff_visc = (
-                    inner_self.shear_viscosity_0
-                    * inner_self.shear_modulus
-                    * inner_self.dt_elastic
-                    / (
-                        inner_self.shear_viscosity_0
-                        + inner_self.dt_elastic * inner_self.shear_modulus
-                    )
-                )
+            # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt)
+            # c0 is a UWexpression routed through PetscDS constants[],
+            # updated each step by _update_bdf_coefficients().
+            eta = inner_self.shear_viscosity_0
+            mu = inner_self.shear_modulus
+            dt_e = inner_self.dt_elastic
+            c0 = inner_self._owning_model._bdf_c0
 
-            # 2nd Order version (need to ask for this one)
-            else:
-                el_eff_visc = (
-                    2
-                    * inner_self.shear_viscosity_0
-                    * inner_self.shear_modulus
-                    * inner_self.dt_elastic
-                    / (
-                        3 * inner_self.shear_viscosity_0
-                        + 2 * inner_self.dt_elastic * inner_self.shear_modulus
-                    )
-                )
+            el_eff_visc = eta * mu * dt_e / (c0 * eta + mu * dt_e)
 
             inner_self._ve_effective_viscosity.sym = el_eff_visc
 
@@ -1252,6 +1244,46 @@ def order(self, value):
         self._reset()
         return
 
+    @property
+    def effective_order(self):
+        """Effective order accounting for DDt history startup.
+
+        During the first few timesteps, the DDt may not have enough history
+        to support the requested order. This property returns the lower of
+        the requested order and the DDt's effective order (which ramps from
+        1 to self.order as history accumulates).
+        """
+        if self.Unknowns is not None and self.Unknowns.DFDt is not None:
+            return min(self._order, self.Unknowns.DFDt.effective_order)
+        return self._order
+
+    def _update_bdf_coefficients(self):
+        """Update BDF coefficient UWexpressions from current dt_elastic and DDt history.
+
+        Call this before each solve so that the constants[] array carries the
+        correct coefficients to the compiled pointwise functions. The coefficient
+        UWexpressions (_bdf_c0..c3) are referenced symbolically in ve_effective_viscosity,
+        E_eff, and stress() — their numeric values flow through PetscDSSetConstants.
+        """
+        if self.Unknowns is not None and self.Unknowns.DFDt is not None:
+            dt_current = self.Parameters.dt_elastic
+            if hasattr(dt_current, 'sym'):
+                dt_current = dt_current.sym
+            coeffs = _bdf_coefficients(
+                self.effective_order, dt_current, self.Unknowns.DFDt._dt_history
+            )
+        else:
+            coeffs = _bdf_coefficients(self.effective_order, None, [])
+
+        # Pad to length 4
+        while len(coeffs) < 4:
+            coeffs.append(sympy.Integer(0))
+
+        self._bdf_c0.sym = coeffs[0]
+        self._bdf_c1.sym = coeffs[1]
+        self._bdf_c2.sym = coeffs[2]
+        self._bdf_c3.sym = coeffs[3]
+
     # The following should have no setters
     @property
     def stress_star(self):
@@ -1287,20 +1319,13 @@ def E_eff(self):
         if self.Unknowns.DFDt is not None:
 
             if self.is_elastic:
-                if self.order != 2:
-                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
-                    E += stress_star / (
-                        2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus
-                    )
+                mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus
+                # BDF history coefficients as UWexpressions (route through constants[])
+                bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3]
 
-                else:
-                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
-                    stress_2star = self.Unknowns.DFDt.psi_star[1].sym
-                    E += stress_star / (
-                        self.Parameters.dt_elastic * self.Parameters.shear_modulus
-                    ) - stress_2star / (
-                        4 * self.Parameters.dt_elastic * self.Parameters.shear_modulus
-                    )
+                # History contribution: -Σ cᵢ·σ_star[i-1] / (2·μ·dt)
+                for i in range(self.Unknowns.DFDt.order):
+                    E += -bdf_cs[i] * self.Unknowns.DFDt.psi_star[i].sym / (2 * mu_dt)
 
         self._E_eff.sym = E
 
@@ -1365,16 +1390,10 @@ def _plastic_effective_viscosity(self):
         if parameters.yield_stress == sympy.oo:
             return sympy.oo
 
-        Edot = self.Unknowns.E
-        if self.Unknowns.DFDt is not None:
-
-            ## First order ...
-            stress_star = self.Unknowns.DFDt.psi_star[0]
-
-            if self.is_elastic:
-                Edot += stress_star.sym / (
-                    2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus
-                )
+        # Use the effective strain rate (including elastic history) for the
+        # yield criterion. This must use the same order-dependent BDF
+        # coefficients as the stress formula.
+        Edot = self.E_eff.sym
 
         strainrate_inv_II = expression(
             R"{\dot\varepsilon_{II}'}",
@@ -1468,8 +1487,6 @@ def flux(self):
         #     plastic_scale_factor = sympy.Max(1, self.plastic_overshoot())
         #     stress /= plastic_scale_factor
 
-        stress = sympy.simplify(stress)
-
         return stress
 
     def stress_projection(self):
@@ -1492,8 +1509,6 @@ def stress_projection(self):
                     / (self.Parameters.dt_elastic * self.Parameters.shear_modulus)
                 )
 
-        stress = sympy.simplify(stress)
-
         return stress
 
     def stress(self):
@@ -1508,34 +1523,15 @@ def stress(self):
         if self.Unknowns.DFDt is not None:
 
             if self.is_elastic:
-                if self.order != 2:
-                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
-                    stress += (
-                        2
-                        * self.viscosity
-                        * (
-                            stress_star
-                            / (2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus)
-                        )
-                    )
+                mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus
+                bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3]
 
-                else:
-                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
-                    stress_2star = self.Unknowns.DFDt.psi_star[1].sym
-
-                    stress += (
-                        2
-                        * self.viscosity
-                        * (
-                            stress_star
-                            / (self.Parameters.dt_elastic * self.Parameters.shear_modulus)
-                            - stress_2star
-                            / (4 * self.Parameters.dt_elastic * self.Parameters.shear_modulus)
-                        )
+                # History contribution: 2·η_eff · (-Σ cᵢ·σ_star[i-1]) / (2·μ·dt)
+                for i in range(self.Unknowns.DFDt.order):
+                    stress += 2 * self.viscosity * (
+                        -bdf_cs[i] * self.Unknowns.DFDt.psi_star[i].sym / (2 * mu_dt)
                     )
 
-        stress = sympy.simplify(stress)
-
         return stress
 
     # def eff_edot(self):
diff --git a/src/underworld3/constitutive_models_new.py b/src/underworld3/constitutive_models_new.py
deleted file mode 100644
index d58be0598..000000000
--- a/src/underworld3/constitutive_models_new.py
+++ /dev/null
@@ -1,1702 +0,0 @@
-r"""
-Alternative constitutive model implementations.
-
-.. deprecated::
-    This module contains a failed attempt to rebuild the constitutive models
-    and is scheduled for removal. Use :mod:`underworld3.constitutive_models`
-    instead for all constitutive model needs.
-
-See Also
---------
-underworld3.constitutive_models : Primary constitutive model module (use this).
-"""
-from typing_extensions import Self
-import sympy
-from sympy import sympify
-from sympy.vector import gradient, divergence
-import numpy as np
-
-from typing import Optional, Callable
-from typing import NamedTuple, Union
-
-from petsc4py import PETSc
-
-import underworld3 as uw
-from underworld3.utilities._api_tools import uw_object
-from underworld3.systems import SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt
-from underworld3.swarm import IndexSwarmVariable
-import underworld3.timing as timing
-
-
-class Constitutive_Model(uw_object):
-    r"""
-    Base class for constitutive models.
-
-    Constitutive laws relate gradients in the unknowns to fluxes of quantities
-    (e.g., heat fluxes related to temperature gradients via thermal conductivity).
-
-    For scalar problems:
-
-    .. math::
-
-        q_i = k_{ij} \frac{\partial T}{\partial x_j}
-
-    where :math:`k_{ij}` are the constitutive parameters. The template assumes
-    :math:`k_{ij} = \delta_{ij}` (identity).
-
-    For vector problems (e.g., Stokes):
-
-    .. math::
-
-        t_{ij} = c_{ijkl} \frac{\partial u_k}{\partial x_l}
-
-    Usually written with symmetrized gradients:
-
-    .. math::
-
-        t_{ij} = c_{ijkl} \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l}
-        + \frac{\partial u_l}{\partial x_k} \right]
-
-    where :math:`c_{ijkl}` are the constitutive parameters. The template assumes
-    :math:`c_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})`,
-    the 4th-rank identity tensor with flux and gradient symmetry.
-    """
-
-    @timing.routine_timer_decorator
-    def __init__(
-        self,
-        dim: int,
-        u_dim: int,
-    ):
-        # Define / identify the various properties in the class but leave
-        # the implementation to child classes. The constitutive tensor is
-        # defined as a template here, but should be instantiated via class
-        # properties as required.
-
-        # We provide a function that converts gradients / gradient history terms
-        # into the relevant flux term.
-
-        self.dim = dim
-        self.u_dim = u_dim
-        self._solver = None
-
-        self.Parameters = self._Parameters()
-        self.Parameters._solver = None
-        self.Parameters._reset = self._reset
-
-        self._material_properties = None
-
-        ## Default consitutive tensor is the identity
-
-        if self.u_dim == 1:
-            self._c = sympy.Matrix.eye(self.dim)
-        else:  # vector problem
-            self._c = uw.maths.tensor.rank4_identity(self.dim)
-
-        self._C = None
-
-        super().__init__()
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-        """
-
-        def __init__(inner_self, k=1):
-            inner_self._k = k
-            inner_self._solver = None
-
-        @property
-        def k(inner_self):
-            return inner_self._k
-
-        @k.setter
-        def k(inner_self, value):
-            inner_self._k = value
-            inner_self._reset()
-
-        """
-        @property
-        def material_properties(self):
-
-            return self._material_properties
-
-        @material_properties.setter
-        def material_properties(self, properties):
-
-            if isinstance(properties, self.Parameters):
-                self._material_properties = properties
-            else:
-                name = self.__class__.__name__
-                raise RuntimeError(f"Use {name}.material_properties = {name}.Parameters(...) ")
-
-            d = self.dim
-            self._build_c_tensor()
-
-            if isinstance(self._solver, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)):
-                self._solver.is_setup = False
-
-            return
-        """
-
-    @property
-    def solver(self):
-        """Each constitutive relationship can, optionally, be associated with one solver object.
-        and a solver object _requires_ a constitive relationship to be defined."""
-        return self._solver
-
-    @solver.setter
-    def solver(self, solver_object):
-        if isinstance(solver_object, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)):
-            self._solver = solver_object
-            self.Parameters._solver = solver_object
-            self._solver.is_setup = False
-
-    ## Properties on all sub-classes
-
-    @property
-    def C(self):
-        """The matrix form of the constitutive model (the `c` property)
-        that relates fluxes to gradients.
-        For scalar problem, this is the matrix representation of the rank 2 tensor.
-        For vector problems, the Mandel form of the rank 4 tensor is returned.
-        NOTE: this is an immutable object that is _a view_ of the underlying tensor
-        """
-
-        d = self.dim
-        rank = len(self.c.shape)
-
-        if rank == 2:
-            return sympy.Matrix(self._c).as_immutable()
-        else:
-            return uw.maths.tensor.rank4_to_mandel(self._c, d).as_immutable()
-
-    @property
-    def c(self):
-        """The tensor form of the constitutive model that relates fluxes to gradients. In scalar
-        problems, `c` and `C` are equivalent (matrices), but in vector problems, `c` is a
-        rank 4 tensor. NOTE: `c` is the canonical form of the constitutive relationship.
-        """
-
-        return self._c.as_immutable()
-
-    def flux(
-        self,
-        ddu: sympy.Matrix = None,
-        ddu_dt: sympy.Matrix = None,
-        u: sympy.Matrix = None,  # may be needed in the case of cylindrical / spherical
-        u_dt: sympy.Matrix = None,
-    ):
-        """Computes the effect of the constitutive tensor on the gradients of the unknowns.
-        (always uses the `c` form of the tensor). In general cases, the history of the gradients
-        may be required to evaluate the flux.
-        """
-
-        c = self.c
-        rank = len(c.shape)
-
-        # tensor multiplication
-
-        if rank == 2:
-            flux = c * ddu.T
-        else:  # rank==4
-            flux = sympy.tensorcontraction(
-                sympy.tensorcontraction(sympy.tensorproduct(c, ddu), (3, 5)), (0, 1)
-            )
-
-        return sympy.Matrix(flux)
-
-    def flux_1d(
-        self,
-        ddu: sympy.Matrix = None,
-        ddu_dt: sympy.Matrix = None,
-        u: sympy.Matrix = None,  # may be needed in the case of cylindrical / spherical
-        u_dt: sympy.Matrix = None,
-    ):
-        """Computes the effect of the constitutive tensor on the gradients of the unknowns.
-        (always uses the `c` form of the tensor). In general cases, the history of the gradients
-        may be required to evaluate the flux. Returns the Voigt form that is flattened so as to
-        match the PETSc field storage pattern for symmetric tensors.
-        """
-
-        flux = self.flux(ddu, ddu_dt, u, u_dt)
-
-        assert (
-            flux.is_symmetric()
-        ), "The conversion to Voigt form is only defined for symmetric tensors in underworld"
-
-        return uw.maths.tensor.rank2_to_voigt(flux, dim=self.dim)
-
-    def _reset(self):
-        d = self.dim
-        self._build_c_tensor()
-
-        if isinstance(self._solver, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)):
-            self._solver.is_setup = False
-
-        return
-
-    def _build_c_tensor(self):
-        """Return the identity tensor of appropriate rank (e.g. for projections)"""
-
-        d = self.dim
-        self._c = self.Parameters.k * uw.maths.tensor.rank4_identity(d)
-
-        return
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-        from textwrap import dedent
-
-        display(
-            Markdown(
-                rf"This consititutive model is formulated for {self.dim} dimensional equations"
-            )
-        )
-
-
-class ViscousFlowModel(Constitutive_Model):
-    r"""
-    Viscous flow constitutive model.
-
-    .. math::
-
-        \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l}
-        + \frac{\partial u_l}{\partial x_k} \right]
-
-    where :math:`\eta` is the viscosity—a scalar constant, sympy function,
-    or mesh variable. This gives an isotropic (but not necessarily homogeneous
-    or linear) relationship between :math:`\tau` and velocity gradients. You can
-    also supply :math:`\eta_{IJ}` (Mandel form) or :math:`\eta_{ijkl}` (rank-4 tensor).
-
-    The Mandel constitutive matrix is in ``viscous_model.C`` and the rank-4
-    tensor form is in ``viscous_model.c``.
-
-    Examples
-    --------
-    >>> viscous_model = ViscousFlowModel(dim)
-    >>> viscous_model.material_properties = viscous_model.Parameters(viscosity=viscosity_fn)
-    >>> solver.constitutive_model = viscous_model
-    >>> tau = viscous_model.flux(gradient_matrix)
-    """
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-        """
-
-        def __init__(
-            inner_self,
-            viscosity: Union[float, sympy.Function] = None,
-        ):
-            if viscosity is None:
-                viscosity = sympy.sympify(1)
-
-            inner_self._viscosity = sympy.sympify(viscosity)
-
-        @property
-        def viscosity(inner_self):
-            return inner_self._viscosity
-
-        @viscosity.setter
-        def viscosity(inner_self, value: Union[float, sympy.Function]):
-            inner_self._viscosity = value
-            inner_self._reset()
-
-    def __init__(self, dim):
-        u_dim = dim
-        super().__init__(dim, u_dim)
-
-        return
-
-    def _build_c_tensor(self):
-        """For this constitutive law, we expect just a viscosity function"""
-
-        d = self.dim
-        viscosity = self.Parameters.viscosity
-
-        try:
-            self._c = 2 * uw.maths.tensor.rank4_identity(d) * viscosity
-        except:
-            d = self.dim
-            dv = uw.maths.tensor.idxmap[d][0]
-            if isinstance(viscosity, sympy.Matrix) and viscosity.shape == (dv, dv):
-                self._c = 2 * uw.maths.tensor.mandel_to_rank4(viscosity, d)
-            elif isinstance(viscosity, sympy.Array) and viscosity.shape == (d, d, d, d):
-                self._c = 2 * viscosity
-            else:
-                raise RuntimeError(
-                    "Viscosity is not a known type (scalar, Mandel matrix, or rank 4 tensor"
-                )
-        return
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-
-        super()._object_viewer()
-
-        ## feedback on this instance
-        display(Latex(r"$\quad\eta = $ " + sympy.sympify(self.Parameters.viscosity)._repr_latex_()))
-
-
-class ViscoPlasticFlowModel(ViscousFlowModel):
-    r"""
-    Viscoplastic flow constitutive model with yield stress.
-
-    .. math::
-
-        \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l}
-        + \frac{\partial u_l}{\partial x_k} \right]
-
-    where :math:`\eta` is the viscosity—a scalar constant, sympy function,
-    or mesh variable. This gives an isotropic relationship between :math:`\tau`
-    and velocity gradients. You can also supply :math:`\eta_{IJ}` (Mandel form)
-    or :math:`\eta_{ijkl}` (rank-4 tensor).
-
-    In a viscoplastic model, the viscosity is defined to cap the overall stress
-    at the *yield stress*. This assumes the yield stress is a scalar limit on
-    the 2nd invariant of the stress. Anisotropic models require careful yield
-    surface definition—only a subset of cases is available.
-
-    The Mandel matrix is in ``viscoplastic_model.C`` and the rank-4 tensor in
-    ``viscoplastic_model.c``.
-
-    Notes
-    -----
-    If ``not~yet~defined`` appears in the effective viscosity, not all required
-    functions have been set. The model defaults to standard viscous behavior
-    if yield terms are not specified.
-
-    Examples
-    --------
-    >>> viscoplastic_model = ViscoPlasticFlowModel(dim)
-    >>> viscoplastic_model.material_properties = viscoplastic_model.Parameters(
-    ...     viscosity=viscosity_fn,
-    ...     yield_stress=yieldstress_fn,
-    ...     min_viscosity=min_viscosity_fn,
-    ...     max_viscosity=max_viscosity_fn,
-    ...     strain_rate_II=strain_rate_inv_fn
-    ... )
-    >>> solver.constitutive_model = viscoplastic_model
-    >>> tau = viscoplastic_model.flux(gradient_matrix)
-    """
-
-    # Init for VP class (not needed ??)
-    def __init__(self, dim):
-        super().__init__(dim)
-
-        return
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-
-        `sympy.oo` (infinity) for default values ensures that sympy.Min simplifies away
-        the conditionals when they are not required
-        """
-
-        def __init__(
-            inner_self,
-            materialIndex: Union[uw.swarm.SwarmVariable, uw.discretisation.MeshVariable] = None,
-            shear_viscosity_0: Union[list, sympy.Function] = [1],
-            shear_viscosity_min: Union[list, sympy.Function] = [-sympy.oo],
-            shear_viscosity_max: Union[list, sympy.Function] = [sympy.oo],
-            yield_stress: Union[list, sympy.Function] = [sympy.oo],
-            yield_stress_min: Union[list, sympy.Function] = [-sympy.oo],
-            strainrate_inv_II: sympy.Function = sympy.oo,
-            strainrate_inv_II_min: float = 0.0,
-            averaging_method: str = "HA",
-        ):
-            if strainrate_inv_II is sympy.oo:
-                strainrate_inv_II = sympy.symbols(
-                    r"\left|\dot\epsilon\right|\rightarrow\textrm{not\ defined}"
-                )
-
-            inner_self._shear_viscosity_0 = sympy.sympify(shear_viscosity_0)
-            inner_self._shear_viscosity_min = sympy.sympify(shear_viscosity_min)
-            inner_self._shear_viscosity_max = sympy.sympify(shear_viscosity_max)
-
-            inner_self._yield_stress = sympy.sympify(yield_stress)
-            inner_self._yield_stress_min = sympy.sympify(yield_stress_min)
-
-            inner_self._strainrate_inv_II = sympy.sympify(strainrate_inv_II)
-            inner_self._strainrate_inv_II_min = sympy.sympify(strainrate_inv_II_min)
-
-            inner_self._averaging_method = averaging_method
-            inner_self._materialIndex = materialIndex
-
-            return
-
-        @property
-        def shear_viscosity_0(inner_self):
-            return inner_self._shear_viscosity_0
-
-        @shear_viscosity_0.setter
-        def shear_viscosity_0(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_0 = value
-            inner_self._reset()
-
-        @property
-        def shear_viscosity_min(inner_self):
-            return inner_self._shear_viscosity_min
-
-        @shear_viscosity_min.setter
-        def shear_viscosity_min(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_min = value
-            inner_self._reset()
-
-        @property
-        def shear_viscosity_max(inner_self):
-            return inner_self._shear_viscosity_max
-
-        @shear_viscosity_max.setter
-        def shear_viscosity_max(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_max = value
-            inner_self._reset()
-
-        @property
-        def yield_stress(inner_self):
-            return inner_self._yield_stress
-
-        @yield_stress.setter
-        def yield_stress(inner_self, value: Union[list, sympy.Function]):
-            inner_self._yield_stress = value
-            inner_self._reset()
-
-        @property
-        def yield_stress_min(inner_self):
-            return inner_self._yield_stress_min
-
-        @yield_stress_min.setter
-        def yield_stress_min(inner_self, value: Union[list, sympy.Function]):
-            inner_self._yield_stress_min = value
-            inner_self._reset()
-
-        @property
-        def strainrate_inv_II(inner_self):
-            return inner_self._strainrate_inv_II
-
-        @strainrate_inv_II.setter
-        def strainrate_inv_II(inner_self, value: sympy.Function):
-            inner_self._strainrate_inv_II = value
-            inner_self._reset()
-
-        @property
-        def strainrate_inv_II_min(inner_self):
-            return inner_self._strainrate_inv_II_min
-
-        @strainrate_inv_II_min.setter
-        def strainrate_inv_II_min(inner_self, value: float):
-            inner_self._strainrate_inv_II_min = sympy.sympify(value)
-            inner_self._reset()
-
-        @property
-        def averaging_method(inner_self):
-            return inner_self._averaging_method
-
-        @averaging_method.setter
-        def averaging_method(inner_self, value: str):
-            inner_self._averaging_method = value
-            inner_self._reset()
-
-        ### Getter and setter for internel mask Variable
-        @property
-        def materialIndex(inner_self):
-            return inner_self._materialIndex
-
-        @materialIndex.setter
-        def materialIndex(inner_self, indexVar):
-            # error checking, only support IndexSwarmVariables for now
-            if isinstance(indexVar, uw.swarm.IndexSwarmVariable):
-                inner_self._materialIndex = indexVar
-                inner_self._reset()
-
-        # This has no setter !!
-        @property
-        def plastic_eff_viscosity(inner_self):
-            import warnings
-
-            if (
-                type(inner_self.yield_stress) != np.ndarray
-                and type(inner_self.yield_stress) != list
-            ):
-                inner_self.yield_stress = list([inner_self.yield_stress])
-            if (
-                type(inner_self.yield_stress_min) != np.ndarray
-                and type(inner_self.yield_stress_min) != list
-            ):
-                inner_self.yield_stress_min = list([inner_self.yield_stress_min])
-
-            if inner_self.materialIndex == None:
-                warnings.warn(
-                    "materialIndex not specified, using the first value for each parameter",
-                    stacklevel=2,
-                )
-
-                yield_stress = inner_self.yield_stress[0]
-                yield_stress_min = inner_self.yield_stress_min[0]
-
-                if yield_stress_min == -sympy.oo:
-                    yield_stress_fn = yield_stress
-                else:
-                    yield_stress_fn = sympy.Max(yield_stress, yield_stress_min)
-
-                if yield_stress_fn == sympy.oo or inner_self.strainrate_inv_II == sympy.oo:
-                    pl_effective_viscosity = sympy.oo
-
-                else:
-                    pl_effective_viscosity = yield_stress_fn / (
-                        (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min
-                    )
-
-            else:
-                ### creates list of values that has the same length as the material index
-                if len(inner_self.yield_stress) != inner_self.materialIndex.indices:
-                    if len(inner_self.yield_stress) > 1:
-                        warnings.warn(
-                            f"Number of values in yield_stress ({len(inner_self.yield_stress)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.yield_stress = list(
-                        np.repeat(inner_self.yield_stress[0], inner_self.materialIndex.indices)
-                    )
-
-                if len(inner_self.yield_stress_min) != inner_self.materialIndex.indices:
-                    if len(inner_self.yield_stress_min) > 1:
-                        warnings.warn(
-                            f"Number of values in yield_stress ({len(inner_self.yield_stress_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.yield_stress_min = list(
-                        np.repeat(
-                            inner_self.yield_stress_min[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                yield_stress = sympy.Matrix(inner_self.yield_stress)
-                yield_stress_min = sympy.Matrix(inner_self.yield_stress_min)
-
-                if yield_stress_min[0] == -sympy.oo:
-                    yield_stress_fn = yield_stress
-                else:
-                    yield_stress_list = []
-
-                    for i in range(len(yield_stress)):
-                        yield_stress_list.append(sympy.Max(yield_stress[i], yield_stress_min[i]))
-
-                    yield_stress_fn = sympy.Matrix(yield_stress_list)
-
-                if yield_stress_fn[0] == sympy.oo or inner_self.strainrate_inv_II == sympy.oo:
-                    pl_effective_viscosity = list(
-                        np.repeat(sympy.oo, inner_self.materialIndex.indices)
-                    )
-
-                else:
-                    pl_effective_viscosity = yield_stress_fn / (
-                        (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min
-                    )
-
-            return pl_effective_viscosity
-
-        # This has no setter !!
-        @property
-        def viscosity(inner_self):
-            import warnings
-
-            if (
-                type(inner_self.shear_viscosity_0) != np.ndarray
-                and type(inner_self.shear_viscosity_0) != list
-            ):
-                inner_self.shear_viscosity_0 = list([inner_self.shear_viscosity_0])
-
-            if (
-                type(inner_self.shear_viscosity_min) != np.ndarray
-                and type(inner_self.shear_viscosity_min) != list
-            ):
-                inner_self.shear_viscosity_min = list([inner_self.shear_viscosity_min])
-            if (
-                type(inner_self.shear_viscosity_max) != np.ndarray
-                and type(inner_self.shear_viscosity_max) != list
-            ):
-                inner_self.shear_viscosity_max = list([inner_self.shear_viscosity_max])
-
-            if inner_self.materialIndex == None:
-                warnings.warn(
-                    "materialIndex not specified, using the first value for each parameter",
-                    stacklevel=2,
-                )
-
-                shear_viscosity_min = inner_self.shear_viscosity_min[0]
-                shear_viscosity_max = inner_self.shear_viscosity_max[0]
-
-                shear_viscosity_0 = inner_self.shear_viscosity_0[0]
-
-                yield_visc = inner_self.plastic_eff_viscosity
-
-                if inner_self.averaging_method.casefold() == "min":
-                    if yield_visc != 0:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(
-                                shear_viscosity_min,
-                                sympy.Min(yield_visc, shear_viscosity_0),
-                            ),
-                        )
-                    else:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(shear_viscosity_min, shear_viscosity_0),
-                        )
-                else:
-                    if yield_visc != 0:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(
-                                shear_viscosity_min,
-                                1.0 / ((1.0 / shear_viscosity_0) + (1.0 / yield_visc)),
-                            ),
-                        )
-                    else:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(shear_viscosity_min, shear_viscosity_0),
-                        )
-
-            else:
-                if len(inner_self.shear_viscosity_0) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_0) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_0 ({len(inner_self.shear_viscosity_0)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_0 = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_0[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                if len(inner_self.shear_viscosity_min) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_min) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_min ({len(inner_self.shear_viscosity_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_min = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_min[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                if len(inner_self.shear_viscosity_max) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_max) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_max ({len(inner_self.shear_viscosity_max)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_max = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_max[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                shear_viscosity_min = sympy.Matrix(inner_self.shear_viscosity_min)
-                shear_viscosity_max = sympy.Matrix(inner_self.shear_viscosity_max)
-
-                shear_viscosity_0 = inner_self.shear_viscosity_0
-
-                yield_visc = inner_self.plastic_eff_viscosity
-
-                viscosity_list = []
-                for i in range(inner_self.materialIndex.indices):
-                    if inner_self.averaging_method.casefold() == "min":
-                        if yield_visc[i] != 0:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(
-                                        shear_viscosity_min[i],
-                                        sympy.Min(yield_visc[i], shear_viscosity_0[i]),
-                                    ),
-                                )
-                            )
-                        else:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(shear_viscosity_min[i], shear_viscosity_0[i]),
-                                )
-                            )
-                    else:
-                        if yield_visc[i] != 0:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(
-                                        shear_viscosity_min[i],
-                                        1 / ((1 / shear_viscosity_0[i]) + (1 / yield_visc[i])),
-                                    ),
-                                )
-                            )
-                        else:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(shear_viscosity_min[i], shear_viscosity_0[i]),
-                                )
-                            )
-
-                    viscosity_fn = sympy.Matrix(viscosity_list)
-
-                effective_viscosity = inner_self.materialIndex.sym.T.dot(viscosity_fn)
-
-            return effective_viscosity
-
-        ## ===== End of parameters sub_class
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-
-        super()._object_viewer()
-
-        ## feedback on this instance
-        display(
-            Latex(
-                r"$\quad\eta_\textrm{0} = $ "
-                + sympy.sympify(self.Parameters.shear_viscosity_0)._repr_latex_()
-            ),
-            Latex(
-                r"$\quad\tau_\textrm{y} = $ "
-                + sympy.sympify(self.Parameters.yield_stress)._repr_latex_(),
-            ),
-            Latex(
-                r"$\quad|\dot\epsilon| = $ "
-                + sympy.sympify(self.Parameters.strainrate_inv_II)._repr_latex_(),
-            ),
-            ## Todo: add all the other properties in here
-        )
-
-
-class ViscoElasticPlasticFlowModel(Constitutive_Model):
-    r"""
-    Viscoelastic-plastic flow constitutive model.
-
-    .. math::
-
-        \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l}
-        + \frac{\partial u_l}{\partial x_k} \right]
-
-    where :math:`\eta` is the viscosity—a scalar constant, sympy function,
-    or mesh variable. This gives an isotropic relationship between :math:`\tau`
-    and velocity gradients. You can also supply :math:`\eta_{IJ}` (Mandel form)
-    or :math:`\eta_{ijkl}` (rank-4 tensor).
-
-    The Mandel matrix is in ``viscoelastic_model.C`` and the rank-4 tensor in
-    ``viscoelastic_model.c``.
-
-    Examples
-    --------
-    >>> viscoelastic_model = ViscoElasticFlowModel(dim)
-    >>> viscoelastic_model.material_properties = viscoelastic_model.Parameters(
-    ...     viscosity=viscosity_fn
-    ... )
-    >>> solver.constitutive_model = viscoelastic_model
-    >>> tau = viscoelastic_model.flux(gradient_matrix)
-    """
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-        """
-
-        def __init__(
-            inner_self,
-            materialIndex: Union[uw.swarm.SwarmVariable, uw.discretisation.MeshVariable] = None,
-            shear_viscosity_0: Union[list, sympy.Function] = [1],
-            shear_modulus: Union[list, sympy.Function] = [sympy.oo],
-            shear_viscosity_min: Union[list, sympy.Function] = [-sympy.oo],
-            shear_viscosity_max: Union[list, sympy.Function] = [sympy.oo],
-            yield_stress: Union[list, sympy.Function] = [sympy.oo],
-            yield_stress_min: Union[list, sympy.Function] = [-sympy.oo],
-            strainrate_inv_II: sympy.Function = sympy.oo,
-            strainrate_inv_II_min: float = 0.0,
-            averaging_method: str = "HA",
-            stress_star: sympy.Function = None,
-            stress_star_star: sympy.Function = None,
-            dt_elastic: Union[float, sympy.Function] = [sympy.oo],
-        ):
-            if strainrate_inv_II is None:
-                strainrate_inv_II = sympy.symbols(
-                    r"\left|\dot\epsilon\right|\rightarrow\textrm{not\ defined}"
-                )
-
-            if stress_star is None:
-                stress_star = sympy.symbols(r"\sigma^*\rightarrow\textrm{not\ defined}")
-
-            inner_self._shear_viscosity_0 = sympy.sympify(shear_viscosity_0)
-            inner_self._shear_modulus = sympy.sympify(shear_modulus)
-            inner_self._dt_elastic = sympy.sympify(dt_elastic)
-            inner_self._yield_stress = sympy.sympify(yield_stress)
-            inner_self._yield_stress_min = sympy.sympify(yield_stress_min)
-            inner_self._shear_viscosity_min = sympy.sympify(shear_viscosity_min)
-            inner_self._shear_viscosity_max = sympy.sympify(shear_viscosity_max)
-            inner_self._strainrate_inv_II = sympy.sympify(strainrate_inv_II)
-            inner_self._stress_star = sympy.sympify(stress_star)
-            inner_self._stress_star_star = sympy.sympify(stress_star_star)
-            inner_self._strainrate_inv_II_min = sympy.sympify(strainrate_inv_II_min)
-
-            inner_self._averaging_method = averaging_method
-            inner_self._materialIndex = materialIndex
-
-            return
-
-        @property
-        def shear_viscosity_0(inner_self):
-            return inner_self._shear_viscosity_0
-
-        @shear_viscosity_0.setter
-        def shear_viscosity_0(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_0 = value
-            inner_self._reset()
-
-        @property
-        def shear_modulus(inner_self):
-            return inner_self._shear_modulus
-
-        @shear_modulus.setter
-        def shear_modulus(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_modulus = value
-            inner_self._reset()
-
-        @property
-        def dt_elastic(inner_self):
-            return inner_self._dt_elastic
-
-        @dt_elastic.setter
-        def dt_elastic(inner_self, value: Union[list, sympy.Function]):
-            inner_self._dt_elastic = value
-            inner_self._reset()
-
-        @property
-        def shear_viscosity_min(inner_self):
-            return inner_self._shear_viscosity_min
-
-        @shear_viscosity_min.setter
-        def shear_viscosity_min(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_min = value
-            inner_self._reset()
-
-        @property
-        def shear_viscosity_max(inner_self):
-            return inner_self._shear_viscosity_max
-
-        @shear_viscosity_max.setter
-        def shear_viscosity_max(inner_self, value: Union[list, sympy.Function]):
-            inner_self._shear_viscosity_max = value
-            inner_self._reset()
-
-        @property
-        def yield_stress(inner_self):
-            return inner_self._yield_stress
-
-        @yield_stress.setter
-        def yield_stress(inner_self, value: Union[list, sympy.Function]):
-            inner_self._yield_stress = value
-            inner_self._reset()
-
-        @property
-        def yield_stress_min(inner_self):
-            return inner_self._yield_stress_min
-
-        @yield_stress_min.setter
-        def yield_stress_min(inner_self, value: Union[list, sympy.Function]):
-            inner_self._yield_stress_min = value
-            inner_self._reset()
-
-        @property
-        def strainrate_inv_II(inner_self):
-            return inner_self._strainrate_inv_II
-
-        @strainrate_inv_II.setter
-        def strainrate_inv_II(inner_self, value: sympy.Function):
-            inner_self._strainrate_inv_II = value
-            inner_self._reset()
-
-        @property
-        def stress_star(inner_self):
-            return inner_self._stress_star
-
-        @stress_star.setter
-        def stress_star(inner_self, value: sympy.Function):
-            inner_self._stress_star = value
-            inner_self._reset()
-
-        @property
-        def stress_star_star(inner_self):
-            return inner_self._stress_star_star
-
-        @stress_star_star.setter
-        def stress_star_star(inner_self, value: sympy.Function):
-            inner_self._stress_star_star = value
-            inner_self._reset()
-
-        @property
-        def strainrate_inv_II_min(inner_self):
-            return inner_self._strainrate_inv_II_min
-
-        @strainrate_inv_II_min.setter
-        def strainrate_inv_II_min(inner_self, value: float):
-            inner_self._strainrate_inv_II_min = sympy.sympify(value)
-            inner_self._reset()
-
-        @property
-        def averaging_method(inner_self):
-            return inner_self._averaging_method
-
-        @averaging_method.setter
-        def averaging_method(inner_self, value: str):
-            inner_self._averaging_method = value
-            inner_self._reset()
-
-        ### Getter and setter for internel mask Variable
-        @property
-        def materialIndex(inner_self):
-            return inner_self._materialIndex
-
-        @materialIndex.setter
-        def materialIndex(inner_self, indexVar):
-            # error checking, only support IndexSwarmVariables for now
-            if isinstance(indexVar, uw.swarm.IndexSwarmVariable):
-                inner_self._materialIndex = indexVar
-                inner_self._reset()
-
-        @property
-        def t_relax(inner_self):
-            # shear modulus defaults to infinity so t_relax goes to zero
-            # in the viscous limit
-
-            if inner_self.materialIndex == None:
-                shear_viscosity_0 = inner_self.shear_viscosity_0[0]
-                shear_modulus = inner_self.shear_modulus[0]
-                t_relax = shear_viscosity_0 / shear_modulus
-            else:
-                shear_viscosity_0 = np.array(inner_self.shear_viscosity_0)
-                shear_modulus = np.array(inner_self.shear_modulus)
-                t_relax = shear_viscosity_0 / shear_modulus
-
-            return t_relax
-
-        @property
-        def ve_effective_viscosity(inner_self):
-            # the dt_elastic defaults to infinity, t_relax to zero,
-            # so this should be well behaved in the viscous limit
-
-            import warnings
-
-            if (
-                type(inner_self.shear_modulus) != np.ndarray
-                and type(inner_self.shear_modulus) != list
-            ):
-                inner_self.shear_modulus = list([inner_self.shear_modulus])
-            if type(inner_self.dt_elastic) != np.ndarray and type(inner_self.dt_elastic) != list:
-                inner_self.dt_elastic = list([inner_self.dt_elastic])
-            if (
-                type(inner_self.shear_viscosity_0) != np.ndarray
-                and type(inner_self.shear_viscosity_0) != list
-            ):
-                inner_self.shear_viscosity_0 = list([inner_self.shear_viscosity_0])
-
-            if inner_self.materialIndex == None:
-                warnings.warn(
-                    "materialIndex not specified, using the first value for each parameter",
-                    stacklevel=2,
-                )
-
-                mu = inner_self.shear_modulus[0]
-                dt = inner_self.dt_elastic[0]
-                eta = inner_self.shear_viscosity_0[0]
-
-                if mu == sympy.oo or dt == sympy.oo:
-                    return eta
-
-                ## The effective viscosity depends on the number of history terms
-                if inner_self.stress_star_star is None:
-                    el_eff_visc = eta * mu * dt / (mu * dt + eta)
-                else:
-                    el_eff_visc = 2 * eta * mu * dt / (2 * mu * dt + 3 * eta)
-
-                return sympy.simplify(el_eff_visc)
-
-            else:
-                ### creates list of values that has the same length as the material index
-                if len(inner_self.shear_modulus) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_modulus) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_modulus ({len(inner_self.shear_modulus)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_modulus = list(
-                        np.repeat(
-                            inner_self.shear_modulus[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                if len(inner_self.dt_elastic) != inner_self.materialIndex.indices:
-                    if len(inner_self.dt_elastic) > 1:
-                        warnings.warn(
-                            f"Number of values in dt_elastic ({len(inner_self.dt_elastic)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.dt_elastic = list(
-                        np.repeat(inner_self.dt_elastic[0], inner_self.materialIndex.indices)
-                    )
-
-                if len(inner_self.shear_viscosity_0) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_0) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_0 ({len(inner_self.shear_viscosity_0)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_0 = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_0[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                mu = np.array(inner_self.shear_modulus)
-                dt = np.array(inner_self.dt_elastic)
-                eta = np.array(inner_self.shear_viscosity_0)
-
-                if mu[0] == sympy.oo or dt[0] == sympy.oo:
-                    return eta
-
-                ## The effective viscosity depends on the number of history terms
-                if inner_self.stress_star_star is None:
-                    el_eff_visc = eta * mu * dt / (mu * dt + eta)
-                else:
-                    el_eff_visc = 2 * eta * mu * dt / (2 * mu * dt + 3 * eta)
-
-                return sympy.simplify(el_eff_visc)
-
-        @property
-        def plastic_eff_viscosity(inner_self):
-            import warnings
-
-            if (
-                type(inner_self.yield_stress) != np.ndarray
-                and type(inner_self.yield_stress) != list
-            ):
-                inner_self.yield_stress = list([inner_self.yield_stress])
-            if (
-                type(inner_self.yield_stress_min) != np.ndarray
-                and type(inner_self.yield_stress_min) != list
-            ):
-                inner_self.yield_stress_min = list([inner_self.yield_stress_min])
-
-            if inner_self.materialIndex == None:
-                warnings.warn(
-                    "materialIndex not specified, using the first value for each parameter",
-                    stacklevel=2,
-                )
-
-                yield_stress = inner_self.yield_stress[0]
-                yield_stress_min = inner_self.yield_stress_min[0]
-
-                if yield_stress_min == -sympy.oo:
-                    yield_stress_fn = yield_stress
-                else:
-                    yield_stress_fn = sympy.Max(yield_stress, yield_stress_min)
-
-                if yield_stress_fn == sympy.oo or inner_self.strainrate_inv_II == sympy.oo:
-                    pl_effective_viscosity = sympy.oo
-
-                else:
-                    pl_effective_viscosity = yield_stress_fn / (
-                        (2 * inner_self.strainrate_inv_II)
-                        # + inner_self.strainrate_inv_II_min
-                    )
-
-            else:
-                ### creates list of values that has the same length as the material index
-                if len(inner_self.yield_stress) != inner_self.materialIndex.indices:
-                    if len(inner_self.yield_stress) > 1:
-                        warnings.warn(
-                            f"Number of values in yield_stress ({len(inner_self.yield_stress)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.yield_stress = list(
-                        np.repeat(inner_self.yield_stress[0], inner_self.materialIndex.indices)
-                    )
-
-                if len(inner_self.yield_stress_min) != inner_self.materialIndex.indices:
-                    if len(inner_self.yield_stress_min) > 1:
-                        warnings.warn(
-                            f"Number of values in yield_stress ({len(inner_self.yield_stress_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.yield_stress_min = list(
-                        np.repeat(
-                            inner_self.yield_stress_min[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                yield_stress = sympy.Matrix(inner_self.yield_stress)
-                yield_stress_min = sympy.Matrix(inner_self.yield_stress_min)
-
-                if yield_stress_min[0] == -sympy.oo:
-                    yield_stress_fn = yield_stress
-                else:
-                    yield_stress_list = []
-
-                    for i in range(len(yield_stress)):
-                        yield_stress_list.append(sympy.Max(yield_stress[i], yield_stress_min[i]))
-
-                    yield_stress_fn = sympy.Matrix(yield_stress_list)
-
-                if yield_stress_fn[0] == sympy.oo or inner_self.strainrate_inv_II == sympy.oo:
-                    pl_effective_viscosity = list(
-                        np.repeat(sympy.oo, inner_self.materialIndex.indices)
-                    )
-
-                else:
-                    pl_effective_viscosity = yield_stress_fn / (
-                        (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min
-                    )
-
-            return pl_effective_viscosity
-
-        # This has no setter !!
-        @property
-        def viscosity(inner_self):
-            # detect if values we need are defined or are placeholder symbols
-
-            ve_eff_visc = inner_self.ve_effective_viscosity
-            yield_visc = inner_self.plastic_eff_viscosity
-
-            import warnings
-
-            if inner_self.materialIndex == None:
-                shear_viscosity_min = inner_self.shear_viscosity_min[0]
-                shear_viscosity_max = inner_self.shear_viscosity_max[0]
-
-                if inner_self.averaging_method.casefold() == "min":
-                    if yield_visc != 0:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(shear_viscosity_min, sympy.Min(yield_visc, ve_eff_visc)),
-                        )
-                    else:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(shear_viscosity_min, ve_eff_visc),
-                        )
-                else:
-                    if yield_visc != 0:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(
-                                shear_viscosity_min,
-                                1.0 / ((1.0 / ve_eff_visc) + (1.0 / yield_visc)),
-                            ),
-                        )
-                    else:
-                        effective_viscosity = sympy.Min(
-                            shear_viscosity_max,
-                            sympy.Max(shear_viscosity_min, ve_eff_visc),
-                        )
-
-            else:
-                if len(inner_self.shear_viscosity_min) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_min) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_min ({len(inner_self.shear_viscosity_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_min = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_min[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                if len(inner_self.shear_viscosity_max) != inner_self.materialIndex.indices:
-                    if len(inner_self.shear_viscosity_max) > 1:
-                        warnings.warn(
-                            f"Number of values in shear_viscosity_max ({len(inner_self.shear_viscosity_max)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.",
-                            stacklevel=2,
-                        )
-                    inner_self.shear_viscosity_max = list(
-                        np.repeat(
-                            inner_self.shear_viscosity_max[0],
-                            inner_self.materialIndex.indices,
-                        )
-                    )
-
-                shear_viscosity_min = sympy.Matrix(inner_self.shear_viscosity_min)
-                shear_viscosity_max = sympy.Matrix(inner_self.shear_viscosity_max)
-
-                viscosity_list = []
-                for i in range(len(yield_visc)):
-                    if inner_self.averaging_method.casefold() == "min":
-                        if yield_visc[i] != 0:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(
-                                        shear_viscosity_min[i],
-                                        sympy.Min(yield_visc[i], ve_eff_visc[i]),
-                                    ),
-                                )
-                            )
-                        else:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(shear_viscosity_min[i], ve_eff_visc[i]),
-                                )
-                            )
-                    else:
-                        if yield_visc[i] != 0:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(
-                                        shear_viscosity_min[i],
-                                        1 / ((1 / ve_eff_visc[i]) + (1 / yield_visc[i])),
-                                    ),
-                                )
-                            )
-                        else:
-                            viscosity_list.append(
-                                sympy.Min(
-                                    shear_viscosity_max[i],
-                                    sympy.Max(shear_viscosity_min[i], ve_eff_visc[i]),
-                                )
-                            )
-
-                    viscosity_fn = sympy.Matrix(viscosity_list)
-
-                effective_viscosity = inner_self.materialIndex.sym.T.dot(viscosity_fn)
-
-            return effective_viscosity
-
-    def __init__(self, dim):
-        u_dim = dim
-        super().__init__(dim, u_dim)
-
-        return
-
-    ## Is this really different from the original ?
-
-    def _build_c_tensor(self):
-        """For this constitutive law, we expect just a viscosity function"""
-
-        d = self.dim
-        viscosity = self.Parameters.viscosity
-        shear_modulus = self.Parameters.shear_modulus
-        dt_elastic = self.Parameters.dt_elastic
-
-        try:
-            self._c = 2 * uw.maths.tensor.rank4_identity(d) * viscosity
-        except:
-            d = self.dim
-            dv = uw.maths.tensor.idxmap[d][0]
-            if isinstance(viscosity, sympy.Matrix) and viscosity.shape == (dv, dv):
-                self._c = 2 * uw.maths.tensor.mandel_to_rank4(viscosity, d)
-            elif isinstance(viscosity, sympy.Array) and viscosity.shape == (d, d, d, d):
-                self._c = 2 * viscosity
-            else:
-                raise RuntimeError(
-                    "Viscosity is not a known type (scalar, Mandel matrix, or rank 4 tensor"
-                )
-        return
-
-    # Modify flux to use the stress history term
-    # This may be preferable to using strain rate which can be discontinuous
-    # and harder to map back and forth between grid and particles without numerical smoothing
-
-    def flux(
-        self,
-        ddu: sympy.Matrix = None,
-        ddu_dt: sympy.Matrix = None,
-        u: sympy.Matrix = None,  # may be needed in the case of cylindrical / spherical
-        u_dt: sympy.Matrix = None,
-    ):
-        """Computes the effect of the constitutive tensor on the gradients of the unknowns.
-        (always uses the `c` form of the tensor). In general cases, the history of the gradients
-        may be required to evaluate the flux. For viscoelasticity, the
-        """
-
-        c = self.c
-        rank = len(c.shape)
-
-        # tensor multiplication
-
-        if rank == 2:
-            flux = c * ddu.T
-        else:  # rank==4
-            flux = sympy.tensorcontraction(
-                sympy.tensorcontraction(sympy.tensorproduct(c, ddu), (3, 5)), (0, 1)
-            )
-
-        # Now add in the stress history. In the
-        # viscous limit, this term is not well behaved
-        # and we need to check that
-
-        if self.is_elastic:
-            if self.Parameters.materialIndex == None:
-                eta = self.Parameters.shear_viscosity_0[0]
-                mu = self.Parameters.shear_modulus[0]
-                dt = self.Parameters.dt_elastic[0]
-                s_star = self.Parameters.stress_star
-                s_star_star = self.Parameters.stress_star_star
-
-                if s_star_star is None:  # 1st order
-                    flux = sympy.Matrix(flux) + eta * s_star / (dt * mu + eta)
-                else:  # 2nd order
-                    flux = (
-                        sympy.Matrix(flux)
-                        + 4 * eta * s_star / (2 * dt * mu + 3 * eta)
-                        - eta * s_star_star / (2 * dt * mu + 3 * eta)
-                    )
-            else:
-                eta = self.Parameters.materialIndex.createMask(self.Parameters.shear_viscosity_0)
-                mu = self.Parameters.materialIndex.createMask(self.Parameters.shear_modulus)
-                dt = np.array(self.Parameters.dt_elastic)
-                s_star = self.Parameters.stress_star
-                s_star_star = self.Parameters.stress_star_star
-
-                if s_star_star is None:  # 1st order
-                    flux = sympy.Matrix(flux) + eta * s_star / (dt * mu + eta)
-                else:  # 2nd order
-                    flux = (
-                        sympy.Matrix(flux)
-                        + 4 * eta * s_star / (2 * dt * mu + 3 * eta)
-                        - eta * s_star_star / (2 * dt * mu + 3 * eta)
-                    )
-
-        return sympy.simplify(sympy.Matrix(flux))
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-
-        super()._object_viewer()
-
-        display(Markdown(r"### Viscous deformation"))
-        display(
-            Latex(
-                r"$\quad\eta_\textrm{0} = $ "
-                + sympy.sympify(self.Parameters.shear_viscosity_0[0])._repr_latex_()
-            ),
-        )
-
-        ## If elasticity is active:
-        display(Markdown(r"#### Elastic deformation"))
-        display(
-            Latex(
-                r"$\quad\mu = $ " + sympy.sympify(self.Parameters.shear_modulus[0])._repr_latex_(),
-            ),
-            Latex(
-                r"$\quad\Delta t_e = $ "
-                + sympy.sympify(self.Parameters.dt_elastic[0])._repr_latex_(),
-            ),
-            Latex(
-                r"$\quad \sigma^* = $ " + sympy.sympify(self.Parameters.stress_star)._repr_latex_(),
-            ),
-        )
-        if self.Parameters.stress_star_star is not None:
-            display(
-                Latex(
-                    r"$\quad \sigma^{**} = $ "
-                    + sympy.sympify(self.Parameters.stress_star_star)._repr_latex_(),
-                ),
-            )
-
-        # If plasticity is active
-        display(Markdown(r"#### Plastic deformation"))
-        display(
-            Latex(
-                r"$\quad\tau_\textrm{y} = $ "
-                + sympy.sympify(self.Parameters.yield_stress[0])._repr_latex_(),
-            ),
-            Latex(
-                r"$\quad|\dot\epsilon| = $ "
-                + sympy.sympify(self.Parameters.strainrate_inv_II)._repr_latex_(),
-            ),
-            ## Todo: add all the other properties in here
-        )
-
-    @property
-    def is_elastic(self):
-        # If any of these is not defined, elasticity is switched off
-
-        if self.Parameters.dt_elastic[0] is sympy.oo:
-            return False
-
-        if self.Parameters.shear_modulus[0] is sympy.oo:
-            return False
-
-        if self.Parameters.stress_star is None:
-            return False
-
-        return True
-
-    @property
-    def is_viscoplastic(self):
-        if self.Parameters.yield_stress[0] == sympy.oo:
-            return False
-
-        if isinstance(self.Parameters.strainrate_inv_II, sympy.core.symbol.Symbol):
-            return False
-
-        return True
-
-
-###
-
-
-class DiffusionModel(Constitutive_Model):
-    r"""
-    Diffusion constitutive model for scalar transport.
-
-    .. math::
-
-        q_{i} = \kappa_{ij} \cdot \frac{\partial \phi}{\partial x_j}
-
-    where :math:`\kappa` is a diffusivity—a scalar constant, sympy function,
-    or mesh variable.
-
-    Examples
-    --------
-    >>> diffusion_model = DiffusionModel(dim)
-    >>> diffusion_model.material_properties = diffusion_model.Parameters(
-    ...     diffusivity=diffusivity_fn
-    ... )
-    >>> scalar_solver.constitutive_model = diffusion_model
-    >>> flux = diffusion_model.flux(gradient_matrix)
-    """
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-        """
-
-        def __init__(
-            inner_self,
-            diffusivity: Union[float, sympy.Function] = 1,
-        ):
-            inner_self._diffusivity = diffusivity
-
-        @property
-        def diffusivity(inner_self):
-            return inner_self._diffusivity
-
-        @diffusivity.setter
-        def diffusivity(inner_self, value: Union[float, sympy.Function]):
-            inner_self._diffusivity = value
-            inner_self._reset()
-
-    def __init__(self, dim):
-        self.u_dim = 1
-        super().__init__(dim, self.u_dim)
-
-        return
-
-    def _build_c_tensor(self):
-        """For this constitutive law, we expect just a diffusivity function"""
-
-        d = self.dim
-        kappa = self.Parameters.diffusivity
-        self._c = sympy.Matrix.eye(d) * kappa
-
-        return
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-
-        super()._object_viewer()
-
-        ## feedback on this instance
-        display(
-            Latex(r"$\quad\kappa = $ " + sympy.sympify(self.Parameters.diffusivity)._repr_latex_())
-        )
-
-        return
-
-
-class TransverseIsotropicFlowModel(Constitutive_Model):
-    r"""
-    Transversely isotropic viscous flow model.
-
-    .. math::
-
-        \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l}
-        + \frac{\partial u_l}{\partial x_k} \right]
-
-    where the viscosity tensor :math:`\eta` is defined as:
-
-    .. math::
-
-        \eta_{ijkl} = \eta_0 I_{ijkl} + (\eta_0 - \eta_1) \left[
-        \frac{1}{2}\left( n_i n_l \delta_{jk} + n_j n_k \delta_{il}
-        + n_i n_l \delta_{jk} + n_j n_l \delta_{ik} \right)
-        - 2 n_i n_j n_k n_l \right]
-
-    and :math:`\hat{\mathbf{n}} \equiv \{n_i\}` is the unit vector defining the
-    local orientation of the weak plane (the director).
-
-    The Mandel matrix is in ``viscous_model.C`` and the rank-4 tensor in
-    ``viscous_model.c``.
-
-    Examples
-    --------
-    >>> viscous_model = TransverseIsotropicFlowModel(dim)
-    >>> viscous_model.material_properties = viscous_model.Parameters(
-    ...     eta_0=viscosity_fn,
-    ...     eta_1=weak_viscosity_fn,
-    ...     director=orientation_vector_fn
-    ... )
-    >>> solver.constitutive_model = viscous_model
-    >>> tau = viscous_model.flux(gradient_matrix)
-    """
-
-    class _Parameters:
-        """Any material properties that are defined by a constitutive relationship are
-        collected in the parameters which can then be defined/accessed by name in
-        individual instances of the class.
-        """
-
-        def __init__(
-            inner_self,
-            eta_0: Union[float, sympy.Function] = 1,
-            eta_1: Union[float, sympy.Function] = 1,
-            director: Union[sympy.Matrix, sympy.Function] = sympy.Matrix([0, 0, 1]),
-        ):
-            inner_self._eta_0 = eta_0
-            inner_self._eta_1 = eta_1
-            inner_self._director = director
-            # inner_self.constitutive_model_class = const_model
-
-        ## Note the inefficiency below if we change all these values one after the other
-
-        @property
-        def eta_0(inner_self):
-            return inner_self._eta_0
-
-        @eta_0.setter
-        def eta_0(
-            inner_self,
-            value: Union[float, sympy.Function],
-        ):
-            inner_self._eta_0 = value
-            inner_self._reset()
-
-        @property
-        def eta_1(inner_self):
-            return inner_self._eta_1
-
-        @eta_1.setter
-        def eta_1(
-            inner_self,
-            value: Union[float, sympy.Function],
-        ):
-            inner_self._eta_1 = value
-            inner_self._reset()
-
-        @property
-        def director(inner_self):
-            return inner_self._director
-
-        @director.setter
-        def director(
-            inner_self,
-            value: Union[sympy.Matrix, sympy.Function],
-        ):
-            inner_self._director = value
-            inner_self._reset()
-
-    def __init__(self, dim):
-        u_dim = dim
-        super().__init__(dim, u_dim)
-
-        # default values ... maybe ??
-        return
-
-    def _build_c_tensor(self):
-        """For this constitutive law, we expect two viscosity functions
-        and a sympy matrix that describes the director components n_{i}"""
-
-        d = self.dim
-        dv = uw.maths.tensor.idxmap[d][0]
-
-        eta_0 = self.Parameters.eta_0
-        eta_1 = self.Parameters.eta_1
-        n = self.Parameters.director
-
-        Delta = eta_1 - eta_0
-
-        lambda_mat = uw.maths.tensor.rank4_identity(d) * eta_0
-
-        for i in range(d):
-            for j in range(d):
-                for k in range(d):
-                    for l in range(d):
-                        lambda_mat[i, j, k, l] += Delta * (
-                            (
-                                n[i] * n[k] * int(j == l)
-                                + n[j] * n[k] * int(l == i)
-                                + n[i] * n[l] * int(j == k)
-                                + n[j] * n[l] * int(k == i)
-                            )
-                            / 2
-                            - 2 * n[i] * n[j] * n[k] * n[l]
-                        )
-
-        lambda_mat = sympy.simplify(uw.maths.tensor.rank4_to_mandel(lambda_mat, d))
-
-        self._c = uw.maths.tensor.mandel_to_rank4(lambda_mat, d)
-
-    def _object_viewer(self):
-        from IPython.display import Latex, Markdown, display
-
-        super()._object_viewer()
-
-        ## feedback on this instance
-        display(Latex(r"$\quad\eta_0 = $ " + sympy.sympify(self.Parameters.eta_0)._repr_latex_()))
-        display(Latex(r"$\quad\eta_1 = $ " + sympy.sympify(self.Parameters.eta_1)._repr_latex_()))
-        display(
-            Latex(
-                r"$\quad\hat{\mathbf{n}} = $ "
-                + sympy.sympify(self.Parameters.director.T)._repr_latex_()
-            )
-        )
-
-
-class MultiMaterial(Constitutive_Model):
-    r"""
-    Manage multiple materials in a constitutive framework.
-
-    Bundles multiple materials into a single consitutive law. The expectation
-    is that these all have compatible flux terms.
-    """
-
-    def __init__(
-        self,
-        material_swarmVariable: Optional[IndexSwarmVariable] = None,
-        constitutive_models: Optional[list] = [],
-    ):
-        self._constitutive_models = constitutive_models
-        self._material_var = material_swarmVariable
-
-        return
diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index d502dd929..f2e5cc103 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -129,6 +129,21 @@ PetscErrorCode UW_PetscDSViewBdWF(PetscDS ds, PetscInt bd)
     return 1;
 }
 
+// Set the time value on a DM. This is passed as `petsc_t` to all
+// pointwise residual and Jacobian functions during assembly.
+// PETSc stores this internally but petsc4py doesn't expose it.
+PetscErrorCode UW_DMSetTime(DM dm, PetscReal time)
+{
+    // DMSetOutputSequenceNumber stores (step, time) on the DM.
+    // The time component is what DMPlexComputeResidual_Internal
+    // passes as petsc_t to the pointwise functions.
+    PetscInt step;
+    PetscReal old_time;
+    PetscCall(DMGetOutputSequenceNumber(dm, &step, &old_time));
+    PetscCall(DMSetOutputSequenceNumber(dm, step, time));
+    return PETSC_SUCCESS;
+}
+
 PetscErrorCode UW_DMPlexSetSNESLocalFEM(DM dm, PetscBool flag, void *ctx)
 {
 
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index c39d9c715..edd73ca15 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -42,6 +42,7 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode UW_PetscDSSetBdTerms   (PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*, void*, void*, void*, void* )
     PetscErrorCode UW_PetscDSViewWF(PetscDS)     
     PetscErrorCode UW_PetscDSViewBdWF(PetscDS, PetscInt)
+    PetscErrorCode UW_DMSetTime( PetscDM, PetscReal )
     PetscErrorCode UW_DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
     PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*)
 
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index e830fdc42..c4b3377e0 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -956,6 +956,122 @@ class SolverBaseClass(uw_object):
 
 
 
+    @property
+    def tau(self):
+        r"""Computed flux from the constitutive model, projected onto the mesh.
+
+        Returns a :class:`~underworld3.discretisation.MeshVariable` containing
+        the flux that the solver actually computed. The variable is created
+        lazily on first access and updated by solving a projection each time
+        the property is read.
+
+        For Stokes-family solvers, this is the deviatoric stress tensor
+        :math:`\boldsymbol{\tau}`. For Poisson/diffusion, it is the
+        diffusive flux :math:`\kappa \nabla u`. The projection handles
+        derivative evaluation internally — users should use this property
+        rather than manually reconstructing the flux from the constitutive
+        formula.
+
+        Subclasses may override this to return pre-computed values (e.g.
+        VE_Stokes returns the stress stored during the solve).
+
+        Returns
+        -------
+        MeshVariable
+            Mesh variable containing the projected flux values.
+            Access numerical data via ``.data`` or ``.array``, symbolic
+            expression via ``.sym``.
+        """
+
+        if self._constitutive_model is None:
+            raise RuntimeError(
+                "No constitutive model set. Assign solver.constitutive_model first."
+            )
+
+        # Lazy initialization of projection infrastructure
+        if not hasattr(self, '_tau_var') or self._tau_var is None:
+            self._setup_tau_projection()
+
+        # Project the constitutive flux onto the mesh variable
+        # Transpose if needed: flux may be (dim,1) column but projector expects row
+        flux = self._constitutive_model.flux
+        if hasattr(flux, 'shape') and flux.shape[1] == 1 and flux.shape[0] > 1:
+            flux = flux.T
+        self._tau_projector.uw_function = flux
+        self._tau_projector.smoothing = 0.0
+        self._tau_projector.solve()
+
+        return self._tau_var
+
+    def _setup_tau_projection(self):
+        """Create the mesh variable and projector for tau (lazy init)."""
+
+        flux = self._constitutive_model.flux
+        dim = self.mesh.dim
+        rows, cols = flux.shape
+
+        # Determine variable type and create appropriate projection solver
+        if rows == cols and rows == dim:
+            # Tensor flux (Stokes stress)
+            self._tau_var = uw.discretisation.MeshVariable(
+                f"tau_{self.instance_number}",
+                self.mesh,
+                (dim, dim),
+                vtype=uw.VarType.SYM_TENSOR,
+                degree=self.u.degree,
+                continuous=True,
+                varsymbol=r"{\tau}",
+            )
+            _work = uw.discretisation.MeshVariable(
+                f"tau_work_{self.instance_number}",
+                self.mesh,
+                1,
+                degree=self.u.degree,
+                continuous=True,
+            )
+            from underworld3.systems.solvers import SNES_Tensor_Projection
+            self._tau_projector = SNES_Tensor_Projection(
+                self.mesh, self._tau_var, _work, verbose=False
+            )
+
+        elif rows == dim and cols == 1:
+            # Vector flux (Poisson heat flux)
+            self._tau_var = uw.discretisation.MeshVariable(
+                f"tau_{self.instance_number}",
+                self.mesh,
+                dim,
+                vtype=uw.VarType.VECTOR,
+                degree=self.u.degree,
+                continuous=True,
+                varsymbol=r"{\mathbf{q}}",
+            )
+            from underworld3.systems.solvers import SNES_Vector_Projection
+            self._tau_projector = SNES_Vector_Projection(
+                self.mesh, self._tau_var, verbose=False
+            )
+
+        elif cols == dim and rows == 1:
+            # Transposed vector flux
+            self._tau_var = uw.discretisation.MeshVariable(
+                f"tau_{self.instance_number}",
+                self.mesh,
+                dim,
+                vtype=uw.VarType.VECTOR,
+                degree=self.u.degree,
+                continuous=True,
+                varsymbol=r"{\mathbf{q}}",
+            )
+            from underworld3.systems.solvers import SNES_Vector_Projection
+            self._tau_projector = SNES_Vector_Projection(
+                self.mesh, self._tau_var, verbose=False
+            )
+
+        else:
+            raise RuntimeError(
+                f"Cannot create tau projection for flux shape {flux.shape}. "
+                f"Expected ({dim},{dim}) tensor, ({dim},1) or (1,{dim}) vector."
+            )
+
     def validate_solver(self):
         """
         Checks to see if the required properties have been set.
@@ -1668,7 +1784,8 @@ class SNES_Scalar(SolverBaseClass):
               _force_setup:    bool =False,
               verbose:         bool=False,
               debug:           bool=False,
-              debug_name:      str=None ):
+              debug_name:      str=None,
+              time=None, ):
         """
         Solve the system of equations.
 
@@ -1691,6 +1808,11 @@ class SNES_Scalar(SolverBaseClass):
             Enable debug output including intermediate residuals.
         debug_name : str, optional
             Name prefix for debug output files.
+        time : float or Quantity, optional
+            Physical time for this solve. Passed as ``petsc_t`` to all
+            pointwise functions. Expressions using ``mesh.t`` evaluate at
+            this time. Non-dimensionalised when scaling is active.
+            Default: None (petsc_t unchanged).
 
         Returns
         -------
@@ -1728,6 +1850,16 @@ class SNES_Scalar(SolverBaseClass):
 
         self._build(verbose, debug, debug_name)
 
+        # Set time on the DM so petsc_t is available in pointwise functions
+        cdef DM _time_dm
+        if time is not None:
+            if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'):
+                t_nd = float(uw.non_dimensionalise(time))
+            else:
+                t_nd = float(time)
+            _time_dm = self.dm
+            UW_DMSetTime(_time_dm.dm, t_nd)
+
         gvec = self.dm.getGlobalVec()
 
         if not zero_init_guess:
@@ -4037,7 +4169,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
               verbose=False,
               debug=False,
               debug_name=None,
-              _force_setup: bool =False, ):
+              _force_setup: bool =False,
+              time=None, ):
         """
         Solve the Stokes system for velocity and pressure.
 
@@ -4064,6 +4197,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             Name prefix for debug output files.
         _force_setup : bool, default=False
             Force rebuild of the solver even if already set up.
+        time : float or Quantity, optional
+            Physical time for this solve. Passed as ``petsc_t`` to all
+            pointwise residual and Jacobian functions. Expressions using
+            ``mesh.t`` will evaluate at this time. Non-dimensionalised
+            automatically when scaling is active. Default: None (petsc_t=0).
 
         Returns
         -------
@@ -4105,6 +4243,17 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
 
         self._build(verbose, debug, debug_name)
 
+        # Set time on the DM so petsc_t is available in pointwise functions.
+        # Non-dimensionalise if the scaling system is active.
+        cdef DM _time_dm_stokes
+        if time is not None:
+            if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'):
+                t_nd = float(uw.non_dimensionalise(time))
+            else:
+                t_nd = float(time)
+            _time_dm_stokes = self.dm
+            UW_DMSetTime(_time_dm_stokes.dm, t_nd)
+
         # Keep a record of these set-up parameters
         tolerance = self.tolerance
         snes_type = self.snes.getType()
diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py
index c111d0fdd..4b11dbd6e 100644
--- a/src/underworld3/discretisation/discretisation_mesh.py
+++ b/src/underworld3/discretisation/discretisation_mesh.py
@@ -629,6 +629,14 @@ def mesh_update_callback(array, change_context):
         self._Gamma.y._ccodestr = "petsc_n[1]"
         self._Gamma.z._ccodestr = "petsc_n[2]"
 
+        # Time coordinate — PETSc passes this as petsc_t to all pointwise
+        # functions. Solvers set dm.time before each solve via solve(time=t).
+        # Users reference it as mesh.t in expressions (e.g. V0 * sympy.sin(omega * mesh.t))
+        from ..utilities.unit_aware_coordinates import TimeSymbol
+
+        self._t = TimeSymbol("t")
+        self._t._units = None  # patched below by _patch_time_units
+
         # Add unit awareness to coordinate symbols if mesh has units or model has scales
         from ..utilities.unit_aware_coordinates import patch_coordinate_units
 
@@ -1497,6 +1505,32 @@ def CoordinateSystem(self) -> CoordinateSystem:
         r"""Alias for :attr:`X` (the coordinate system object)."""
         return self._CoordinateSystem
 
+    @property
+    def t(self):
+        r"""Symbolic time coordinate.
+
+        PETSc passes a time value (``petsc_t``) to all pointwise residual
+        and Jacobian functions. Use ``mesh.t`` in expressions to reference
+        this time without forcing JIT recompilation each timestep.
+
+        The low-level PETSc solver accepts ``time=t`` to set the value
+        of ``petsc_t`` for pointwise functions. If not provided, ``petsc_t``
+        defaults to 0. Note: the high-level Python ``solve()`` wrappers
+        do not yet pass ``time=`` through — set it directly via
+        ``UW_DMSetTime`` at the Cython level if needed.
+
+        When the scaling system is active, ``mesh.t`` carries time units
+        (derived from the model's time scale) so that dimensional analysis
+        works correctly in expressions.
+
+        Examples
+        --------
+        >>> omega = 2 * np.pi / period
+        >>> stokes.add_dirichlet_bc((V0 * sympy.sin(omega * mesh.t), 0.0), "Top")
+        >>> stokes.solve(time=current_time)   # sets petsc_t before SNES
+        """
+        return self._t
+
     @property
     def r(self) -> Tuple[sympy.vector.BaseScalar]:
         r"""Tuple of coordinate scalars :math:`(x, y)` or :math:`(x, y, z)`.
diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx
index d9a0b825a..98f218e3b 100644
--- a/src/underworld3/function/_function.pyx
+++ b/src/underworld3/function/_function.pyx
@@ -548,6 +548,12 @@ def _project_to_work_variable(expr, mesh, smoothing=1e-6):
     import underworld3 as uw
 
     # Handle matrix expressions - need multi-component work variable
+    # TODO(BUG): This fails for dim×dim matrices (e.g. 2×2 stress tensor)
+    # because MeshVariable can't infer vtype from num_components=4 (flat int).
+    # Needs: pass num_components=(rows,cols) with vtype=TENSOR, and use
+    # Tensor_Projection instead of per-component scalar Projection.
+    # Currently, callers like SemiLagrangian.update_pre_solve fall back to
+    # their own projection solver via except-clause when this fails.
     if hasattr(expr, 'shape') and expr.shape != (1, 1):
         rows, cols = expr.shape
         n_components = rows * cols
diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py
index 25e2e3e0f..64c10b255 100644
--- a/src/underworld3/systems/ddt.py
+++ b/src/underworld3/systems/ddt.py
@@ -445,7 +445,9 @@ def effective_order(self):
         startup, ``effective_order`` ramps from 1 to ``self.order`` as
         successive solves populate the history slots with distinct values.
         """
-        return min(self.order, max(1, self._n_solves_completed + 1))
+        # BDF-k requires k completed solves to have k distinct history values.
+        # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves.
+        return min(self.order, max(1, self._n_solves_completed))
 
     def update_history_fn(self):
         r"""Copy current :math:`\psi` to the first history slot ``psi_star[0]``."""
@@ -524,6 +526,11 @@ def update_post_solve(
 
         return
 
+    @property
+    def bdf_coefficients(self):
+        """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps."""
+        return _bdf_coefficients(self.effective_order, self._dt, self._dt_history)
+
     def bdf(self, order: Optional[int] = None):
         r"""Backward differentiation approximation of the time-derivative of ψ.
 
@@ -763,7 +770,9 @@ def effective_order(self):
         startup, ``effective_order`` ramps from 1 to ``self.order`` as
         successive solves populate the history slots with distinct values.
         """
-        return min(self.order, max(1, self._n_solves_completed + 1))
+        # BDF-k requires k completed solves to have k distinct history values.
+        # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves.
+        return min(self.order, max(1, self._n_solves_completed))
 
     def update_history_fn(self):
         r"""Copy current :math:`\psi` to ``psi_star[0]`` via evaluation or projection."""
@@ -890,6 +899,11 @@ def update_post_solve(
 
         return
 
+    @property
+    def bdf_coefficients(self):
+        """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps."""
+        return _bdf_coefficients(self.effective_order, self._dt, self._dt_history)
+
     def bdf(self, order=None):
         r"""Backward differentiation approximation of the time-derivative of :math:`\psi`.
 
@@ -1158,7 +1172,9 @@ def effective_order(self):
         startup, ``effective_order`` ramps from 1 to ``self.order`` as
         successive solves populate the history slots with distinct values.
         """
-        return min(self.order, max(1, self._n_solves_completed + 1))
+        # BDF-k requires k completed solves to have k distinct history values.
+        # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves.
+        return min(self.order, max(1, self._n_solves_completed))
 
     def initialise_history(self):
         r"""Initialize all history slots to the current value of :math:`\psi`.
@@ -1229,6 +1245,7 @@ def update_post_solve(
 
         if self._n_solves_completed < self.order:
             self._n_solves_completed += 1
+
         return
 
     def update_pre_solve(
@@ -1237,11 +1254,22 @@ def update_pre_solve(
         evalf: Optional[bool] = False,
         verbose: Optional[bool] = False,
         dt_physical: Optional[float] = None,
+        store_result: Optional[bool] = True,
     ):
         """Sample upstream values along characteristics before solve.
 
         On the first call, automatically initialises history from the
         current field values so that bdf() returns zero on the first step.
+
+        Parameters
+        ----------
+        store_result : bool, optional
+            If True (default), evaluate psi_fn at current positions and store
+            in psi_star[0] before advection — the standard DDt behaviour.
+            If False, skip this step and the history shift: only advect the
+            existing psi_star levels upstream. Used by VE_Stokes where
+            psi_star[0] already contains the projected actual stress from
+            the previous solve.
         """
 
         self._dt = dt
@@ -1275,15 +1303,19 @@ def update_pre_solve(
         else:
             phi = sympy.sympify(1)
 
-        for i in range(self.order - 1, 0, -1):
-            self.psi_star[i].array[...] = (
-                phi * self.psi_star[i - 1].array[...] + (1 - phi) * self.psi_star[i].array[...]
-            )
+        if store_result:
+            for i in range(self.order - 1, 0, -1):
+                self.psi_star[i].array[...] = (
+                    phi * self.psi_star[i - 1].array[...] + (1 - phi) * self.psi_star[i].array[...]
+                )
 
         # 2. Compute the current value of psi_fn which we store in psi_star[0]
         #    Note the need to do a try/except to handle unsupported evaluations
         #    (e.g. of derivatives)
         #
+        #    When store_result=False (e.g. VE stress history), skip this step —
+        #    psi_star[0] already contains the projected actual stress from
+        #    the previous solve and we want to advect *that*, not the flux.
 
         # CRITICAL FIX (2025-11-28): Handle coordinates correctly for unit-aware mode.
         # Previous bug: extracting .magnitude gives METERS (e.g., 1000000), but:
@@ -1322,28 +1354,29 @@ def update_pre_solve(
             :, :
         ]
 
-        try:
-            # Use shifted ND coords to avoid quad mesh boundary issues
-            # node_coords_nd is slightly shifted toward cell centroids (lines 703-709)
-            # evaluate() treats plain numpy as ND [0-1] coordinates
-            eval_result = uw.function.evaluate(
-                self.psi_fn,
-                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)
+        if store_result:
+            try:
+                # Use shifted ND coords to avoid quad mesh boundary issues
+                # node_coords_nd is slightly shifted toward cell centroids
+                # evaluate() treats plain numpy as ND [0-1] coordinates
+                eval_result = uw.function.evaluate(
+                    self.psi_fn,
+                    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
+                self.psi_star[0].array[...] = eval_result
 
-        except Exception:
-            # Fallback to projection solver for expressions that can't be directly evaluated
-            # (e.g., containing derivatives)
-            self._psi_star_projection_solver.uw_function = self.psi_fn
-            self._psi_star_projection_solver.smoothing = 0.0
-            self._psi_star_projection_solver.solve(verbose=verbose)
+            except Exception:
+                # Fallback to projection solver for expressions that can't be directly evaluated
+                # (e.g., containing derivatives)
+                self._psi_star_projection_solver.uw_function = self.psi_fn
+                self._psi_star_projection_solver.smoothing = 0.0
+                self._psi_star_projection_solver.solve(verbose=verbose)
 
         # 3. Compute the upstream values from the psi_fn
 
@@ -1595,6 +1628,11 @@ def update_pre_solve(
 
         return
 
+    @property
+    def bdf_coefficients(self):
+        """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps."""
+        return _bdf_coefficients(self.effective_order, self._dt, self._dt_history)
+
     def bdf(self, order=None):
         r"""Backward differentiation approximation of the time-derivative of :math:`\psi`.
 
@@ -1773,7 +1811,9 @@ def _object_viewer(self):
     @property
     def effective_order(self):
         """Current effective BDF order, accounting for history startup."""
-        return min(self.order, max(1, self._n_solves_completed + 1))
+        # BDF-k requires k completed solves to have k distinct history values.
+        # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves.
+        return min(self.order, max(1, self._n_solves_completed))
 
     def initialise_history(self):
         r"""Initialize all history slots to the current value of :math:`\psi`.
@@ -1886,6 +1926,11 @@ def update_post_solve(
         if self._n_solves_completed < self.order:
             self._n_solves_completed += 1
 
+    @property
+    def bdf_coefficients(self):
+        """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps."""
+        return _bdf_coefficients(self.effective_order, self._dt, self._dt_history)
+
     def bdf(self, order=None):
         r"""Backward differentiation approximation of the time-derivative of :math:`\psi`.
 
@@ -2067,7 +2112,9 @@ def _object_viewer(self):
     @property
     def effective_order(self):
         """Current effective BDF order, accounting for history startup."""
-        return min(self.order, max(1, self._n_solves_completed + 1))
+        # BDF-k requires k completed solves to have k distinct history values.
+        # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves.
+        return min(self.order, max(1, self._n_solves_completed))
 
     def initialise_history(self):
         r"""Initialize all history slots to the current value of :math:`\psi`.
@@ -2181,6 +2228,11 @@ def update_post_solve(
 
         return
 
+    @property
+    def bdf_coefficients(self):
+        """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps."""
+        return _bdf_coefficients(self.effective_order, self._dt, self._dt_history)
+
     def bdf(self, order=None):
         r"""Backward differentiation approximation of the time-derivative of :math:`\psi`.
 
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 550632306..62387707a 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -1246,6 +1246,22 @@ def delta_t(self):
         """Elastic timestep from the constitutive model."""
         return self.constitutive_model.Parameters.dt_elastic
 
+    @property
+    def tau(self):
+        r"""Deviatoric stress from the most recent solve (stored in history).
+
+        For VE_Stokes, the stress is projected into ``psi_star[0]`` after
+        each solve. This override returns that variable directly — no
+        additional projection needed.
+
+        Returns
+        -------
+        MeshVariable
+            The stress history variable containing the actual deviatoric
+            stress from the most recent solve.
+        """
+        return self.DFDt.psi_star[0]
+
     ## Solver needs to update the stress history terms as well as call the SNES solve:
 
     @timing.routine_timer_decorator
@@ -1274,12 +1290,24 @@ def solve(
         if timestep is None:
             timestep = self.delta_t.sym
 
-        if timestep != self.delta_t:
-            self._constitutive_model.Parameters.dt_elastic = timestep  # this will force an initialisation because the functions need to be updated
+        # dt_elastic is a constitutive parameter (relaxation timescale) —
+        # never overwritten by solve(). The advection timestep (for departure
+        # point tracing) and the elastic relaxation timescale are independent.
 
         if _force_setup:
             self.is_setup = False
 
+        # Re-setup when effective_order changes (e.g. DDt history ramp-up
+        # from order 1 to order 2). The JIT-compiled pointwise functions
+        # depend on the order used in the constitutive model.
+        _current_eff_order = self.constitutive_model.effective_order
+        if not hasattr(self, '_prev_effective_order'):
+            self._prev_effective_order = None
+        if _current_eff_order != self._prev_effective_order:
+            self.is_setup = False
+            self.constitutive_model._solver_is_setup = False
+        self._prev_effective_order = _current_eff_order
+
         if not self.constitutive_model._solver_is_setup:
             self.is_setup = False
             self.DFDt.psi_fn = self.constitutive_model.flux.T
@@ -1289,11 +1317,29 @@ def solve(
             self._setup_discretisation(verbose)
             self._setup_solver(verbose)
 
+        # --- Stress history management via standard DDt pathway ---
+        #
+        # update_pre_solve(advect_only=True) performs:
+        #   1. History shift: psi_star[i] ← psi_star[i-1]
+        #   2. Skip psi_fn evaluation (psi_star[0] already has projected stress)
+        #   3. Advect all history levels to upstream positions along characteristics
+        #
+        # After this: psi_star[0] = σ* (previous stress at upstream),
+        #             psi_star[1] = σ** (stress from 2 steps ago, double-traced).
+
         if uw.mpi.rank == 0 and verbose:
-            print(f"VE Stokes solver - pre-solve DFDt update", flush=True)
+            print(f"VE Stokes solver - advect stress history", flush=True)
 
-        # Update SemiLagrange Flux terms
-        self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf)
+        self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf,
+                                   store_result=False)
+
+        # Update BDF coefficients from current dt_elastic and DDt history.
+        # These are UWexpressions that route through PetscDS constants[],
+        # so the compiled pointwise functions pick up the new values without
+        # JIT recompilation.
+        self.constitutive_model._update_bdf_coefficients()
+
+        # 2. SOLVE: PETSc uses the advected σ*, σ** via the constitutive model
 
         if uw.mpi.rank == 0 and verbose:
             print(f"VE Stokes solver - solve Stokes flow", flush=True)
@@ -1305,8 +1351,47 @@ def solve(
             picard=0,
         )
 
+        # 3. STORE ACTUAL STRESS and SHIFT HISTORY.
+        #
+        #    After advection + solve:
+        #      psi_star[0] = advected σ*  (used by the solver)
+        #      psi_star[1] = advected σ** (used by the solver)
+        #
+        #    We need to:
+        #    a) Project actual stress → psi_star[0] (while σ*, σ** are intact)
+        #    b) Save the advected σ* into psi_star[1] for next step's σ**
+        #       (chained characteristic tracing)
+        #
+        #    The projection reads psi_star[0..1] via stress_deviator, so we
+        #    must project BEFORE shifting. Then save the advected σ* and
+        #    overwrite psi_star[0] with the projected stress.
+
         if uw.mpi.rank == 0 and verbose:
-            print(f"VEP Stokes solver - post-solve DFDt update", flush=True)
+            print(f"VE Stokes solver - store stress and shift history", flush=True)
+
+        import numpy as np
+
+        # Save advected σ* before projection modifies psi_star[0]
+        _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...])
+
+        # Project actual stress into psi_star[0].
+        # Uses the constitutive formula so that σ* and σ** from psi_star
+        # are read correctly during projection.
+        self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux
+        self.DFDt._psi_star_projection_solver.smoothing = 0.0
+        self.DFDt._psi_star_projection_solver.solve(verbose=verbose)
+
+        # Now psi_star[0] = projected τ (actual stress from this solve).
+        # Shift history: psi_star[i] ← previous psi_star[i-1] (advected values).
+        # psi_star[1] gets the advected σ* (saved before projection).
+        # Higher levels shift down the chain.
+        for i in range(self.DFDt.order - 1, 0, -1):
+            if i == 1:
+                self.DFDt.psi_star[i].array[...] = _advected_sigma_star
+            else:
+                self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...]
+
+        # 5. BOOKKEEPING
 
         self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf)
 
diff --git a/src/underworld3/utilities/unit_aware_coordinates.py b/src/underworld3/utilities/unit_aware_coordinates.py
index 304aadebf..0073f2830 100644
--- a/src/underworld3/utilities/unit_aware_coordinates.py
+++ b/src/underworld3/utilities/unit_aware_coordinates.py
@@ -83,6 +83,37 @@ def get_units(self):
         return self._units
 
 
+class TimeSymbol(sympy.Symbol):
+    """A sympy Symbol subclass that carries _ccodestr and _units for JIT.
+
+    Standard sympy Symbols are immutable and don't allow setting arbitrary
+    attributes. This subclass permits _ccodestr (for C code generation)
+    and _units (for dimensional analysis), following the same pattern as
+    UnitAwareBaseScalar for spatial coordinates.
+    """
+
+    _ccodestr = "petsc_t"
+    _units = None
+
+    def __new__(cls, name="t", **kwargs):
+        obj = super().__new__(cls, name, real=True, nonnegative=True, **kwargs)
+        return obj
+
+    @property
+    def units(self):
+        return self._units
+
+    @units.setter
+    def units(self, value):
+        self._units = value
+
+    def get_units(self):
+        return self._units
+
+    def _ccode(self, printer):
+        return self._ccodestr
+
+
 def create_unit_aware_coordinate_system(name, units=None):
     """
     Create a coordinate system with unit-aware coordinates.
@@ -207,6 +238,40 @@ def _patch_coordinate(coord, units):
         # If _Gamma needs any units in the future, they should be explicitly
         # dimensionless (None), not inherited from mesh coordinates.
 
+        # Patch time coordinate with time units from the model
+        _patch_time_units(mesh)
+
+
+def _patch_time_units(mesh):
+    """
+    Patch mesh.t with time units from the model's scaling system.
+
+    If the model has a time scale defined, mesh.t gets those units.
+    Otherwise mesh.t remains dimensionless (._units = None).
+    """
+    if not hasattr(mesh, "_t"):
+        return
+
+    time_units = None
+    try:
+        import underworld3 as uw
+
+        model = uw.get_default_model()
+        if hasattr(model, "_fundamental_scales") and model._fundamental_scales:
+            scales = model._fundamental_scales
+            if "time" in scales:
+                time_scale = scales["time"]
+                if hasattr(time_scale, "units"):
+                    time_units = time_scale.units
+                elif hasattr(time_scale, "_pint_qty"):
+                    time_units = time_scale._pint_qty.units
+    except Exception:
+        pass
+
+    mesh._t._units = time_units
+    if not hasattr(mesh._t, "get_units"):
+        mesh._t.get_units = lambda: mesh._t._units
+
 
 def get_coordinate_units(coord):
     """
diff --git a/tests/plot_ve_oscillatory_validation.py b/tests/plot_ve_oscillatory_validation.py
new file mode 100644
index 000000000..f3ac103f7
--- /dev/null
+++ b/tests/plot_ve_oscillatory_validation.py
@@ -0,0 +1,239 @@
+"""Viscoelastic oscillatory shear — validation plot for PR.
+
+Produces a figure showing:
+  - Top panel: applied shear rate γ̇(t) = γ̇₀ sin(ωt)
+  - Bottom panel: stress σ_xy(t) for order 1 and order 2 vs analytical
+
+Maxwell analytical solution including startup transient:
+    σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)]
+
+Results are saved as .npz checkpoint files for re-analysis.
+
+Usage:
+    python tests/plot_ve_oscillatory_validation.py
+    python tests/plot_ve_oscillatory_validation.py --replot   # replot from saved data
+"""
+
+import numpy as np
+import sympy
+import os
+import sys
+import underworld3 as uw
+from underworld3.function import expression
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient)."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def run_oscillatory(order, n_steps, dt, V0, H, ETA, MU, omega, save_prefix=None):
+    """Run oscillatory VE shear box, return time series and save checkpoints."""
+    W = 2.0
+    gamma_dot_0 = 2.0 * V0 / H
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W / 2, -H / 2), maxCoords=(W / 2, H / 2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity")
+
+    stokes.add_dirichlet_bc((V_bc, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_bc, 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-6
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num = [], []
+    time_phys = 0.0
+
+    checkpoint_interval = max(1, n_steps // 20)
+
+    for step in range(n_steps):
+        time_phys += dt
+
+        V_t = V0 * np.sin(omega * time_phys)
+        V_bc.sym = V_t
+
+        stokes.solve(zero_init_guess=False, timestep=dt, evalf=False)
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+
+        if save_prefix and (step + 1) % checkpoint_interval == 0:
+            np.savez(
+                f"{save_prefix}_order{order}_checkpoint.npz",
+                times=np.array(times),
+                stress=np.array(stress_num),
+                step=step + 1,
+                dt=dt,
+                omega=omega,
+                ETA=ETA,
+                MU=MU,
+                V0=V0,
+                H=H,
+                order=order,
+            )
+
+    del stokes, mesh
+
+    t_arr = np.array(times)
+    s_arr = np.array(stress_num)
+
+    if save_prefix:
+        np.savez(
+            f"{save_prefix}_order{order}_final.npz",
+            times=t_arr,
+            stress=s_arr,
+            dt=dt,
+            omega=omega,
+            ETA=ETA,
+            MU=MU,
+            V0=V0,
+            H=H,
+            order=order,
+            gamma_dot_0=gamma_dot_0,
+        )
+
+    return t_arr, s_arr
+
+
+def make_plot(results, params, output_path):
+    """Generate the validation figure from results dict."""
+
+    ETA = params["ETA"]
+    MU = params["MU"]
+    V0 = params["V0"]
+    H = params["H"]
+    omega = params["omega"]
+    dt = params["dt"]
+    t_r = ETA / MU
+    De = omega * t_r
+    gamma_dot_0 = 2.0 * V0 / H
+
+    # Determine time range from results
+    t_max = max(r[0][-1] for r in results.values())
+
+    # Analytical solution (fine sampling)
+    t_fine = np.linspace(0, t_max, 1000)
+    sigma_analytical = maxwell_oscillatory(t_fine, ETA, MU, gamma_dot_0, omega)
+    shear_rate = gamma_dot_0 * np.sin(omega * t_fine)
+
+    phase_lag_deg = np.degrees(np.arctan(De))
+    steady_amp = ETA * gamma_dot_0 / np.sqrt(1 + De**2)
+    period = 2.0 * np.pi / omega
+
+    fig, (ax1, ax2) = plt.subplots(
+        2, 1, figsize=(10, 6), sharex=True,
+        gridspec_kw={"height_ratios": [1, 2.5]},
+    )
+
+    # Top panel: applied shear rate
+    ax1.plot(t_fine / t_r, shear_rate, "k-", linewidth=1.2,
+             label=r"$\dot{\gamma}(t) = \dot{\gamma}_0 \sin(\omega t)$")
+    ax1.set_ylabel(r"Shear rate $\dot{\gamma}$")
+    ax1.axhline(0, color="grey", linewidth=0.5)
+    ax1.legend(loc="upper right", fontsize=9)
+    ax1.set_xlim(0, t_max / t_r)
+
+    # Bottom panel: stress response
+    ax2.plot(t_fine / t_r, sigma_analytical, "k-", linewidth=1.5,
+             label="Analytical (Maxwell)", zorder=3)
+
+    colors = {1: "#2196F3", 2: "#E91E63", 3: "#4CAF50"}
+    markers = {1: "o", 2: "s", 3: "^"}
+
+    for order in sorted(results.keys()):
+        t, stress = results[order]
+        ax2.plot(t / t_r, stress, markers.get(order, "o"),
+                 color=colors.get(order, "grey"),
+                 markersize=2.5, alpha=0.7,
+                 label=f"Order {order} (BDF-{order})", zorder=2)
+
+    ax2.set_xlabel(r"Time $t / t_r$")
+    ax2.set_ylabel(r"Stress $\sigma_{xy}$")
+    ax2.axhline(0, color="grey", linewidth=0.5)
+    ax2.legend(loc="lower right", fontsize=9)
+
+    ax2.set_title(
+        f"Maxwell viscoelastic shear: De = {De:.1f}, "
+        r"$\delta t / t_r$" + f" = {dt / t_r:.2f}, "
+        f"phase lag = {phase_lag_deg:.0f}\u00b0, "
+        f"steady amplitude = {steady_amp:.3f}",
+        fontsize=10,
+    )
+
+    plt.tight_layout()
+    plt.savefig(output_path, dpi=150, bbox_inches="tight")
+    print(f"Saved: {output_path}")
+
+
+if __name__ == "__main__":
+
+    output_dir = "tests"
+    save_prefix = os.path.join(output_dir, "ve_oscillatory")
+
+    replot = "--replot" in sys.argv
+
+    # Physical parameters
+    ETA = 1.0
+    MU = 1.0
+    V0 = 0.5
+    H = 1.0
+    t_r = ETA / MU
+
+    De = 5.0  # High Deborah number — strong elastic effects, 79° phase lag
+    omega = De / t_r
+    period = 2.0 * np.pi / omega
+    n_periods = 4
+    dt = 0.02 * t_r  # Fine enough to resolve the fast oscillations
+    n_steps = int(n_periods * period / dt)
+
+    params = dict(ETA=ETA, MU=MU, V0=V0, H=H, omega=omega, dt=dt)
+
+    print(f"Maxwell oscillatory shear: De={De}, omega={omega:.3f}, period={period:.3f}")
+    print(f"dt={dt}, dt/t_r={dt/t_r}, n_steps={n_steps}, t_end={n_steps * dt:.2f}")
+    print(f"Phase lag = {np.degrees(np.arctan(De)):.1f} deg")
+    print()
+
+    if replot:
+        # Load saved data
+        results = {}
+        for order in [1, 2]:
+            fpath = f"{save_prefix}_order{order}_final.npz"
+            if os.path.exists(fpath):
+                data = np.load(fpath)
+                results[order] = (data["times"], data["stress"])
+                print(f"Loaded order {order} from {fpath}")
+            else:
+                print(f"Missing {fpath} — run without --replot first")
+                sys.exit(1)
+    else:
+        results = {}
+        for order in [1, 2]:
+            print(f"Running order {order} ({n_steps} steps)...", end=" ", flush=True)
+            t, stress = run_oscillatory(
+                order, n_steps, dt, V0, H, ETA, MU, omega,
+                save_prefix=save_prefix,
+            )
+            results[order] = (t, stress)
+            print("done")
+
+    make_plot(results, params, os.path.join(output_dir, "ve_oscillatory_validation.png"))
diff --git a/tests/plot_ve_vep_oscillatory.py b/tests/plot_ve_vep_oscillatory.py
new file mode 100644
index 000000000..abc234ca6
--- /dev/null
+++ b/tests/plot_ve_vep_oscillatory.py
@@ -0,0 +1,41 @@
+"""Plot VE vs VEP oscillatory shear from saved checkpoint."""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+data = np.load("output/ve_vep_oscillatory.npz")
+t_ve, s_ve = data["t_ve"], data["s_ve"]
+t_vep, s_vep = data["t_vep"], data["s_vep"]
+t_ana, s_ana = data["t_ana"], data["s_ana"]
+De = float(data["De"])
+tau_y = float(data["tau_y"])
+ve_amp = float(data["ve_amp"])
+omega = float(data["omega"])
+
+fig, ax = plt.subplots(1, 1, figsize=(10, 5))
+
+# Analytical VE (smooth curve)
+ax.plot(t_ana, s_ana, "k-", linewidth=0.8, alpha=0.5, label="VE analytical")
+
+# Numerical VE
+ax.plot(t_ve, s_ve, "b-o", markersize=3, linewidth=1.2, label="VE numerical (order 1)")
+
+# Numerical VEP
+ax.plot(t_vep, s_vep, "r-s", markersize=3, linewidth=1.2, label=f"VEP numerical ($\\tau_y$={tau_y})")
+
+# Yield stress lines
+ax.axhline(tau_y, color="r", linestyle="--", alpha=0.4, linewidth=0.8)
+ax.axhline(-tau_y, color="r", linestyle="--", alpha=0.4, linewidth=0.8)
+ax.text(0.3, tau_y + 0.02, f"$\\tau_y$ = {tau_y}", color="r", fontsize=9, alpha=0.6)
+
+ax.set_xlabel("Time ($t / t_r$)")
+ax.set_ylabel("$\\sigma_{xy}$")
+ax.set_title(f"Oscillatory Maxwell VE vs VEP shear (De = {De})")
+ax.legend(loc="upper left", fontsize=9)
+ax.set_xlim(0, t_ve[-1])
+ax.grid(True, alpha=0.3)
+
+plt.tight_layout()
+plt.savefig("output/ve_vep_oscillatory.png", dpi=150)
+plt.savefig("output/ve_vep_oscillatory.pdf")
+print("Saved output/ve_vep_oscillatory.png and .pdf")
diff --git a/tests/run_ve_order2_debug.py b/tests/run_ve_order2_debug.py
new file mode 100644
index 000000000..286b73595
--- /dev/null
+++ b/tests/run_ve_order2_debug.py
@@ -0,0 +1,81 @@
+"""Debug order-2 VE: trace psi_star values at each step."""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0
+dt = 0.1
+gamma_dot = 2.0 * V0 / H
+
+mesh = uw.meshing.StructuredQuadBox(
+    elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+)
+v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=2)
+stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+stokes.constitutive_model.Parameters.shear_modulus = MU
+stokes.constitutive_model.Parameters.dt_elastic = dt
+
+stokes.add_dirichlet_bc((V0, 0.0), "Top")
+stokes.add_dirichlet_bc((-V0, 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-6
+
+Stress = uw.discretisation.MeshVariable(
+    "Stress", mesh, (2, 2), vtype=uw.VarType.SYM_TENSOR, degree=2, continuous=True,
+)
+work = uw.discretisation.MeshVariable("W", mesh, 1, degree=2)
+sigma_proj = uw.systems.Tensor_Projection(mesh, tensor_Field=Stress, scalar_Field=work)
+
+ddt = stokes.DFDt
+centre = np.array([[0.0, 0.0]])
+
+time_phys = 0.0
+for step in range(5):
+    eff_order = stokes.constitutive_model.effective_order
+
+    # Read psi_star values at centre BEFORE solve
+    psi0_pre = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre)
+    psi1_pre = uw.function.evaluate(ddt.psi_star[1].sym[0, 1], centre)
+    psi0_val = float(psi0_pre.flatten()[0])
+    psi1_val = float(psi1_pre.flatten()[0])
+
+    # Check psi_fn formula
+    psi_fn_01 = str(stokes.DFDt.psi_fn[0, 1])
+    has_2star = "**" in psi_fn_01 or "psi_star" in psi_fn_01
+
+    t0 = timer.time()
+    stokes.solve(zero_init_guess=False, evalf=False)
+    solve_t = timer.time() - t0
+
+    psi_fn_01_post = str(stokes.DFDt.psi_fn[0, 1])
+    print(f"  psi_fn[0,1] BEFORE solve: {psi_fn_01[:100]}...")
+    print(f"  psi_fn[0,1] AFTER  solve: {psi_fn_01_post[:100]}...")
+    time_phys += dt
+
+    # Read psi_star AFTER solve (after update_post_solve)
+    psi0_post = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre)
+    psi1_post = uw.function.evaluate(ddt.psi_star[1].sym[0, 1], centre)
+    psi0_post_val = float(psi0_post.flatten()[0])
+    psi1_post_val = float(psi1_post.flatten()[0])
+
+    sigma_proj.uw_function = stokes.stress_deviator
+    sigma_proj.solve()
+    val = uw.function.evaluate(Stress.sym[0, 1], centre)
+    sigma_xy = float(val.flatten()[0])
+    ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA))
+
+    print(f"step {step}  eff_order={eff_order}  t={time_phys:.2f}  solve={solve_t:.1f}s")
+    print(f"  PRE:  psi_star[0]_xy={psi0_val:.6f}  psi_star[1]_xy={psi1_val:.6f}")
+    print(f"  POST: psi_star[0]_xy={psi0_post_val:.6f}  psi_star[1]_xy={psi1_post_val:.6f}")
+    print(f"  sigma_xy={sigma_xy:.6f}  analytical={ana:.6f}  "
+          f"rel_err={abs(sigma_xy - ana) / ana:.3e}")
+    print()
+
+print("Done")
diff --git a/tests/run_ve_oscillatory.py b/tests/run_ve_oscillatory.py
new file mode 100644
index 000000000..2cf7e21da
--- /dev/null
+++ b/tests/run_ve_oscillatory.py
@@ -0,0 +1,122 @@
+"""Oscillatory VE shear box — time-harmonic analytical validation.
+
+Maxwell material under sinusoidal boundary velocity V(t) = V0 sin(ωt).
+Full solution including startup transient:
+
+    σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)]
+
+where De = ω t_r is the Deborah number.
+
+At steady state (t >> t_r), the transient dies out and the stress
+oscillates with amplitude η γ̇₀/√(1+De²) and phase lag arctan(De).
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+from underworld3.function import expression
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient)."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def run_oscillatory(order, n_steps, dt_over_tr, De):
+    """Run oscillatory VE shear box."""
+
+    ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0
+    t_r = ETA / MU
+    omega = De / t_r
+    V0 = 0.5
+    gamma_dot_0 = 2.0 * V0 / H
+    dt = dt_over_tr * t_r
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    # Use a UWexpression for time-dependent BCs — avoids JIT recompilation
+    V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity")
+
+    stokes.add_dirichlet_bc((V_bc, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_bc, 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-6
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num, stress_ana = [], [], []
+    time_phys = 0.0
+
+    for step in range(n_steps):
+        time_phys += dt
+
+        # Update boundary velocity for this timestep
+        V_t = V0 * np.sin(omega * time_phys)
+        V_bc.sym = V_t
+
+        stokes.solve(zero_init_guess=False, evalf=False)
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+
+        ana = maxwell_oscillatory(time_phys, ETA, MU, gamma_dot_0, omega)
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+        stress_ana.append(ana)
+
+    del stokes, mesh
+    return np.array(times), np.array(stress_num), np.array(stress_ana)
+
+
+if __name__ == "__main__":
+    De = 1.0        # Deborah number (equal viscous and elastic timescales)
+    dt_ratio = 0.1   # coarser for speed; finer (0.05) for publication quality
+    n_periods = 3
+    t_r = 1.0
+    omega = De / t_r
+    period = 2 * np.pi / omega
+    n_steps = int(n_periods * period / (dt_ratio * t_r))
+
+    print(f"De = {De}, omega = {omega:.3f}, period = {period:.3f}")
+    print(f"dt/t_r = {dt_ratio}, n_steps = {n_steps}")
+    print()
+
+    for order in [1, 2]:
+        t0 = timer.time()
+        t, num, ana = run_oscillatory(order, n_steps, dt_ratio, De)
+        wall = timer.time() - t0
+
+        # Error over the last period (steady state)
+        last_period = t > (n_periods - 1) * period
+        if np.any(last_period):
+            rms = np.sqrt(np.mean((num[last_period] - ana[last_period])**2))
+            amp_ana = np.max(np.abs(ana[last_period]))
+            rel_rms = rms / amp_ana if amp_ana > 0 else float('nan')
+        else:
+            rel_rms = float('nan')
+
+        print(f"Order {order}: wall={wall:.0f}s, "
+              f"last-period relative RMS = {rel_rms:.4e}")
+
+        # Samples from the last period
+        if np.any(last_period):
+            idx = np.where(last_period)[0]
+            sample = idx[::max(1, len(idx)//10)]
+            for i in sample:
+                err = abs(num[i] - ana[i])
+                print(f"  t={t[i]:.3f}  num={num[i]:+.6f}  ana={ana[i]:+.6f}  err={err:.4e}")
+        print()
diff --git a/tests/run_ve_shear_order2_quick.py b/tests/run_ve_shear_order2_quick.py
new file mode 100644
index 000000000..fb40d71f9
--- /dev/null
+++ b/tests/run_ve_shear_order2_quick.py
@@ -0,0 +1,48 @@
+"""Quick order-2 VE shear box test — read stress from psi_star[0]."""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0
+dt = 0.1
+gamma_dot = 2.0 * V0 / H
+
+mesh = uw.meshing.StructuredQuadBox(
+    elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+)
+v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=2)
+stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+stokes.constitutive_model.Parameters.shear_modulus = MU
+stokes.constitutive_model.Parameters.dt_elastic = dt
+
+stokes.add_dirichlet_bc((V0, 0.0), "Top")
+stokes.add_dirichlet_bc((-V0, 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-6
+
+centre = np.array([[0.0, 0.0]])
+ddt = stokes.DFDt
+
+time_phys = 0.0
+for step in range(20):
+    eff = stokes.constitutive_model.effective_order
+    t0 = timer.time()
+    stokes.solve(zero_init_guess=False, evalf=False)
+    solve_t = timer.time() - t0
+    time_phys += dt
+
+    val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+    sigma_xy = float(val.flatten()[0])
+    ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA))
+    print(f"step {step:2d}  eff={eff}  t={time_phys:.2f}  solve={solve_t:.1f}s  "
+          f"sigma_xy={sigma_xy:.6f}  analytical={ana:.6f}  "
+          f"rel_err={abs(sigma_xy - ana) / ana:.3e}")
+
+print("Done")
diff --git a/tests/run_ve_shear_quick.py b/tests/run_ve_shear_quick.py
new file mode 100644
index 000000000..debced41d
--- /dev/null
+++ b/tests/run_ve_shear_quick.py
@@ -0,0 +1,48 @@
+"""Quick order-1 VE shear box validation (10 steps)."""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0
+dt = 0.1
+gamma_dot = 2.0 * V0 / H
+
+mesh = uw.meshing.StructuredQuadBox(
+    elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+)
+v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=1)
+stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+stokes.constitutive_model.Parameters.shear_modulus = MU
+stokes.constitutive_model.Parameters.dt_elastic = dt
+
+stokes.add_dirichlet_bc((V0, 0.0), "Top")
+stokes.add_dirichlet_bc((-V0, 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-6
+
+centre = np.array([[0.0, 0.0]])
+ddt = stokes.DFDt
+
+time_phys = 0.0
+for step in range(10):
+    t0 = timer.time()
+    stokes.solve(zero_init_guess=False, evalf=False)
+    solve_t = timer.time() - t0
+    time_phys += dt
+
+    # Read stress from solver.tau (the projected actual stress)
+    val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+    sigma_xy = float(val.flatten()[0])
+    ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA))
+    print(f"step {step}  t={time_phys:.2f}  solve={solve_t:.1f}s  "
+          f"sigma_xy={sigma_xy:.6f}  analytical={ana:.6f}  "
+          f"rel_err={abs(sigma_xy - ana) / ana:.3e}")
+
+print("Done")
diff --git a/tests/run_ve_shear_validation.py b/tests/run_ve_shear_validation.py
new file mode 100644
index 000000000..759616edb
--- /dev/null
+++ b/tests/run_ve_shear_validation.py
@@ -0,0 +1,112 @@
+"""Quick validation script for the VE shear box test.
+
+Run with: pixi run -e amr-dev python tests/run_ve_shear_validation.py
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+
+def maxwell_xy(t, eta, mu, gamma_dot):
+    """Analytical σ_xy for Maxwell material under constant shear rate."""
+    return eta * gamma_dot * (1.0 - np.exp(-t * mu / eta))
+
+
+def run(order, n_steps, dt_ratio):
+    ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0
+    t_relax = ETA / MU
+    dt = dt_ratio * t_relax
+    gamma_dot = 2.0 * V0 / H
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8),
+        minCoords=(-W / 2, -H / 2),
+        maxCoords=(W / 2, H / 2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    stokes.add_dirichlet_bc((V0, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V0, 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-6
+
+    Stress = uw.discretisation.MeshVariable(
+        "Stress", mesh, (2, 2),
+        vtype=uw.VarType.SYM_TENSOR, degree=2, continuous=True,
+    )
+    work = uw.discretisation.MeshVariable("W", mesh, 1, degree=2)
+    sigma_proj = uw.systems.Tensor_Projection(mesh, tensor_Field=Stress, scalar_Field=work)
+
+    times, num, ana = [], [], []
+    time_phys = 0.0
+
+    for step in range(n_steps):
+        t0 = timer.time()
+        stokes.solve(zero_init_guess=False, evalf=False)
+        solve_time = timer.time() - t0
+        time_phys += dt
+
+        sigma_proj.uw_function = stokes.stress_deviator
+        sigma_proj.solve()
+
+        val = uw.function.evaluate(Stress.sym[0, 1], np.array([[0.0, 0.0]]))
+        sigma_xy = float(val.flatten()[0])
+        ana_val = maxwell_xy(time_phys, ETA, MU, gamma_dot)
+
+        times.append(time_phys)
+        num.append(sigma_xy)
+        ana.append(ana_val)
+
+        rel_err = abs(sigma_xy - ana_val) / max(ana_val, 1e-10)
+        print(f"  step {step:3d}  t={time_phys:.2f}  solve={solve_time:.1f}s  "
+              f"sigma_xy={sigma_xy:.6f}  analytical={ana_val:.6f}  rel_err={rel_err:.3e}")
+
+    del stokes, mesh
+    return np.array(times), np.array(num), np.array(ana)
+
+
+if __name__ == "__main__":
+    print("=== Order 1, dt/tr=0.1, 20 steps ===")
+    t0 = timer.time()
+    t, n, a = run(1, 20, 0.1)
+    wall = timer.time() - t0
+    mask = a > 0.01 * a[-1]
+    rms = np.sqrt(np.mean(((n[mask] - a[mask]) / a[mask]) ** 2))
+    final_err = abs(n[-1] - a[-1]) / a[-1]
+    print(f"  Wall time: {wall:.0f}s")
+    print(f"  Final rel error: {final_err:.4e}")
+    print(f"  RMS rel error:   {rms:.4e}")
+    print()
+
+    print("=== Order 2, dt/tr=0.1, 20 steps ===")
+    t0 = timer.time()
+    t2, n2, a2 = run(2, 20, 0.1)
+    wall2 = timer.time() - t0
+    mask2 = a2 > 0.01 * a2[-1]
+    rms2 = np.sqrt(np.mean(((n2[mask2] - a2[mask2]) / a2[mask2]) ** 2))
+    final_err2 = abs(n2[-1] - a2[-1]) / a2[-1]
+    print(f"  Wall time: {wall2:.0f}s")
+    print(f"  Final rel error: {final_err2:.4e}")
+    print(f"  RMS rel error:   {rms2:.4e}")
+    print()
+
+    viscous_limit = 1.0
+    print(f"Steady-state check (order 1, t=2.0):")
+    print(f"  sigma_xy = {n[-1]:.6f}, viscous limit = {viscous_limit:.6f}, "
+          f"rel error = {abs(n[-1] - viscous_limit) / viscous_limit:.4e}")
+    print()
+
+    if rms2 < rms:
+        print("PASS: Order 2 has smaller RMS error than order 1")
+    else:
+        print("FAIL: Order 2 should have smaller error than order 1")
diff --git a/tests/run_ve_vep_oscillatory_checkpoint.py b/tests/run_ve_vep_oscillatory_checkpoint.py
new file mode 100644
index 000000000..8ae07f015
--- /dev/null
+++ b/tests/run_ve_vep_oscillatory_checkpoint.py
@@ -0,0 +1,124 @@
+"""Run VE and VEP oscillatory shear, save time series to .npz for plotting."""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+from underworld3.function import expression
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for oscillatory Maxwell shear."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def run_oscillatory(order, n_steps, dt, omega, V0, tau_y=None):
+    """Run oscillatory shear (VE or VEP depending on tau_y)."""
+
+    ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0
+    gamma_dot_0 = 2.0 * V0 / H
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    if tau_y is not None:
+        stokes.constitutive_model.Parameters.yield_stress = tau_y
+        stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6
+
+    v.data[:, 0] = v.coords[:, 1] * gamma_dot_0
+    v.data[:, 1] = 0.0
+
+    V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity")
+    stokes.add_dirichlet_bc((V_bc, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_bc, 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-5
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    stokes.petsc_options["snes_max_it"] = 50
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num = [], []
+    time_phys = 0.0
+
+    for step in range(n_steps):
+        time_phys += dt
+        V_bc.sym = V0 * np.sin(omega * time_phys)
+
+        stokes.solve(zero_init_guess=False, evalf=False)
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+
+        label = "VEP" if tau_y is not None else "VE"
+        print(f"  {label} step {step:3d}  t={time_phys:.2f}  sigma={sigma_xy:+.6f}", flush=True)
+
+    del stokes, mesh
+    return np.array(times), np.array(stress_num)
+
+
+if __name__ == "__main__":
+    De = 1.0
+    ETA, MU = 1.0, 1.0
+    t_r = ETA / MU
+    omega = De / t_r
+    V0 = 0.5
+    gamma_dot_0 = 2.0 * V0 / 1.0
+    period = 2.0 * np.pi / omega
+    TAU_Y = 0.4
+
+    dt_ratio = 0.1
+    dt = dt_ratio * t_r
+    n_periods = 2
+    n_steps = int(n_periods * period / dt)
+
+    ve_amp = ETA * gamma_dot_0 / np.sqrt(1.0 + De**2)
+
+    print(f"De={De}, omega={omega:.3f}, period={period:.3f}")
+    print(f"VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}")
+    print(f"dt={dt}, n_steps={n_steps}")
+    print()
+
+    # VE (no yield)
+    print("=== VE (order 1) ===")
+    t0 = timer.time()
+    t_ve, s_ve = run_oscillatory(1, n_steps, dt, omega, V0, tau_y=None)
+    print(f"  Wall: {timer.time()-t0:.0f}s\n")
+
+    # VEP (with yield)
+    print("=== VEP (order 1) ===")
+    t0 = timer.time()
+    t_vep, s_vep = run_oscillatory(1, n_steps, dt, omega, V0, tau_y=TAU_Y)
+    print(f"  Wall: {timer.time()-t0:.0f}s\n")
+
+    # Analytical
+    t_ana = np.linspace(dt, n_periods * period, 500)
+    s_ana = maxwell_oscillatory(t_ana, ETA, MU, gamma_dot_0, omega)
+
+    # Save
+    outfile = "output/ve_vep_oscillatory.npz"
+    import os
+    os.makedirs("output", exist_ok=True)
+    np.savez(
+        outfile,
+        t_ve=t_ve, s_ve=s_ve,
+        t_vep=t_vep, s_vep=s_vep,
+        t_ana=t_ana, s_ana=s_ana,
+        De=De, tau_y=TAU_Y, ve_amp=ve_amp, omega=omega,
+    )
+    print(f"Saved to {outfile}")
diff --git a/tests/run_ve_vep_oscillatory_plot.py b/tests/run_ve_vep_oscillatory_plot.py
new file mode 100644
index 000000000..b5ee4907c
--- /dev/null
+++ b/tests/run_ve_vep_oscillatory_plot.py
@@ -0,0 +1,161 @@
+"""Run VE and VEP oscillatory shear from consistent initial conditions.
+
+Both cases start from the analytical Maxwell solution at t0=0.01 to avoid
+zero-velocity singularity and ensure identical starting stress.
+
+Saves time series to stdout and generates plot directly.
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+from underworld3.function import expression
+import os
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for oscillatory Maxwell shear."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def run_oscillatory(order, n_steps, dt, omega, V0, t0, tau_y=None):
+    """Run oscillatory shear from analytical initial condition at t0."""
+
+    ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0
+    gamma_dot_0 = 2.0 * V0 / H
+    t_r = ETA / MU
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    if tau_y is not None:
+        stokes.constitutive_model.Parameters.yield_stress = tau_y
+        stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6
+
+    # Initialise velocity to match analytical solution at t0
+    V_t0 = V0 * np.sin(omega * t0)
+    gamma_dot_t0 = 2.0 * V_t0 / H
+    v.data[:, 0] = v.coords[:, 1] * gamma_dot_t0
+    v.data[:, 1] = 0.0
+
+    # Initialise stress history to analytical at t0
+    sigma_t0 = maxwell_oscillatory(t0, ETA, MU, gamma_dot_0, omega)
+    # psi_star[0] is (N, dim, dim) sym tensor — set the xy component
+    stokes.DFDt.psi_star[0].array[:, 0, 1] = sigma_t0
+    stokes.DFDt.psi_star[0].array[:, 1, 0] = sigma_t0
+    stokes.DFDt._history_initialised = True
+
+    V_bc = expression(r"{V_{bc}}", V_t0, "Time-dependent boundary velocity")
+    stokes.add_dirichlet_bc((V_bc, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_bc, 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-5
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    stokes.petsc_options["snes_max_it"] = 50
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num = [], []
+    time_phys = t0
+
+    label = "VEP" if tau_y is not None else "VE"
+
+    for step in range(n_steps):
+        time_phys += dt
+        V_bc.sym = V0 * np.sin(omega * time_phys)
+
+        stokes.solve(zero_init_guess=False, evalf=False)
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+
+        print(f"  {label} step {step:3d}  t={time_phys:.3f}  sigma={sigma_xy:+.6f}", flush=True)
+
+    del stokes, mesh
+    return np.array(times), np.array(stress_num)
+
+
+if __name__ == "__main__":
+    De = 1.0
+    ETA, MU = 1.0, 1.0
+    t_r = ETA / MU
+    omega = De / t_r
+    V0 = 0.5
+    gamma_dot_0 = 2.0 * V0 / 1.0
+    period = 2.0 * np.pi / omega
+    TAU_Y = 0.4
+
+    dt_ratio = 0.1
+    dt = dt_ratio * t_r
+    t0 = 0.01
+    n_periods = 2
+    n_steps = int(n_periods * period / dt)
+
+    ve_amp = ETA * gamma_dot_0 / np.sqrt(1.0 + De**2)
+    sigma_t0 = maxwell_oscillatory(t0, ETA, MU, gamma_dot_0, omega)
+
+    print(f"De={De}, omega={omega:.3f}, period={period:.3f}")
+    print(f"VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}")
+    print(f"t0={t0}, sigma(t0)={sigma_t0:.6f}")
+    print(f"dt={dt}, n_steps={n_steps}")
+    print()
+
+    # VE
+    print("=== VE (order 1) ===")
+    t0w = timer.time()
+    t_ve, s_ve = run_oscillatory(1, n_steps, dt, omega, V0, t0, tau_y=None)
+    print(f"  Wall: {timer.time()-t0w:.0f}s\n")
+
+    # VEP
+    print("=== VEP (order 1) ===")
+    t0w = timer.time()
+    t_vep, s_vep = run_oscillatory(1, n_steps, dt, omega, V0, t0, tau_y=TAU_Y)
+    print(f"  Wall: {timer.time()-t0w:.0f}s\n")
+
+    # Analytical
+    t_ana = np.linspace(t0, t_ve[-1], 500)
+    s_ana = maxwell_oscillatory(t_ana, ETA, MU, gamma_dot_0, omega)
+
+    # Plot
+    import matplotlib
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(1, 1, figsize=(10, 5))
+
+    ax.plot(t_ana, s_ana, "k-", linewidth=0.8, alpha=0.5, label="VE analytical")
+    ax.plot(t_ve, s_ve, "b-o", markersize=3, linewidth=1.2, label="VE numerical (order 1)")
+    ax.plot(t_vep, s_vep, "r-s", markersize=3, linewidth=1.2, label=f"VEP numerical ($\\tau_y$={TAU_Y})")
+
+    ax.axhline(TAU_Y, color="r", linestyle="--", alpha=0.4, linewidth=0.8)
+    ax.axhline(-TAU_Y, color="r", linestyle="--", alpha=0.4, linewidth=0.8)
+    ax.text(0.3, TAU_Y + 0.02, f"$\\tau_y$ = {TAU_Y}", color="r", fontsize=9, alpha=0.6)
+
+    ax.set_xlabel("Time ($t / t_r$)")
+    ax.set_ylabel("$\\sigma_{xy}$")
+    ax.set_title(f"Oscillatory Maxwell VE vs VEP shear (De = {De})")
+    ax.legend(loc="upper left", fontsize=9)
+    ax.set_xlim(0, t_ve[-1])
+    ax.grid(True, alpha=0.3)
+
+    plt.tight_layout()
+    os.makedirs("output", exist_ok=True)
+    plt.savefig("output/ve_vep_oscillatory.png", dpi=150)
+    plt.savefig("output/ve_vep_oscillatory.pdf")
+    print("Saved output/ve_vep_oscillatory.png and .pdf")
diff --git a/tests/run_vep_oscillatory.py b/tests/run_vep_oscillatory.py
new file mode 100644
index 000000000..e76e212e5
--- /dev/null
+++ b/tests/run_vep_oscillatory.py
@@ -0,0 +1,134 @@
+"""Oscillatory VEP shear box — does the yield stress cap the oscillation?
+
+Same setup as the VE oscillatory test but with a yield stress.
+The VE amplitude at steady state is η γ̇₀/√(1+De²).
+If τ_y is below this amplitude, the stress should be clipped.
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+from underworld3.function import expression
+
+
+def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega):
+    """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient)."""
+    t_r = eta / mu
+    De = omega * t_r
+    prefactor = eta * gamma_dot_0 / (1.0 + De**2)
+    return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r))
+
+
+def run_vep_oscillatory(order, n_steps, dt_over_tr, De, tau_y):
+    """Run oscillatory VEP shear box."""
+
+    ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0
+    t_r = ETA / MU
+    omega = De / t_r
+    V0 = 0.5
+    gamma_dot_0 = 2.0 * V0 / H
+    dt = dt_over_tr * t_r
+
+    # VE steady-state amplitude
+    ve_amplitude = ETA * gamma_dot_0 / np.sqrt(1 + De**2)
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+    stokes.constitutive_model.Parameters.yield_stress = tau_y
+    stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6
+
+    # Seed with linear shear
+    v.data[:, 0] = v.coords[:, 1] * gamma_dot_0
+    v.data[:, 1] = 0.0
+
+    V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity")
+    stokes.add_dirichlet_bc((V_bc, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V_bc, 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-5
+
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    stokes.petsc_options["snes_max_it"] = 50
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num, stress_ve = [], [], []
+    time_phys = 0.0
+
+    for step in range(n_steps):
+        time_phys += dt
+        V_t = V0 * np.sin(omega * time_phys)
+        V_bc.sym = V_t
+
+        t0s = __import__('time').time()
+        stokes.solve(zero_init_guess=False, evalf=False)
+        solve_t = __import__('time').time() - t0s
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+        ve_stress = maxwell_oscillatory(time_phys, ETA, MU, gamma_dot_0, omega)
+
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+        stress_ve.append(ve_stress)
+        print(f"  step {step:3d}  t={time_phys:.2f}  solve={solve_t:.1f}s  "
+              f"vep={sigma_xy:+.6f}  ve={ve_stress:+.6f}", flush=True)
+
+    del stokes, mesh
+    return np.array(times), np.array(stress_num), np.array(stress_ve), ve_amplitude
+
+
+if __name__ == "__main__":
+    De = 1.0
+    dt_ratio = 0.1
+    n_periods = 2
+    t_r = 1.0
+    omega = De / t_r
+    period = 2 * np.pi / omega
+    n_steps = int(n_periods * period / (dt_ratio * t_r))
+
+    ETA, gamma_dot_0 = 1.0, 1.0
+    ve_amp = ETA * gamma_dot_0 / np.sqrt(1 + De**2)
+
+    # Set τ_y below the VE amplitude so we see clipping
+    TAU_Y = 0.4
+
+    print(f"De={De}, VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}")
+    print(f"Expect clipping at ±{TAU_Y}")
+    print(f"dt/t_r={dt_ratio}, n_steps={n_steps}")
+    print()
+
+    for order in [1]:
+        t0 = timer.time()
+        t, num, ve, amp = run_vep_oscillatory(order, n_steps, dt_ratio, De, TAU_Y)
+        wall = timer.time() - t0
+
+        print(f"Order {order} (wall={wall:.0f}s):")
+        print(f"  VE amplitude: {amp:.4f}")
+        print(f"  Max |σ_xy|:   {np.max(np.abs(num)):.6f}")
+        print(f"  Min σ_xy:     {np.min(num):.6f}")
+        print(f"  Max σ_xy:     {np.max(num):.6f}")
+        print(f"  Clipped at τ_y={TAU_Y}? {np.max(np.abs(num)) <= TAU_Y + 0.01}")
+        print()
+
+        # Print last period
+        last_period = t > (n_periods - 1) * period
+        if np.any(last_period):
+            idx = np.where(last_period)[0]
+            sample = idx[::max(1, len(idx)//12)]
+            print(f"  Last period samples:")
+            for i in sample:
+                clipped = " CLIPPED" if abs(num[i]) > TAU_Y - 0.01 else ""
+                print(f"    t={t[i]:.3f}  vep={num[i]:+.6f}  ve={ve[i]:+.6f}{clipped}")
+        print()
diff --git a/tests/run_vep_shear_box.py b/tests/run_vep_shear_box.py
new file mode 100644
index 000000000..7d4fabe3f
--- /dev/null
+++ b/tests/run_vep_shear_box.py
@@ -0,0 +1,121 @@
+"""Viscoelastic-plastic shear box — stress buildup with yield cap.
+
+Maxwell VE material with Drucker-Prager yield stress under constant shear.
+
+Phase 1 (elastic buildup): σ_xy follows the Maxwell curve
+Phase 2 (yielded): σ_xy = τ_y (plastic cap)
+
+Transition time: t_yield = -t_r ln(1 - τ_y/(η γ̇))
+(only exists if η γ̇ > τ_y)
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+
+def maxwell_stress_xy(t, eta, mu, gamma_dot):
+    """Uncapped VE stress (for comparison)."""
+    t_r = eta / mu
+    return eta * gamma_dot * (1.0 - np.exp(-t / t_r))
+
+
+def vep_stress_xy(t, eta, mu, gamma_dot, tau_y):
+    """VEP stress: Maxwell buildup capped at τ_y."""
+    ve = maxwell_stress_xy(t, eta, mu, gamma_dot)
+    return np.minimum(ve, tau_y)
+
+
+def run_vep_shear(order, n_steps, dt_over_tr, tau_y):
+    """Run VEP shear box, return (times, num_stress, ana_stress)."""
+
+    ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0
+    t_r = ETA / MU
+    dt = dt_over_tr * t_r
+    gamma_dot = 2.0 * V0 / H  # = 1.0
+
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order)
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+    stokes.constitutive_model.Parameters.yield_stress = tau_y
+    stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6
+
+    stokes.add_dirichlet_bc((V0, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V0, 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-5
+
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    stokes.petsc_options["snes_max_it"] = 50
+
+    # Seed with linear shear
+    v.data[:, 0] = v.coords[:, 1] * gamma_dot
+    v.data[:, 1] = 0.0
+
+    centre = np.array([[0.0, 0.0]])
+    times, stress_num, stress_ana, stress_ve = [], [], [], []
+    time_phys = 0.0
+
+    for step in range(n_steps):
+        t0 = timer.time()
+        stokes.solve(zero_init_guess=False, evalf=False)
+        solve_t = timer.time() - t0
+        time_phys += dt
+
+        val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+        sigma_xy = float(val.flatten()[0])
+
+        ana = vep_stress_xy(time_phys, ETA, MU, gamma_dot, tau_y)
+        ve_only = maxwell_stress_xy(time_phys, ETA, MU, gamma_dot)
+
+        times.append(time_phys)
+        stress_num.append(sigma_xy)
+        stress_ana.append(ana)
+        stress_ve.append(ve_only)
+
+    del stokes, mesh
+    return np.array(times), np.array(stress_num), np.array(stress_ana), np.array(stress_ve)
+
+
+if __name__ == "__main__":
+    ETA, MU = 1.0, 1.0
+    gamma_dot = 1.0  # 2*V0/H
+    t_r = ETA / MU
+
+    # Choose τ_y below the viscous steady state (η·γ̇ = 1.0)
+    TAU_Y = 0.5
+    t_yield = -t_r * np.log(1 - TAU_Y / (ETA * gamma_dot))
+
+    print(f"eta={ETA}, mu={MU}, gamma_dot={gamma_dot}")
+    print(f"tau_y={TAU_Y}, viscous_limit={ETA*gamma_dot}")
+    print(f"t_yield={t_yield:.4f} (stress reaches tau_y)")
+    print()
+
+    for order in [1, 2]:
+        t0 = timer.time()
+        t, num, ana, ve = run_vep_shear(order, n_steps=30, dt_over_tr=0.1, tau_y=TAU_Y)
+        wall = timer.time() - t0
+
+        print(f"Order {order} (wall={wall:.0f}s):")
+        for i in range(len(t)):
+            marker = " *" if t[i] > t_yield - 0.05 and t[i] < t_yield + 0.15 else ""
+            print(f"  t={t[i]:.2f}  num={num[i]:.6f}  ana={ana[i]:.6f}  "
+                  f"ve={ve[i]:.6f}  err={abs(num[i]-ana[i]):.4e}{marker}")
+
+        # Error in post-yield regime
+        yielded = t > t_yield + 0.2
+        if np.any(yielded):
+            post_err = np.max(np.abs(num[yielded] - ana[yielded]))
+            print(f"  Max post-yield error: {post_err:.4e}")
+        print()
diff --git a/tests/run_vp_shear_box.py b/tests/run_vp_shear_box.py
new file mode 100644
index 000000000..01a0c1651
--- /dev/null
+++ b/tests/run_vp_shear_box.py
@@ -0,0 +1,90 @@
+"""Viscoplastic shear box — yield stress validation.
+
+For simple shear with constant strain rate:
+  - Below yield: σ_xy = η · γ̇  (viscous)
+  - At/above yield: σ_xy = τ_y  (plastic cap)
+
+The transition occurs when η · γ̇ = τ_y, i.e. γ̇_crit = τ_y / η.
+
+We test with a range of driving velocities and check the stress
+transition matches the analytical prediction.
+"""
+
+import time as timer
+import numpy as np
+import sympy
+import underworld3 as uw
+
+
+def run_vp_shear(V0, eta, yield_stress):
+    """Single-step VP shear solve, return σ_xy at centre."""
+
+    H, W = 1.0, 2.0
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2),
+    )
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = eta
+    stokes.constitutive_model.Parameters.yield_stress = yield_stress
+    stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6
+
+    stokes.add_dirichlet_bc((V0, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V0, 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-4
+
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+    stokes.petsc_options["snes_max_it"] = 50
+
+    # Initialise with uniform shear to avoid singular Jacobian at zero strain rate
+    gamma_dot = 2.0 * V0 / H
+    v.data[:, 0] = v.coords[:, 1] * gamma_dot
+    v.data[:, 1] = 0.0
+
+    stokes.solve(zero_init_guess=False)
+
+    # Read stress via tau
+    centre = np.array([[0.0, 0.0]])
+    val = uw.function.evaluate(stokes.tau.sym[0, 1], centre)
+    sigma_xy = float(val.flatten()[0])
+
+    converged = stokes.snes.getConvergedReason() > 0
+    del stokes, mesh
+    return sigma_xy, converged
+
+
+if __name__ == "__main__":
+    ETA = 1.0
+    TAU_Y = 0.5
+    H = 1.0
+
+    # γ̇ = 2V0/H, σ_xy = η·γ̇ in viscous regime, σ_xy = τ_y when η·γ̇ ≥ τ_y
+    gamma_dot_crit = TAU_Y / ETA
+    print(f"eta={ETA}, tau_y={TAU_Y}, gamma_dot_crit={gamma_dot_crit}")
+    print()
+
+    V0_values = np.array([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1.0, 2.0])
+
+    print(f"{'V0':>6} {'gamma_dot':>10} {'sigma_xy':>10} {'analytical':>10} {'rel_err':>10} {'status':>8}")
+    print("-" * 62)
+
+    for V0 in V0_values:
+        gamma_dot = 2.0 * V0 / H
+        ana = min(ETA * gamma_dot, TAU_Y)
+
+        t0 = timer.time()
+        sigma_xy, converged = run_vp_shear(V0, ETA, TAU_Y)
+        dt = timer.time() - t0
+
+        rel_err = abs(sigma_xy - ana) / max(ana, 1e-10)
+        status = "OK" if converged else "FAIL"
+        print(f"{V0:6.3f} {gamma_dot:10.4f} {sigma_xy:10.6f} {ana:10.6f} {rel_err:10.3e} {status:>8}")
+
+    print()
+    print("Done")
diff --git a/tests/test_1051_VE_shear_box.py b/tests/test_1051_VE_shear_box.py
new file mode 100644
index 000000000..1608a7291
--- /dev/null
+++ b/tests/test_1051_VE_shear_box.py
@@ -0,0 +1,181 @@
+"""
+Viscoelastic shear box — analytical validation.
+
+A Maxwell viscoelastic material under uniform simple shear has the exact
+solution for the shear stress component:
+
+    σ_xy(t) = η · γ̇ · (1 - exp(-t / t_r))
+
+where:
+    t_r = η/μ           Maxwell relaxation time
+    γ̇  = dv_x/dy       engineering shear rate (velocity gradient)
+    ε̇_xy = γ̇/2         tensor strain rate
+
+The steady-state stress is σ_xy = η·γ̇ = 2η·ε̇_xy, the purely viscous limit.
+
+Boundary conditions:
+    Top / Bottom:  prescribed horizontal velocity ±V₀, zero vertical
+    Left / Right:  free horizontal (outflow), zero vertical velocity
+
+This avoids periodic BCs (untested) and corner singularities from free-slip.
+The result is a spatially uniform simple shear — effectively a 1D problem.
+
+The stress is read from psi_star[0] (the projected actual stress stored by
+VE_Stokes after each solve), not recomputed from the constitutive formula.
+"""
+
+import pytest
+import numpy as np
+import sympy
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_3, pytest.mark.tier_b]
+
+
+# ---------------------------------------------------------------------------
+# Analytical solution
+# ---------------------------------------------------------------------------
+
+def maxwell_stress_xy(t, eta, mu, gamma_dot):
+    """Exact σ_xy for constant-rate simple shear of a Maxwell material."""
+    t_relax = eta / mu
+    return eta * gamma_dot * (1.0 - np.exp(-t / t_relax))
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _run_ve_shear(order, n_steps, dt_over_tr):
+    """Run a VE shear-box, return (times, numerical_stress, analytical_stress)."""
+
+    ETA = 1.0
+    MU = 1.0
+    V0 = 0.5
+    H = 1.0
+    W = 2.0
+
+    t_relax = ETA / MU
+    dt = dt_over_tr * t_relax
+    gamma_dot = 2.0 * V0 / H
+
+    res = 8
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(2 * res, res),
+        minCoords=(-W / 2.0, -H / 2.0),
+        maxCoords=(W / 2.0, H / 2.0),
+    )
+
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)
+
+    stokes = uw.systems.VE_Stokes(
+        mesh, velocityField=v, pressureField=p, order=order, verbose=False,
+    )
+
+    stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA
+    stokes.constitutive_model.Parameters.shear_modulus = MU
+    stokes.constitutive_model.Parameters.dt_elastic = dt
+
+    stokes.add_dirichlet_bc((V0, 0.0), "Top")
+    stokes.add_dirichlet_bc((-V0, 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-6
+    stokes.petsc_options["snes_type"] = "newtonls"
+    stokes.petsc_options["ksp_type"] = "fgmres"
+
+    times = []
+    stress_num = []
+    stress_ana = []
+    centre = np.array([[0.0, 0.0]])
+
+    time = 0.0
+    for step in range(n_steps):
+        stokes.solve(zero_init_guess=False, evalf=False)
+        time += dt
+
+        # Read stress from solver.tau — the actual projected stress
+        sigma_xy_val = uw.function.evaluate(
+            stokes.tau.sym[0, 1], centre
+        )
+        sigma_xy = float(sigma_xy_val.flatten()[0])
+
+        times.append(time)
+        stress_num.append(sigma_xy)
+        stress_ana.append(maxwell_stress_xy(time, ETA, MU, gamma_dot))
+
+    del stokes, mesh
+
+    return np.array(times), np.array(stress_num), np.array(stress_ana)
+
+
+# ---------------------------------------------------------------------------
+# Tests
+# ---------------------------------------------------------------------------
+
+class TestVEShearBox:
+    """Viscoelastic shear box with analytical Maxwell solution."""
+
+    def test_order1_converges(self):
+        """Order-1 VE stress should track the analytical Maxwell curve."""
+
+        times, num, ana = _run_ve_shear(order=1, n_steps=20, dt_over_tr=0.1)
+
+        rel_error_final = abs(num[-1] - ana[-1]) / abs(ana[-1])
+        mask = ana > 0.01 * ana[-1]
+        rel_error_rms = np.sqrt(np.mean(((num[mask] - ana[mask]) / ana[mask]) ** 2))
+
+        print(f"Order 1: final rel error = {rel_error_final:.4e}, "
+              f"RMS rel error = {rel_error_rms:.4e}")
+
+        assert rel_error_final < 0.05
+        assert rel_error_rms < 0.06
+
+    def test_order2_converges(self):
+        """Order-2 VE stress should converge more accurately than order-1."""
+
+        times, num, ana = _run_ve_shear(order=2, n_steps=20, dt_over_tr=0.1)
+
+        rel_error_final = abs(num[-1] - ana[-1]) / abs(ana[-1])
+        mask = ana > 0.01 * ana[-1]
+        rel_error_rms = np.sqrt(np.mean(((num[mask] - ana[mask]) / ana[mask]) ** 2))
+
+        print(f"Order 2: final rel error = {rel_error_final:.4e}, "
+              f"RMS rel error = {rel_error_rms:.4e}")
+
+        # Order-2 should be significantly better than order-1
+        assert rel_error_final < 0.005
+        assert rel_error_rms < 0.02
+
+    def test_order2_better_than_order1(self):
+        """Order 2 should have smaller error than order 1 at the same Δt."""
+
+        _, num1, ana1 = _run_ve_shear(order=1, n_steps=20, dt_over_tr=0.1)
+        _, num2, ana2 = _run_ve_shear(order=2, n_steps=20, dt_over_tr=0.1)
+
+        mask = ana1 > 0.01 * ana1[-1]
+        err1 = abs(num1[-1] - ana1[-1]) / ana1[-1]
+        err2 = abs(num2[-1] - ana2[-1]) / ana2[-1]
+
+        print(f"Final errors: Order 1 = {err1:.4e}, Order 2 = {err2:.4e}")
+
+        assert err2 < err1, (
+            f"Order 2 error ({err2:.4e}) should be less than order 1 ({err1:.4e})"
+        )
+
+    def test_steady_state_approach(self):
+        """Stress should monotonically approach the viscous limit η·γ̇."""
+
+        times, num, ana = _run_ve_shear(order=1, n_steps=30, dt_over_tr=0.1)
+
+        final_err = abs(num[-1] - ana[-1]) / ana[-1]
+
+        print(f"Steady state approach: σ_xy = {num[-1]:.6f}, "
+              f"analytical = {ana[-1]:.6f}")
+
+        diffs = np.diff(num)
+        assert np.all(diffs >= -1.0e-10), "Stress should be monotonically increasing"
+        assert final_err < 0.02