diff --git a/docs/advanced/figures/ti_vep_benchmark_final.png b/docs/advanced/figures/ti_vep_benchmark_final.png new file mode 100644 index 000000000..f23362660 Binary files /dev/null and b/docs/advanced/figures/ti_vep_benchmark_final.png differ diff --git a/docs/advanced/index.md b/docs/advanced/index.md index c9e9a5e6b..94826f71d 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -23,6 +23,12 @@ Implement advanced material models and constitutive laws. **[→ Complex Rheologies](complex-rheologies.md)** +### VEP with Transverse Isotropy for Fault Mechanics +Viscoelastic-plastic rheology with anisotropic weak planes and resolved +fault-plane yield for modelling fault zones. + +**[→ VEP + Transverse Isotropy](vep-transverse-isotropy-faults.md)** + ### Custom Meshes Create complex geometries with gmsh for research problems. @@ -74,6 +80,7 @@ Ready to contribute to Underworld3? parallel-computing performance complex-rheologies +vep-transverse-isotropy-faults custom-meshes curved-boundary-conditions mesh-adaptation diff --git a/docs/advanced/vep-transverse-isotropy-faults.md b/docs/advanced/vep-transverse-isotropy-faults.md new file mode 100644 index 000000000..8bf3e56a4 --- /dev/null +++ b/docs/advanced/vep-transverse-isotropy-faults.md @@ -0,0 +1,265 @@ +--- +title: "Viscoelastic-Plastic Rheology with Transverse Isotropy for Fault Mechanics" +--- + +# Viscoelastic-Plastic Rheology with Transverse Isotropy + +Fault zones in the lithosphere are thin regions of localised deformation where the mechanical response differs from the surrounding rock. They are weaker in shear along the fault plane than in the bulk, they accumulate elastic stress between slip events, and they yield when that stress exceeds a threshold. Capturing all three behaviours -- anisotropic weakness, elastic memory, and plastic yield -- requires a constitutive model that combines transverse isotropy (TI) with viscoelastic-plastic (VEP) rheology. + +This document develops the mathematical formulation used in the `TransverseIsotropicVEPFlowModel` class, starting from the isotropic VEP model and the TI viscosity tensor, then showing how they combine through a resolved fault-plane yield criterion. + +## Isotropic Viscoelastic-Plastic Rheology + +### Maxwell Viscoelasticity + +A Maxwell viscoelastic material partitions the total strain rate into viscous and elastic contributions: + +$$ +\dot\varepsilon_{ij}^{\text{total}} = \dot\varepsilon_{ij}^{\text{viscous}} + \dot\varepsilon_{ij}^{\text{elastic}} += \frac{\sigma_{ij}}{2\eta} + \frac{1}{2\mu}\frac{D\sigma_{ij}}{Dt} +$$ + +where $\eta$ is the shear viscosity, $\mu$ is the shear modulus, and $D/Dt$ denotes the Jaumann (or other objective) derivative. + +Discretising the time derivative using a BDF-$k$ scheme with leading coefficient $c_0$ and history coefficients $c_1, c_2, \ldots$ gives: + +$$ +\dot\varepsilon_{ij}^{\text{total}} = \frac{\sigma_{ij}}{2\eta} + \frac{c_0 \sigma_{ij} + c_1 \sigma_{ij}^{*} + c_2 \sigma_{ij}^{**} + \cdots}{2\mu\Delta t} +$$ + +where $\sigma^{*}$ and $\sigma^{**}$ are the stress at the previous and second-previous timesteps, advected to the current particle positions (the Lagrangian stress history). Solving for the current stress: + +$$ +\sigma_{ij} = 2\eta_{\text{ve}}\,\dot\varepsilon_{ij}^{\text{eff}} +$$ + +with the **viscoelastic effective viscosity** and **effective strain rate**: + +$$ +\eta_{\text{ve}} = \frac{\eta\,\mu\,\Delta t}{c_0\,\eta + \mu\,\Delta t}, +\qquad +\dot\varepsilon_{ij}^{\text{eff}} = \dot\varepsilon_{ij}^{\text{total}} +- \frac{c_1 \sigma_{ij}^{*} + c_2 \sigma_{ij}^{**} + \cdots}{2\mu\,\Delta t} +$$ + +The effective strain rate incorporates stress history: it is the strain rate that a purely viscous material with viscosity $\eta_{\text{ve}}$ would need to produce the current stress. For BDF-1, $c_0 = 1$ and $c_1 = -1$; higher orders improve temporal accuracy. + +The Maxwell relaxation time $t_r = \eta / \mu$ controls the elastic-to-viscous transition. When $\Delta t \gg t_r$, the material behaves viscously ($\eta_{\text{ve}} \to \eta$); when $\Delta t \ll t_r$, it behaves elastically ($\eta_{\text{ve}} \to \mu\Delta t / c_0$). + +### Plastic Yield + +When stress exceeds a yield threshold $\tau_y$, the material yields plastically. In the isotropic case, yield is tested against the second invariant of the effective strain rate: + +$$ +\dot\varepsilon_{II} = \sqrt{\tfrac{1}{2}\dot\varepsilon_{ij}^{\text{eff}}\,\dot\varepsilon_{ij}^{\text{eff}}} +$$ + +The effective viscosity is capped so that the resulting stress does not exceed the yield stress: + +$$ +\eta_{\text{vep}} = \min\!\left(\eta_{\text{ve}},\;\frac{\tau_y}{2\,\dot\varepsilon_{II}}\right) +$$ + +This is the Drucker-Prager yield criterion expressed as a viscosity cap. The stress is then $\sigma_{ij} = 2\eta_{\text{vep}}\,\dot\varepsilon_{ij}^{\text{eff}}$. + +## Transverse Isotropy + +### The Muhlhaus-Moresi Viscosity Tensor + +A transversely isotropic material has a single weak plane defined by its unit normal $\hat{n}$ (the "director"). The fourth-rank viscosity tensor is: + +$$ +\eta_{ijkl} = 2\eta_0\,I_{ijkl} +- (\eta_0 - \eta_1)\left[ +\frac{1}{2}\!\left(n_i n_k \delta_{jl} + 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] +$$ + +where $I_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$ is the symmetric identity tensor, $\eta_0$ is the bulk viscosity, and $\eta_1$ is the viscosity for shear along the weak plane. When $\eta_1 = \eta_0$, the anisotropic correction vanishes and the tensor reduces to the isotropic case. + +The stress is: + +$$ +\sigma_{ij} = \eta_{ijkl}\,\dot\varepsilon_{kl} +$$ + +The director $\hat{n}$ typically comes from the fault surface normals, transferred to the mesh via nearest-neighbour interpolation. Far from the fault, $\eta_1$ is set equal to $\eta_0$ (via an influence function), so the material reverts to isotropic. + +### Fault Representation + +In Underworld3, fault zones are represented as embedded surfaces. The `Surface` class provides: + +- A signed distance field from the fault +- Normal vectors at each surface vertex +- Influence functions that smoothly transition material properties from fault-zone values (near the surface) to background values (far away) + +Available influence profiles include Gaussian ($e^{-(d/w)^2}$), smoothstep ($3t^2 - 2t^3$), linear ramp, and step function. The width parameter $w$ controls the fault zone thickness. + +A typical setup transfers fault normals to the mesh as a director field, and uses an influence function to interpolate yield stress or viscosity ratio between fault-zone and background values. + +## Combined TI-VEP: Resolved Fault-Plane Yield + +The key insight in combining TI with VEP is that yield should be tested against the **resolved shear stress on the fault plane**, not the global stress invariant. A fault with normal $\hat{n}$ oriented at an angle to the imposed deformation may have a low global strain rate invariant while experiencing high shear along the fault plane itself. + +### Viscoelastic Effective Viscosities + +Both viscosity parameters receive the VE treatment: + +$$ +\eta_{0,\text{ve}} = \frac{\eta_0\,\mu\,\Delta t}{c_0\,\eta_0 + \mu\,\Delta t}, +\qquad +\eta_{1,\text{ve}} = \frac{\eta_1\,\mu\,\Delta t}{c_0\,\eta_1 + \mu\,\Delta t} +$$ + +Using $\eta_{0,\text{ve}}$ in the tensor (rather than the raw $\eta_0$) ensures that the anisotropic correction $\Delta = \eta_{0,\text{ve}} - \eta_{1,\text{eff}}$ vanishes when $\eta_1 = \eta_0$ and yield is inactive. Without this, the tensor would have a spurious anisotropic component even for isotropic materials. + +### Resolved Shear on the Fault Plane + +Given the effective strain rate tensor $\dot\varepsilon_{ij}^{\text{eff}}$ (which includes the stress history), the traction-like vector on the fault plane is: + +$$ +T_i = \dot\varepsilon_{ij}^{\text{eff}}\,n_j +$$ + +This has a component normal to the fault and a component tangent to it. The normal component is: + +$$ +\dot\varepsilon_n = T_i\,n_i +$$ + +The in-plane (tangential) shear magnitude follows from Pythagoras: + +$$ +|\dot\gamma| = \sqrt{|T|^2 - \dot\varepsilon_n^2} += \sqrt{T_i T_i - (T_j n_j)^2} +$$ + +This formulation works in both 2D and 3D without constructing an explicit tangent vector (which is not unique in 3D). The quantity $|\dot\gamma|$ is the magnitude of the shear strain rate resolved onto the fault plane. + +### Fault-Plane Yield Criterion + +The plastic viscosity is determined by the resolved fault-plane shear: + +$$ +\eta_{1,\text{pl}} = \frac{\tau_y}{2\,|\dot\gamma|} +$$ + +This is the same Drucker-Prager pattern as the isotropic case, but projected onto the fault plane. The yield-limited fault-plane viscosity is: + +$$ +\eta_{1,\text{eff}} = \min\!\left(\eta_{1,\text{ve}},\;\eta_{1,\text{pl}}\right) +$$ + +In practice, a smooth approximation replaces the $\min$ to aid solver convergence. The default "smooth" yield mode uses: + +$$ +\eta_{1,\text{eff}} = \eta_{1,\text{ve}}\,\frac{1 + f}{1 + f + f^2}, +\quad f = \frac{\eta_{1,\text{ve}}}{\eta_{1,\text{pl}}} +$$ + +which transitions smoothly from $\eta_{1,\text{ve}}$ (when $f \ll 1$, below yield) to $\eta_{1,\text{pl}}$ (when $f \gg 1$, above yield). + +### The Full Stress Formula + +The stress is computed from the anisotropic tensor with the yield-limited viscosities: + +$$ +\sigma_{ij} = C_{ijkl}(\eta_{0,\text{ve}},\;\eta_{1,\text{eff}},\;\hat{n}) +\;\dot\varepsilon_{kl}^{\text{eff}} +$$ + +where $C_{ijkl}$ is the Muhlhaus-Moresi tensor with $\eta_0 \to \eta_{0,\text{ve}}$ and $\eta_1 \to \eta_{1,\text{eff}}$. This naturally separates the response into: + +- **Normal to the fault**: governed by $\eta_{0,\text{ve}}$ (pure VE, no yield) +- **Shear along the fault**: governed by $\eta_{1,\text{eff}}$ (VEP, yield-limited) + +When the fault-plane shear stress reaches $\tau_y$, only the fault-plane component yields. Stress normal to the fault continues to build elastically. This is physically correct: faults slip in shear, not in compression. + +## Using `TransverseIsotropicVEPFlowModel` + +```python +import underworld3 as uw +import numpy as np +import sympy + +# Mesh and variables +mesh = uw.meshing.StructuredQuadBox(elementRes=(64, 64)) +v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2, vtype=uw.VarType.VECTOR) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, + continuous=True, vtype=uw.VarType.SCALAR) + +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + +# Create TI-VEP model (order=1 for BDF-1 time integration) +cm = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=1 +) +stokes.constitutive_model = cm + +# Set parameters +cm.Parameters.shear_viscosity_0 = 1.0 # bulk viscosity +cm.Parameters.shear_viscosity_1 = 0.1 # fault-plane viscosity +cm.Parameters.shear_modulus = 1.0 # elastic shear modulus +cm.Parameters.yield_stress = 0.15 # fault-plane yield stress + +# Director from fault normal (e.g., fault at 15 degrees from horizontal) +theta = np.radians(15) +cm.Parameters.director = sympy.Matrix([-np.sin(theta), np.cos(theta)]) +``` + +The director can also be a spatially varying field (e.g., from a `Surface` object's normals transferred to a mesh variable), and the yield stress can vary spatially using an influence function to localise yielding near the fault. + +## Smooth Yield Approximations + +The `"softmin"` yield mode (default) uses a smooth approximation to $\min(\eta_{\text{ve}}, \eta_{\text{pl}})$ to avoid the non-differentiable kink that causes problems for the SNES solver. The approximation is: + +$$g(f) = 1 + \text{softplus}(f-1) - \text{softplus}(-1), \qquad \eta_{\text{eff}} = \eta_{\text{ve}} / g(f)$$ + +where $\text{softplus}(x) = (x + \sqrt{x^2 + \delta^2})/2$ and $f = \eta_{\text{ve}}/\eta_{\text{pl}}$. The offset correction ensures $g(0) = 1$ exactly, so there is no spurious yield correction when the material is below yield. + +The sharpness parameter $\delta$ (default 0.1) controls the width of the smooth transition around the yield point. Smaller $\delta$ gives a sharper cap (closer to the true $\min$) but a stiffer nonlinearity for the solver. + +### Choosing $\delta$ + +The accuracy of the smooth approximation depends on the ratio $f_{ss} = \eta_{\text{ve}} / \eta_{\text{pl}}$ at steady state. For a simple shear problem, this simplifies to: + +$$f_{ss} = \frac{\sigma_{\text{viscous}}}{\tau_y} = \frac{\eta\,\dot\gamma}{\tau_y}$$ + +The softmin is accurate when $\delta \ll f_{ss}$, i.e., when the viscous stress substantially exceeds the yield stress. Practical guidance: + +| $\delta$ | Accuracy at $f_{ss} = 1.5$ | Accuracy at $f_{ss} = 3$ | Solver cost | +|----------|---------------------------|--------------------------|-------------| +| 0.5 | ~85% of $\tau_y$ | ~99% | lowest | +| 0.1 | ~99% of $\tau_y$ | ~100% | low | +| 0.01 | ~100% | ~100% | moderate | + +The default $\delta = 0.1$ is accurate for all cases where the viscous stress exceeds the yield stress by at least 50% ($f_{ss} > 1.5$). For problems where SNES convergence is difficult at yield onset, increase $\delta$ toward 0.3--0.5 as a relaxation parameter. Set it via `cm.yield_softness = 0.1`. + +## Benchmark Results + +The figure below shows the TI-VEP model validated against the analytical Maxwell viscoelastic solution with plastic yield cap, for a simple shear box with an embedded fault. Two yield stresses are tested ($\tau_y = 0.15$ and $\tau_y = 0.30$) at both 0 and 15 degrees fault angle. Solid curves show the analytical VE solution capped at $\tau_y$; markers show the numerical results. + +```{figure} figures/ti_vep_benchmark_final.png +:name: fig-tivep-benchmark + +TI-VEP shear box benchmark. **Left**: horizontal fault ($\theta = 0°$), where resolved shear equals $\sigma_{xy}$. **Right**: angled fault ($\theta = 15°$), showing resolved fault-plane shear (circles) capping at $\tau_y$ while the global $\sigma_{xy}$ (crosses) continues to build as the bulk VE component grows. With the corrected softmin ($\delta = 0.1$), all cases reach within 1--2% of the analytical yield cap. +``` + +At 0 degrees, the resolved shear is simply $\sigma_{xy}$ and the yield cap is exact. At 15 degrees, the anisotropic tensor creates a mechanical coupling between normal and shear components on the fault plane: the resolved shear caps at $\tau_y$ while the global stress tensor reflects contributions from both the yielded fault-plane component (governed by $\eta_{1,\text{eff}}$) and the non-yielding bulk component (governed by $\eta_{0,\text{ve}}$). + +## Summary of Constitutive Models + +| Model | Viscosity | Elasticity | Yield | Anisotropy | +|-------|-----------|------------|-------|------------| +| `ViscousFlowModel` | $\eta$ | -- | -- | -- | +| `ViscoPlasticFlowModel` | $\eta$ | -- | $\dot\varepsilon_{II}$ | -- | +| `ViscoElasticPlasticFlowModel` | $\eta$ | $\mu$, BDF-$k$ | $\dot\varepsilon_{II}$ | -- | +| `TransverseIsotropicFlowModel` | $\eta_0, \eta_1, \hat{n}$ | -- | -- | TI tensor | +| `TransverseIsotropicVEPFlowModel` | $\eta_0, \eta_1, \hat{n}$ | $\mu$, BDF-$k$ | $\|\dot\gamma\|$ (fault-plane) | TI tensor | + +## References + +- Moresi, L., Muhlhaus, H.-B., 2006. Anisotropic viscous models of large-deformation Mohr-Coulomb failure. *Phil. Mag.*, 86, 3287-3305. +- Muhlhaus, H.-B., Moresi, L., Hobbs, B., Dufour, F., 2002. Large amplitude folding in finely layered viscoelastic rock structures. *Pure Appl. Geophys.*, 159, 2311-2333. diff --git a/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py b/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py new file mode 100644 index 000000000..4c045ba0f --- /dev/null +++ b/docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py @@ -0,0 +1,330 @@ +# %% [markdown] +r""" +# Viscoelastic-Plastic Shear Box with Angled Fault + +**PHYSICS:** solid_mechanics +**DIFFICULTY:** intermediate +**RUNTIME:** ~2 minutes + +## Description + +A 2D shear box with an embedded fault at 15 degrees from horizontal, using the +`TransverseIsotropicVEPFlowModel` constitutive model. This combines: + +- **Transverse isotropy**: anisotropic viscosity with the fault normal as director +- **Viscoelasticity**: Maxwell stress buildup with BDF-1 time integration +- **Plastic yield**: resolved fault-plane shear limits the stress + +The key physical result: stress builds elastically until the **resolved shear +stress on the fault plane** reaches the yield stress $\tau_y$. At that point, +fault-plane shear yields while the stress component normal to the fault +continues to build as pure viscoelastic. + +## Physical Setup + +| Parameter | Value | +|-----------|-------| +| Domain | 1 x 1 | +| Fault | centred, 15 deg from horizontal | +| $\eta_0$ (bulk) | 1 | +| $\eta_1$ (fault-plane) | 1 | +| $\mu$ (shear modulus) | 1 | +| $\tau_y$ (fault yield) | 0.15 | +| Fault width | 0.08 | +| Top velocity | 0.5 | +""" + +# %% +#| echo: false +import nest_asyncio +nest_asyncio.apply() + +# %% +import os +import numpy as np +import sympy +import underworld3 as uw + +import matplotlib +if not os.environ.get('DISPLAY') and not os.environ.get('WAYLAND_DISPLAY'): + matplotlib.use('Agg') +import matplotlib.pyplot as plt + +os.makedirs("output", exist_ok=True) + +# %% [markdown] +""" +## Parameters +""" + +# %% +# Physical parameters +ETA_0 = 1.0 # bulk viscosity +ETA_1 = 1.0 # fault-plane viscosity (same as bulk for this test) +MU = 1.0 # shear modulus +TAU_Y = 0.15 # fault-plane yield stress +V_TOP = 0.5 # top boundary velocity +DT = 0.025 # timestep +N_STEPS = 80 # number of steps + +# Mesh +RES = 64 # mesh resolution (RES x RES) + +# Fault geometry +FAULT_ANGLE_DEG = 15.0 # angle from horizontal (degrees) +FAULT_WIDTH = 0.08 # influence function width +FAULT_LENGTH = 0.6 # fault length (centered in domain) + +# %% [markdown] +r""" +## Mesh and Variables +""" + +# %% +mesh = uw.meshing.StructuredQuadBox( + elementRes=(RES, RES), + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + qdegree=3, +) + +v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2, vtype=uw.VarType.VECTOR) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, + continuous=True, vtype=uw.VarType.SCALAR) + +# %% [markdown] +r""" +## Fault Surface + +The fault is a 1D polyline at 15 degrees from horizontal, centered in the domain. +The `Surface` class computes the signed distance field and normals automatically. +""" + +# %% +theta = np.radians(FAULT_ANGLE_DEG) +cx, cy = 0.5, 0.5 # centre of domain + +# Fault endpoints +dx = FAULT_LENGTH / 2 * np.cos(theta) +dy = FAULT_LENGTH / 2 * np.sin(theta) +fault_points = np.array([ + [cx - dx, cy - dy], + [cx + dx, cy + dy], +]) + +fault = uw.meshing.Surface("fault", mesh, fault_points, symbol="F") +fault.discretize() + +# Director: fault normal (perpendicular to fault, pointing "up") +n_x = -np.sin(theta) +n_y = np.cos(theta) +director = sympy.Matrix([n_x, n_y]) + +print(f"Fault angle: {FAULT_ANGLE_DEG} deg") +print(f"Director (fault normal): [{n_x:.4f}, {n_y:.4f}]") + +# %% [markdown] +r""" +## Yield Stress Field + +The yield stress varies spatially: low near the fault, high in the bulk. +We interpolate the weakness (1/$\tau_y$) to avoid steep gradients. +""" + +# %% +TAU_Y_BULK = 200.0 # effectively infinite for the bulk + +weakness = fault.influence_function( + width=FAULT_WIDTH, + value_near=1 / TAU_Y, + value_far=1 / TAU_Y_BULK, + profile="gaussian", +) +tau_y_field = 1 / weakness + +# %% [markdown] +r""" +## Solver Setup + +The `TransverseIsotropicVEPFlowModel` combines the anisotropic viscosity tensor +(Muhlhaus-Moresi) with viscoelastic stress history and plastic yield on the +fault plane. +""" + +# %% +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + +# Create model with BDF-1 time integration +cm = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=1 +) +stokes.constitutive_model = cm + +# Set parameters +cm.Parameters.shear_viscosity_0 = ETA_0 +cm.Parameters.shear_viscosity_1 = ETA_1 +cm.Parameters.shear_modulus = MU +cm.Parameters.yield_stress = tau_y_field +cm.Parameters.director = director +cm.Parameters.shear_viscosity_min = ETA_0 * 1.0e-3 +cm.Parameters.strainrate_inv_II_min = 1.0e-6 +cm.yield_mode = "softmin" # smooth approximation to min (default delta=0.1) + +# Solver settings +stokes.saddle_preconditioner = 1 / cm.K +stokes.tolerance = 1.0e-4 +stokes.petsc_options["ksp_type"] = "fgmres" + +# Boundary conditions: simple shear +stokes.add_essential_bc(sympy.Matrix([V_TOP, 0.0]), "Top") +stokes.add_essential_bc(sympy.Matrix([0.0, 0.0]), "Bottom") +stokes.add_essential_bc((sympy.oo, 0.0), "Left") +stokes.add_essential_bc((sympy.oo, 0.0), "Right") +stokes.bodyforce = sympy.Matrix([0.0, 0.0]) + +# %% [markdown] +r""" +## Time Stepping + +Track the stress components at a point on the fault to show elastic buildup +and yield cap behavior. +""" + +# %% +# Monitoring point: centre of the fault +monitor_coord = np.array([[0.5, 0.5]]) + +times = [] +sigma_xy_history = [] +sigma_resolved_history = [] + +# Analytical VE reference (no yield): sigma_xy = eta * gamma_dot * (1 - exp(-t / t_r)) +gamma_dot = V_TOP # approximate global shear rate +t_relax = ETA_1 / MU + +for step in range(N_STEPS): + stokes.solve(timestep=DT, zero_init_guess=(step == 0)) + + t = (step + 1) * DT + reason = stokes.snes.getConvergedReason() + its = stokes.snes.getIterationNumber() + + # Sample stress at monitoring point + tau = stokes.tau + tau_data = tau.data + tau_coords = tau.coords + + # Find nearest stress evaluation point to monitor location + dists = np.linalg.norm(tau_coords - monitor_coord, axis=1) + idx = np.argmin(dists) + s_xx, s_yy, s_xy = tau_data[idx, 0], tau_data[idx, 1], tau_data[idx, 2] + + # Resolved shear on fault plane: tau_resolved = n^T sigma n_perp + # For a fault with normal (n_x, n_y), the tangent is (n_y, -n_x) + # Resolved shear = t^T sigma n = sigma_ij * t_i * n_j + t_x, t_y = n_y, -n_x # tangent vector + resolved_shear = (s_xx * t_x * n_x + s_xy * (t_x * n_y + t_y * n_x) + + s_yy * t_y * n_y) + + times.append(t) + sigma_xy_history.append(s_xy) + sigma_resolved_history.append(resolved_shear) + + flag = " ***" if reason < 0 else "" + if step % 10 == 0 or reason < 0: + print(f"Step {step+1:3d}, t={t:.3f}: " + f"sigma_xy={s_xy:.4f}, resolved={resolved_shear:.4f}, " + f"SNES={reason}, its={its}{flag}") + +# %% [markdown] +r""" +## Stress Evolution + +The plot shows stress at the fault centre over time. The dashed line shows +the analytical VE solution (no yield). The resolved shear stress on the +fault plane should cap at $\tau_y = 0.15$. +""" + +# %% +times = np.array(times) +sigma_xy_history = np.array(sigma_xy_history) +sigma_resolved_history = np.array(sigma_resolved_history) + +# Analytical VE solution (no yield) +t_analytical = np.linspace(0, times[-1], 200) +sigma_ve_analytical = ETA_1 * gamma_dot * (1 - np.exp(-t_analytical / t_relax)) + +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + +# Left panel: sigma_xy and resolved shear vs time +ax = axes[0] +ax.plot(times, sigma_xy_history, 'b-o', markersize=2, label=r'$\sigma_{xy}$ (global)') +ax.plot(times, sigma_resolved_history, 'r-s', markersize=2, + label=r'$\tau_{\mathrm{resolved}}$ (fault plane)') +ax.plot(t_analytical, sigma_ve_analytical, 'k--', alpha=0.5, label='VE analytical (no yield)') +ax.axhline(y=TAU_Y, color='gray', linestyle=':', linewidth=2, + label=rf'$\tau_y = {TAU_Y}$') +ax.set_xlabel('Time') +ax.set_ylabel('Stress') +ax.set_title(f'Stress at fault centre (angle={FAULT_ANGLE_DEG} deg)') +ax.legend(fontsize=9) +ax.grid(True, alpha=0.3) + +# Right panel: stress profile across fault at end of simulation +ax = axes[1] +n_samples = 100 +# Profile perpendicular to the fault, through the centre +profile_dist = np.linspace(-0.4, 0.4, n_samples) +profile_x = 0.5 + profile_dist * n_x # along fault normal direction +profile_y = 0.5 + profile_dist * n_y +# Clip to domain +valid = (profile_x > 0.02) & (profile_x < 0.98) & (profile_y > 0.02) & (profile_y < 0.98) +profile_coords = np.column_stack([profile_x[valid], profile_y[valid]]) + +# Evaluate tau_y field along profile +tau_y_profile = uw.function.evaluate(tau_y_field, profile_coords).flatten() + +# For stress, find nearest tau evaluation points +tau_coords = stokes.tau.coords +tau_data = stokes.tau.data +stress_profile = np.zeros(len(profile_coords)) +for i, pc in enumerate(profile_coords): + dists = np.linalg.norm(tau_coords - pc, axis=1) + idx = np.argmin(dists) + s_xx, s_yy, s_xy = tau_data[idx, 0], tau_data[idx, 1], tau_data[idx, 2] + t_x, t_y = n_y, -n_x + stress_profile[i] = (s_xx * t_x * n_x + s_xy * (t_x * n_y + t_y * n_x) + + s_yy * t_y * n_y) + +ax.plot(profile_dist[valid], np.abs(stress_profile), 'r-', linewidth=2, + label=r'$|\tau_{\mathrm{resolved}}|$') +ax.plot(profile_dist[valid], tau_y_profile, 'k--', linewidth=1, + label=r'$\tau_y$ (yield stress)') +ax.axvline(0, color='gray', linestyle=':', alpha=0.5, label='Fault centre') +ax.set_xlabel('Distance from fault (along normal)') +ax.set_ylabel('Stress') +ax.set_title('Final stress profile across fault') +ax.legend(fontsize=9) +ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.savefig("output/ti_vep_angled_fault.png", dpi=150) +plt.show() + +# %% [markdown] +r""" +## Summary + +This example demonstrates the `TransverseIsotropicVEPFlowModel`: + +1. **Elastic stress buildup**: stress grows from zero following the Maxwell solution +2. **Fault-plane yield**: the resolved shear stress on the fault plane caps at $\tau_y$ +3. **Anisotropic yield**: only the fault-parallel shear component yields; the normal + component continues to build elastically +4. **Smooth spatial transition**: the Gaussian influence function localises yield to + the fault zone, with background material remaining elastic + +The resolved shear criterion ($|\dot\gamma| = \sqrt{|T|^2 - \dot\varepsilon_n^2}$) +correctly identifies fault-plane shear regardless of the fault orientation, avoiding +the orientation-dependent errors that arise from using the global strain rate invariant. +""" diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 9a74d5222..b6c7c2b9f 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1110,8 +1110,8 @@ 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._yield_mode = "softmin" # "min", "harmonic", "smooth", or "softmin" + self._yield_softness = 0.1 # δ parameter for "softmin" mode 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. @@ -1484,12 +1484,15 @@ def viscosity(self): elif self._yield_mode == "softmin": # Smooth approximation to Min(η_ve, η_pl): # η_eff = η_ve / g(f) - # g(f) = (1+f)/2 + √((f-1)² + δ²)/2 ≈ max(1, f) - # where f = η_ve/η_pl and δ = yield_softness. + # g(f) = 1 + softplus(f-1) - softplus(-1) ≈ max(1, f) + # where softplus(x) = (x + √(x² + δ²))/2 and f = η_ve/η_pl. + # Corrected so g(0) = 1 exactly (no spurious yield below onset). # Approaches exact Min as δ→0. No Min/Max in expression. delta = self._yield_softness f = effective_viscosity / vp_effective_viscosity - g = (1 + f) / 2 + sympy.sqrt((f - 1)**2 + delta**2) / 2 + import math # float offset avoids sympy expression blowup in tensor + offset = (-1 + math.sqrt(1 + delta**2)) / 2 + g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset effective_viscosity = effective_viscosity / g else: effective_viscosity = sympy.Min(effective_viscosity, vp_effective_viscosity) @@ -1745,7 +1748,8 @@ def yield_softness(self): Smaller values → sharper yield (closer to Min, less robust). Larger values → smoother transition (more robust, lower stress). - Default 0.5. Only used when ``yield_mode == "softmin"``. + Default 0.1. Only used when ``yield_mode == "softmin"``. + Increase toward 0.5 if SNES convergence is difficult at yield onset. """ return self._yield_softness @@ -2283,7 +2287,7 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias): director = api_tools.Parameter( r"\hat{n}", - lambda inner_self: 1, + lambda inner_self: sympy.Matrix([0] * (inner_self._owning_model.dim - 1) + [1]), "Director orientation", units=None, # Dimensionless unit vector ) @@ -2404,6 +2408,578 @@ def _object_viewer(self): ) +class TransverseIsotropicVEPFlowModel(TransverseIsotropicFlowModel): + r"""Transversely isotropic viscoelastic-plastic flow model for fault mechanics. + + Combines the anisotropic viscosity tensor from :class:`TransverseIsotropicFlowModel` + with viscoelastic stress history and plastic yield limiting on the fault plane. + + The anisotropic viscosity tensor uses two viscosities (η₀ for the bulk, + η₁ for fault-plane shear) and a director n̂ defining the weak plane. + The yield stress τ_y limits the shear stress resolved on the fault plane. + + Parameters + ---------- + unknowns : Unknowns + Solver unknowns (velocity, pressure). + order : int, default=1 + Time integration order for stress history (1 or 2). + material_name : str, optional + Name for disambiguation in multi-material setups. + + See Also + -------- + TransverseIsotropicFlowModel : Anisotropic viscous model (no yield/elasticity). + ViscoElasticPlasticFlowModel : Isotropic VEP model. + """ + + def __init__(self, unknowns, order=1, material_name: str = None): + + self._material_name = material_name + + # Stress history expressions + self._stress_star = expression( + r"{\tau^{*}}", None, + r"Lagrangian Stress at $t - \delta_t$", + ) + self._stress_2star = expression( + r"{\tau^{**}}", None, + r"Lagrangian Stress at $t - 2\delta_t$", + ) + self._E_eff = expression( + r"{\dot{\varepsilon}_{\textrm{eff}}}", None, + "Equivalent value of strain rate (accounting for stress history)", + ) + self._E_eff_inv_II = expression( + r"{\dot{\varepsilon}_{II,\textrm{eff}}}", None, + "Equivalent value of strain rate 2nd invariant (accounting for stress history)", + ) + + self._order = order + self._yield_mode = "softmin" + self._yield_softness = 0.1 + self._bdf_blend = 0.5 + self._max_dt_ratio_for_higher_order = 2.0 + + # Timestep (set by solver) + self._dt = expression(r"{\Delta t}", sympy.oo, "Timestep (set by solver)") + + # BDF coefficients (initialised to BDF-1) + 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) + + return + + class _Parameters(_ParameterBase, _ViscousParameterAlias): + """Parameters for transverse isotropic VEP model. + + Combines anisotropic parameters (η₀, η₁, director) with VEP + parameters (shear_modulus, yield_stress, etc.). + """ + + import underworld3.utilities._api_tools as api_tools + + # Anisotropic parameters + shear_viscosity_0 = api_tools.Parameter( + r"\eta_0", lambda inner_self: 1, + "Bulk shear viscosity", units="Pa*s", + ) + shear_viscosity_1 = api_tools.Parameter( + r"\eta_1", lambda inner_self: 1, + "Fault-plane shear viscosity", units="Pa*s", + ) + director = api_tools.Parameter( + r"\hat{n}", + lambda inner_self: sympy.Matrix([0] * (inner_self._owning_model.dim - 1) + [1]), + "Director orientation (fault normal)", units=None, + ) + + # Elastic parameter + shear_modulus = api_tools.Parameter( + R"{\mu}", lambda inner_self: sympy.oo, + "Shear modulus", units="Pa", + ) + + # Timestep (managed by solver) + @property + def dt_elastic(inner_self): + """Timestep for VE formulas. Set by the solver.""" + return inner_self._owning_model._dt + + @dt_elastic.setter + def dt_elastic(inner_self, value): + if hasattr(value, 'sym'): + inner_self._owning_model._dt.sym = value.sym + else: + inner_self._owning_model._dt.sym = value + + # Viscosity limits + shear_viscosity_min = api_tools.Parameter( + R"{\eta_{\textrm{min}}}", + lambda inner_self: -sympy.oo, + "Shear viscosity, minimum cutoff", units="Pa*s", + ) + + # Yield parameters (applied to fault-plane shear) + yield_stress = api_tools.Parameter( + R"{\tau_{y}}", lambda inner_self: sympy.oo, + "Yield stress (fault-plane shear)", units="Pa", + ) + yield_stress_min = api_tools.Parameter( + R"{\tau_{y, \mathrm{min}}}", + lambda inner_self: -sympy.oo, + "Yield stress minimum cutoff", units="Pa", + ) + strainrate_inv_II_min = api_tools.Parameter( + R"{\dot\varepsilon_{II,\mathrm{min}}}", + lambda inner_self: 0, + "Strain rate invariant minimum value", units="1/s", + ) + + def __init__(inner_self, _owning_model): + inner_self._owning_model = _owning_model + + inner_self._ve_effective_viscosity = expression( + R"{\eta_{\mathrm{eff}}}", None, + "Effective viscosity (elastic, fault-plane)", + ) + inner_self._t_relax = expression( + R"{t_{\mathrm{relax}}}", None, + "Maxwell relaxation time", + ) + + @property + def ve_effective_viscosity(inner_self): + r"""VE effective viscosity using η₁ (fault-plane viscosity).""" + mu_val = inner_self.shear_modulus.sym if hasattr(inner_self.shear_modulus, 'sym') else inner_self.shear_modulus + if mu_val is sympy.oo: + return inner_self.shear_viscosity_1 + + eta = inner_self.shear_viscosity_1 + mu = inner_self.shear_modulus + dt_e = inner_self.dt_elastic + c0 = inner_self._owning_model._bdf_c0 + + el_eff_visc = eta * mu * dt_e / (c0 * eta + mu * dt_e) + inner_self._ve_effective_viscosity.sym = el_eff_visc + return inner_self._ve_effective_viscosity + + @property + def t_relax(inner_self): + r"""Maxwell relaxation time: η₁ / μ.""" + inner_self._t_relax.sym = inner_self.shear_viscosity_1 / inner_self.shear_modulus + return inner_self._t_relax + + ## End of parameters + + @property + def is_elastic(self): + """True if elastic behavior is active (finite shear_modulus).""" + if self.Parameters.shear_modulus.sym is sympy.oo: + return False + return True + + @property + def is_viscoplastic(self): + """True if plastic yielding is active (finite yield_stress).""" + if self.Parameters.yield_stress.sym is sympy.oo: + return False + return True + + @property + def order(self): + """Time integration order (1 or 2).""" + return self._order + + @order.setter + def order(self, value): + """Set time integration order (warns if DFDt already created).""" + self._order = value + self._reset() + solver = getattr(self.Parameters, '_solver', None) + if solver is not None: + ddt = getattr(solver.Unknowns, 'DFDt', None) + if ddt is not None and ddt.order < value: + import warnings + warnings.warn( + f"Setting order={value} but DFDt was created with order={ddt.order}. " + f"Create the model with the desired order before assigning to the solver.", + UserWarning, stacklevel=2, + ) + elif ddt is not None: + solver._order = value + return + + @property + def effective_order(self): + """Effective order accounting for DDt history startup.""" + if self.Unknowns is not None and self.Unknowns.DFDt is not None: + ddt_eff = self.Unknowns.DFDt.effective_order + return min(self._order, ddt_eff) + return self._order + + def _update_bdf_coefficients(self): + """Update BDF coefficient UWexpressions with blending.""" + order = self.effective_order + + 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 + + dt_history = self.Unknowns.DFDt._dt_history + if order >= 2 and len(dt_history) > 0 and dt_history[0] is not None: + try: + ratio = float(dt_current) / float(dt_history[0]) + if ratio > self._max_dt_ratio_for_higher_order: + order = 1 + except (TypeError, ZeroDivisionError): + pass + + coeffs = _bdf_coefficients(order, dt_current, dt_history) + + alpha = self._bdf_blend + if 0 < alpha < 1 and order >= 2: + coeffs_o1 = _bdf_coefficients(1, dt_current, dt_history) + while len(coeffs_o1) < len(coeffs): + coeffs_o1.append(sympy.Integer(0)) + coeffs = [ + (1 - alpha) * c1 + alpha * ck + for c1, ck in zip(coeffs_o1, coeffs) + ] + else: + coeffs = _bdf_coefficients(order, None, []) + + 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] + + @property + def stress_star(self): + r"""Previous timestep stress from history.""" + if self.Unknowns.DFDt is not None: + self._stress_star.sym = self.Unknowns.DFDt.psi_star[0].sym + return self._stress_star + + @property + def E_eff(self): + r"""Effective strain rate including elastic history.""" + E = self.Unknowns.E + + if self.Unknowns.DFDt is not None and self.is_elastic: + mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus + bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] + 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 + return self._E_eff + + @property + def E_eff_inv_II(self): + r"""Second invariant of effective strain rate.""" + E_eff = self.E_eff.sym + self._E_eff_inv_II.sym = sympy.sqrt((E_eff**2).trace() / 2) + return self._E_eff_inv_II + + @property + def viscosity(self): + r"""Effective viscosity for the fault-plane shear component. + + Applies the yield mode (smooth/softmin/min/harmonic) to η₁, + leaving η₀ (bulk) unchanged. The anisotropic tensor handles + the directional dependence. + """ + inner_self = self.Parameters + + if inner_self.yield_stress.sym == sympy.oo: + return inner_self.shear_viscosity_0 + + # η₁ is the fault-plane viscosity that gets yield-limited + eta_1_eff = inner_self.ve_effective_viscosity + + if self.is_viscoplastic: + vp_eff = self._plastic_effective_viscosity + if self._yield_mode == "harmonic": + eta_1_eff = 1 / (1 / eta_1_eff + 1 / vp_eff) + elif self._yield_mode == "smooth": + f = eta_1_eff / vp_eff + eta_1_eff = eta_1_eff * (1 + f) / (1 + f + f**2) + elif self._yield_mode == "softmin": + delta = self._yield_softness + f = eta_1_eff / vp_eff + import math # float offset avoids sympy expression blowup in tensor + offset = (-1 + math.sqrt(1 + delta**2)) / 2 + g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset + eta_1_eff = eta_1_eff / g + else: + eta_1_eff = sympy.Min(eta_1_eff, vp_eff) + + return inner_self.shear_viscosity_0 + + @property + def K(self): + """Effective stiffness for preconditioner.""" + return self.Parameters.shear_viscosity_0 + + @property + def _plastic_effective_viscosity(self): + """Plastic viscosity from resolved fault-plane shear strain rate. + + Computes the in-plane shear magnitude using Pythagoras: + T = ε̇_eff · n (traction-like vector on fault) + ε̇_n = T · n (normal component) + |γ̇| = √(|T|² - ε̇_n²) (in-plane shear magnitude) + + This works in both 2D and 3D — no explicit tangent vector needed. + The formula 2η₁_pl = τ_y / |γ̇| is the same pattern as isotropic + Drucker-Prager but projected onto the fault plane. + """ + parameters = self.Parameters + + ty_val = parameters.yield_stress.sym if hasattr(parameters.yield_stress, 'sym') else parameters.yield_stress + if ty_val is sympy.oo: + return sympy.oo + + Edot = self.E_eff.sym + + # Resolve strain rate onto fault plane via Pythagoras + n = parameters.director.sym + T = Edot * n # "traction" vector on fault + edot_n = (n.T * T)[0, 0] # normal component + T_sq = (T.T * T)[0, 0] # |T|² + gamma_dot_sq = T_sq - edot_n**2 # in-plane shear² + gamma_dot_abs = sympy.sqrt(sympy.Max(gamma_dot_sq, 0)) + + tau_y = parameters.yield_stress + if parameters.yield_stress_min.sym != 0: + tau_y = sympy.Max(parameters.yield_stress_min, tau_y) + + if parameters.strainrate_inv_II_min.sym != 0: + viscosity_yield = tau_y / ( + 2 * (gamma_dot_abs + parameters.strainrate_inv_II_min) + ) + else: + viscosity_yield = tau_y / (2 * gamma_dot_abs) + + return viscosity_yield + + def _build_c_tensor(self): + """Build the anisotropic tensor with VE effective viscosities. + + Both η₀ and η₁ are replaced by their VE effective values: + η₀_ve = η₀·μ·dt / (c₀·η₀ + μ·dt) + η₁_ve = η₁·μ·dt / (c₀·η₁ + μ·dt) + Then η₁_ve is further yield-limited to η₁_eff. This ensures + Δ = η₀_ve - η₁_eff = 0 when η₁ = η₀ and yield is inactive. + """ + + if self._is_setup: + return + + d = self.dim + + # η₀: VE effective (no yield) + eta_0_raw = self.Parameters.shear_viscosity_0 + mu = self.Parameters.shear_modulus + dt_e = self.Parameters.dt_elastic + c0 = self._bdf_c0 + + mu_val = mu.sym if hasattr(mu, 'sym') else mu + if mu_val is sympy.oo: + eta_0 = eta_0_raw.sym if hasattr(eta_0_raw, 'sym') else eta_0_raw + else: + eta_0 = eta_0_raw * mu * dt_e / (c0 * eta_0_raw + mu * dt_e) + + # η₁: VE effective + yield limited + eta_1_eff = self.Parameters.ve_effective_viscosity + + if self.is_viscoplastic: + vp_eff = self._plastic_effective_viscosity + if self._yield_mode == "harmonic": + eta_1_eff = 1 / (1 / eta_1_eff + 1 / vp_eff) + elif self._yield_mode == "smooth": + f = eta_1_eff / vp_eff + eta_1_eff = eta_1_eff * (1 + f) / (1 + f + f**2) + elif self._yield_mode == "softmin": + delta = self._yield_softness + f = eta_1_eff / vp_eff + import math # float offset avoids sympy expression blowup in tensor + offset = (-1 + math.sqrt(1 + delta**2)) / 2 + g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset + eta_1_eff = eta_1_eff / g + else: + eta_1_eff = sympy.Min(eta_1_eff, vp_eff) + + n = self.Parameters.director.sym + Delta = eta_0 - eta_1_eff + + identity = uw.maths.tensor.rank4_identity(d) + lambda_mat = sympy.MutableDenseNDimArray.zeros(d, d, d, d) + + for i in range(d): + for j in range(d): + for k in range(d): + for l in range(d): + base_val = 2 * identity[i, j, k, l] * eta_0 + aniso_correction = ( + 2 * 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] + ) + ) + val = base_val - aniso_correction + if hasattr(val, '__getitem__') and not isinstance(val, (sympy.MatrixBase, sympy.NDimArray)): + val = sympy.Mul(sympy.S.One, val, evaluate=False) + lambda_mat[i, j, k, l] = val + + lambda_mat = uw.maths.tensor.rank4_to_mandel(lambda_mat, d) + self._c = uw.maths.tensor.mandel_to_rank4(lambda_mat, d) + + self._is_setup = True + self._solver_is_setup = False + return + + @property + def flux(self): + """Stress flux for the weak form.""" + return self.stress() + + def stress_projection(self): + """VE stress without plastic correction (for history storage). + + Uses the anisotropic tensor with VE effective viscosities but + no yield limiting (η₁_ve, not η₁_eff). This is the stress that + should be stored in the DFDt history for the next timestep. + """ + edot = self.grad_u + self._build_c_tensor_ve() + # Contract with the VE-only tensor (not self._c which has yield) + c_ve = self._c_ve + if len(c_ve.shape) == 2: + flux = c_ve * edot + else: + flux = sympy.tensorcontraction( + sympy.tensorcontraction(sympy.tensorproduct(c_ve, edot), (1, 5)), (0, 3) + ) + return sympy.Matrix(flux) + + def _build_c_tensor_ve(self): + """Build anisotropic tensor with VE η₁ only (no yield).""" + d = self.dim + eta_0 = self.Parameters.shear_viscosity_0.sym + eta_1_ve = self.Parameters.ve_effective_viscosity + n = self.Parameters.director.sym + Delta = eta_0 - eta_1_ve + + identity = uw.maths.tensor.rank4_identity(d) + lambda_mat = sympy.MutableDenseNDimArray.zeros(d, d, d, d) + + for i in range(d): + for j in range(d): + for k in range(d): + for l in range(d): + base_val = 2 * identity[i, j, k, l] * eta_0 + aniso_correction = ( + 2 * 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] + ) + ) + val = base_val - aniso_correction + if hasattr(val, '__getitem__') and not isinstance(val, (sympy.MatrixBase, sympy.NDimArray)): + val = sympy.Mul(sympy.S.One, val, evaluate=False) + lambda_mat[i, j, k, l] = val + + lambda_mat = uw.maths.tensor.rank4_to_mandel(lambda_mat, d) + self._c_ve = uw.maths.tensor.mandel_to_rank4(lambda_mat, d) + + def stress(self): + """Viscoelastic-plastic anisotropic stress for the weak form. + + Matches isotropic VEP pattern: tensor contraction for current strain + rate, scalar yield-limited viscosity for VE history terms. The tensor + C(η₁_eff) handles anisotropy; the history uses the same η₁_eff as + a scalar multiplier (consistent with how isotropic VEP uses + self.viscosity for both). + """ + self._build_c_tensor() + + # Apply the anisotropic tensor to the effective strain rate + # (current + VE history): σ = C(η₀_ve, η₁_eff) : ε̇_eff + # This is the correct VE formula — the tensor handles anisotropy + # for both current and history contributions uniformly. + edot_eff = self.E_eff.sym if hasattr(self.E_eff, 'sym') else self.E_eff + stress = self._q(edot_eff) + + return stress + + @property + def yield_mode(self): + r"""How to apply yield limiting to the fault-plane viscosity. + + Same options as :class:`ViscoElasticPlasticFlowModel`: + ``"smooth"`` (default), ``"softmin"``, ``"harmonic"``, ``"min"``. + """ + return self._yield_mode + + @yield_mode.setter + def yield_mode(self, value): + if value not in ("min", "harmonic", "smooth", "softmin"): + raise ValueError(f"yield_mode must be 'min', 'harmonic', 'smooth', or 'softmin', got '{value}'") + self._yield_mode = value + self._reset() + + @property + def yield_softness(self): + """Regularisation parameter δ for softmin mode.""" + return self._yield_softness + + @yield_softness.setter + def yield_softness(self, value): + self._yield_softness = value + self._reset() + + @property + def bdf_blend(self): + """BDF coefficient blending: 0=pure O1, 0.5=default, 1=pure O2.""" + return self._bdf_blend + + @bdf_blend.setter + def bdf_blend(self, value): + self._bdf_blend = value + + @property + def requires_stress_history(self): + """Transverse isotropic VEP requires stress history tracking.""" + return True + + @property + def plastic_fraction(self): + """Fraction of strain rate that is plastic.""" + eta_1_ve = self.Parameters.ve_effective_viscosity + eta_1_eff = self.viscosity + # viscosity property returns η₀, need to compare η₁ effective vs η₁ ve + # This is approximate for the anisotropic case + return sympy.Max(0, 1 - eta_1_eff / eta_1_ve.sym if hasattr(eta_1_ve, 'sym') else 0) + + class MultiMaterialConstitutiveModel(Constitutive_Model): r""" Multi-material constitutive model using level-set weighted flux averaging.