Auto-detect bdf_blend: 1.0 for pure VE, 0.75 for VEP

The previous default (0.5) degraded order-2 VE accuracy. Now:
- bdf_blend defaults to None (auto-detect)
- Pure VE (yield_stress = oo): uses 1.0 (pure BDF-k, full accuracy)
- VEP (finite yield_stress): uses 0.75 (stable, near-optimal accuracy)
- Users can still override explicitly with stokes.constitutive_model.bdf_blend = value

Also fixes _update_bdf_coefficients to read the property (self.bdf_blend)
instead of the backing field (self._bdf_blend) so auto-detect works.

Test test_1051 no longer needs explicit bdf_blend — auto-detect gives 1.0
for its pure VE setup.

Underworld development team with AI support from Claude Code

Author: lmoresi

src/underworld3/constitutive_models.py

@@ -1112,7 +1112,7 @@ def __init__(self, unknowns, order=1, material_name: str = None):
         self._order = order
         self._yield_mode = "smooth"  # "min", "harmonic", "smooth", or "softmin"
         self._yield_softness = 0.5  # δ parameter for "softmin" mode
-        self._bdf_blend = 0.5  # blend O1/O2 coefficients: 0=pure O1, 1=pure O2
+        self._bdf_blend = None  # auto: 1.0 for VE, 0.75 for VEP
 
         # Timestep — set by the solver before each solve(). Not a user parameter.
         # Initialised to oo (viscous limit). The solver overwrites this with the
@@ -1363,7 +1363,7 @@ def _update_bdf_coefficients(self):
 
             # Blend with O1 coefficients for stability
             # 0 = pure O1, 0.5 = balanced (default), 1 = pure requested order
-            alpha = self._bdf_blend
+            alpha = self.bdf_blend  # property resolves None → auto-detect
             if 0 < alpha < 1 and order >= 2:
                 coeffs_o1 = _bdf_coefficients(1, dt_current, dt_history)
                 while len(coeffs_o1) < len(coeffs):
@@ -1761,9 +1761,12 @@ def bdf_blend(self):
         Blends O1 and O2 BDF coefficients: ``c = (1-α)·c_O1 + α·c_O2``.
 
         - ``α = 0``: pure BDF-1 (most stable, first-order accurate)
-        - ``α = 0.5`` (default): balanced blend (stable, improved accuracy)
-        - ``α = 1``: pure BDF-2 (second-order, can be unstable for VEP)
+        - ``α = 0.75``: default for VEP (stable, near-optimal accuracy)
+        - ``α = 1``: pure BDF-2 (default for pure VE, second-order accurate)
+        - ``None`` (default): auto-detect — 1.0 for VE, 0.75 for VEP
         """
+        if self._bdf_blend is None:
+            return 0.75 if self.is_viscoplastic else 1.0
         return self._bdf_blend
 
     @bdf_blend.setter

tests/test_1051_VE_shear_box.py

@@ -77,7 +77,7 @@ def _run_ve_shear(order, n_steps, dt_over_tr):
     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.bdf_blend = 1.0  # pure BDF-k (no O1/O2 blending)
+    # bdf_blend auto-detects: 1.0 for pure VE, 0.75 for VEP
 
     stokes.add_dirichlet_bc((V0, 0.0), "Top")
     stokes.add_dirichlet_bc((-V0, 0.0), "Bottom")