Unify Stokes/VE_Stokes and VEP infrastructure

[lmoresi]

Summary

Eliminates the need for users to choose between Stokes and VE_Stokes solvers.
The Stokes solver now auto-creates stress history infrastructure when assigned
a viscoelastic-plastic constitutive model. Also fixes several VEP correctness
and usability issues.

This is a working baseline — the embedded fault demonstration is the next step.

Changes

Solver unification

• Stokes solver auto-creates DFDt when requires_stress_history=True
• _create_stress_history_ddt() extracted from VE_Stokes
• solve() includes VE pre/post hooks when DFDt is active
• tau property returns psi_star[0] when DFDt exists
• VE_Stokes retained as thin backward-compat wrapper
• Constitutive model order preserved (not overwritten when solver _order=0)

dt_elastic elimination

• No longer a user-facing parameter — solver sets it from solve(timestep=dt)
• Prevents the silent NaN bug when dt_elastic defaulted to infinity
• Parameters.dt_elastic is now a property that reads the solver-managed _dt

Constitutive model fixes

• requires_stress_history property on base class (False) and VEP (True)
• is_viscoplastic comparison fix (.sym is sympy.oo not == sympy.oo)
• plastic_fraction property: max(0, 1 - eta_vep/eta_ve) for VEP,
  max(0, 1 - eta_vp/eta_viscous) for VP, 0 for viscous
• BDF order fallback when timestep ratio > 2x (prevents negative stress)
• Default order=2 for VEP (BDF-2, the validated recommendation)

Performance

• simplify=False default in uw.function.evaluate

Validated benchmarks

• Pure VE loading: stress matches Maxwell analytical to ~96%
• VEP uniform yield: elastic loading → yield → stress clips at tau_y
• Time-dependent yield: elastic loading → yield → stress tracks weakening
  tau_y → elastic rebound → re-yield at residual
• Strain weakening: elastic loading → yield → plastic strain accumulates
  → tau_y weakens → stress drops (full stick-slip cycle)

Design note

docs/developer/design/SOLVER_UNIFICATION_DESIGN.md documents the architecture
and the path to VE-NavierStokes (future work).

Test plan

• VE shear box (test_1051) passes
• Stokes solver tests (test_1010) pass
• VEP diagnostic: VE, VP, VEP side-by-side comparison
• Time-dependent yield with elastic rebound
• Strain weakening with plastic fraction feedback
• CI pipeline
• Embedded fault demonstration (next step)

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

Unifies Stokes and VE_Stokes behavior by making stress-history infrastructure (DFDt) auto-managed based on the constitutive model, while also improving VEP timestep handling and performance defaults.

Changes:

• Auto-create/manage DFDt stress history when the constitutive model requires it; integrate VE/VEP pre/post hooks into Stokes.solve() and expose tau consistently.
• Replace user-facing dt_elastic with a solver-managed timestep property and add BDF safety fallback for large timestep jumps.
• Default uw.function.evaluate(..., simplify=...) to simplify=False and add new VEP benchmark scripts.

Reviewed changes

Copilot reviewed 6 out of 8 changed files in this pull request and generated 8 comments.

Show a summary per file

File	Description
tests/vep_timedep_yield.py	Adds a standalone benchmark script for time-dependent yield stress behavior.
tests/vep_strain_weakening.py	Adds a standalone benchmark script for strain-weakening / plastic-fraction behavior.
src/underworld3/systems/solvers.py	Adds stress-history creation helper and integrates VE/VEP pre/post logic into Stokes solve + tau accessor.
src/underworld3/function/functions_unit_system.py	Changes default evaluation simplification to improve performance.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Preserves constitutive model order when solver order is unset; lazily creates DFDt based on model requirement.
src/underworld3/constitutive_models.py	Adds requires_stress_history / plastic_fraction, changes VEP default order, and makes dt_elastic solver-managed.

Comments suppressed due to low confidence (1)

src/underworld3/function/functions_unit_system.py:1

• Changing the default simplify from True to False is an API behavior change that may affect downstream expectations (symbolic form, performance vs. expression reduction, and potentially numerical stability for some expressions). Consider documenting this change prominently (docstring/release notes), or preserving backward compatibility by keeping the old default and introducing a separate performance-oriented entry point/config flag.

