diff --git a/docs/examples/WIP/VEP_Fault_Investigation.py b/docs/examples/WIP/VEP_Fault_Investigation.py new file mode 100644 index 000000000..97cab167d --- /dev/null +++ b/docs/examples/WIP/VEP_Fault_Investigation.py @@ -0,0 +1,280 @@ +# %% [markdown] +""" +# VEP Embedded Fault — Investigation Notebook + +Horizontal fault at y=0.5 with gaussian influence function for tau_y. +Step through the model and visualise to find the source of SNES divergence. +""" + +# %% +#| echo: false # Hide in html version + +# This is required to fix pyvista +# (visualisation) crashes in interactive notebooks (including on binder) + +import nest_asyncio +nest_asyncio.apply() + +# %% +import numpy as np +import sympy +import underworld3 as uw +from underworld3.systems import Stokes + +# %% [markdown] +""" +## Parameters +""" + +# %% +ETA = 1.0 +MU = 1.0 +TAU_Y_FAULT = 0.2 +TAU_Y_BULK = 200.0 +FAULT_WIDTH = 0.04 +DT = 0.05 +V_TOP = 0.5 + +# %% [markdown] +""" +## Mesh and variables +""" + +# %% +mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 32), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=2) + +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] +""" +## Fault surface and yield stress +""" + +# %% +fault_points = np.array([[0.0, 0.5, 0.0], [1.0, 0.5, 0.0]]) +fault = uw.meshing.Surface("fault", mesh, fault_points) +fault.discretize() + +# Interpolate weakness (1/tau_y) — avoids steep gaussian ramp in tau_y +weakness = fault.influence_function( + width=FAULT_WIDTH, + value_near=1 / TAU_Y_FAULT, + value_far=1 / TAU_Y_BULK, + profile="gaussian", +) +tau_y_field = 1 / weakness + +print(f"Fault: {fault.n_vertices} vertices") + +# %% [markdown] +""" +## Solver setup +""" + +# %% + +# %% +stokes = Stokes(mesh, velocityField=v, pressureField=p) +# Create model instance BEFORE assigning — DFDt order is set at assignment time +cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=1) +stokes.constitutive_model = cm + +# %% + +cm.Parameters.shear_viscosity_0 = ETA +cm.Parameters.shear_modulus = MU +cm.Parameters.yield_stress = tau_y_field +# yield_mode="smooth" is the default (corrected harmonic, no Min/Max) +cm.Parameters.shear_viscosity_min = ETA * 1.0e-2 +cm.Parameters.strainrate_inv_II_min = 1.0e-5 + +stokes.saddle_preconditioner = 1 / cm.K +stokes.tolerance = 1.0e-4 + +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]) +stokes.petsc_options["ksp_type"] = "fgmres" +stokes.petsc_options["snes_force_iteration"] = None # always do at least 1 SNES iteration + +# %% [markdown] +""" +## Visualisation helpers +""" + +# %% +def plot_state(step_num, title_extra=""): + """Plot velocity, pressure, and stress for the current state.""" + + if not uw.is_notebook(): + print("Skipping visualisation (not a notebook)") + return + + import pyvista as pv + import underworld3.visualisation as vis + + pvmesh = vis.mesh_to_pv_mesh(mesh) + pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v.sym) + pvmesh.point_data["P"] = vis.scalar_fn_to_pv_points(pvmesh, p.sym) + pvmesh.point_data["Vmag"] = np.linalg.norm(pvmesh.point_data["V"], axis=1) + + # Stress from tau + tau_data = stokes.tau.data + tau_coords = stokes.tau.coords + + # Velocity arrows + velocity_points = vis.meshVariable_to_pv_cloud(v) + velocity_points.point_data["V"] = vis.vector_fn_to_pv_points(velocity_points, v.sym) + + pl = pv.Plotter(shape=(1, 3), window_size=(1500, 500)) + + # Panel 1: Velocity magnitude + pl.subplot(0, 0) + pl.add_mesh(pvmesh, scalars="Vmag", cmap="viridis", show_edges=False) + pl.add_arrows(velocity_points.points[::3], velocity_points.point_data["V"][::3], + mag=0.15, color="white") + pl.add_title(f"Velocity (step {step_num}){title_extra}") + pl.camera_position = "xy" + + # Panel 2: Pressure + pl.subplot(0, 1) + pl.add_mesh(pvmesh, scalars="P", cmap="coolwarm", show_edges=False) + pl.add_title("Pressure") + pl.camera_position = "xy" + + # Panel 3: Stress sigma_xy from tau + tau_cloud = pv.PolyData(np.column_stack([tau_coords, np.zeros(len(tau_coords))])) + tau_cloud.point_data["sigma_xy"] = tau_data[:, 2] + pl.subplot(0, 2) + pl.add_mesh(tau_cloud, scalars="sigma_xy", cmap="RdBu_r", + point_size=8, render_points_as_spheres=True) + pl.add_title("sigma_xy (from tau)") + pl.camera_position = "xy" + + pl.show() + +# %% [markdown] +""" +## Check tau_y field before solving +""" + +# %% +# Evaluate tau_y at mesh nodes to verify the fault zone +tau_y_vals = uw.function.evaluate(tau_y_field, mesh.X.coords) +print(f"tau_y range: [{tau_y_vals.min():.4f}, {tau_y_vals.max():.4f}]") + +if uw.is_notebook(): + import pyvista as pv + import underworld3.visualisation as vis + + pvmesh = vis.mesh_to_pv_mesh(mesh) + pvmesh.point_data["tau_y"] = tau_y_vals.flatten() + + pl = pv.Plotter(window_size=(600, 500)) + pl.add_mesh(pvmesh, scalars="tau_y", cmap="coolwarm", show_edges=True) + pl.add_title("Yield stress field") + pl.camera_position = "xy" + pl.show() + +# %% [markdown] +""" +## Step through the model + +Run a few steps in the elastic regime, then step carefully through yield onset. +""" + +# %% +import numpy as _np + +for step in range(1 + int(2/DT)): + # Inspect psi_star before solve + s0 = stokes.DFDt.psi_star[0].data + s1 = stokes.DFDt.psi_star[1].data if stokes.DFDt.order >= 2 else None + s0_xy = f"[{s0[:,2].min():.4f}, {s0[:,2].max():.4f}]" + s1_xy = f"[{s1[:,2].min():.4f}, {s1[:,2].max():.4f}]" if s1 is not None else "n/a" + + stokes.solve(timestep=DT, zero_init_guess=(step == 0)) + t = (step + 1) * DT + reason = stokes.snes.getConvergedReason() + its = stokes.snes.getIterationNumber() + sigma_xy = stokes.tau.data[:, 2] + + flag = " ***" if reason < 0 else "" + print(f"Step {step+1:3d}, t={t:.2f}: σ[{sigma_xy.min():.4f},{sigma_xy.max():.4f}] " + f"σ*{s0_xy} σ**{s1_xy} SNES={reason} its={its}{flag}") + + +# %% [markdown] +""" +## Yield onset — step carefully +""" + +# %% +# Continue stepping — yield should start around step 9-10 +# for step in range(9, 15): +# stokes.solve(timestep=DT, zero_init_guess=False) +# t = (step + 1) * DT +# reason = stokes.snes.getConvergedReason() +# its = stokes.snes.getIterationNumber() +# sigma_xy = stokes.tau.data[:, 2] +# print(f"Step {step+1}, t={t:.2f}: sigma_xy [{sigma_xy.min():.4f}, {sigma_xy.max():.4f}], " +# f"SNES={reason}, its={its}") + +# if reason < 0: +# print(f" *** DIVERGED at step {step+1} ***") +# plot_state(step + 1, f" (DIVERGED, SNES={reason})") +# break + +# %% [markdown] +""" +## Continue past divergence +""" + +# %% +# # Keep going to see if the solution stabilises +# for step in range(15, 25): +# stokes.solve(timestep=DT, zero_init_guess=False) +# t = (step + 1) * DT +# reason = stokes.snes.getConvergedReason() +# its = stokes.snes.getIterationNumber() +# sigma_xy = stokes.tau.data[:, 2] +# print(f"Step {step+1}, t={t:.2f}: sigma_xy [{sigma_xy.min():.4f}, {sigma_xy.max():.4f}], SNES={reason}, its={its}") + +plot_state(25, " (final)") + +# %% [markdown] +""" +## Stress profile across fault +""" + +# %% +import matplotlib +if not uw.is_notebook(): + matplotlib.use("Agg") +import matplotlib.pyplot as plt + +coords = stokes.tau.coords +sd = stokes.tau.data +x_coords = coords[:, 0] +y_coords = coords[:, 1] + +# Get nodes near x=0.5 +near_centre = np.abs(x_coords - 0.5) < 0.05 +y_profile = y_coords[near_centre] +sigma_profile = sd[near_centre, 2] +sort_idx = y_profile.argsort() + +fig, ax = plt.subplots(figsize=(6, 4)) +ax.plot(y_profile[sort_idx], sigma_profile[sort_idx], 'r-o', markersize=3) +ax.axvline(0.5, color='gray', linestyle=':', alpha=0.5, label="fault") +ax.set_xlabel("y") +ax.set_ylabel(r"$\sigma_{xy}$") +ax.set_title("Stress profile at x=0.5") +ax.legend() +ax.grid(True, alpha=0.3) +plt.show() diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 4629c91ea..889d3ee76 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -589,6 +589,15 @@ def requires_stress_history(self): """ return False + @property + def plastic_fraction(self): + """Fraction of strain rate that is plastic (0 for non-plastic models). + + Returns a sympy expression that can be evaluated post-solve via + ``uw.function.evaluate(cm.plastic_fraction, coords)``. + """ + return sympy.Integer(0) + def _build_c_tensor(self): """Return the identity tensor of appropriate rank (e.g. for projections)""" @@ -763,6 +772,11 @@ def grad_u(self): # edot = (ddu + ddu.T) / 2 # return edot + @property + def plastic_fraction(self): + """Fraction of strain rate that is plastic: 1 - η_vp / η_viscous.""" + return sympy.Max(0, 1 - self.viscosity / self.Parameters.shear_viscosity_0) + def _build_c_tensor(self): """For this constitutive law, we expect just a viscosity function""" @@ -1096,6 +1110,14 @@ 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 + + # Timestep — set by the solver before each solve(). Not a user parameter. + # Initialised to oo (viscous limit). The solver overwrites this with the + # actual timestep on every call to solve(timestep=dt). + self._dt = expression(r"{\Delta t}", sympy.oo, "Timestep (set by solver)") # BDF coefficients as UWexpressions — route through PetscDS constants[]. # Updated each step by _update_bdf_coefficients() before solve. @@ -1138,12 +1160,23 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias): units="Pa", ) - dt_elastic = api_tools.Parameter( - R"{\Delta t_{e}}", - lambda inner_self: sympy.oo, - "Elastic timestep", - units="s", - ) + @property + def dt_elastic(inner_self): + """Timestep for VE formulas. Set by the solver, not a user parameter. + + Returns the UWexpression that the solver updates before each solve. + This flows through PetscDS constants[] so the JIT-compiled pointwise + functions always see the current timestep. + """ + return inner_self._owning_model._dt + + @dt_elastic.setter + def dt_elastic(inner_self, value): + """Allow the solver to set dt via Parameters.dt_elastic = timestep.""" + if hasattr(value, 'sym'): + inner_self._owning_model._dt.sym = value.sym + else: + inner_self._owning_model._dt.sym = value shear_viscosity_min = api_tools.Parameter( R"{\eta_{\textrm{min}}}", @@ -1249,9 +1282,33 @@ def order(self): @order.setter def order(self, value): - """Set the time integration order.""" + """Set the time integration order. + + If the model is already attached to a solver with a DFDt, this will + warn if the DFDt was created with a lower order (since it can't be + changed after creation — the DFDt allocates history buffers at init). + """ self._order = value self._reset() + + # Propagate to connected solver if present + 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 the solver's DFDt was already " + f"created with order={ddt.order}. The DFDt order cannot be " + f"changed after creation. To use order={value}, create the " + f"model with the desired order before assigning to the solver:\n" + f" cm = ViscoElasticPlasticFlowModel(stokes.Unknowns, order={value})\n" + f" stokes.constitutive_model = cm", + UserWarning, + stacklevel=2, + ) + elif ddt is not None: + solver._order = value return @property @@ -1267,6 +1324,11 @@ def effective_order(self): return min(self._order, self.Unknowns.DFDt.effective_order) return self._order + # Maximum timestep ratio (dt_new / dt_old) for which BDF-2+ is safe. + # Beyond this, fall back to BDF-1 to avoid negative-stress extrapolation + # when stress history is non-smooth (e.g. yield events). + _max_dt_ratio_for_higher_order = 2.0 + def _update_bdf_coefficients(self): """Update BDF coefficient UWexpressions from current dt_elastic and DDt history. @@ -1274,16 +1336,44 @@ def _update_bdf_coefficients(self): correct coefficients to the compiled pointwise functions. The coefficient UWexpressions (_bdf_c0..c3) are referenced symbolically in ve_effective_viscosity, E_eff, and stress() — their numeric values flow through PetscDSSetConstants. + + When the timestep ratio exceeds ``_max_dt_ratio_for_higher_order``, + BDF-2+ coefficients can cause negative stress extrapolation if the + stress history is non-smooth (e.g. after a yield event). In this case + we fall back to BDF-1 coefficients for safety. """ + 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 - coeffs = _bdf_coefficients( - self.effective_order, dt_current, self.Unknowns.DFDt._dt_history - ) + + # Guard: fall back to BDF-1 when timestep increases too rapidly + 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 # symbolic dt — can't evaluate, keep requested order + + coeffs = _bdf_coefficients(order, dt_current, dt_history) + + # Blend with O1 coefficients for stability + # 0 = pure O1, 0.5 = balanced (default), 1 = pure requested order + 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(self.effective_order, None, []) + coeffs = _bdf_coefficients(order, None, []) # Pad to length 4 while len(coeffs) < 4: @@ -1358,12 +1448,16 @@ def K(self): def viscosity(self): r"""Effective viscosity combining visco-elastic and plastic limits. - Returns :math:`\min(\eta_{\mathrm{ve}}, \tau_y / 2\dot{\varepsilon}_{II})`. - """ - # detect if values we need are defined or are placeholder symbols + The yield mode controls how η_ve and η_pl are combined: - ## Do we want this to be an expression of its own ? If so, define above in __init__() and - ## make sure it is updated in this call, rather than being replaced. + - ``"smooth"`` (default): corrected harmonic ``η_ve·(1+f)/(1+f+f²)`` + where ``f = η_ve/η_pl``. Converges to η_pl at deep yielding, + no Min/Max discontinuities. + - ``"harmonic"``: ``1/(1/η_ve + 1/η_pl)``. Smooth but undershoots τ_y + when η_ve is small relative to η_pl. + - ``"min"``: sharp ``Min(η_ve, η_pl)``. Exact yield stress but can + cause SNES divergence with higher-order BDF time integration. + """ inner_self = self.Parameters @@ -1374,21 +1468,44 @@ def viscosity(self): if self.is_viscoplastic: vp_effective_viscosity = self._plastic_effective_viscosity - effective_viscosity = sympy.Min(effective_viscosity, vp_effective_viscosity) - - ## Why is it p**2 here ? - # p = self.plastic_correction() - # effective_viscosity *= 2 * p**2 / (1 + p**2) - - # effective_viscosity *= self.plastic_correction() + if self._yield_mode == "harmonic": + effective_viscosity = 1 / (1 / effective_viscosity + 1 / vp_effective_viscosity) + elif self._yield_mode == "smooth": + # Corrected harmonic: cancels the excess 1/η_ve contribution + # at deep yielding while staying smooth everywhere. + # η_eff = η_ve · (1+f) / (1 + f + f²) + # where f = η_ve/η_pl measures yield overshoot. + # + # f → 0 (elastic): η_eff → η_ve (no correction) + # f → ∞ (yielding): η_eff → η_pl (exact yield) + # No Min/Max — just arithmetic. Continuous derivatives. + f = effective_viscosity / vp_effective_viscosity + effective_viscosity = effective_viscosity * (1 + f) / (1 + f + f**2) + 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. + # 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 + effective_viscosity = effective_viscosity / g + else: + effective_viscosity = sympy.Min(effective_viscosity, vp_effective_viscosity) - # If we want to apply limits to the viscosity but see caveat above + # Apply viscosity floor — but skip for smooth/harmonic yield modes + # where the outer Max creates a nested Min/Max that breaks the + # BDF-2 Jacobian. Those modes are already smooth and bounded. if inner_self.shear_viscosity_min.sym != -sympy.oo: - return sympy.Max( - effective_viscosity, - inner_self.shear_viscosity_min, - ) + if self.is_viscoplastic and self._yield_mode in ("harmonic", "smooth", "softmin"): + return effective_viscosity + else: + return sympy.Max( + effective_viscosity, + inner_self.shear_viscosity_min, + ) else: return effective_viscosity @@ -1522,12 +1639,10 @@ def stress_projection(self): return stress def stress(self): - """viscoelastic stress projection (no plastic response)""" + """Viscoelastic(-plastic) deviatoric stress for the weak form.""" edot = self.grad_u - # This is a scalar viscosity ... - stress = 2 * self.viscosity * edot if self.Unknowns.DFDt is not None: @@ -1536,7 +1651,6 @@ def stress(self): mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] - # History contribution: 2·η_eff · (-Σ cᵢ·σ_star[i-1]) / (2·μ·dt) for i in range(self.Unknowns.DFDt.order): stress += 2 * self.viscosity * ( -bdf_cs[i] * self.Unknowns.DFDt.psi_star[i].sym / (2 * mu_dt) @@ -1598,11 +1712,74 @@ def _object_viewer(self): ## Todo: add all the other properties in here ) + @property + def yield_mode(self): + r"""How to combine VE and plastic viscosities. + + ``"smooth"`` (default): corrected harmonic — + ``η_ve · (1+f) / (1+f+f²)`` where ``f = η_ve/η_pl``. + Smooth, no Min/Max. Best balance of accuracy and robustness. + ``"softmin"``: smooth approximation to Min — + ``η_ve / g(f)`` where ``g(f) ≈ max(1, f)`` with smoothing + parameter δ (``yield_softness``, default 0.5). + Closer to exact yield than ``"smooth"`` but less robust. + ``"harmonic"``: parallel blending — ``1/(1/η_ve + 1/η_pl)``. + Smooth but undershoots τ_y for soft materials. + ``"min"``: sharp cutoff — ``Min(η_ve, η_pl)``. + Exact yield but can cause SNES divergence with BDF-2. + """ + 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): + r"""Regularisation parameter δ for ``"softmin"`` yield mode. + + Controls how closely the soft minimum approximates the sharp Min. + 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"``. + """ + return self._yield_softness + + @yield_softness.setter + def yield_softness(self, value): + self._yield_softness = value + self._reset() + + @property + def bdf_blend(self): + r"""Blending parameter α for BDF history coefficients. + + 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) + """ + return self._bdf_blend + + @bdf_blend.setter + def bdf_blend(self, value): + self._bdf_blend = value + @property def requires_stress_history(self): """VEP models always require stress history tracking.""" return True + @property + def plastic_fraction(self): + """Fraction of strain rate that is plastic: 1 - η_vep / η_ve.""" + return sympy.Max(0, 1 - self.viscosity / self.Parameters.ve_effective_viscosity.sym) + @property def is_elastic(self): """True if elastic behavior is active (finite dt_elastic and shear_modulus).""" diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 82bbd8fb2..5316c87fb 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -930,7 +930,10 @@ class SolverBaseClass(uw_object): self._constitutive_model = model_or_class self._constitutive_model.Unknowns = self.Unknowns self._constitutive_model._solver_is_setup = False - self._constitutive_model.order = self._order + # Only override the constitutive model's order if the solver has + # an explicit VE order (> 0). Otherwise preserve the model's default. + if self._order > 0: + self._constitutive_model.order = self._order # Establish bidirectional reference so parameter changes can propagate to solver self._constitutive_model.Parameters._solver = self @@ -938,7 +941,8 @@ class SolverBaseClass(uw_object): ### checking if it's a class elif type(model_or_class) == type(uw.constitutive_models.Constitutive_Model): self._constitutive_model = model_or_class(self.Unknowns) - self._constitutive_model.order = self._order + if self._order > 0: + self._constitutive_model.order = self._order # Establish bidirectional reference so parameter changes can propagate to solver self._constitutive_model.Parameters._solver = self @@ -950,15 +954,11 @@ class SolverBaseClass(uw_object): "constitutive_model must be a valid class or instance of a valid class" ) - # Check that the solver can support this constitutive model's requirements. - # Models with stress history (VEP) need a solver that manages DFDt — e.g. VE_Stokes. - # Using them on a plain Stokes solver silently drops the history terms. + # If the constitutive model requires stress history (e.g. VEP), create the + # DFDt infrastructure lazily. This means users don't need to choose between + # Stokes and VE_Stokes — the solver adapts to the constitutive model. if self._constitutive_model.requires_stress_history and self.Unknowns.DFDt is None: - raise TypeError( - f"{type(self._constitutive_model).__name__} requires stress history tracking " - f"(DFDt). Use uw.systems.VE_Stokes instead of uw.systems.Stokes, or provide " - f"a DFDt object when constructing the solver." - ) + self._create_stress_history_ddt(order=self._constitutive_model.order) # May not work due to flux being incomplete if self.Unknowns.DFDt is not None: diff --git a/src/underworld3/function/functions_unit_system.py b/src/underworld3/function/functions_unit_system.py index fddc8d873..c4d275f0f 100644 --- a/src/underworld3/function/functions_unit_system.py +++ b/src/underworld3/function/functions_unit_system.py @@ -35,7 +35,7 @@ def evaluate( coords, coord_sys=None, other_arguments=None, - simplify=True, + simplify=False, verbose=False, evalf=False, mode="default", @@ -364,7 +364,7 @@ def global_evaluate( coords=None, coord_sys=None, other_arguments=None, - simplify=True, + simplify=False, verbose=False, evalf=False, mode="default", diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 067e4686f..1c8ed0f60 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -687,6 +687,169 @@ def __init__( return + def _create_stress_history_ddt(self, order=2): + """Create DFDt for stress history tracking (VE/VEP models). + + Called automatically when a constitutive model with + ``requires_stress_history = True`` is assigned. Can also be called + explicitly to pre-create the DFDt with a specific order. + """ + if self.Unknowns.DFDt is not None: + return # already created + + self._order = order + self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( + self.mesh, + sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), + self.u.sym, + vtype=uw.VarType.SYM_TENSOR, + degree=self.u.degree - 1, + continuous=True, + varsymbol=rf"{{F[ {self.u.symbol} ] }}", + verbose=self.verbose, + bcs=None, + order=order, + smoothing=0.0001, + ) + + @timing.routine_timer_decorator + def solve( + self, + zero_init_guess: bool = True, + timestep: float = None, + _force_setup: bool = False, + verbose=False, + evalf=False, + order=None, + picard: int = 0, + ): + """Solve the Stokes system, with optional viscoelastic stress history. + + When a constitutive model with stress history is active (DFDt is not None), + the solve includes pre/post hooks for advecting stress history, updating + BDF coefficients, and projecting the actual stress after the solve. + + Parameters + ---------- + zero_init_guess : bool + If True, use zero initial guess. Otherwise use current field values. + timestep : float, optional + Advection timestep. Required when stress history is active. + _force_setup : bool + Force rebuild of pointwise functions. + verbose : bool + Enable verbose output. + evalf : bool + Force numerical evaluation during history updates. + order : int, optional + Override the VE time integration order. + picard : int, default=0 + Number of Picard iterations before switching to Newton. + Picard uses a simplified Jacobian and can help convergence + for strongly nonlinear problems like VEP at yield onset. + """ + + has_stress_history = self.Unknowns.DFDt is not None + + if has_stress_history: + if timestep is None: + raise ValueError( + "timestep is required for viscoelastic solve. " + "Call stokes.solve(timestep=dt)" + ) + + # dt_elastic must always equal the solve timestep. The constitutive + # model's VE formulas (eta_eff, stress history terms) all reference + # Parameters.dt_elastic. If it differs from the actual timestep, + # the stress computation is inconsistent with the time integration. + self.constitutive_model.Parameters.dt_elastic = timestep + + if order is None or order > self._order: + order = self._order + + if _force_setup: + self.is_setup = False + + # Re-setup when effective_order changes (DDt history ramp-up) + _current_eff_order = self.constitutive_model.effective_order + if not hasattr(self, '_prev_effective_order'): + self._prev_effective_order = None + if _current_eff_order != self._prev_effective_order: + self.is_setup = False + self.constitutive_model._solver_is_setup = False + self._prev_effective_order = _current_eff_order + + if not self.constitutive_model._solver_is_setup: + self.is_setup = False + self.DFDt.psi_fn = self.constitutive_model.flux.T + + if not self.is_setup: + self._setup_pointwise_functions(verbose) + self._setup_discretisation(verbose) + self._setup_solver(verbose) + + # 1. ADVECT stress history along characteristics + if uw.mpi.rank == 0 and verbose: + print(f"Stokes solver - advect stress history", flush=True) + + self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, + store_result=False) + self.constitutive_model._update_bdf_coefficients() + + # 2. SOLVE + if uw.mpi.rank == 0 and verbose: + print(f"Stokes solver - solve", flush=True) + + super().solve( + zero_init_guess, + _force_setup=_force_setup, + verbose=verbose, + picard=picard, + ) + + # 3. PROJECT actual stress and SHIFT history + if uw.mpi.rank == 0 and verbose: + print(f"Stokes solver - store stress and shift history", flush=True) + + import numpy as np + + _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) + + self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux + self.DFDt._psi_star_projection_solver.smoothing = 0.0 + self.DFDt._psi_star_projection_solver.solve(verbose=verbose) + + for i in range(self.DFDt.order - 1, 0, -1): + if i == 1: + self.DFDt.psi_star[i].array[...] = _advected_sigma_star + else: + self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] + + self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf) + + self.is_setup = True + self.constitutive_model._solver_is_setup = True + + else: + # Plain Stokes — no stress history + super().solve( + zero_init_guess, + _force_setup=_force_setup, + verbose=verbose, + ) + + @property + def tau(self): + r"""Deviatoric stress from the most recent solve. + + When stress history is active (VEP), returns ``psi_star[0]`` which + contains the actual projected stress. Otherwise falls through to the + base class lazy projection. + """ + if self.Unknowns.DFDt is not None: + return self.DFDt.psi_star[0] + return super().tau + # ========================================================================= # PETSc Residual Templates # These define the weak form terms assembled by PETSc's finite element system. @@ -1126,72 +1289,26 @@ def estimate_dt(self): class SNES_VE_Stokes(SNES_Stokes): - r""" - Viscoelastic Stokes equation solver. - - Provides a discrete representation of the Stokes flow equations with - incompressibility (or near-incompressibility) constraint and a flux - history term for viscoelastic modelling. Inherits from :class:`SNES_Stokes`. - - Momentum equation: - - .. math:: - - -\nabla \cdot \underbrace{\left[ \boldsymbol{\tau} - p \mathbf{I} - \right]}_{\mathbf{F}} = \underbrace{\mathbf{f}}_{\mathbf{h}} - - Continuity equation: - - .. math:: + r"""Viscoelastic Stokes solver (backward-compatibility wrapper). - \underbrace{\nabla \cdot \mathbf{u}}_{\mathbf{h}_p} = 0 + .. deprecated:: + Use ``uw.systems.Stokes`` directly with a + ``ViscoElasticPlasticFlowModel`` constitutive model. The Stokes + solver now creates stress history infrastructure automatically + when the constitutive model requires it. - The flux term is a deviatoric stress :math:`\boldsymbol{\tau}` related - to velocity gradients :math:`\nabla \mathbf{u}` through a viscosity - tensor :math:`\eta`, plus a volumetric (pressure) part :math:`p`: - - .. math:: - - \mathbf{F}: \quad \boldsymbol{\tau} = \frac{\eta}{2} - \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) - - The constraint equation :math:`\mathbf{h}_p = 0` is incompressible flow - by default but can be set to any function of :math:`\mathbf{u}` and - :math:`\nabla \cdot \mathbf{u}`. + This wrapper pre-creates the DFDt with a specific ``order`` parameter, + which is useful when you want to control the BDF order before assigning + the constitutive model. Parameters ---------- mesh : Mesh The computational mesh. - velocityField : MeshVariable, optional - Mesh variable for velocity. Created automatically if not provided. - pressureField : MeshVariable, optional - Mesh variable for pressure. Created automatically if not provided. - degree : int, default=2 - Polynomial degree for velocity elements. order : int, default=2 - Order parameter (typically same as degree). - p_continuous : bool, default=True - If False, use discontinuous pressure elements. - verbose : bool, default=False - Enable verbose output. - DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional - Time derivative operator (may be used in child classes). - - Notes - ----- - - The viscosity tensor :math:`\boldsymbol{\eta}` is set via the - ``constitutive_model`` property - - For viscoelastic problems, the flux term contains stress history - tracked on a particle swarm - - Augmented Lagrangian approach adds :math:`\lambda \nabla \cdot \mathbf{u}` - to penalize incompressibility - - Pressure element order determines mixed FEM integration order - - See Also - -------- - SNES_Stokes : Base Stokes solver. - uw.constitutive_models.ViscoElasticPlasticFlowModel : Constitutive model for VE flow. + BDF time integration order for stress history. + **kwargs + All other arguments are passed to :class:`SNES_Stokes`. """ instances = 0 @@ -1205,12 +1322,9 @@ def __init__( order: Optional[int] = 2, p_continuous: Optional[bool] = True, verbose: Optional[bool] = False, - # DuDt Not used in VE, but may be in child classes DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None, ): - - # Stokes is parent (will not build DuDt or DFDt) super().__init__( mesh, velocityField, @@ -1222,22 +1336,8 @@ def __init__( DFDt=DFDt, ) - self._order = order # VE time-order - - if self.Unknowns.DFDt is None: - self.Unknowns.DFDt = uw.systems.ddt.SemiLagrangian( - self.mesh, - sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim), - self.u.sym, - vtype=uw.VarType.SYM_TENSOR, - degree=self.u.degree - 1, - continuous=True, - varsymbol=rf"{{F[ {self.u.symbol} ] }}", - verbose=self.verbose, - bcs=None, - order=self._order, - smoothing=0.0001, - ) + # Pre-create DFDt so it's available before constitutive model is set + self._create_stress_history_ddt(order=order) return @@ -1246,160 +1346,6 @@ def delta_t(self): """Elastic timestep from the constitutive model.""" return self.constitutive_model.Parameters.dt_elastic - @property - def tau(self): - r"""Deviatoric stress from the most recent solve (stored in history). - - For VE_Stokes, the stress is projected into ``psi_star[0]`` after - each solve. This override returns that variable directly — no - additional projection needed. - - Returns - ------- - MeshVariable - The stress history variable containing the actual deviatoric - stress from the most recent solve. - """ - return self.DFDt.psi_star[0] - - ## Solver needs to update the stress history terms as well as call the SNES solve: - - @timing.routine_timer_decorator - def solve( - self, - zero_init_guess: bool = True, - timestep: float = None, - _force_setup: bool = False, - verbose=False, - evalf=False, - order=None, - ): - """ - Generates solution to constructed system. - - Params - ------ - zero_init_guess: - If `True`, a zero initial guess will be used for the - system solution. Otherwise, the current values of `self.u` will be used. - """ - - if order is None or order > self._order: - order = self._order - - if timestep is None: - timestep = self.delta_t.sym - - # dt_elastic is a constitutive parameter (relaxation timescale) — - # never overwritten by solve(). The advection timestep (for departure - # point tracing) and the elastic relaxation timescale are independent. - - if _force_setup: - self.is_setup = False - - # Re-setup when effective_order changes (e.g. DDt history ramp-up - # from order 1 to order 2). The JIT-compiled pointwise functions - # depend on the order used in the constitutive model. - _current_eff_order = self.constitutive_model.effective_order - if not hasattr(self, '_prev_effective_order'): - self._prev_effective_order = None - if _current_eff_order != self._prev_effective_order: - self.is_setup = False - self.constitutive_model._solver_is_setup = False - self._prev_effective_order = _current_eff_order - - if not self.constitutive_model._solver_is_setup: - self.is_setup = False - self.DFDt.psi_fn = self.constitutive_model.flux.T - - if not self.is_setup: - self._setup_pointwise_functions(verbose) - self._setup_discretisation(verbose) - self._setup_solver(verbose) - - # --- Stress history management via standard DDt pathway --- - # - # update_pre_solve(advect_only=True) performs: - # 1. History shift: psi_star[i] ← psi_star[i-1] - # 2. Skip psi_fn evaluation (psi_star[0] already has projected stress) - # 3. Advect all history levels to upstream positions along characteristics - # - # After this: psi_star[0] = σ* (previous stress at upstream), - # psi_star[1] = σ** (stress from 2 steps ago, double-traced). - - if uw.mpi.rank == 0 and verbose: - print(f"VE Stokes solver - advect stress history", flush=True) - - self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, - store_result=False) - - # Update BDF coefficients from current dt_elastic and DDt history. - # These are UWexpressions that route through PetscDS constants[], - # so the compiled pointwise functions pick up the new values without - # JIT recompilation. - self.constitutive_model._update_bdf_coefficients() - - # 2. SOLVE: PETSc uses the advected σ*, σ** via the constitutive model - - if uw.mpi.rank == 0 and verbose: - print(f"VE Stokes solver - solve Stokes flow", flush=True) - - super().solve( - zero_init_guess, - _force_setup=_force_setup, - verbose=verbose, - picard=0, - ) - - # 3. STORE ACTUAL STRESS and SHIFT HISTORY. - # - # After advection + solve: - # psi_star[0] = advected σ* (used by the solver) - # psi_star[1] = advected σ** (used by the solver) - # - # We need to: - # a) Project actual stress → psi_star[0] (while σ*, σ** are intact) - # b) Save the advected σ* into psi_star[1] for next step's σ** - # (chained characteristic tracing) - # - # The projection reads psi_star[0..1] via stress_deviator, so we - # must project BEFORE shifting. Then save the advected σ* and - # overwrite psi_star[0] with the projected stress. - - if uw.mpi.rank == 0 and verbose: - print(f"VE Stokes solver - store stress and shift history", flush=True) - - import numpy as np - - # Save advected σ* before projection modifies psi_star[0] - _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) - - # Project actual stress into psi_star[0]. - # Uses the constitutive formula so that σ* and σ** from psi_star - # are read correctly during projection. - self.DFDt._psi_star_projection_solver.uw_function = self.constitutive_model.flux - self.DFDt._psi_star_projection_solver.smoothing = 0.0 - self.DFDt._psi_star_projection_solver.solve(verbose=verbose) - - # Now psi_star[0] = projected τ (actual stress from this solve). - # Shift history: psi_star[i] ← previous psi_star[i-1] (advected values). - # psi_star[1] gets the advected σ* (saved before projection). - # Higher levels shift down the chain. - for i in range(self.DFDt.order - 1, 0, -1): - if i == 1: - self.DFDt.psi_star[i].array[...] = _advected_sigma_star - else: - self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] - - # 5. BOOKKEEPING - - self.DFDt.update_post_solve(timestep, verbose=verbose, evalf=evalf) - - self.is_setup = True - self.constitutive_model._solver_is_setup = True - - return - class SNES_Projection(SNES_Scalar): r""" diff --git a/tests/vep_fault_weakening.py b/tests/vep_fault_weakening.py new file mode 100644 index 000000000..e186ce77b --- /dev/null +++ b/tests/vep_fault_weakening.py @@ -0,0 +1,158 @@ +"""VEP shear box with embedded fault — convergence study. + +Horizontal fault at y=0.5 using Surface gaussian influence function. +Runs at two vertical resolutions to check convergence. + +Run: pixi run -e amr-dev python tests/vep_fault_weakening.py +""" + +import time +import numpy as np +import sympy +import underworld3 as uw + +ETA = 1.0 +MU = 1.0 +TAU_Y_FAULT = 0.2 +TAU_Y_BULK = 2.0 +FAULT_WIDTH = 0.08 +DT = 0.1 +NSTEPS = 25 +V_TOP = 0.5 + + +def run_fault_model(res_x, res_y): + """Run the fault model at given resolution, return time series.""" + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res_x, res_y), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=2) + 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) + + fault_points = np.array([[0.0, 0.5, 0.0], [1.0, 0.5, 0.0]]) + fault = uw.meshing.Surface("fault", mesh, fault_points) + fault.discretize() + + tau_y_field = fault.influence_function( + width=FAULT_WIDTH, value_near=TAU_Y_FAULT, value_far=TAU_Y_BULK, profile="gaussian") + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + cm = stokes.constitutive_model + cm.Parameters.shear_viscosity_0 = ETA + cm.Parameters.shear_modulus = MU + cm.Parameters.yield_stress = tau_y_field + cm.Parameters.shear_viscosity_min = ETA * 1.0e-2 + cm.Parameters.strainrate_inv_II_min = 1.0e-10 + # saddle_preconditioner left at default (uses constitutive stiffness) + stokes.tolerance = 1.0e-4 + + 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]) + stokes.petsc_options["ksp_type"] = "fgmres" + + fault_sample = np.array([[0.5, 0.5]]) + bulk_sample = np.array([[0.5, 0.25]]) + + times = [] + fault_stresses = [] + bulk_stresses = [] + converged = [] + + for step in range(NSTEPS): + stokes.solve(timestep=DT, zero_init_guess=False) + t = (step + 1) * DT + times.append(t) + + reason = stokes.snes.getConvergedReason() + converged.append(reason) + + sd = stokes.tau.data + coords = stokes.tau.coords + fault_idx = np.argmin(np.sum((coords - fault_sample)**2, axis=1)) + bulk_idx = np.argmin(np.sum((coords - bulk_sample)**2, axis=1)) + fault_stresses.append(sd[fault_idx, 2]) + bulk_stresses.append(sd[bulk_idx, 2]) + + flag = "" if reason > 0 else f" SNES={reason}" + if (step + 1) % 5 == 0 or step < 3 or reason < 0: + uw.pprint(0, f" [{res_x}x{res_y}] step {step+1}, t={t:.2f}, " + f"fault={fault_stresses[-1]:.4f}, bulk={bulk_stresses[-1]:.4f}{flag}") + + # Cross-section at final step + y_coords = coords[:, 1] + x_coords = coords[:, 0] + near_centre = np.abs(x_coords - 0.5) < 0.1 + y_profile = y_coords[near_centre] + sigma_profile = sd[near_centre, 2] + sort_idx = y_profile.argsort() + + return { + "times": np.array(times), + "fault": np.array(fault_stresses), + "bulk": np.array(bulk_stresses), + "converged": np.array(converged), + "profile_y": y_profile[sort_idx], + "profile_sigma": sigma_profile[sort_idx], + } + + +# --- Run both resolutions --- + +t0 = time.time() +results = {} +for res_x, res_y in [(16, 16), (16, 32)]: + uw.pprint(0, f"\n=== Resolution {res_x}x{res_y} ===") + t1 = time.time() + results[(res_x, res_y)] = run_fault_model(res_x, res_y) + uw.pprint(0, f" done ({time.time()-t1:.0f}s)") + +uw.pprint(0, f"\nTotal: {time.time()-t0:.0f}s") + +# --- Plot --- + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + +t_anal = np.linspace(0.001, NSTEPS * DT, 200) +t_r = ETA / MU +maxwell = ETA * V_TOP * (1 - np.exp(-t_anal / t_r)) + +axes[0].plot(t_anal, maxwell, 'k--', linewidth=1, alpha=0.4, label="Maxwell") +for (rx, ry), r in results.items(): + axes[0].plot(r["times"], r["fault"], linewidth=2, label=f"fault {rx}x{ry}") + axes[0].plot(r["times"], r["bulk"], linewidth=2, linestyle='--', label=f"bulk {rx}x{ry}") + # Mark diverged steps + div_mask = r["converged"] < 0 + if div_mask.any(): + axes[0].plot(r["times"][div_mask], r["fault"][div_mask], 'x', color='red', markersize=6) +axes[0].axhline(TAU_Y_FAULT, color='gray', linestyle=':', alpha=0.5) +axes[0].set_xlabel("Time") +axes[0].set_ylabel(r"$\sigma_{xy}$") +axes[0].set_title("Stress history") +axes[0].legend(fontsize=8) +axes[0].grid(True, alpha=0.3) + +for (rx, ry), r in results.items(): + axes[1].plot(r["profile_y"], r["profile_sigma"], '-o', markersize=2, + linewidth=2, label=f"{rx}x{ry}") +axes[1].axvline(0.5, color='gray', linestyle=':', alpha=0.5, label="fault") +axes[1].set_xlabel("y") +axes[1].set_ylabel(r"$\sigma_{xy}$ (final step)") +axes[1].set_title("Stress profile across fault") +axes[1].legend(fontsize=9) +axes[1].grid(True, alpha=0.3) + +fig.suptitle(f"VEP embedded fault convergence: $\\tau_y$={TAU_Y_FAULT}/{TAU_Y_BULK}, width={FAULT_WIDTH}", + fontsize=12, y=1.02) +fig.tight_layout() +out_path = "/Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/solver-unification/vep_fault.png" +fig.savefig(out_path, dpi=150, bbox_inches='tight') +uw.pprint(0, f"Saved {out_path}") diff --git a/tests/vep_strain_weakening.py b/tests/vep_strain_weakening.py new file mode 100644 index 000000000..db6c3ae26 --- /dev/null +++ b/tests/vep_strain_weakening.py @@ -0,0 +1,150 @@ +"""VEP shear box with strain-weakening yield stress. + +Plastic fraction computed directly from stored data: + eta_ve = eta * mu * dt / (eta + mu * dt) (known scalar) + edot_II_eff = edot_kin + sigma_star / (2 * mu * dt) (from psi_star data) + eta_vep = min(eta_ve, tau_y / (2 * edot_II_eff)) (yield viscosity) + plastic_fraction = max(0, 1 - eta_vep / eta_ve) + +No evaluate, no projection — pure numpy on stored data. + +Run: pixi run -e amr-dev python tests/vep_strain_weakening.py +""" + +import time +import numpy as np +import sympy +import underworld3 as uw + +ETA = 1.0 +MU = 1.0 +TAU_Y0 = 0.3 +TAU_RESIDUAL = 0.1 +EPS_CRIT = 0.5 +DT = 0.1 +NSTEPS = 30 +V_TOP = 0.5 + +t0 = time.time() + +mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=2) +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) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel +cm = stokes.constitutive_model +cm.Parameters.shear_viscosity_0 = ETA +cm.Parameters.shear_modulus = MU +cm.Parameters.yield_stress = TAU_Y0 +cm.Parameters.shear_viscosity_min = ETA * 1.0e-3 +cm.Parameters.strainrate_inv_II_min = 1.0e-10 +stokes.tolerance = 1.0e-4 + +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]) +stokes.petsc_options["ksp_type"] = "fgmres" + +print(f"Setup: {time.time()-t0:.1f}s") +print(f"eta={ETA}, mu={MU}, t_relax={ETA/MU}") +print(f"tau_y0={TAU_Y0}, tau_residual={TAU_RESIDUAL}, eps_crit={EPS_CRIT}") +print() + +times = [] +max_stresses = [] +tau_y_values = [] +eps_p_cum = 0.0 +eps_p_values = [] +edot_p_values = [] + +eta_ve = ETA * MU * DT / (ETA + MU * DT) +edot_kin = V_TOP / 2.0 # tensor shear rate for simple shear + +print(f"{'step':>4} {'t':>5} {'sigma_xy':>10} {'tau_y':>8} {'eps_p':>8} {'edot_p':>8}") +print("-" * 60) + +current_tau_y = TAU_Y0 + +for step in range(NSTEPS): + cm.Parameters.yield_stress.sym = current_tau_y + + t1 = time.time() + stokes.solve(timestep=DT, zero_init_guess=False) + solve_time = time.time() - t1 + + t = (step + 1) * DT + times.append(t) + + sigma_xy = stokes.tau.data[:, 2].max() + max_stresses.append(sigma_xy) + + # Compute plastic fraction from stored data + sigma_star_xy = stokes.DFDt.psi_star[0].data[:, 2].mean() + edot_history = sigma_star_xy / (2 * MU * DT) + edot_II_eff = edot_kin + edot_history + + if edot_II_eff > 0: + eta_vep = min(eta_ve, current_tau_y / (2 * edot_II_eff)) + else: + eta_vep = eta_ve + + plastic_fraction = max(0.0, 1.0 - eta_vep / eta_ve) + edot_p = edot_II_eff * plastic_fraction + edot_p_values.append(edot_p) + + # Accumulate and weaken + eps_p_cum += edot_p * DT + weakening = min(eps_p_cum / EPS_CRIT, 1.0) + current_tau_y = TAU_Y0 + (TAU_RESIDUAL - TAU_Y0) * weakening + + tau_y_values.append(current_tau_y) + eps_p_values.append(eps_p_cum) + + if (step + 1) % 5 == 0 or step < 12: + print(f"{step+1:4d} {t:5.2f} {sigma_xy:10.4f} {current_tau_y:8.4f} {eps_p_cum:8.4f} {edot_p:8.4f} ({solve_time:.1f}s)") + +print() +print(f"Total: {time.time()-t0:.0f}s") + +# --- Plot --- + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True) + +t_anal = np.linspace(0.001, max(times), 200) +t_r = ETA / MU +maxwell = ETA * V_TOP * (1 - np.exp(-t_anal / t_r)) + +axes[0].plot(t_anal, maxwell, 'k--', linewidth=1, alpha=0.4, label="Maxwell (no yield)") +axes[0].plot(times, max_stresses, 'r-', linewidth=2, label=r"$\sigma_{xy}$") +axes[0].plot(times, tau_y_values, 'b--', linewidth=1.5, label=r"$\tau_y(\varepsilon_p)$") +axes[0].axhline(TAU_Y0, color='gray', linestyle=':', alpha=0.5) +axes[0].axhline(TAU_RESIDUAL, color='gray', linestyle='-.', alpha=0.5) +axes[0].set_ylabel("Stress") +axes[0].legend(fontsize=9) +axes[0].grid(True, alpha=0.3) +axes[0].set_title(f"VEP strain weakening: $\\eta$={ETA}, $\\mu$={MU}") + +axes[1].plot(times, edot_p_values, 'm-', linewidth=2) +axes[1].set_ylabel(r"Plastic $\dot{\varepsilon}$") +axes[1].grid(True, alpha=0.3) + +axes[2].plot(times, eps_p_values, 'g-', linewidth=2) +axes[2].axhline(EPS_CRIT, color='gray', linestyle=':', alpha=0.5, label=f"$\\varepsilon_{{crit}}$={EPS_CRIT}") +axes[2].set_xlabel("Time") +axes[2].set_ylabel(r"Accumulated $\varepsilon_p$") +axes[2].legend(fontsize=9) +axes[2].grid(True, alpha=0.3) + +fig.tight_layout() +out_path = "/Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/solver-unification/vep_strain_weakening.png" +fig.savefig(out_path, dpi=150) +print(f"Saved {out_path}") diff --git a/tests/vep_timedep_yield.py b/tests/vep_timedep_yield.py new file mode 100644 index 000000000..4b9cc5af7 --- /dev/null +++ b/tests/vep_timedep_yield.py @@ -0,0 +1,137 @@ +"""VEP with time-dependent yield stress (no feedback). + +Prescribed tau_y(t) that decreases linearly over time. +VE stress builds, hits tau_y, then tau_y drops and stress follows. + +No strain accumulation, no projection of viscosity ratios — +just the solver responding to a changing yield stress parameter. + +Run: pixi run -e amr-dev python tests/vep_timedep_yield.py +""" + +import time +import numpy as np +import sympy +import underworld3 as uw + +ETA = 1.0 +MU = 1.0 +V_TOP = 0.5 +DT = 0.1 +NSTEPS = 25 + +# tau_y schedule: starts above Maxwell steady state, drops below it +TAU_Y_START = 1.5 # well above Maxwell steady state (2*eta*edot = 1.0) +TAU_Y_END = 0.2 +TAU_Y_DROP_START = 0.8 # start dropping at t=0.8 +TAU_Y_DROP_END = 1.8 # reach minimum at t=1.8 + +t0 = time.time() + +mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4), + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + qdegree=2, +) + +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) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel +cm = stokes.constitutive_model +cm.Parameters.shear_viscosity_0 = ETA +cm.Parameters.shear_modulus = MU +cm.Parameters.shear_viscosity_min = ETA * 1.0e-3 +cm.Parameters.strainrate_inv_II_min = 1.0e-10 +stokes.saddle_preconditioner = 1.0 +stokes.tolerance = 1.0e-4 + +# Start with high yield stress +cm.Parameters.yield_stress = TAU_Y_START + +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]) +stokes.petsc_options["ksp_type"] = "fgmres" + +print(f"Setup: {time.time()-t0:.1f}s") +print(f"eta={ETA}, mu={MU}, t_relax={ETA/MU}") +print(f"tau_y: {TAU_Y_START} -> {TAU_Y_END} (t={TAU_Y_DROP_START}..{TAU_Y_DROP_END})") +print(f"Maxwell steady state: eta*gamma_dot = {ETA*V_TOP}") +print() + +times = [] +max_stresses = [] +tau_y_values = [] + +print(f"{'step':>4} {'t':>5} {'sigma_xy':>10} {'tau_y':>8} {'phase':>10}") +print("-" * 50) + +for step in range(NSTEPS): + t = (step + 1) * DT + + # Update tau_y for this step + if t < TAU_Y_DROP_START: + current_tau_y = TAU_Y_START + elif t > TAU_Y_DROP_END: + current_tau_y = TAU_Y_END + else: + frac = (t - TAU_Y_DROP_START) / (TAU_Y_DROP_END - TAU_Y_DROP_START) + current_tau_y = TAU_Y_START + (TAU_Y_END - TAU_Y_START) * frac + + cm.Parameters.yield_stress = current_tau_y + + t1 = time.time() + stokes.solve(timestep=DT, zero_init_guess=False) + solve_time = time.time() - t1 + + times.append(t) + tau_y_values.append(current_tau_y) + + sd = stokes.tau.data + sigma_xy = sd[:, 2].max() + max_stresses.append(sigma_xy) + + # Determine phase + if sigma_xy > 0.95 * current_tau_y: + phase = "yield" + else: + phase = "elastic" + + print(f"{step+1:4d} {t:5.2f} {sigma_xy:10.4f} {current_tau_y:8.3f} {phase:>10} ({solve_time:.1f}s)") + +print() +print(f"Total: {time.time()-t0:.0f}s") + +# --- Plot --- + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +fig, ax = plt.subplots(figsize=(8, 5)) + +t_anal = np.linspace(0.001, max(times), 200) +t_r = ETA / MU +# Maxwell: sigma_xy = eta * gamma_dot * (1 - exp(-t/t_r)) +# gamma_dot = V_TOP / H = V_TOP (H=1) +maxwell = ETA * V_TOP * (1 - np.exp(-t_anal / t_r)) + +ax.plot(t_anal, maxwell, 'k--', linewidth=1, alpha=0.4, label="Maxwell (no yield)") +ax.plot(times, max_stresses, 'r-o', linewidth=2, markersize=3, label=r"$\sigma_{xy}$") +ax.plot(times, tau_y_values, 'b--', linewidth=1.5, label=r"$\tau_y(t)$") +ax.set_xlabel("Time") +ax.set_ylabel("Stress") +ax.set_title(f"VEP with prescribed weakening: $\\eta$={ETA}, $\\mu$={MU}") +ax.legend(fontsize=10) +ax.grid(True, alpha=0.3) + +fig.tight_layout() +out_path = "/Users/lmoresi/+Underworld/underworld3-pixi/.claude/worktrees/solver-unification/vep_timedep_yield.png" +fig.savefig(out_path, dpi=150) +print(f"Saved {out_path}") diff --git a/vep_fault.png b/vep_fault.png new file mode 100644 index 000000000..690ffba4b Binary files /dev/null and b/vep_fault.png differ diff --git a/vep_strain_weakening.png b/vep_strain_weakening.png new file mode 100644 index 000000000..66cf748cb Binary files /dev/null and b/vep_strain_weakening.png differ diff --git a/vep_timedep_yield.png b/vep_timedep_yield.png new file mode 100644 index 000000000..0203648f4 Binary files /dev/null and b/vep_timedep_yield.png differ