Explicit stress history for viscoelastic Stokes solver

[lmoresi]

Summary

Explicit stress history architecture for VE_Stokes: store actual computed stress
in psi_star (rather than flux), advect stored values along characteristics, and
use BDF-1/2/3 discretisation of the Maxwell constitutive law.

Key changes

• Explicit stress history: after each solve, the actual deviatoric stress is
  projected into psi_star[0] and shifted down the history chain. This replaces
  the previous approach of storing/advecting the flux expression.
• solver.tau property: access the computed stress tensor (projected, not the
  symbolic formula)
• BDF-3 constitutive formulae: ve_effective_viscosity, E_eff, and stress()
  now have three-way branching for BDF-1/2/3 with correct coefficients
• effective_order property: ramps from 1 to requested order as DDt history
  accumulates, preventing use of higher-order coefficients before distinct
  history is available
• dt_elastic is independent of timestep: solve(timestep=) no longer
  overwrites the constitutive relaxation parameter
• Analytical formula fix: removed spurious De factor in Maxwell oscillatory
  solution (was invisible at De=1, caused apparent errors at other De)
• Benchmark documentation: oscillatory shear validation with convergence
  tables and resolution study

Convergence (constant shear, Maxwell analytical)

dt/t_r	BDF-1	BDF-2	BDF-3
0.200	3.0e-2	4.0e-3	6.4e-3
0.100	1.5e-2	9.3e-4	1.7e-3
0.050	7.8e-3	2.3e-4	4.3e-4
0.020	3.1e-3	3.7e-5	--

BDF-2 is the recommended default (second-order convergence, 85x better than
BDF-1 at small dt).

Oscillatory validation

A shear layer with a harmonic (in time) velocity boundary condition. For a Maxwell material, the analytic solution can be derived - it shows a phase lag in stress compared to the driving velocity, and the amplitude evolves with time.

The Maxwell constitutive law for simple shear:

$$\dot{\sigma}{xy} + \frac{\sigma{xy}}{t_r} = \mu \dot{\gamma}_0 \sin(\omega t)$$

where $t_r = \eta/\mu$ is the relaxation time and $\text{De} = \omega t_r$ is the Deborah number.

Full solution with $\sigma(0) = 0$:

$$\sigma_{xy}(t) = \frac{\eta \dot{\gamma}_0}{1 + \text{De}^2}\left[\sin(\omega t) - \text{De} \cdot \cos(\omega t) + \text{De},e^{-t/t_r}\right]$$

Steady-state amplitude: $A = \eta \dot{\gamma}_0 / \sqrt{1 + \text{De}^2}$

Phase lag: $\delta = \arctan(\text{De})$

Notes on dt_elastic separation

A running-average approach for history accumulation when dt << dt_elastic was
investigated but found to be extremely diffusive for semi-Lagrangian transport
and is not implemented. To prevent unstable behaviour when timesteps become
small, we advise limiting the minimum effective viscosity in line with the
physics of the problem.

Test plan

• Constant shear benchmark (order 1, 2, 3) against Maxwell analytical
• Oscillatory shear validation (De=1, 1.5, 5) with phase lag and amplitude
• Resolution study (De=5, three timestep sizes)
• VP and VEP shear box tests
• CI pipeline

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR refactors the viscoelastic Stokes (VE_Stokes) stress-history handling to store and advect the computed deviatoric stress explicitly, adds BDF-1/2/3 branching for Maxwell constitutive terms, and introduces a symbolic mesh time coordinate (mesh.t) wired through PETSc’s petsc_t plumbing.

Changes:

• Replace “advect flux expression” with explicit stress history storage in psi_star[0], plus history shifting after each solve.
• Add solver.tau as a projected/computed stress (base solver projection) and override for VE_Stokes to return stored history directly.
• Add mesh.t / petsc_t support and benchmark/validation scripts & docs for oscillatory shear.

Reviewed changes

Copilot reviewed 25 out of 25 changed files in this pull request and generated 16 comments.

