Explicit stress history for viscoelastic Stokes solver

docs/advanced/benchmarks/index.md

@@ -0,0 +1,13 @@
+---
+title: "Benchmarks"
+---
+
+# Solver Benchmarks
+
+Validation benchmarks comparing Underworld3 solvers against analytical solutions.
+
+```{toctree}
+:maxdepth: 1
+
+ve-oscillatory-shear
+```

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.

src/underworld3/constitutive_models.py

@@ -971,12 +971,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
@@ -1198,32 +1198,19 @@ 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
-                    )
-                )
-
-            # 2nd Order version (need to ask for this one)
+            # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt)
+            # c0: 1 (BDF-1), 3/2 (BDF-2), 11/6 (BDF-3)
+            eta = inner_self.shear_viscosity_0
+            mu = inner_self.shear_modulus
+            dt_e = inner_self.dt_elastic
+            eff_order = inner_self._owning_model.effective_order
+
+            if eff_order == 2:
+                el_eff_visc = 2 * eta * mu * dt_e / (3 * eta + 2 * mu * dt_e)
+            elif eff_order >= 3:
+                el_eff_visc = 6 * eta * mu * dt_e / (11 * eta + 6 * mu * dt_e)
             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 / (eta + mu * dt_e)
 
             inner_self._ve_effective_viscosity.sym = el_eff_visc
 
@@ -1252,6 +1239,19 @@ 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
+
     # The following should have no setters
     @property
     def stress_star(self):
@@ -1287,20 +1287,23 @@ 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
+                eff_order = self.effective_order
 
-                else:
+                if eff_order == 2:
                     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
-                    )
+                    E += stress_star / mu_dt - stress_2star / (4 * mu_dt)
+                elif eff_order >= 3:
+                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
+                    stress_2star = self.Unknowns.DFDt.psi_star[1].sym
+                    stress_3star = self.Unknowns.DFDt.psi_star[2].sym
+                    E += (3 * stress_star / (2 * mu_dt)
+                          - 3 * stress_2star / (4 * mu_dt)
+                          + stress_3star / (6 * mu_dt))
+                else:
+                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
+                    E += stress_star / (2 * mu_dt)
 
         self._E_eff.sym = E
 
@@ -1365,16 +1368,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 +1465,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 +1487,6 @@ def stress_projection(self):
                     / (self.Parameters.dt_elastic * self.Parameters.shear_modulus)
                 )
 
-        stress = sympy.simplify(stress)
-
         return stress
 
     def stress(self):
@@ -1508,33 +1501,29 @@ def stress(self):
         if self.Unknowns.DFDt is not None:
 
             if self.is_elastic:
-                if self.order != 2:
+                mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus
+                eff_order = self.effective_order
+
+                if eff_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)
-                        )
+                    stress_2star = self.Unknowns.DFDt.psi_star[1].sym
+                    stress += 2 * self.viscosity * (
+                        stress_star / mu_dt - stress_2star / (4 * mu_dt)
                     )
-
-                else:
+                elif eff_order >= 3:
                     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)
-                        )
+                    stress_3star = self.Unknowns.DFDt.psi_star[2].sym
+                    stress += 2 * self.viscosity * (
+                        3 * stress_star / (2 * mu_dt)
+                        - 3 * stress_2star / (4 * mu_dt)
+                        + stress_3star / (6 * mu_dt)
+                    )
+                else:
+                    stress_star = self.Unknowns.DFDt.psi_star[0].sym
+                    stress += 2 * self.viscosity * (
+                        stress_star / (2 * mu_dt)
                     )
-
-        stress = sympy.simplify(stress)
 
         return stress
 

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

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