"""

---

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

[Copilot]

The output path is hard-coded to a developer-specific absolute location, which will fail on other machines and in CI environments. Use a repo-relative path (e.g., alongside the script), a configurable CLI argument/env var, or write to a temporary directory and print the resolved path.

[Copilot]

This hard-coded absolute output path is user-specific and will break for other developers and CI runs. Prefer a relative path (e.g., Path(__file__).with_suffix('.png')), or make the output path configurable (argument/env var) and default to a safe writable location.

[Copilot]

This mutates the underlying sympy value directly and bypasses the parameter interface used elsewhere (e.g., cm.Parameters.yield_stress = ...). Use the public setter to keep behavior consistent (and to ensure any side effects/validation are applied uniformly).

Suggested change

	cm.Parameters.yield_stress.sym = current_tau_y
	cm.Parameters.yield_stress = current_tau_y

[Copilot]

The order argument is clamped but never applied afterward (it doesn’t update the solver/DFDt order nor the constitutive model order), so callers cannot actually override the VE integration order even though the API suggests they can. Either (a) remove order from the signature, or (b) plumb it through by updating the relevant state (e.g., solver/model/DFDt order) and triggering the necessary re-setup when it changes.

[Copilot]

The order argument is clamped but never applied afterward (it doesn’t update the solver/DFDt order nor the constitutive model order), so callers cannot actually override the VE integration order even though the API suggests they can. Either (a) remove order from the signature, or (b) plumb it through by updating the relevant state (e.g., solver/model/DFDt order) and triggering the necessary re-setup when it changes.

Suggested change

	if order is None or order > self._order:
	order = self._order
	# Clamp requested VE history order to the solver's maximum order.
	if order is None or order > self._order:
	order = self._order
	# If the requested order differs from the current DFDt / constitutive
	# model configuration, update it and force a re-setup.
	if order is not None:
	# Update DFDt order if available.
	if hasattr(self, "DFDt") and hasattr(self.DFDt, "order"):
	if self.DFDt.order != order:
	self.DFDt.order = order
	self.is_setup = False
	if hasattr(self.constitutive_model, "_solver_is_setup"):
	self.constitutive_model._solver_is_setup = False
	# Propagate order to constitutive model if it has an 'order' attribute.
	if hasattr(self.constitutive_model, "order"):
	if self.constitutive_model.order != order:
	self.constitutive_model.order = order
	self.is_setup = False
	if hasattr(self.constitutive_model, "_solver_is_setup"):
	self.constitutive_model._solver_is_setup = False

[Copilot]

SolverBaseClass unconditionally calls _create_stress_history_ddt(...) when requires_stress_history=True. If a user assigns a stress-history model to a solver that doesn’t implement that method, this will raise an AttributeError at runtime. Add an explicit guard (e.g., hasattr(self, '_create_stress_history_ddt')) and raise a clear ValueError saying this solver type doesn’t support stress history, or move the lazy-creation logic into the solver(s) that actually support it.

Suggested change

	if self._constitutive_model.requires_stress_history and self.Unknowns.DFDt is None:
	if self._constitutive_model.requires_stress_history and self.Unknowns.DFDt is None:
	if not hasattr(self, "_create_stress_history_ddt"):
	raise ValueError(
	"This solver type does not support stress history: "
	"_create_stress_history_ddt is not implemented."
	)

[Copilot]

This docstring is now inconsistent with the new solver unification behavior (the Stokes solver auto-creates DFDt when required). Update the wording to reflect the lazy creation pathway (and/or the actual failure mode if a non-supporting solver is used).

Suggested change

	Models that return True require a solver with stress history
	management (e.g. VE_Stokes). Assigning such a model to a plain
	Stokes solver will raise an error.
	Models that return True require a solver with stress-history
	management. For the unified Stokes solver this causes the DFDt
	history field to be created and managed lazily when the model is
	attached. Solvers that do not support stress-history management
	may either ignore this flag or raise an error; see the specific
	solver documentation for details.

[Copilot]

Importing NumPy inside solve() adds repeated overhead in a hot path. Prefer a module-level import (if not already present) or use the array’s .copy() method directly to avoid the local import.

Suggested change

	import numpy as np
	_advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...])
	_advected_sigma_star = self.DFDt.psi_star[0].array[...].copy()