Show a summary per file

File	Description
uw	Build-script changes to aggressively clean build artifacts and attempt to detect/correct stale installed .py files.
tests/test_1051_VE_shear_box.py	New pytest validation against Maxwell constant-shear analytical solution using stokes.tau.
tests/run_vp_shear_box.py	New manual VP shear-box validation script using stokes.tau.
tests/run_vep_shear_box.py	New manual VEP shear-box validation script.
tests/run_vep_oscillatory.py	New manual oscillatory VEP validation script (includes analytical helper).
tests/run_ve_vep_oscillatory_plot.py	New VE vs VEP oscillatory comparison/plot script with analytical init.
tests/run_ve_vep_oscillatory_checkpoint.py	New checkpoint writer for VE vs VEP oscillatory runs.
tests/run_ve_shear_validation.py	New VE shear-box validation helper script.
tests/run_ve_shear_quick.py	New quick VE order-1 shear script.
tests/run_ve_shear_order2_quick.py	New quick VE order-2 shear script.
tests/run_ve_oscillatory.py	New VE oscillatory validation script.
tests/run_ve_order2_debug.py	New debugging script to inspect psi_star evolution.
tests/plot_ve_vep_oscillatory.py	New plotting script for VE/VEP oscillatory checkpoints.
tests/plot_ve_oscillatory_validation.py	New benchmark plotting script and checkpointing for oscillatory VE.
src/underworld3/utilities/unit_aware_coordinates.py	Add TimeSymbol and time-unit patching for unit-aware mesh.t.
src/underworld3/systems/solvers.py	VE_Stokes solve sequence refactor: advect_history() + post-solve projection into psi_star[0] + history shift; add VE override tau.
src/underworld3/systems/ddt.py	Change effective-order ramp logic and add SemiLagrangian.advect_history() for pure-history advection.
src/underworld3/function/_function.pyx	Add TODO note about matrix projection limitations.
src/underworld3/discretisation/discretisation_mesh.py	Add mesh.t property and initialize _t as TimeSymbol.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Add base SolverBaseClass.tau projection property and plumb time= into SNES solves via UW_DMSetTime.
src/underworld3/cython/petsc_extras.pxi	Add Cython binding for UW_DMSetTime.
src/underworld3/cython/petsc_compat.h	Implement UW_DMSetTime() via DM(Set/Get)OutputSequenceNumber.
src/underworld3/constitutive_models.py	Implement BDF-1/2/3 branching for VE effective viscosity / E_eff / stress; introduce constitutive effective_order.
docs/advanced/benchmarks/ve-oscillatory-shear.md	New benchmark write-up for oscillatory Maxwell shear + timestep study + dt_elastic note.
docs/advanced/benchmarks/index.md	Add benchmarks index and toctree entry.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Same effective_order off-by-one issue here as in the symbolic DDt: using max(1, _n_solves_completed) delays BDF-2/3 usage by one extra timestep even when distinct history values exist. This can degrade the advertised order of accuracy during startup.

Suggested change

	# 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))
	# Distinct history values = initial condition + completed solves.
	# With 0 completed solves → 1 value (order 1); with 1 solve → 2 values (order 2).
	return min(self.order, max(1, self._n_solves_completed + 1))

[Copilot]

This always copies the full psi_star[0] array even when self.DFDt.order == 1 (where no history shift is performed). For large meshes this extra allocation/copy can be expensive; consider guarding the copy/shift logic with if self.DFDt.order > 1: and only snapshotting when it’s actually needed.

[Copilot]

Copying differing .py sources directly into site-packages can hide packaging/install issues and create an environment that doesn’t correspond to any wheel/sdist artifact (harder to reproduce/debug). If the goal is always-up-to-date sources for dev, consider using an editable install (pip install -e .) or bumping the package version/build tag so pip’s wheel caching can’t serve stale files, rather than mutating the installed tree post-install.

Suggested change

	echo " Fixing: copying source files to installed location..."
	for f in "${stale_files[@]}"; do
	cp "src/underworld3/$f" "$site_packages/$f"
	done
	echo -e " ${GREEN}✓${NC} ${#stale_files[@]} file(s) updated"
	echo " Installed package appears out of sync with the source tree."
	echo " Please reinstall the 'underworld3' package to refresh installed files."
	echo " For development, consider using an editable install: 'pip install -e .'."
	echo " Alternatively, bump the package version so pip does not reuse a stale wheel."
	echo -e " ${YELLOW}No files were modified automatically by this script.${NC}"

[Copilot]

The docstring’s analytical solution still shows a De factor in the prefactor (η·γ̇0·De/(1+De²)), but maxwell_oscillatory() below uses the corrected η·γ̇0/(1+De²) form. Updating the docstring keeps the benchmark description consistent with the corrected formula elsewhere in the PR.

Suggested change

	σ_xy(t) = η γ̇₀ De/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)]
	σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)]

[Copilot]

Same effective_order off-by-one issue here: BDF-2 should be available after 1 completed solve (with initial condition as the older history). Using _n_solves_completed directly postpones higher-order coefficients longer than necessary and can reduce accuracy.

Suggested change

	# 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))
	# BDF-k requires k distinct history values.
	# With the initial condition stored as the oldest history entry,
	# we have (_n_solves_completed + 1) history states available.
	return min(self.order, max(1, self._n_solves_completed + 1))

[Copilot]

The docstring and ve_amplitude computation include an extra De factor in the steady-state amplitude. That contradicts the corrected analytical solution used elsewhere in this PR (and the new benchmark doc), where the steady-state amplitude is η·γ̇0 / sqrt(1+De^2). Please update the amplitude formula and the prefactor/return expression accordingly so the validation script isn’t checking against a wrong analytic baseline.

[Copilot]

maxwell_oscillatory() here still uses prefactor = η·γ̇0·De/(1+De²) (and the corresponding De * exp(-t/t_r) term). Elsewhere in this PR the spurious De factor was removed; if this helper is meant to match the corrected Maxwell solution it should use prefactor = η·γ̇0/(1+De²) with the same bracketed terms.

Suggested change

	prefactor = eta * gamma_dot_0 * De / (1.0 + De**2)
	prefactor = eta * gamma_dot_0 / (1.0 + De**2)

[Copilot]

effective_order appears off by one for BDF startup. After 1 completed solve you already have two distinct history states (t0 and t1), so BDF-2 can be used on the second solve. Returning max(1, _n_solves_completed) delays BDF-2 until the third solve and can reduce the expected convergence order; consider reverting to max(1, _n_solves_completed + 1) (and likewise for BDF-3).

Suggested change

	# 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))
	# After n completed solves, we have n+1 distinct time levels (t0..tn),
	# so BDF-k can be used once k ≤ n+1. This corresponds to an effective
	# order of min(self.order, max(1, _n_solves_completed + 1)).
	return min(self.order, max(1, self._n_solves_completed + 1))

[Copilot]

This monotonicity assertion is very strict (diffs > 0) and may be flaky due to solver tolerances / roundoff (e.g. tiny negative differences around 1e-12). Consider allowing a small tolerance (e.g. diffs >= -eps) so the test remains robust across platforms and PETSc/KSP configurations.

Suggested change

	assert np.all(diffs > 0), "Stress should be monotonically increasing"
	eps = 1.0e-10
	assert np.all(diffs >= -eps), "Stress should be monotonically increasing"

[Copilot]

The analytical formula in the module docstring still includes an extra De factor in the prefactor and in the steady-state amplitude. The implementation below uses prefactor = η·γ̇0/(1+De²), so the docstring should be updated to match the corrected Maxwell solution (and avoid confusing users reading this validation script).

Suggested change

	σ_xy(t) = η γ̇₀ De/(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 η γ̇₀ De/√(1+De²) and phase lag arctan(De).
	σ_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).