diff --git a/docs/advanced/benchmarks/index.md b/docs/advanced/benchmarks/index.md new file mode 100644 index 00000000..774102ea --- /dev/null +++ b/docs/advanced/benchmarks/index.md @@ -0,0 +1,13 @@ +--- +title: "Benchmarks" +--- + +# Solver Benchmarks + +Validation benchmarks comparing Underworld3 solvers against analytical solutions. + +```{toctree} +:maxdepth: 1 + +ve-oscillatory-shear +``` diff --git a/docs/advanced/benchmarks/ve-oscillatory-shear.md b/docs/advanced/benchmarks/ve-oscillatory-shear.md new file mode 100644 index 00000000..03e95c43 --- /dev/null +++ b/docs/advanced/benchmarks/ve-oscillatory-shear.md @@ -0,0 +1,108 @@ +--- +title: "Viscoelastic Oscillatory Shear Benchmark" +--- + +# Maxwell Oscillatory Shear + +This benchmark validates the viscoelastic Stokes solver against the analytical +solution for a Maxwell material under oscillatory simple shear. + +## Problem Setup + +A box with height $H$ and width $2H$ is sheared by imposing time-dependent +velocities on the top and bottom boundaries: + +$$v_x(y=\pm H/2, t) = \pm V_0 \sin(\omega t)$$ + +The left and right boundaries are free-slip (no vertical velocity). The shear +rate is $\dot\gamma(t) = \dot\gamma_0 \sin(\omega t)$ where $\dot\gamma_0 = 2V_0/H$. + +## Analytical Solution + +The Maxwell constitutive law gives: + +$$\dot\sigma_{xy} + \frac{\sigma_{xy}}{t_r} = \mu \dot\gamma_0 \sin(\omega t)$$ + +where $t_r = \eta/\mu$ is the relaxation time. With $\sigma(0) = 0$, the full +solution (including the startup transient) is: + +$$\sigma_{xy}(t) = \frac{\eta \dot\gamma_0}{1 + \text{De}^2} +\left[\sin(\omega t) - \text{De}\cos(\omega t) + \text{De}\,e^{-t/t_r}\right]$$ + +where $\text{De} = \omega t_r$ is the Deborah number. + +**Steady-state properties** (after transient decays): + +- Amplitude: $A = \eta \dot\gamma_0 / \sqrt{1 + \text{De}^2}$ +- Phase lag: $\delta = \arctan(\text{De})$ + +At $\text{De} = 0$ (viscous limit): $A = \eta\dot\gamma_0$, $\delta = 0$. +At $\text{De} \to \infty$ (elastic limit): $A \to 0$, $\delta \to 90°$. + +## Convergence with BDF Order + +The VE solver uses BDF-$k$ time integration ($k = 1, 2, 3$). The convergence +study (constant shear, $\text{De} = 1$) shows: + +| $\Delta t / t_r$ | BDF-1 error | BDF-2 error | BDF-3 error | +|-------------------|-------------|-------------|-------------| +| 0.200 | 3.0e-02 | 4.0e-03 | 6.4e-03 | +| 0.100 | 1.5e-02 | 9.3e-04 | 1.7e-03 | +| 0.050 | 7.8e-03 | 2.3e-04 | 4.3e-04 | +| 0.020 | 3.1e-03 | 3.7e-05 | — | + +BDF-2 achieves second-order convergence (~4x error reduction per halving) and +is the recommended default. BDF-1 is first-order. BDF-3 converges at nearly +second order but with a larger error constant. + +## Resolution Study (Oscillatory, De = 5) + +At high Deborah number, the oscillation period is short relative to the +relaxation time, requiring fine time resolution. The plot below shows the +effect of timestep size at $\text{De} = 5$ ($\omega t_r = 5$, phase lag = 79°): + +- **63 pts/period** ($\Delta t/t_r = 0.02$): both orders match analytical +- **31 pts/period** ($\Delta t/t_r = 0.04$): O1 shows slight amplitude reduction, O2 still accurate +- **16 pts/period** ($\Delta t/t_r = 0.08$): O1 amplitude visibly damped, O2 remains good + +```{note} +The amplitude reduction at coarse timesteps is numerical dissipation from +the BDF-1 discrete transfer function, not a cumulative error. The discrete +steady-state amplitude is a fixed fraction of the analytical amplitude, +determined by $\omega \Delta t$. +``` + +## Running the Benchmarks + +```bash +# Oscillatory validation (De=1.5, order 1 and 2) +python tests/plot_ve_oscillatory_validation.py + +# Resolution study (De=5, three timestep sizes) +# Saves .npz data files for re-analysis +python tests/plot_ve_oscillatory_validation.py + +# Replot from saved data (no re-running) +python tests/plot_ve_oscillatory_validation.py --replot +``` + +## Notes on `dt_elastic` + +The parameter `dt_elastic` on the constitutive model is the elastic relaxation +timescale used in the BDF discretisation. It controls the effective viscosity +$\eta_{\text{eff}}$ and the stress history weighting: + +- BDF-1: $\eta_{\text{eff}} = \eta\mu\Delta t_e / (\eta + \mu\Delta t_e)$ +- BDF-2: $\eta_{\text{eff}} = 2\eta\mu\Delta t_e / (3\eta + 2\mu\Delta t_e)$ +- BDF-3: $\eta_{\text{eff}} = 6\eta\mu\Delta t_e / (11\eta + 6\mu\Delta t_e)$ + +When `timestep` is passed to `VE_Stokes.solve()`, it controls the advection +step for semi-Lagrangian history transport. It does **not** overwrite +`dt_elastic` — these are independent parameters. + +A running-average approach for accumulating history when $\Delta t \ll \Delta t_e$ +was investigated but found to be extremely diffusive for semi-Lagrangian +transport and is not implemented. To prevent runaway or unstable behaviour +when timesteps become small (e.g. due to CFL constraints or failure events), +we advise limiting the minimum effective viscosity, in line with the physics +of the problem. diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 8d02acde..47293a68 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -48,6 +48,7 @@ from underworld3.swarm import IndexSwarmVariable from underworld3.discretisation import MeshVariable from underworld3.systems.ddt import SemiLagrangian as SemiLagrangian_DDt +from underworld3.systems.ddt import _bdf_coefficients from underworld3.function.quantities import UWQuantity from underworld3.systems.ddt import Lagrangian as Lagrangian_DDt @@ -971,12 +972,12 @@ def viscosity(self): # Keep this as an sub-expression for clarity if inner_self.shear_viscosity_min.sym != -sympy.oo: - self._plastic_eff_viscosity._sym = sympy.simplify( - sympy.Max(effective_viscosity, inner_self.shear_viscosity_min) + self._plastic_eff_viscosity._sym = sympy.Max( + effective_viscosity, inner_self.shear_viscosity_min ) else: - self._plastic_eff_viscosity._sym = sympy.simplify(effective_viscosity) + self._plastic_eff_viscosity._sym = effective_viscosity # Returns an expression that has a different description return self._plastic_eff_viscosity @@ -1086,6 +1087,14 @@ def __init__(self, unknowns, order=1, material_name: str = None): self._order = order + # BDF coefficients as UWexpressions — route through PetscDS constants[]. + # Updated each step by _update_bdf_coefficients() before solve. + # Initialised to BDF-1 values: [1, -1, 0, 0]. + 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) @@ -1198,32 +1207,15 @@ def ve_effective_viscosity(inner_self): if inner_self.shear_modulus == sympy.oo: return inner_self.shear_viscosity_0 - # Note, 1st order only here but we should add higher order versions of this - - # 1st Order version (default) - if inner_self._owning_model.order != 2: - el_eff_visc = ( - inner_self.shear_viscosity_0 - * inner_self.shear_modulus - * inner_self.dt_elastic - / ( - inner_self.shear_viscosity_0 - + inner_self.dt_elastic * inner_self.shear_modulus - ) - ) + # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt) + # c0 is a UWexpression routed through PetscDS constants[], + # updated each step by _update_bdf_coefficients(). + eta = inner_self.shear_viscosity_0 + mu = inner_self.shear_modulus + dt_e = inner_self.dt_elastic + c0 = inner_self._owning_model._bdf_c0 - # 2nd Order version (need to ask for this one) - else: - el_eff_visc = ( - 2 - * inner_self.shear_viscosity_0 - * inner_self.shear_modulus - * inner_self.dt_elastic - / ( - 3 * inner_self.shear_viscosity_0 - + 2 * inner_self.dt_elastic * inner_self.shear_modulus - ) - ) + el_eff_visc = eta * mu * dt_e / (c0 * eta + mu * dt_e) inner_self._ve_effective_viscosity.sym = el_eff_visc @@ -1252,6 +1244,46 @@ def order(self, value): self._reset() return + @property + def effective_order(self): + """Effective order accounting for DDt history startup. + + During the first few timesteps, the DDt may not have enough history + to support the requested order. This property returns the lower of + the requested order and the DDt's effective order (which ramps from + 1 to self.order as history accumulates). + """ + if self.Unknowns is not None and self.Unknowns.DFDt is not None: + return min(self._order, self.Unknowns.DFDt.effective_order) + return self._order + + def _update_bdf_coefficients(self): + """Update BDF coefficient UWexpressions from current dt_elastic and DDt history. + + Call this before each solve so that the constants[] array carries the + 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. + """ + 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 + ) + else: + coeffs = _bdf_coefficients(self.effective_order, None, []) + + # Pad to length 4 + 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] + # The following should have no setters @property def stress_star(self): @@ -1287,20 +1319,13 @@ def E_eff(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.order != 2: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - E += stress_star / ( - 2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus - ) + mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus + # BDF history coefficients as UWexpressions (route through constants[]) + bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] - else: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - E += stress_star / ( - self.Parameters.dt_elastic * self.Parameters.shear_modulus - ) - stress_2star / ( - 4 * self.Parameters.dt_elastic * self.Parameters.shear_modulus - ) + # History contribution: -Σ cᵢ·σ_star[i-1] / (2·μ·dt) + 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 @@ -1365,16 +1390,10 @@ def _plastic_effective_viscosity(self): if parameters.yield_stress == sympy.oo: return sympy.oo - Edot = self.Unknowns.E - if self.Unknowns.DFDt is not None: - - ## First order ... - stress_star = self.Unknowns.DFDt.psi_star[0] - - if self.is_elastic: - Edot += stress_star.sym / ( - 2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus - ) + # Use the effective strain rate (including elastic history) for the + # yield criterion. This must use the same order-dependent BDF + # coefficients as the stress formula. + Edot = self.E_eff.sym strainrate_inv_II = expression( R"{\dot\varepsilon_{II}'}", @@ -1468,8 +1487,6 @@ def flux(self): # plastic_scale_factor = sympy.Max(1, self.plastic_overshoot()) # stress /= plastic_scale_factor - stress = sympy.simplify(stress) - return stress def stress_projection(self): @@ -1492,8 +1509,6 @@ def stress_projection(self): / (self.Parameters.dt_elastic * self.Parameters.shear_modulus) ) - stress = sympy.simplify(stress) - return stress def stress(self): @@ -1508,34 +1523,15 @@ def stress(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.order != 2: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress += ( - 2 - * self.viscosity - * ( - stress_star - / (2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus) - ) - ) + mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus + bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] - else: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - - stress += ( - 2 - * self.viscosity - * ( - stress_star - / (self.Parameters.dt_elastic * self.Parameters.shear_modulus) - - stress_2star - / (4 * self.Parameters.dt_elastic * self.Parameters.shear_modulus) - ) + # 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) ) - stress = sympy.simplify(stress) - return stress # def eff_edot(self): diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h index d502dd92..f2e5cc10 100644 --- a/src/underworld3/cython/petsc_compat.h +++ b/src/underworld3/cython/petsc_compat.h @@ -129,6 +129,21 @@ PetscErrorCode UW_PetscDSViewBdWF(PetscDS ds, PetscInt bd) return 1; } +// Set the time value on a DM. This is passed as `petsc_t` to all +// pointwise residual and Jacobian functions during assembly. +// PETSc stores this internally but petsc4py doesn't expose it. +PetscErrorCode UW_DMSetTime(DM dm, PetscReal time) +{ + // DMSetOutputSequenceNumber stores (step, time) on the DM. + // The time component is what DMPlexComputeResidual_Internal + // passes as petsc_t to the pointwise functions. + PetscInt step; + PetscReal old_time; + PetscCall(DMGetOutputSequenceNumber(dm, &step, &old_time)); + PetscCall(DMSetOutputSequenceNumber(dm, step, time)); + return PETSC_SUCCESS; +} + PetscErrorCode UW_DMPlexSetSNESLocalFEM(DM dm, PetscBool flag, void *ctx) { diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi index c39d9c71..edd73ca1 100644 --- a/src/underworld3/cython/petsc_extras.pxi +++ b/src/underworld3/cython/petsc_extras.pxi @@ -42,6 +42,7 @@ cdef extern from "petsc_compat.h": PetscErrorCode UW_PetscDSSetBdTerms (PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*, void*, void*, void*, void* ) PetscErrorCode UW_PetscDSViewWF(PetscDS) PetscErrorCode UW_PetscDSViewBdWF(PetscDS, PetscInt) + PetscErrorCode UW_DMSetTime( PetscDM, PetscReal ) PetscErrorCode UW_DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *) PetscErrorCode UW_DMPlexComputeBdIntegral( PetscDM, PetscVec, PetscDMLabel, PetscInt, const PetscInt*, void*, PetscScalar*, void*) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 1f9cebea..4c0da8a8 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -951,6 +951,122 @@ class SolverBaseClass(uw_object): + @property + def tau(self): + r"""Computed flux from the constitutive model, projected onto the mesh. + + Returns a :class:`~underworld3.discretisation.MeshVariable` containing + the flux that the solver actually computed. The variable is created + lazily on first access and updated by solving a projection each time + the property is read. + + For Stokes-family solvers, this is the deviatoric stress tensor + :math:`\boldsymbol{\tau}`. For Poisson/diffusion, it is the + diffusive flux :math:`\kappa \nabla u`. The projection handles + derivative evaluation internally — users should use this property + rather than manually reconstructing the flux from the constitutive + formula. + + Subclasses may override this to return pre-computed values (e.g. + VE_Stokes returns the stress stored during the solve). + + Returns + ------- + MeshVariable + Mesh variable containing the projected flux values. + Access numerical data via ``.data`` or ``.array``, symbolic + expression via ``.sym``. + """ + + if self._constitutive_model is None: + raise RuntimeError( + "No constitutive model set. Assign solver.constitutive_model first." + ) + + # Lazy initialization of projection infrastructure + if not hasattr(self, '_tau_var') or self._tau_var is None: + self._setup_tau_projection() + + # Project the constitutive flux onto the mesh variable + # Transpose if needed: flux may be (dim,1) column but projector expects row + flux = self._constitutive_model.flux + if hasattr(flux, 'shape') and flux.shape[1] == 1 and flux.shape[0] > 1: + flux = flux.T + self._tau_projector.uw_function = flux + self._tau_projector.smoothing = 0.0 + self._tau_projector.solve() + + return self._tau_var + + def _setup_tau_projection(self): + """Create the mesh variable and projector for tau (lazy init).""" + + flux = self._constitutive_model.flux + dim = self.mesh.dim + rows, cols = flux.shape + + # Determine variable type and create appropriate projection solver + if rows == cols and rows == dim: + # Tensor flux (Stokes stress) + self._tau_var = uw.discretisation.MeshVariable( + f"tau_{self.instance_number}", + self.mesh, + (dim, dim), + vtype=uw.VarType.SYM_TENSOR, + degree=self.u.degree, + continuous=True, + varsymbol=r"{\tau}", + ) + _work = uw.discretisation.MeshVariable( + f"tau_work_{self.instance_number}", + self.mesh, + 1, + degree=self.u.degree, + continuous=True, + ) + from underworld3.systems.solvers import SNES_Tensor_Projection + self._tau_projector = SNES_Tensor_Projection( + self.mesh, self._tau_var, _work, verbose=False + ) + + elif rows == dim and cols == 1: + # Vector flux (Poisson heat flux) + self._tau_var = uw.discretisation.MeshVariable( + f"tau_{self.instance_number}", + self.mesh, + dim, + vtype=uw.VarType.VECTOR, + degree=self.u.degree, + continuous=True, + varsymbol=r"{\mathbf{q}}", + ) + from underworld3.systems.solvers import SNES_Vector_Projection + self._tau_projector = SNES_Vector_Projection( + self.mesh, self._tau_var, verbose=False + ) + + elif cols == dim and rows == 1: + # Transposed vector flux + self._tau_var = uw.discretisation.MeshVariable( + f"tau_{self.instance_number}", + self.mesh, + dim, + vtype=uw.VarType.VECTOR, + degree=self.u.degree, + continuous=True, + varsymbol=r"{\mathbf{q}}", + ) + from underworld3.systems.solvers import SNES_Vector_Projection + self._tau_projector = SNES_Vector_Projection( + self.mesh, self._tau_var, verbose=False + ) + + else: + raise RuntimeError( + f"Cannot create tau projection for flux shape {flux.shape}. " + f"Expected ({dim},{dim}) tensor, ({dim},1) or (1,{dim}) vector." + ) + def validate_solver(self): """ Checks to see if the required properties have been set. @@ -1659,7 +1775,8 @@ class SNES_Scalar(SolverBaseClass): _force_setup: bool =False, verbose: bool=False, debug: bool=False, - debug_name: str=None ): + debug_name: str=None, + time=None, ): """ Solve the system of equations. @@ -1682,6 +1799,11 @@ class SNES_Scalar(SolverBaseClass): Enable debug output including intermediate residuals. debug_name : str, optional Name prefix for debug output files. + time : float or Quantity, optional + Physical time for this solve. Passed as ``petsc_t`` to all + pointwise functions. Expressions using ``mesh.t`` evaluate at + this time. Non-dimensionalised when scaling is active. + Default: None (petsc_t unchanged). Returns ------- @@ -1719,6 +1841,16 @@ class SNES_Scalar(SolverBaseClass): self._build(verbose, debug, debug_name) + # Set time on the DM so petsc_t is available in pointwise functions + cdef DM _time_dm + if time is not None: + if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'): + t_nd = float(uw.non_dimensionalise(time)) + else: + t_nd = float(time) + _time_dm = self.dm + UW_DMSetTime(_time_dm.dm, t_nd) + gvec = self.dm.getGlobalVec() if not zero_init_guess: @@ -3755,7 +3887,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): verbose=False, debug=False, debug_name=None, - _force_setup: bool =False, ): + _force_setup: bool =False, + time=None, ): """ Solve the Stokes system for velocity and pressure. @@ -3782,6 +3915,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Name prefix for debug output files. _force_setup : bool, default=False Force rebuild of the solver even if already set up. + time : float or Quantity, optional + Physical time for this solve. Passed as ``petsc_t`` to all + pointwise residual and Jacobian functions. Expressions using + ``mesh.t`` will evaluate at this time. Non-dimensionalised + automatically when scaling is active. Default: None (petsc_t=0). Returns ------- @@ -3823,6 +3961,17 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._build(verbose, debug, debug_name) + # Set time on the DM so petsc_t is available in pointwise functions. + # Non-dimensionalise if the scaling system is active. + cdef DM _time_dm_stokes + if time is not None: + if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'): + t_nd = float(uw.non_dimensionalise(time)) + else: + t_nd = float(time) + _time_dm_stokes = self.dm + UW_DMSetTime(_time_dm_stokes.dm, t_nd) + # Keep a record of these set-up parameters tolerance = self.tolerance snes_type = self.snes.getType() diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index c111d0fd..4b11dbd6 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -629,6 +629,14 @@ def mesh_update_callback(array, change_context): self._Gamma.y._ccodestr = "petsc_n[1]" self._Gamma.z._ccodestr = "petsc_n[2]" + # Time coordinate — PETSc passes this as petsc_t to all pointwise + # functions. Solvers set dm.time before each solve via solve(time=t). + # Users reference it as mesh.t in expressions (e.g. V0 * sympy.sin(omega * mesh.t)) + from ..utilities.unit_aware_coordinates import TimeSymbol + + self._t = TimeSymbol("t") + self._t._units = None # patched below by _patch_time_units + # Add unit awareness to coordinate symbols if mesh has units or model has scales from ..utilities.unit_aware_coordinates import patch_coordinate_units @@ -1497,6 +1505,32 @@ def CoordinateSystem(self) -> CoordinateSystem: r"""Alias for :attr:`X` (the coordinate system object).""" return self._CoordinateSystem + @property + def t(self): + r"""Symbolic time coordinate. + + PETSc passes a time value (``petsc_t``) to all pointwise residual + and Jacobian functions. Use ``mesh.t`` in expressions to reference + this time without forcing JIT recompilation each timestep. + + The low-level PETSc solver accepts ``time=t`` to set the value + of ``petsc_t`` for pointwise functions. If not provided, ``petsc_t`` + defaults to 0. Note: the high-level Python ``solve()`` wrappers + do not yet pass ``time=`` through — set it directly via + ``UW_DMSetTime`` at the Cython level if needed. + + When the scaling system is active, ``mesh.t`` carries time units + (derived from the model's time scale) so that dimensional analysis + works correctly in expressions. + + Examples + -------- + >>> omega = 2 * np.pi / period + >>> stokes.add_dirichlet_bc((V0 * sympy.sin(omega * mesh.t), 0.0), "Top") + >>> stokes.solve(time=current_time) # sets petsc_t before SNES + """ + return self._t + @property def r(self) -> Tuple[sympy.vector.BaseScalar]: r"""Tuple of coordinate scalars :math:`(x, y)` or :math:`(x, y, z)`. diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index d9a0b825..98f218e3 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -548,6 +548,12 @@ def _project_to_work_variable(expr, mesh, smoothing=1e-6): import underworld3 as uw # Handle matrix expressions - need multi-component work variable + # TODO(BUG): This fails for dim×dim matrices (e.g. 2×2 stress tensor) + # because MeshVariable can't infer vtype from num_components=4 (flat int). + # Needs: pass num_components=(rows,cols) with vtype=TENSOR, and use + # Tensor_Projection instead of per-component scalar Projection. + # Currently, callers like SemiLagrangian.update_pre_solve fall back to + # their own projection solver via except-clause when this fails. if hasattr(expr, 'shape') and expr.shape != (1, 1): rows, cols = expr.shape n_components = rows * cols diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 89e44827..400c33b8 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -331,7 +331,9 @@ def effective_order(self): startup, ``effective_order`` ramps from 1 to ``self.order`` as successive solves populate the history slots with distinct values. """ - return min(self.order, max(1, self._n_solves_completed + 1)) + # BDF-k requires k completed solves to have k distinct history values. + # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves. + return min(self.order, max(1, self._n_solves_completed)) def update_history_fn(self): r"""Copy current :math:`\psi` to the first history slot ``psi_star[0]``.""" @@ -406,6 +408,11 @@ def update_post_solve( return + @property + def bdf_coefficients(self): + """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" + return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + def bdf(self, order: Optional[int] = None): r"""Compute the backward differentiation approximation of the time-derivative of ψ. For order 1: bdf ≡ ψ - psi_star[0] @@ -679,7 +686,9 @@ def effective_order(self): startup, ``effective_order`` ramps from 1 to ``self.order`` as successive solves populate the history slots with distinct values. """ - return min(self.order, max(1, self._n_solves_completed + 1)) + # BDF-k requires k completed solves to have k distinct history values. + # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves. + return min(self.order, max(1, self._n_solves_completed)) def update_history_fn(self): r"""Copy current :math:`\psi` to ``psi_star[0]`` via evaluation or projection.""" @@ -802,6 +811,11 @@ def update_post_solve( return + @property + def bdf_coefficients(self): + """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" + return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + def bdf(self, order=None): r"""Backwards differentiation form for calculating DuDt. @@ -1113,7 +1127,9 @@ def effective_order(self): startup, ``effective_order`` ramps from 1 to ``self.order`` as successive solves populate the history slots with distinct values. """ - return min(self.order, max(1, self._n_solves_completed + 1)) + # BDF-k requires k completed solves to have k distinct history values. + # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves. + return min(self.order, max(1, self._n_solves_completed)) def initialise_history(self): r"""Initialize all history slots to the current value of :math:`\psi`. @@ -1184,6 +1200,7 @@ def update_post_solve( if self._n_solves_completed < self.order: self._n_solves_completed += 1 + return def update_pre_solve( @@ -1192,11 +1209,22 @@ def update_pre_solve( evalf: Optional[bool] = False, verbose: Optional[bool] = False, dt_physical: Optional[float] = None, + store_result: Optional[bool] = True, ): """Sample upstream values along characteristics before solve. On the first call, automatically initialises history from the current field values so that bdf() returns zero on the first step. + + Parameters + ---------- + store_result : bool, optional + If True (default), evaluate psi_fn at current positions and store + in psi_star[0] before advection — the standard DDt behaviour. + If False, skip this step and the history shift: only advect the + existing psi_star levels upstream. Used by VE_Stokes where + psi_star[0] already contains the projected actual stress from + the previous solve. """ if not self._history_initialised: @@ -1224,15 +1252,19 @@ def update_pre_solve( else: phi = sympy.sympify(1) - for i in range(self.order - 1, 0, -1): - self.psi_star[i].array[...] = ( - phi * self.psi_star[i - 1].array[...] + (1 - phi) * self.psi_star[i].array[...] - ) + if store_result: + for i in range(self.order - 1, 0, -1): + self.psi_star[i].array[...] = ( + phi * self.psi_star[i - 1].array[...] + (1 - phi) * self.psi_star[i].array[...] + ) # 2. Compute the current value of psi_fn which we store in psi_star[0] # Note the need to do a try/except to handle unsupported evaluations # (e.g. of derivatives) # + # When store_result=False (e.g. VE stress history), skip this step — + # psi_star[0] already contains the projected actual stress from + # the previous solve and we want to advect *that*, not the flux. # CRITICAL FIX (2025-11-28): Handle coordinates correctly for unit-aware mode. # Previous bug: extracting .magnitude gives METERS (e.g., 1000000), but: @@ -1271,28 +1303,29 @@ def update_pre_solve( :, : ] - try: - # Use shifted ND coords to avoid quad mesh boundary issues - # node_coords_nd is slightly shifted toward cell centroids (lines 703-709) - # evaluate() treats plain numpy as ND [0-1] coordinates - eval_result = uw.function.evaluate( - self.psi_fn, - node_coords_nd, - evalf=evalf, - ) - # Wrap result with units if psi_star has units but eval didn't return UnitAwareArray - psi_star_units = self.psi_star[0].units - if psi_star_units is not None and not isinstance(eval_result, UnitAwareArray): - eval_result = UnitAwareArray(eval_result, units=psi_star_units) + if store_result: + try: + # Use shifted ND coords to avoid quad mesh boundary issues + # node_coords_nd is slightly shifted toward cell centroids + # evaluate() treats plain numpy as ND [0-1] coordinates + eval_result = uw.function.evaluate( + self.psi_fn, + node_coords_nd, + evalf=evalf, + ) + # Wrap result with units if psi_star has units but eval didn't return UnitAwareArray + psi_star_units = self.psi_star[0].units + if psi_star_units is not None and not isinstance(eval_result, UnitAwareArray): + eval_result = UnitAwareArray(eval_result, units=psi_star_units) - self.psi_star[0].array[...] = eval_result + self.psi_star[0].array[...] = eval_result - except Exception: - # Fallback to projection solver for expressions that can't be directly evaluated - # (e.g., containing derivatives) - self._psi_star_projection_solver.uw_function = self.psi_fn - self._psi_star_projection_solver.smoothing = 0.0 - self._psi_star_projection_solver.solve(verbose=verbose) + except Exception: + # Fallback to projection solver for expressions that can't be directly evaluated + # (e.g., containing derivatives) + self._psi_star_projection_solver.uw_function = self.psi_fn + self._psi_star_projection_solver.smoothing = 0.0 + self._psi_star_projection_solver.solve(verbose=verbose) # 3. Compute the upstream values from the psi_fn @@ -1536,6 +1569,128 @@ def update_pre_solve( return + def _advect_history_PARKED(self, dt, evalf=False, verbose=False): + """Advect all psi_star history levels to upstream positions. + + Pure advection only — no copy-down, no psi_fn evaluation. + Each psi_star[i] is sampled at positions traced back by Δt along + the velocity characteristics, and the result is written back in place. + + Used by VE_Stokes for explicit stress history management where the + actual stress is stored in psi_star[0] after each solve. This method + advects those stored values to upstream positions before the next solve. + + Parameters + ---------- + dt : float or sympy expression + Timestep for characteristic tracing. + evalf : bool + Force numerical evaluation. + verbose : bool + Verbose output. + """ + + if not self._history_initialised: + self.initialise_history() + + # --- Coordinate setup (shared with update_pre_solve) --- + from underworld3.utilities.unit_aware_array import UnitAwareArray + + psi_star_0_coords = self.psi_star[0].coords + if hasattr(psi_star_0_coords, "magnitude"): + psi_star_0_coords_nd = uw.non_dimensionalise(psi_star_0_coords) + if isinstance(psi_star_0_coords_nd, UnitAwareArray): + psi_star_0_coords_nd = np.array(psi_star_0_coords_nd) + elif hasattr(psi_star_0_coords_nd, 'magnitude'): + psi_star_0_coords_nd = psi_star_0_coords_nd.magnitude + else: + psi_star_0_coords_nd = psi_star_0_coords + + cellid = self.mesh.get_closest_cells(psi_star_0_coords_nd) + centroid_coords = self.mesh._centroids[cellid] + shift = 0.001 + node_coords_nd = (1.0 - shift) * psi_star_0_coords_nd + shift * centroid_coords + + # --- Unit handling for dt --- + model = uw.get_default_model() + coords_template = self.psi_star[0].coords + has_units = hasattr(coords_template, "magnitude") or hasattr(coords_template, "_magnitude") + + if has_units: + dt_for_calc = dt.to("second") if hasattr(dt, "to") else dt + else: + if hasattr(dt, "magnitude") or hasattr(dt, "value"): + dt_nondim = uw.non_dimensionalise(dt, model) + if hasattr(dt_nondim, "magnitude"): + dt_for_calc = float(dt_nondim.magnitude) + elif hasattr(dt_nondim, "value"): + dt_for_calc = float(dt_nondim.value) + else: + dt_for_calc = float(dt_nondim) + else: + dt_for_calc = dt + + # --- Advect each history level (RK2 midpoint method) --- + for i in range(self.order - 1, -1, -1): + # Evaluate velocity at node positions + v_result = uw.function.evaluate(self.V_fn, node_coords_nd) + if isinstance(v_result, UnitAwareArray): + v_at_node_pts = v_result[:, 0, :] + if not isinstance(v_at_node_pts, UnitAwareArray): + v_at_node_pts = UnitAwareArray(v_at_node_pts, units=v_result.units) + else: + v_at_node_pts = v_result[:, 0, :] + + # Non-dimensionalize velocities if needed + if not has_units and isinstance(v_at_node_pts, UnitAwareArray): + v_nondim = uw.non_dimensionalise(v_at_node_pts, model) + v_at_node_pts = np.array(v_nondim) if isinstance(v_nondim, UnitAwareArray) else v_nondim + + coords = self.psi_star[i].coords + if not has_units and isinstance(coords, UnitAwareArray): + coords = np.array(coords) + + # RK2: midpoint velocity + mid_pt_coords = coords - v_at_node_pts * (0.5 * dt_for_calc) + mid_pt_coords = self.mesh.return_coords_to_bounds(mid_pt_coords) + v_mid_result = uw.function.global_evaluate(self.V_fn, mid_pt_coords) + if isinstance(v_mid_result, UnitAwareArray): + v_at_mid_pts = v_mid_result[:, 0, :] + if not isinstance(v_at_mid_pts, UnitAwareArray): + v_at_mid_pts = UnitAwareArray(v_at_mid_pts, units=v_mid_result.units) + else: + v_at_mid_pts = v_mid_result[:, 0, :] + + if not has_units and isinstance(v_at_mid_pts, UnitAwareArray): + v_nondim = uw.non_dimensionalise(v_at_mid_pts, model) + v_at_mid_pts = np.array(v_nondim) if isinstance(v_nondim, UnitAwareArray) else v_nondim + + # Upstream position + end_pt_coords = coords - v_at_mid_pts * dt_for_calc + end_pt_coords = self.mesh.return_coords_to_bounds(end_pt_coords) + + # Sample psi_star[i] at upstream position + expr_to_evaluate = self.psi_star[i].sym + if hasattr(expr_to_evaluate, 'shape') and expr_to_evaluate.shape == (1, 1): + expr_to_evaluate = expr_to_evaluate[0, 0] + + value_at_end_points = uw.function.global_evaluate( + expr_to_evaluate, end_pt_coords, + ) + + psi_star_units = self.psi_star[i].units + if psi_star_units is not None and not isinstance(value_at_end_points, UnitAwareArray): + value_at_end_points = UnitAwareArray(value_at_end_points, units=psi_star_units) + + self.psi_star[i].array[...] = value_at_end_points + + return + + @property + def bdf_coefficients(self): + """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" + return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + def bdf(self, order=None): r"""Backwards differentiation form for calculating DuDt. @@ -1748,7 +1903,9 @@ def _object_viewer(self): @property def effective_order(self): """Current effective BDF order, accounting for history startup.""" - return min(self.order, max(1, self._n_solves_completed + 1)) + # BDF-k requires k completed solves to have k distinct history values. + # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves. + return min(self.order, max(1, self._n_solves_completed)) def initialise_history(self): r"""Initialize all history slots to the current value of :math:`\psi`. @@ -1854,6 +2011,11 @@ def update_post_solve( if self._n_solves_completed < self.order: self._n_solves_completed += 1 + @property + def bdf_coefficients(self): + """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" + return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + def bdf(self, order=None): r"""Backwards differentiation form for calculating DuDt. @@ -2067,7 +2229,9 @@ def _object_viewer(self): @property def effective_order(self): """Current effective BDF order, accounting for history startup.""" - return min(self.order, max(1, self._n_solves_completed + 1)) + # BDF-k requires k completed solves to have k distinct history values. + # With 0 or 1 completed solves → order 1. Order 2 needs ≥2 solves. + return min(self.order, max(1, self._n_solves_completed)) def initialise_history(self): r"""Initialize all history slots to the current value of :math:`\psi`. @@ -2174,6 +2338,11 @@ def update_post_solve( return + @property + def bdf_coefficients(self): + """Current BDF coefficients [c0, c1, ...] accounting for variable timesteps.""" + return _bdf_coefficients(self.effective_order, self._dt, self._dt_history) + def bdf(self, order=None): r"""Backwards differentiation form for calculating DuDt. diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 25c85b16..0bf156ed 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1246,6 +1246,22 @@ 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 @@ -1274,12 +1290,24 @@ def solve( if timestep is None: timestep = self.delta_t.sym - if timestep != self.delta_t: - self._constitutive_model.Parameters.elastic_dt = timestep # this will force an initialisation because the functions need to be updated + # 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 @@ -1289,11 +1317,29 @@ def solve( 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 - pre-solve DFDt update", flush=True) + print(f"VE Stokes solver - advect stress history", flush=True) - # Update SemiLagrange Flux terms - self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf) + 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) @@ -1305,8 +1351,47 @@ def solve( 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"VEP Stokes solver - post-solve DFDt update", flush=True) + 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) diff --git a/src/underworld3/utilities/unit_aware_coordinates.py b/src/underworld3/utilities/unit_aware_coordinates.py index 304aadeb..0073f283 100644 --- a/src/underworld3/utilities/unit_aware_coordinates.py +++ b/src/underworld3/utilities/unit_aware_coordinates.py @@ -83,6 +83,37 @@ def get_units(self): return self._units +class TimeSymbol(sympy.Symbol): + """A sympy Symbol subclass that carries _ccodestr and _units for JIT. + + Standard sympy Symbols are immutable and don't allow setting arbitrary + attributes. This subclass permits _ccodestr (for C code generation) + and _units (for dimensional analysis), following the same pattern as + UnitAwareBaseScalar for spatial coordinates. + """ + + _ccodestr = "petsc_t" + _units = None + + def __new__(cls, name="t", **kwargs): + obj = super().__new__(cls, name, real=True, nonnegative=True, **kwargs) + return obj + + @property + def units(self): + return self._units + + @units.setter + def units(self, value): + self._units = value + + def get_units(self): + return self._units + + def _ccode(self, printer): + return self._ccodestr + + def create_unit_aware_coordinate_system(name, units=None): """ Create a coordinate system with unit-aware coordinates. @@ -207,6 +238,40 @@ def _patch_coordinate(coord, units): # If _Gamma needs any units in the future, they should be explicitly # dimensionless (None), not inherited from mesh coordinates. + # Patch time coordinate with time units from the model + _patch_time_units(mesh) + + +def _patch_time_units(mesh): + """ + Patch mesh.t with time units from the model's scaling system. + + If the model has a time scale defined, mesh.t gets those units. + Otherwise mesh.t remains dimensionless (._units = None). + """ + if not hasattr(mesh, "_t"): + return + + time_units = None + try: + import underworld3 as uw + + model = uw.get_default_model() + if hasattr(model, "_fundamental_scales") and model._fundamental_scales: + scales = model._fundamental_scales + if "time" in scales: + time_scale = scales["time"] + if hasattr(time_scale, "units"): + time_units = time_scale.units + elif hasattr(time_scale, "_pint_qty"): + time_units = time_scale._pint_qty.units + except Exception: + pass + + mesh._t._units = time_units + if not hasattr(mesh._t, "get_units"): + mesh._t.get_units = lambda: mesh._t._units + def get_coordinate_units(coord): """ diff --git a/tests/plot_ve_oscillatory_validation.py b/tests/plot_ve_oscillatory_validation.py new file mode 100644 index 00000000..f3ac103f --- /dev/null +++ b/tests/plot_ve_oscillatory_validation.py @@ -0,0 +1,239 @@ +"""Viscoelastic oscillatory shear — validation plot for PR. + +Produces a figure showing: + - Top panel: applied shear rate γ̇(t) = γ̇₀ sin(ωt) + - Bottom panel: stress σ_xy(t) for order 1 and order 2 vs analytical + +Maxwell analytical solution including startup transient: + σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)] + +Results are saved as .npz checkpoint files for re-analysis. + +Usage: + python tests/plot_ve_oscillatory_validation.py + python tests/plot_ve_oscillatory_validation.py --replot # replot from saved data +""" + +import numpy as np +import sympy +import os +import sys +import underworld3 as uw +from underworld3.function import expression +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient).""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def run_oscillatory(order, n_steps, dt, V0, H, ETA, MU, omega, save_prefix=None): + """Run oscillatory VE shear box, return time series and save checkpoints.""" + W = 2.0 + gamma_dot_0 = 2.0 * V0 / H + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W / 2, -H / 2), maxCoords=(W / 2, H / 2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity") + + stokes.add_dirichlet_bc((V_bc, 0.0), "Top") + stokes.add_dirichlet_bc((-V_bc, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + + centre = np.array([[0.0, 0.0]]) + times, stress_num = [], [] + time_phys = 0.0 + + checkpoint_interval = max(1, n_steps // 20) + + for step in range(n_steps): + time_phys += dt + + V_t = V0 * np.sin(omega * time_phys) + V_bc.sym = V_t + + stokes.solve(zero_init_guess=False, timestep=dt, evalf=False) + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + times.append(time_phys) + stress_num.append(sigma_xy) + + if save_prefix and (step + 1) % checkpoint_interval == 0: + np.savez( + f"{save_prefix}_order{order}_checkpoint.npz", + times=np.array(times), + stress=np.array(stress_num), + step=step + 1, + dt=dt, + omega=omega, + ETA=ETA, + MU=MU, + V0=V0, + H=H, + order=order, + ) + + del stokes, mesh + + t_arr = np.array(times) + s_arr = np.array(stress_num) + + if save_prefix: + np.savez( + f"{save_prefix}_order{order}_final.npz", + times=t_arr, + stress=s_arr, + dt=dt, + omega=omega, + ETA=ETA, + MU=MU, + V0=V0, + H=H, + order=order, + gamma_dot_0=gamma_dot_0, + ) + + return t_arr, s_arr + + +def make_plot(results, params, output_path): + """Generate the validation figure from results dict.""" + + ETA = params["ETA"] + MU = params["MU"] + V0 = params["V0"] + H = params["H"] + omega = params["omega"] + dt = params["dt"] + t_r = ETA / MU + De = omega * t_r + gamma_dot_0 = 2.0 * V0 / H + + # Determine time range from results + t_max = max(r[0][-1] for r in results.values()) + + # Analytical solution (fine sampling) + t_fine = np.linspace(0, t_max, 1000) + sigma_analytical = maxwell_oscillatory(t_fine, ETA, MU, gamma_dot_0, omega) + shear_rate = gamma_dot_0 * np.sin(omega * t_fine) + + phase_lag_deg = np.degrees(np.arctan(De)) + steady_amp = ETA * gamma_dot_0 / np.sqrt(1 + De**2) + period = 2.0 * np.pi / omega + + fig, (ax1, ax2) = plt.subplots( + 2, 1, figsize=(10, 6), sharex=True, + gridspec_kw={"height_ratios": [1, 2.5]}, + ) + + # Top panel: applied shear rate + ax1.plot(t_fine / t_r, shear_rate, "k-", linewidth=1.2, + label=r"$\dot{\gamma}(t) = \dot{\gamma}_0 \sin(\omega t)$") + ax1.set_ylabel(r"Shear rate $\dot{\gamma}$") + ax1.axhline(0, color="grey", linewidth=0.5) + ax1.legend(loc="upper right", fontsize=9) + ax1.set_xlim(0, t_max / t_r) + + # Bottom panel: stress response + ax2.plot(t_fine / t_r, sigma_analytical, "k-", linewidth=1.5, + label="Analytical (Maxwell)", zorder=3) + + colors = {1: "#2196F3", 2: "#E91E63", 3: "#4CAF50"} + markers = {1: "o", 2: "s", 3: "^"} + + for order in sorted(results.keys()): + t, stress = results[order] + ax2.plot(t / t_r, stress, markers.get(order, "o"), + color=colors.get(order, "grey"), + markersize=2.5, alpha=0.7, + label=f"Order {order} (BDF-{order})", zorder=2) + + ax2.set_xlabel(r"Time $t / t_r$") + ax2.set_ylabel(r"Stress $\sigma_{xy}$") + ax2.axhline(0, color="grey", linewidth=0.5) + ax2.legend(loc="lower right", fontsize=9) + + ax2.set_title( + f"Maxwell viscoelastic shear: De = {De:.1f}, " + r"$\delta t / t_r$" + f" = {dt / t_r:.2f}, " + f"phase lag = {phase_lag_deg:.0f}\u00b0, " + f"steady amplitude = {steady_amp:.3f}", + fontsize=10, + ) + + plt.tight_layout() + plt.savefig(output_path, dpi=150, bbox_inches="tight") + print(f"Saved: {output_path}") + + +if __name__ == "__main__": + + output_dir = "tests" + save_prefix = os.path.join(output_dir, "ve_oscillatory") + + replot = "--replot" in sys.argv + + # Physical parameters + ETA = 1.0 + MU = 1.0 + V0 = 0.5 + H = 1.0 + t_r = ETA / MU + + De = 5.0 # High Deborah number — strong elastic effects, 79° phase lag + omega = De / t_r + period = 2.0 * np.pi / omega + n_periods = 4 + dt = 0.02 * t_r # Fine enough to resolve the fast oscillations + n_steps = int(n_periods * period / dt) + + params = dict(ETA=ETA, MU=MU, V0=V0, H=H, omega=omega, dt=dt) + + print(f"Maxwell oscillatory shear: De={De}, omega={omega:.3f}, period={period:.3f}") + print(f"dt={dt}, dt/t_r={dt/t_r}, n_steps={n_steps}, t_end={n_steps * dt:.2f}") + print(f"Phase lag = {np.degrees(np.arctan(De)):.1f} deg") + print() + + if replot: + # Load saved data + results = {} + for order in [1, 2]: + fpath = f"{save_prefix}_order{order}_final.npz" + if os.path.exists(fpath): + data = np.load(fpath) + results[order] = (data["times"], data["stress"]) + print(f"Loaded order {order} from {fpath}") + else: + print(f"Missing {fpath} — run without --replot first") + sys.exit(1) + else: + results = {} + for order in [1, 2]: + print(f"Running order {order} ({n_steps} steps)...", end=" ", flush=True) + t, stress = run_oscillatory( + order, n_steps, dt, V0, H, ETA, MU, omega, + save_prefix=save_prefix, + ) + results[order] = (t, stress) + print("done") + + make_plot(results, params, os.path.join(output_dir, "ve_oscillatory_validation.png")) diff --git a/tests/plot_ve_vep_oscillatory.py b/tests/plot_ve_vep_oscillatory.py new file mode 100644 index 00000000..abc234ca --- /dev/null +++ b/tests/plot_ve_vep_oscillatory.py @@ -0,0 +1,41 @@ +"""Plot VE vs VEP oscillatory shear from saved checkpoint.""" + +import numpy as np +import matplotlib.pyplot as plt + +data = np.load("output/ve_vep_oscillatory.npz") +t_ve, s_ve = data["t_ve"], data["s_ve"] +t_vep, s_vep = data["t_vep"], data["s_vep"] +t_ana, s_ana = data["t_ana"], data["s_ana"] +De = float(data["De"]) +tau_y = float(data["tau_y"]) +ve_amp = float(data["ve_amp"]) +omega = float(data["omega"]) + +fig, ax = plt.subplots(1, 1, figsize=(10, 5)) + +# Analytical VE (smooth curve) +ax.plot(t_ana, s_ana, "k-", linewidth=0.8, alpha=0.5, label="VE analytical") + +# Numerical VE +ax.plot(t_ve, s_ve, "b-o", markersize=3, linewidth=1.2, label="VE numerical (order 1)") + +# Numerical VEP +ax.plot(t_vep, s_vep, "r-s", markersize=3, linewidth=1.2, label=f"VEP numerical ($\\tau_y$={tau_y})") + +# Yield stress lines +ax.axhline(tau_y, color="r", linestyle="--", alpha=0.4, linewidth=0.8) +ax.axhline(-tau_y, color="r", linestyle="--", alpha=0.4, linewidth=0.8) +ax.text(0.3, tau_y + 0.02, f"$\\tau_y$ = {tau_y}", color="r", fontsize=9, alpha=0.6) + +ax.set_xlabel("Time ($t / t_r$)") +ax.set_ylabel("$\\sigma_{xy}$") +ax.set_title(f"Oscillatory Maxwell VE vs VEP shear (De = {De})") +ax.legend(loc="upper left", fontsize=9) +ax.set_xlim(0, t_ve[-1]) +ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.savefig("output/ve_vep_oscillatory.png", dpi=150) +plt.savefig("output/ve_vep_oscillatory.pdf") +print("Saved output/ve_vep_oscillatory.png and .pdf") diff --git a/tests/run_ve_order2_debug.py b/tests/run_ve_order2_debug.py new file mode 100644 index 00000000..286b7359 --- /dev/null +++ b/tests/run_ve_order2_debug.py @@ -0,0 +1,81 @@ +"""Debug order-2 VE: trace psi_star values at each step.""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + +ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 +dt = 0.1 +gamma_dot = 2.0 * V0 / H + +mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), +) +v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + +stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=2) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA +stokes.constitutive_model.Parameters.shear_modulus = MU +stokes.constitutive_model.Parameters.dt_elastic = dt + +stokes.add_dirichlet_bc((V0, 0.0), "Top") +stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") +stokes.tolerance = 1.0e-6 + +Stress = uw.discretisation.MeshVariable( + "Stress", mesh, (2, 2), vtype=uw.VarType.SYM_TENSOR, degree=2, continuous=True, +) +work = uw.discretisation.MeshVariable("W", mesh, 1, degree=2) +sigma_proj = uw.systems.Tensor_Projection(mesh, tensor_Field=Stress, scalar_Field=work) + +ddt = stokes.DFDt +centre = np.array([[0.0, 0.0]]) + +time_phys = 0.0 +for step in range(5): + eff_order = stokes.constitutive_model.effective_order + + # Read psi_star values at centre BEFORE solve + psi0_pre = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre) + psi1_pre = uw.function.evaluate(ddt.psi_star[1].sym[0, 1], centre) + psi0_val = float(psi0_pre.flatten()[0]) + psi1_val = float(psi1_pre.flatten()[0]) + + # Check psi_fn formula + psi_fn_01 = str(stokes.DFDt.psi_fn[0, 1]) + has_2star = "**" in psi_fn_01 or "psi_star" in psi_fn_01 + + t0 = timer.time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_t = timer.time() - t0 + + psi_fn_01_post = str(stokes.DFDt.psi_fn[0, 1]) + print(f" psi_fn[0,1] BEFORE solve: {psi_fn_01[:100]}...") + print(f" psi_fn[0,1] AFTER solve: {psi_fn_01_post[:100]}...") + time_phys += dt + + # Read psi_star AFTER solve (after update_post_solve) + psi0_post = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre) + psi1_post = uw.function.evaluate(ddt.psi_star[1].sym[0, 1], centre) + psi0_post_val = float(psi0_post.flatten()[0]) + psi1_post_val = float(psi1_post.flatten()[0]) + + sigma_proj.uw_function = stokes.stress_deviator + sigma_proj.solve() + val = uw.function.evaluate(Stress.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA)) + + print(f"step {step} eff_order={eff_order} t={time_phys:.2f} solve={solve_t:.1f}s") + print(f" PRE: psi_star[0]_xy={psi0_val:.6f} psi_star[1]_xy={psi1_val:.6f}") + print(f" POST: psi_star[0]_xy={psi0_post_val:.6f} psi_star[1]_xy={psi1_post_val:.6f}") + print(f" sigma_xy={sigma_xy:.6f} analytical={ana:.6f} " + f"rel_err={abs(sigma_xy - ana) / ana:.3e}") + print() + +print("Done") diff --git a/tests/run_ve_oscillatory.py b/tests/run_ve_oscillatory.py new file mode 100644 index 00000000..2cf7e21d --- /dev/null +++ b/tests/run_ve_oscillatory.py @@ -0,0 +1,122 @@ +"""Oscillatory VE shear box — time-harmonic analytical validation. + +Maxwell material under sinusoidal boundary velocity V(t) = V0 sin(ωt). +Full solution including startup transient: + + σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)] + +where De = ω t_r is the Deborah number. + +At steady state (t >> t_r), the transient dies out and the stress +oscillates with amplitude η γ̇₀/√(1+De²) and phase lag arctan(De). +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw +from underworld3.function import expression + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient).""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def run_oscillatory(order, n_steps, dt_over_tr, De): + """Run oscillatory VE shear box.""" + + ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0 + t_r = ETA / MU + omega = De / t_r + V0 = 0.5 + gamma_dot_0 = 2.0 * V0 / H + dt = dt_over_tr * t_r + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + # Use a UWexpression for time-dependent BCs — avoids JIT recompilation + V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity") + + stokes.add_dirichlet_bc((V_bc, 0.0), "Top") + stokes.add_dirichlet_bc((-V_bc, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + + centre = np.array([[0.0, 0.0]]) + times, stress_num, stress_ana = [], [], [] + time_phys = 0.0 + + for step in range(n_steps): + time_phys += dt + + # Update boundary velocity for this timestep + V_t = V0 * np.sin(omega * time_phys) + V_bc.sym = V_t + + stokes.solve(zero_init_guess=False, evalf=False) + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + + ana = maxwell_oscillatory(time_phys, ETA, MU, gamma_dot_0, omega) + times.append(time_phys) + stress_num.append(sigma_xy) + stress_ana.append(ana) + + del stokes, mesh + return np.array(times), np.array(stress_num), np.array(stress_ana) + + +if __name__ == "__main__": + De = 1.0 # Deborah number (equal viscous and elastic timescales) + dt_ratio = 0.1 # coarser for speed; finer (0.05) for publication quality + n_periods = 3 + t_r = 1.0 + omega = De / t_r + period = 2 * np.pi / omega + n_steps = int(n_periods * period / (dt_ratio * t_r)) + + print(f"De = {De}, omega = {omega:.3f}, period = {period:.3f}") + print(f"dt/t_r = {dt_ratio}, n_steps = {n_steps}") + print() + + for order in [1, 2]: + t0 = timer.time() + t, num, ana = run_oscillatory(order, n_steps, dt_ratio, De) + wall = timer.time() - t0 + + # Error over the last period (steady state) + last_period = t > (n_periods - 1) * period + if np.any(last_period): + rms = np.sqrt(np.mean((num[last_period] - ana[last_period])**2)) + amp_ana = np.max(np.abs(ana[last_period])) + rel_rms = rms / amp_ana if amp_ana > 0 else float('nan') + else: + rel_rms = float('nan') + + print(f"Order {order}: wall={wall:.0f}s, " + f"last-period relative RMS = {rel_rms:.4e}") + + # Samples from the last period + if np.any(last_period): + idx = np.where(last_period)[0] + sample = idx[::max(1, len(idx)//10)] + for i in sample: + err = abs(num[i] - ana[i]) + print(f" t={t[i]:.3f} num={num[i]:+.6f} ana={ana[i]:+.6f} err={err:.4e}") + print() diff --git a/tests/run_ve_shear_order2_quick.py b/tests/run_ve_shear_order2_quick.py new file mode 100644 index 00000000..fb40d71f --- /dev/null +++ b/tests/run_ve_shear_order2_quick.py @@ -0,0 +1,48 @@ +"""Quick order-2 VE shear box test — read stress from psi_star[0].""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + +ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 +dt = 0.1 +gamma_dot = 2.0 * V0 / H + +mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), +) +v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + +stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=2) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA +stokes.constitutive_model.Parameters.shear_modulus = MU +stokes.constitutive_model.Parameters.dt_elastic = dt + +stokes.add_dirichlet_bc((V0, 0.0), "Top") +stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") +stokes.tolerance = 1.0e-6 + +centre = np.array([[0.0, 0.0]]) +ddt = stokes.DFDt + +time_phys = 0.0 +for step in range(20): + eff = stokes.constitutive_model.effective_order + t0 = timer.time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_t = timer.time() - t0 + time_phys += dt + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA)) + print(f"step {step:2d} eff={eff} t={time_phys:.2f} solve={solve_t:.1f}s " + f"sigma_xy={sigma_xy:.6f} analytical={ana:.6f} " + f"rel_err={abs(sigma_xy - ana) / ana:.3e}") + +print("Done") diff --git a/tests/run_ve_shear_quick.py b/tests/run_ve_shear_quick.py new file mode 100644 index 00000000..debced41 --- /dev/null +++ b/tests/run_ve_shear_quick.py @@ -0,0 +1,48 @@ +"""Quick order-1 VE shear box validation (10 steps).""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + +ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 +dt = 0.1 +gamma_dot = 2.0 * V0 / H + +mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), +) +v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + +stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=1) +stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA +stokes.constitutive_model.Parameters.shear_modulus = MU +stokes.constitutive_model.Parameters.dt_elastic = dt + +stokes.add_dirichlet_bc((V0, 0.0), "Top") +stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") +stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") +stokes.tolerance = 1.0e-6 + +centre = np.array([[0.0, 0.0]]) +ddt = stokes.DFDt + +time_phys = 0.0 +for step in range(10): + t0 = timer.time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_t = timer.time() - t0 + time_phys += dt + + # Read stress from solver.tau (the projected actual stress) + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + ana = ETA * gamma_dot * (1.0 - np.exp(-time_phys * MU / ETA)) + print(f"step {step} t={time_phys:.2f} solve={solve_t:.1f}s " + f"sigma_xy={sigma_xy:.6f} analytical={ana:.6f} " + f"rel_err={abs(sigma_xy - ana) / ana:.3e}") + +print("Done") diff --git a/tests/run_ve_shear_validation.py b/tests/run_ve_shear_validation.py new file mode 100644 index 00000000..759616ed --- /dev/null +++ b/tests/run_ve_shear_validation.py @@ -0,0 +1,112 @@ +"""Quick validation script for the VE shear box test. + +Run with: pixi run -e amr-dev python tests/run_ve_shear_validation.py +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + + +def maxwell_xy(t, eta, mu, gamma_dot): + """Analytical σ_xy for Maxwell material under constant shear rate.""" + return eta * gamma_dot * (1.0 - np.exp(-t * mu / eta)) + + +def run(order, n_steps, dt_ratio): + ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 + t_relax = ETA / MU + dt = dt_ratio * t_relax + gamma_dot = 2.0 * V0 / H + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), + minCoords=(-W / 2, -H / 2), + maxCoords=(W / 2, H / 2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + stokes.add_dirichlet_bc((V0, 0.0), "Top") + stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + + Stress = uw.discretisation.MeshVariable( + "Stress", mesh, (2, 2), + vtype=uw.VarType.SYM_TENSOR, degree=2, continuous=True, + ) + work = uw.discretisation.MeshVariable("W", mesh, 1, degree=2) + sigma_proj = uw.systems.Tensor_Projection(mesh, tensor_Field=Stress, scalar_Field=work) + + times, num, ana = [], [], [] + time_phys = 0.0 + + for step in range(n_steps): + t0 = timer.time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_time = timer.time() - t0 + time_phys += dt + + sigma_proj.uw_function = stokes.stress_deviator + sigma_proj.solve() + + val = uw.function.evaluate(Stress.sym[0, 1], np.array([[0.0, 0.0]])) + sigma_xy = float(val.flatten()[0]) + ana_val = maxwell_xy(time_phys, ETA, MU, gamma_dot) + + times.append(time_phys) + num.append(sigma_xy) + ana.append(ana_val) + + rel_err = abs(sigma_xy - ana_val) / max(ana_val, 1e-10) + print(f" step {step:3d} t={time_phys:.2f} solve={solve_time:.1f}s " + f"sigma_xy={sigma_xy:.6f} analytical={ana_val:.6f} rel_err={rel_err:.3e}") + + del stokes, mesh + return np.array(times), np.array(num), np.array(ana) + + +if __name__ == "__main__": + print("=== Order 1, dt/tr=0.1, 20 steps ===") + t0 = timer.time() + t, n, a = run(1, 20, 0.1) + wall = timer.time() - t0 + mask = a > 0.01 * a[-1] + rms = np.sqrt(np.mean(((n[mask] - a[mask]) / a[mask]) ** 2)) + final_err = abs(n[-1] - a[-1]) / a[-1] + print(f" Wall time: {wall:.0f}s") + print(f" Final rel error: {final_err:.4e}") + print(f" RMS rel error: {rms:.4e}") + print() + + print("=== Order 2, dt/tr=0.1, 20 steps ===") + t0 = timer.time() + t2, n2, a2 = run(2, 20, 0.1) + wall2 = timer.time() - t0 + mask2 = a2 > 0.01 * a2[-1] + rms2 = np.sqrt(np.mean(((n2[mask2] - a2[mask2]) / a2[mask2]) ** 2)) + final_err2 = abs(n2[-1] - a2[-1]) / a2[-1] + print(f" Wall time: {wall2:.0f}s") + print(f" Final rel error: {final_err2:.4e}") + print(f" RMS rel error: {rms2:.4e}") + print() + + viscous_limit = 1.0 + print(f"Steady-state check (order 1, t=2.0):") + print(f" sigma_xy = {n[-1]:.6f}, viscous limit = {viscous_limit:.6f}, " + f"rel error = {abs(n[-1] - viscous_limit) / viscous_limit:.4e}") + print() + + if rms2 < rms: + print("PASS: Order 2 has smaller RMS error than order 1") + else: + print("FAIL: Order 2 should have smaller error than order 1") diff --git a/tests/run_ve_square_wave.py b/tests/run_ve_square_wave.py new file mode 100644 index 00000000..d56ffc73 --- /dev/null +++ b/tests/run_ve_square_wave.py @@ -0,0 +1,244 @@ +"""Variable-timestep VE benchmark: square-wave forcing. + +Maxwell material under square-wave shear rate (Fourier series, N harmonics): + + γ̇(t) = (4γ̇₀/π) Σ_{k=1..N} sin((2k-1)ωt) / (2k-1) + +Since the Maxwell equation is linear, the analytical stress is the +superposition of single-frequency solutions: + + σ(t) = Σ_{k=1..N} maxwell_oscillatory(t, η, μ, aₖ, ωₖ) + +where aₖ = 4γ̇₀/(π(2k-1)) and ωₖ = (2k-1)ω. + +The sharp transitions demand small dt at the edges, large dt in the +flat regions — testing variable-dt BDF correctness. + +Usage: + python tests/run_ve_square_wave.py +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw +from underworld3.function import expression + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for single-frequency oscillatory Maxwell shear.""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def square_wave_analytical(t, eta, mu, gamma_dot_0, omega, n_harmonics=20): + """Analytical stress for square-wave forcing via Fourier superposition.""" + sigma = np.zeros_like(t) + for k in range(1, n_harmonics + 1): + n = 2 * k - 1 # odd harmonics: 1, 3, 5, ... + a_k = 4.0 * gamma_dot_0 / (np.pi * n) + omega_k = n * omega + sigma += maxwell_oscillatory(t, eta, mu, a_k, omega_k) + return sigma + + +def square_wave_shear_rate(t, gamma_dot_0, omega, n_harmonics=20): + """Square-wave shear rate via truncated Fourier series.""" + rate = np.zeros_like(t) + for k in range(1, n_harmonics + 1): + n = 2 * k - 1 + rate += 4.0 * gamma_dot_0 / (np.pi * n) * np.sin(n * omega * t) + return rate + + +def adaptive_dt(t_current, omega, dt_min, dt_max): + """Adaptive timestep: small near square-wave transitions, large on plateaux. + + Transitions occur at t = (2m+1)·π/(2ω) for integer m, i.e. at odd + multiples of quarter-period. We use distance to nearest transition + to interpolate between dt_min and dt_max. + """ + half_period = np.pi / omega + # Phase within half-period [0, half_period) + phase = t_current % half_period + # Distance to nearest transition (0 or half_period boundary) + dist = min(phase, half_period - phase) + # Normalise to [0, 1] where 0 = at transition, 1 = mid-plateau + frac = dist / (half_period / 2.0) + # Smooth interpolation + return dt_min + (dt_max - dt_min) * frac**2 + + +def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, + n_harmonics=20, uniform=False): + """Run VE square-wave shear box with adaptive or uniform timestep.""" + + ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0 + t_r = ETA / MU + omega = De / t_r + V0 = 0.5 + gamma_dot_0 = 2.0 * V0 / H + dt_min = dt_min_over_tr * t_r + dt_max = dt_max_over_tr * t_r + T = 2.0 * np.pi / omega + t_end = n_periods * T + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W / 2, -H / 2), maxCoords=(W / 2, H / 2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + + # Boundary conditions: simple shear driven by top/bottom velocity + # V_top updated numerically each timestep to produce square-wave γ̇ + V_top = expression(R"V_{\mathrm{top}}", sympy.Float(0.0), "Top boundary velocity") + + stokes.add_dirichlet_bc((V_top, 0.0), "Top") + stokes.add_dirichlet_bc((-V_top, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + + # Time loop + times = [] + numerical_stress = [] + timesteps_used = [] + t_current = 0.0 + step = 0 + + while t_current < t_end: + # Determine timestep + if uniform: + dt = dt_min + else: + dt = adaptive_dt(t_current, omega, dt_min, dt_max) + + t_next = t_current + dt + + # Update boundary velocity: V_top(t) such that γ̇ = 2·V_top/H = square wave + # square_wave_shear_rate returns the actual γ̇, so V_top = γ̇ · H/2 + gamma_dot_t = square_wave_shear_rate( + np.array([t_next]), gamma_dot_0, omega, n_harmonics + )[0] + V_top.sym = sympy.Float(gamma_dot_t * H / 2.0) + + stokes.constitutive_model.Parameters.dt_elastic = dt + + stokes.solve(zero_init_guess=False, timestep=dt) + + t_current += dt + step += 1 + + # Extract stress at mesh centre + centre = np.array([[0.0, 0.0]]) + tau_xy = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(tau_xy.flatten()[0]) + + times.append(t_current) + numerical_stress.append(sigma_xy) + timesteps_used.append(dt) + + if step % 50 == 0: + ana = square_wave_analytical( + np.array([t_current]), ETA, MU, gamma_dot_0, omega, n_harmonics + )[0] + print(f" Step {step:4d}: t/t_r = {t_current / t_r:.3f}, " + f"dt/t_r = {dt / t_r:.4f}, " + f"σ_xy = {sigma_xy:.6f}, ana = {ana:.6f}") + + times = np.array(times) + numerical_stress = np.array(numerical_stress) + timesteps_used = np.array(timesteps_used) + + # Analytical solution + analytical_stress = square_wave_analytical(times, ETA, MU, gamma_dot_0, omega, n_harmonics) + + # Error metrics (skip first period for startup transient) + mask = times > T + if mask.sum() > 0: + l2_err = np.sqrt(np.mean((numerical_stress[mask] - analytical_stress[mask]) ** 2)) + linf_err = np.max(np.abs(numerical_stress[mask] - analytical_stress[mask])) + else: + l2_err = linf_err = np.nan + + return { + "times": times, + "numerical": numerical_stress, + "analytical": analytical_stress, + "timesteps": timesteps_used, + "l2_error": l2_err, + "linf_error": linf_err, + "n_steps": step, + } + + +if __name__ == "__main__": + De = 1.5 + order = 2 + n_periods = 3 + n_harmonics = 3 + + print("=" * 60) + print(f"Square-wave VE benchmark: De={De}, order={order}") + print(f" {n_harmonics} Fourier harmonics, {n_periods} periods") + print("=" * 60) + + # Run with adaptive dt + print("\n--- Adaptive timestep ---") + t0 = timer.time() + result_adaptive = run_square_wave( + order=order, De=De, n_periods=n_periods, + dt_min_over_tr=0.02, dt_max_over_tr=0.15, + n_harmonics=n_harmonics, uniform=False, + ) + t_adaptive = timer.time() - t0 + print(f" {result_adaptive['n_steps']} steps in {t_adaptive:.1f}s") + print(f" L2 error: {result_adaptive['l2_error']:.6e}") + print(f" Linf error: {result_adaptive['linf_error']:.6e}") + print(f" dt range: [{result_adaptive['timesteps'].min():.4f}, " + f"{result_adaptive['timesteps'].max():.4f}]") + + # Run with uniform dt (reference) + print("\n--- Uniform timestep (dt_min) ---") + t0 = timer.time() + result_uniform = run_square_wave( + order=order, De=De, n_periods=n_periods, + dt_min_over_tr=0.02, dt_max_over_tr=0.02, + n_harmonics=n_harmonics, uniform=True, + ) + t_uniform = timer.time() - t0 + print(f" {result_uniform['n_steps']} steps in {t_uniform:.1f}s") + print(f" L2 error: {result_uniform['l2_error']:.6e}") + print(f" Linf error: {result_uniform['linf_error']:.6e}") + + # Summary + print("\n" + "=" * 60) + print("Summary:") + print(f" Adaptive: {result_adaptive['n_steps']} steps, " + f"L2 = {result_adaptive['l2_error']:.2e}") + print(f" Uniform: {result_uniform['n_steps']} steps, " + f"L2 = {result_uniform['l2_error']:.2e}") + ratio = result_adaptive['l2_error'] / result_uniform['l2_error'] + print(f" Error ratio (adaptive/uniform): {ratio:.2f}") + print(f" Step savings: {1 - result_adaptive['n_steps'] / result_uniform['n_steps']:.0%}") + + # Save results + np.savez( + "tests/ve_square_wave_benchmark.npz", + adaptive_times=result_adaptive["times"], + adaptive_numerical=result_adaptive["numerical"], + adaptive_analytical=result_adaptive["analytical"], + adaptive_timesteps=result_adaptive["timesteps"], + uniform_times=result_uniform["times"], + uniform_numerical=result_uniform["numerical"], + uniform_analytical=result_uniform["analytical"], + uniform_timesteps=result_uniform["timesteps"], + ) + print("\nResults saved to tests/ve_square_wave_benchmark.npz") diff --git a/tests/run_ve_vep_oscillatory_checkpoint.py b/tests/run_ve_vep_oscillatory_checkpoint.py new file mode 100644 index 00000000..8ae07f01 --- /dev/null +++ b/tests/run_ve_vep_oscillatory_checkpoint.py @@ -0,0 +1,124 @@ +"""Run VE and VEP oscillatory shear, save time series to .npz for plotting.""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw +from underworld3.function import expression + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for oscillatory Maxwell shear.""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def run_oscillatory(order, n_steps, dt, omega, V0, tau_y=None): + """Run oscillatory shear (VE or VEP depending on tau_y).""" + + ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0 + gamma_dot_0 = 2.0 * V0 / H + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + if tau_y is not None: + stokes.constitutive_model.Parameters.yield_stress = tau_y + stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6 + + v.data[:, 0] = v.coords[:, 1] * gamma_dot_0 + v.data[:, 1] = 0.0 + + V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity") + stokes.add_dirichlet_bc((V_bc, 0.0), "Top") + stokes.add_dirichlet_bc((-V_bc, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-5 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["snes_max_it"] = 50 + + centre = np.array([[0.0, 0.0]]) + times, stress_num = [], [] + time_phys = 0.0 + + for step in range(n_steps): + time_phys += dt + V_bc.sym = V0 * np.sin(omega * time_phys) + + stokes.solve(zero_init_guess=False, evalf=False) + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + times.append(time_phys) + stress_num.append(sigma_xy) + + label = "VEP" if tau_y is not None else "VE" + print(f" {label} step {step:3d} t={time_phys:.2f} sigma={sigma_xy:+.6f}", flush=True) + + del stokes, mesh + return np.array(times), np.array(stress_num) + + +if __name__ == "__main__": + De = 1.0 + ETA, MU = 1.0, 1.0 + t_r = ETA / MU + omega = De / t_r + V0 = 0.5 + gamma_dot_0 = 2.0 * V0 / 1.0 + period = 2.0 * np.pi / omega + TAU_Y = 0.4 + + dt_ratio = 0.1 + dt = dt_ratio * t_r + n_periods = 2 + n_steps = int(n_periods * period / dt) + + ve_amp = ETA * gamma_dot_0 / np.sqrt(1.0 + De**2) + + print(f"De={De}, omega={omega:.3f}, period={period:.3f}") + print(f"VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}") + print(f"dt={dt}, n_steps={n_steps}") + print() + + # VE (no yield) + print("=== VE (order 1) ===") + t0 = timer.time() + t_ve, s_ve = run_oscillatory(1, n_steps, dt, omega, V0, tau_y=None) + print(f" Wall: {timer.time()-t0:.0f}s\n") + + # VEP (with yield) + print("=== VEP (order 1) ===") + t0 = timer.time() + t_vep, s_vep = run_oscillatory(1, n_steps, dt, omega, V0, tau_y=TAU_Y) + print(f" Wall: {timer.time()-t0:.0f}s\n") + + # Analytical + t_ana = np.linspace(dt, n_periods * period, 500) + s_ana = maxwell_oscillatory(t_ana, ETA, MU, gamma_dot_0, omega) + + # Save + outfile = "output/ve_vep_oscillatory.npz" + import os + os.makedirs("output", exist_ok=True) + np.savez( + outfile, + t_ve=t_ve, s_ve=s_ve, + t_vep=t_vep, s_vep=s_vep, + t_ana=t_ana, s_ana=s_ana, + De=De, tau_y=TAU_Y, ve_amp=ve_amp, omega=omega, + ) + print(f"Saved to {outfile}") diff --git a/tests/run_ve_vep_oscillatory_plot.py b/tests/run_ve_vep_oscillatory_plot.py new file mode 100644 index 00000000..b5ee4907 --- /dev/null +++ b/tests/run_ve_vep_oscillatory_plot.py @@ -0,0 +1,161 @@ +"""Run VE and VEP oscillatory shear from consistent initial conditions. + +Both cases start from the analytical Maxwell solution at t0=0.01 to avoid +zero-velocity singularity and ensure identical starting stress. + +Saves time series to stdout and generates plot directly. +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw +from underworld3.function import expression +import os + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for oscillatory Maxwell shear.""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def run_oscillatory(order, n_steps, dt, omega, V0, t0, tau_y=None): + """Run oscillatory shear from analytical initial condition at t0.""" + + ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0 + gamma_dot_0 = 2.0 * V0 / H + t_r = ETA / MU + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + if tau_y is not None: + stokes.constitutive_model.Parameters.yield_stress = tau_y + stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6 + + # Initialise velocity to match analytical solution at t0 + V_t0 = V0 * np.sin(omega * t0) + gamma_dot_t0 = 2.0 * V_t0 / H + v.data[:, 0] = v.coords[:, 1] * gamma_dot_t0 + v.data[:, 1] = 0.0 + + # Initialise stress history to analytical at t0 + sigma_t0 = maxwell_oscillatory(t0, ETA, MU, gamma_dot_0, omega) + # psi_star[0] is (N, dim, dim) sym tensor — set the xy component + stokes.DFDt.psi_star[0].array[:, 0, 1] = sigma_t0 + stokes.DFDt.psi_star[0].array[:, 1, 0] = sigma_t0 + stokes.DFDt._history_initialised = True + + V_bc = expression(r"{V_{bc}}", V_t0, "Time-dependent boundary velocity") + stokes.add_dirichlet_bc((V_bc, 0.0), "Top") + stokes.add_dirichlet_bc((-V_bc, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-5 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["snes_max_it"] = 50 + + centre = np.array([[0.0, 0.0]]) + times, stress_num = [], [] + time_phys = t0 + + label = "VEP" if tau_y is not None else "VE" + + for step in range(n_steps): + time_phys += dt + V_bc.sym = V0 * np.sin(omega * time_phys) + + stokes.solve(zero_init_guess=False, evalf=False) + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + times.append(time_phys) + stress_num.append(sigma_xy) + + print(f" {label} step {step:3d} t={time_phys:.3f} sigma={sigma_xy:+.6f}", flush=True) + + del stokes, mesh + return np.array(times), np.array(stress_num) + + +if __name__ == "__main__": + De = 1.0 + ETA, MU = 1.0, 1.0 + t_r = ETA / MU + omega = De / t_r + V0 = 0.5 + gamma_dot_0 = 2.0 * V0 / 1.0 + period = 2.0 * np.pi / omega + TAU_Y = 0.4 + + dt_ratio = 0.1 + dt = dt_ratio * t_r + t0 = 0.01 + n_periods = 2 + n_steps = int(n_periods * period / dt) + + ve_amp = ETA * gamma_dot_0 / np.sqrt(1.0 + De**2) + sigma_t0 = maxwell_oscillatory(t0, ETA, MU, gamma_dot_0, omega) + + print(f"De={De}, omega={omega:.3f}, period={period:.3f}") + print(f"VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}") + print(f"t0={t0}, sigma(t0)={sigma_t0:.6f}") + print(f"dt={dt}, n_steps={n_steps}") + print() + + # VE + print("=== VE (order 1) ===") + t0w = timer.time() + t_ve, s_ve = run_oscillatory(1, n_steps, dt, omega, V0, t0, tau_y=None) + print(f" Wall: {timer.time()-t0w:.0f}s\n") + + # VEP + print("=== VEP (order 1) ===") + t0w = timer.time() + t_vep, s_vep = run_oscillatory(1, n_steps, dt, omega, V0, t0, tau_y=TAU_Y) + print(f" Wall: {timer.time()-t0w:.0f}s\n") + + # Analytical + t_ana = np.linspace(t0, t_ve[-1], 500) + s_ana = maxwell_oscillatory(t_ana, ETA, MU, gamma_dot_0, omega) + + # Plot + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(1, 1, figsize=(10, 5)) + + ax.plot(t_ana, s_ana, "k-", linewidth=0.8, alpha=0.5, label="VE analytical") + ax.plot(t_ve, s_ve, "b-o", markersize=3, linewidth=1.2, label="VE numerical (order 1)") + ax.plot(t_vep, s_vep, "r-s", markersize=3, linewidth=1.2, label=f"VEP numerical ($\\tau_y$={TAU_Y})") + + ax.axhline(TAU_Y, color="r", linestyle="--", alpha=0.4, linewidth=0.8) + ax.axhline(-TAU_Y, color="r", linestyle="--", alpha=0.4, linewidth=0.8) + ax.text(0.3, TAU_Y + 0.02, f"$\\tau_y$ = {TAU_Y}", color="r", fontsize=9, alpha=0.6) + + ax.set_xlabel("Time ($t / t_r$)") + ax.set_ylabel("$\\sigma_{xy}$") + ax.set_title(f"Oscillatory Maxwell VE vs VEP shear (De = {De})") + ax.legend(loc="upper left", fontsize=9) + ax.set_xlim(0, t_ve[-1]) + ax.grid(True, alpha=0.3) + + plt.tight_layout() + os.makedirs("output", exist_ok=True) + plt.savefig("output/ve_vep_oscillatory.png", dpi=150) + plt.savefig("output/ve_vep_oscillatory.pdf") + print("Saved output/ve_vep_oscillatory.png and .pdf") diff --git a/tests/run_vep_oscillatory.py b/tests/run_vep_oscillatory.py new file mode 100644 index 00000000..e76e212e --- /dev/null +++ b/tests/run_vep_oscillatory.py @@ -0,0 +1,134 @@ +"""Oscillatory VEP shear box — does the yield stress cap the oscillation? + +Same setup as the VE oscillatory test but with a yield stress. +The VE amplitude at steady state is η γ̇₀/√(1+De²). +If τ_y is below this amplitude, the stress should be clipped. +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw +from underworld3.function import expression + + +def maxwell_oscillatory(t, eta, mu, gamma_dot_0, omega): + """Full analytical σ_xy for oscillatory Maxwell shear (incl. transient).""" + t_r = eta / mu + De = omega * t_r + prefactor = eta * gamma_dot_0 / (1.0 + De**2) + return prefactor * (np.sin(omega * t) - De * np.cos(omega * t) + De * np.exp(-t / t_r)) + + +def run_vep_oscillatory(order, n_steps, dt_over_tr, De, tau_y): + """Run oscillatory VEP shear box.""" + + ETA, MU, H, W = 1.0, 1.0, 1.0, 2.0 + t_r = ETA / MU + omega = De / t_r + V0 = 0.5 + gamma_dot_0 = 2.0 * V0 / H + dt = dt_over_tr * t_r + + # VE steady-state amplitude + ve_amplitude = ETA * gamma_dot_0 / np.sqrt(1 + De**2) + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.constitutive_model.Parameters.yield_stress = tau_y + stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6 + + # Seed with linear shear + v.data[:, 0] = v.coords[:, 1] * gamma_dot_0 + v.data[:, 1] = 0.0 + + V_bc = expression(r"{V_{bc}}", 0.0, "Time-dependent boundary velocity") + stokes.add_dirichlet_bc((V_bc, 0.0), "Top") + stokes.add_dirichlet_bc((-V_bc, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-5 + + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["snes_max_it"] = 50 + + centre = np.array([[0.0, 0.0]]) + times, stress_num, stress_ve = [], [], [] + time_phys = 0.0 + + for step in range(n_steps): + time_phys += dt + V_t = V0 * np.sin(omega * time_phys) + V_bc.sym = V_t + + t0s = __import__('time').time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_t = __import__('time').time() - t0s + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + ve_stress = maxwell_oscillatory(time_phys, ETA, MU, gamma_dot_0, omega) + + times.append(time_phys) + stress_num.append(sigma_xy) + stress_ve.append(ve_stress) + print(f" step {step:3d} t={time_phys:.2f} solve={solve_t:.1f}s " + f"vep={sigma_xy:+.6f} ve={ve_stress:+.6f}", flush=True) + + del stokes, mesh + return np.array(times), np.array(stress_num), np.array(stress_ve), ve_amplitude + + +if __name__ == "__main__": + De = 1.0 + dt_ratio = 0.1 + n_periods = 2 + t_r = 1.0 + omega = De / t_r + period = 2 * np.pi / omega + n_steps = int(n_periods * period / (dt_ratio * t_r)) + + ETA, gamma_dot_0 = 1.0, 1.0 + ve_amp = ETA * gamma_dot_0 / np.sqrt(1 + De**2) + + # Set τ_y below the VE amplitude so we see clipping + TAU_Y = 0.4 + + print(f"De={De}, VE amplitude={ve_amp:.4f}, tau_y={TAU_Y}") + print(f"Expect clipping at ±{TAU_Y}") + print(f"dt/t_r={dt_ratio}, n_steps={n_steps}") + print() + + for order in [1]: + t0 = timer.time() + t, num, ve, amp = run_vep_oscillatory(order, n_steps, dt_ratio, De, TAU_Y) + wall = timer.time() - t0 + + print(f"Order {order} (wall={wall:.0f}s):") + print(f" VE amplitude: {amp:.4f}") + print(f" Max |σ_xy|: {np.max(np.abs(num)):.6f}") + print(f" Min σ_xy: {np.min(num):.6f}") + print(f" Max σ_xy: {np.max(num):.6f}") + print(f" Clipped at τ_y={TAU_Y}? {np.max(np.abs(num)) <= TAU_Y + 0.01}") + print() + + # Print last period + last_period = t > (n_periods - 1) * period + if np.any(last_period): + idx = np.where(last_period)[0] + sample = idx[::max(1, len(idx)//12)] + print(f" Last period samples:") + for i in sample: + clipped = " CLIPPED" if abs(num[i]) > TAU_Y - 0.01 else "" + print(f" t={t[i]:.3f} vep={num[i]:+.6f} ve={ve[i]:+.6f}{clipped}") + print() diff --git a/tests/run_vep_shear_box.py b/tests/run_vep_shear_box.py new file mode 100644 index 00000000..7d4fabe3 --- /dev/null +++ b/tests/run_vep_shear_box.py @@ -0,0 +1,121 @@ +"""Viscoelastic-plastic shear box — stress buildup with yield cap. + +Maxwell VE material with Drucker-Prager yield stress under constant shear. + +Phase 1 (elastic buildup): σ_xy follows the Maxwell curve +Phase 2 (yielded): σ_xy = τ_y (plastic cap) + +Transition time: t_yield = -t_r ln(1 - τ_y/(η γ̇)) +(only exists if η γ̇ > τ_y) +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + + +def maxwell_stress_xy(t, eta, mu, gamma_dot): + """Uncapped VE stress (for comparison).""" + t_r = eta / mu + return eta * gamma_dot * (1.0 - np.exp(-t / t_r)) + + +def vep_stress_xy(t, eta, mu, gamma_dot, tau_y): + """VEP stress: Maxwell buildup capped at τ_y.""" + ve = maxwell_stress_xy(t, eta, mu, gamma_dot) + return np.minimum(ve, tau_y) + + +def run_vep_shear(order, n_steps, dt_over_tr, tau_y): + """Run VEP shear box, return (times, num_stress, ana_stress).""" + + ETA, MU, V0, H, W = 1.0, 1.0, 0.5, 1.0, 2.0 + t_r = ETA / MU + dt = dt_over_tr * t_r + gamma_dot = 2.0 * V0 / H # = 1.0 + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=order) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.constitutive_model.Parameters.yield_stress = tau_y + stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6 + + stokes.add_dirichlet_bc((V0, 0.0), "Top") + stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-5 + + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["snes_max_it"] = 50 + + # Seed with linear shear + v.data[:, 0] = v.coords[:, 1] * gamma_dot + v.data[:, 1] = 0.0 + + centre = np.array([[0.0, 0.0]]) + times, stress_num, stress_ana, stress_ve = [], [], [], [] + time_phys = 0.0 + + for step in range(n_steps): + t0 = timer.time() + stokes.solve(zero_init_guess=False, evalf=False) + solve_t = timer.time() - t0 + time_phys += dt + + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + + ana = vep_stress_xy(time_phys, ETA, MU, gamma_dot, tau_y) + ve_only = maxwell_stress_xy(time_phys, ETA, MU, gamma_dot) + + times.append(time_phys) + stress_num.append(sigma_xy) + stress_ana.append(ana) + stress_ve.append(ve_only) + + del stokes, mesh + return np.array(times), np.array(stress_num), np.array(stress_ana), np.array(stress_ve) + + +if __name__ == "__main__": + ETA, MU = 1.0, 1.0 + gamma_dot = 1.0 # 2*V0/H + t_r = ETA / MU + + # Choose τ_y below the viscous steady state (η·γ̇ = 1.0) + TAU_Y = 0.5 + t_yield = -t_r * np.log(1 - TAU_Y / (ETA * gamma_dot)) + + print(f"eta={ETA}, mu={MU}, gamma_dot={gamma_dot}") + print(f"tau_y={TAU_Y}, viscous_limit={ETA*gamma_dot}") + print(f"t_yield={t_yield:.4f} (stress reaches tau_y)") + print() + + for order in [1, 2]: + t0 = timer.time() + t, num, ana, ve = run_vep_shear(order, n_steps=30, dt_over_tr=0.1, tau_y=TAU_Y) + wall = timer.time() - t0 + + print(f"Order {order} (wall={wall:.0f}s):") + for i in range(len(t)): + marker = " *" if t[i] > t_yield - 0.05 and t[i] < t_yield + 0.15 else "" + print(f" t={t[i]:.2f} num={num[i]:.6f} ana={ana[i]:.6f} " + f"ve={ve[i]:.6f} err={abs(num[i]-ana[i]):.4e}{marker}") + + # Error in post-yield regime + yielded = t > t_yield + 0.2 + if np.any(yielded): + post_err = np.max(np.abs(num[yielded] - ana[yielded])) + print(f" Max post-yield error: {post_err:.4e}") + print() diff --git a/tests/run_vp_shear_box.py b/tests/run_vp_shear_box.py new file mode 100644 index 00000000..01a0c165 --- /dev/null +++ b/tests/run_vp_shear_box.py @@ -0,0 +1,90 @@ +"""Viscoplastic shear box — yield stress validation. + +For simple shear with constant strain rate: + - Below yield: σ_xy = η · γ̇ (viscous) + - At/above yield: σ_xy = τ_y (plastic cap) + +The transition occurs when η · γ̇ = τ_y, i.e. γ̇_crit = τ_y / η. + +We test with a range of driving velocities and check the stress +transition matches the analytical prediction. +""" + +import time as timer +import numpy as np +import sympy +import underworld3 as uw + + +def run_vp_shear(V0, eta, yield_stress): + """Single-step VP shear solve, return σ_xy at centre.""" + + H, W = 1.0, 2.0 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-W/2, -H/2), maxCoords=(W/2, H/2), + ) + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.constitutive_model.Parameters.yield_stress = yield_stress + stokes.constitutive_model.Parameters.strainrate_inv_II_min = 1.0e-6 + + stokes.add_dirichlet_bc((V0, 0.0), "Top") + stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-4 + + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["snes_max_it"] = 50 + + # Initialise with uniform shear to avoid singular Jacobian at zero strain rate + gamma_dot = 2.0 * V0 / H + v.data[:, 0] = v.coords[:, 1] * gamma_dot + v.data[:, 1] = 0.0 + + stokes.solve(zero_init_guess=False) + + # Read stress via tau + centre = np.array([[0.0, 0.0]]) + val = uw.function.evaluate(stokes.tau.sym[0, 1], centre) + sigma_xy = float(val.flatten()[0]) + + converged = stokes.snes.getConvergedReason() > 0 + del stokes, mesh + return sigma_xy, converged + + +if __name__ == "__main__": + ETA = 1.0 + TAU_Y = 0.5 + H = 1.0 + + # γ̇ = 2V0/H, σ_xy = η·γ̇ in viscous regime, σ_xy = τ_y when η·γ̇ ≥ τ_y + gamma_dot_crit = TAU_Y / ETA + print(f"eta={ETA}, tau_y={TAU_Y}, gamma_dot_crit={gamma_dot_crit}") + print() + + V0_values = np.array([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1.0, 2.0]) + + print(f"{'V0':>6} {'gamma_dot':>10} {'sigma_xy':>10} {'analytical':>10} {'rel_err':>10} {'status':>8}") + print("-" * 62) + + for V0 in V0_values: + gamma_dot = 2.0 * V0 / H + ana = min(ETA * gamma_dot, TAU_Y) + + t0 = timer.time() + sigma_xy, converged = run_vp_shear(V0, ETA, TAU_Y) + dt = timer.time() - t0 + + rel_err = abs(sigma_xy - ana) / max(ana, 1e-10) + status = "OK" if converged else "FAIL" + print(f"{V0:6.3f} {gamma_dot:10.4f} {sigma_xy:10.6f} {ana:10.6f} {rel_err:10.3e} {status:>8}") + + print() + print("Done") diff --git a/tests/test_1051_VE_shear_box.py b/tests/test_1051_VE_shear_box.py new file mode 100644 index 00000000..1608a729 --- /dev/null +++ b/tests/test_1051_VE_shear_box.py @@ -0,0 +1,181 @@ +""" +Viscoelastic shear box — analytical validation. + +A Maxwell viscoelastic material under uniform simple shear has the exact +solution for the shear stress component: + + σ_xy(t) = η · γ̇ · (1 - exp(-t / t_r)) + +where: + t_r = η/μ Maxwell relaxation time + γ̇ = dv_x/dy engineering shear rate (velocity gradient) + ε̇_xy = γ̇/2 tensor strain rate + +The steady-state stress is σ_xy = η·γ̇ = 2η·ε̇_xy, the purely viscous limit. + +Boundary conditions: + Top / Bottom: prescribed horizontal velocity ±V₀, zero vertical + Left / Right: free horizontal (outflow), zero vertical velocity + +This avoids periodic BCs (untested) and corner singularities from free-slip. +The result is a spatially uniform simple shear — effectively a 1D problem. + +The stress is read from psi_star[0] (the projected actual stress stored by +VE_Stokes after each solve), not recomputed from the constitutive formula. +""" + +import pytest +import numpy as np +import sympy +import underworld3 as uw + +pytestmark = [pytest.mark.level_3, pytest.mark.tier_b] + + +# --------------------------------------------------------------------------- +# Analytical solution +# --------------------------------------------------------------------------- + +def maxwell_stress_xy(t, eta, mu, gamma_dot): + """Exact σ_xy for constant-rate simple shear of a Maxwell material.""" + t_relax = eta / mu + return eta * gamma_dot * (1.0 - np.exp(-t / t_relax)) + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _run_ve_shear(order, n_steps, dt_over_tr): + """Run a VE shear-box, return (times, numerical_stress, analytical_stress).""" + + ETA = 1.0 + MU = 1.0 + V0 = 0.5 + H = 1.0 + W = 2.0 + + t_relax = ETA / MU + dt = dt_over_tr * t_relax + gamma_dot = 2.0 * V0 / H + + res = 8 + mesh = uw.meshing.StructuredQuadBox( + elementRes=(2 * res, res), + minCoords=(-W / 2.0, -H / 2.0), + maxCoords=(W / 2.0, H / 2.0), + ) + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + + stokes = uw.systems.VE_Stokes( + mesh, velocityField=v, pressureField=p, order=order, verbose=False, + ) + + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.constitutive_model.Parameters.shear_modulus = MU + stokes.constitutive_model.Parameters.dt_elastic = dt + + stokes.add_dirichlet_bc((V0, 0.0), "Top") + stokes.add_dirichlet_bc((-V0, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + + times = [] + stress_num = [] + stress_ana = [] + centre = np.array([[0.0, 0.0]]) + + time = 0.0 + for step in range(n_steps): + stokes.solve(zero_init_guess=False, evalf=False) + time += dt + + # Read stress from solver.tau — the actual projected stress + sigma_xy_val = uw.function.evaluate( + stokes.tau.sym[0, 1], centre + ) + sigma_xy = float(sigma_xy_val.flatten()[0]) + + times.append(time) + stress_num.append(sigma_xy) + stress_ana.append(maxwell_stress_xy(time, ETA, MU, gamma_dot)) + + del stokes, mesh + + return np.array(times), np.array(stress_num), np.array(stress_ana) + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + +class TestVEShearBox: + """Viscoelastic shear box with analytical Maxwell solution.""" + + def test_order1_converges(self): + """Order-1 VE stress should track the analytical Maxwell curve.""" + + times, num, ana = _run_ve_shear(order=1, n_steps=20, dt_over_tr=0.1) + + rel_error_final = abs(num[-1] - ana[-1]) / abs(ana[-1]) + mask = ana > 0.01 * ana[-1] + rel_error_rms = np.sqrt(np.mean(((num[mask] - ana[mask]) / ana[mask]) ** 2)) + + print(f"Order 1: final rel error = {rel_error_final:.4e}, " + f"RMS rel error = {rel_error_rms:.4e}") + + assert rel_error_final < 0.05 + assert rel_error_rms < 0.06 + + def test_order2_converges(self): + """Order-2 VE stress should converge more accurately than order-1.""" + + times, num, ana = _run_ve_shear(order=2, n_steps=20, dt_over_tr=0.1) + + rel_error_final = abs(num[-1] - ana[-1]) / abs(ana[-1]) + mask = ana > 0.01 * ana[-1] + rel_error_rms = np.sqrt(np.mean(((num[mask] - ana[mask]) / ana[mask]) ** 2)) + + print(f"Order 2: final rel error = {rel_error_final:.4e}, " + f"RMS rel error = {rel_error_rms:.4e}") + + # Order-2 should be significantly better than order-1 + assert rel_error_final < 0.005 + assert rel_error_rms < 0.02 + + def test_order2_better_than_order1(self): + """Order 2 should have smaller error than order 1 at the same Δt.""" + + _, num1, ana1 = _run_ve_shear(order=1, n_steps=20, dt_over_tr=0.1) + _, num2, ana2 = _run_ve_shear(order=2, n_steps=20, dt_over_tr=0.1) + + mask = ana1 > 0.01 * ana1[-1] + err1 = abs(num1[-1] - ana1[-1]) / ana1[-1] + err2 = abs(num2[-1] - ana2[-1]) / ana2[-1] + + print(f"Final errors: Order 1 = {err1:.4e}, Order 2 = {err2:.4e}") + + assert err2 < err1, ( + f"Order 2 error ({err2:.4e}) should be less than order 1 ({err1:.4e})" + ) + + def test_steady_state_approach(self): + """Stress should monotonically approach the viscous limit η·γ̇.""" + + times, num, ana = _run_ve_shear(order=1, n_steps=30, dt_over_tr=0.1) + + final_err = abs(num[-1] - ana[-1]) / ana[-1] + + print(f"Steady state approach: σ_xy = {num[-1]:.6f}, " + f"analytical = {ana[-1]:.6f}") + + diffs = np.diff(num) + assert np.all(diffs >= -1.0e-10), "Stress should be monotonically increasing" + assert final_err < 0.02 diff --git a/uw b/uw index 77929c6b..efd4db2a 100755 --- a/uw +++ b/uw @@ -159,8 +159,12 @@ run_build() { fi if [ "$need_clean" = true ]; then echo " PETSc target changed — cleaning build cache..." - rm -rf "$SCRIPT_DIR"/build/lib.* "$SCRIPT_DIR"/build/temp.* fi + # Always clean the setuptools build directory. --no-cache-dir only + # bypasses pip's wheel cache; the build/ tree (Cython .c files, + # compiled .so objects) persists across runs and silently reuses + # stale artefacts when only .py sources change. + rm -rf "$SCRIPT_DIR"/build/lib.* "$SCRIPT_DIR"/build/temp.* "$SCRIPT_DIR"/build/bdist.* mkdir -p "$SCRIPT_DIR/build" echo "$current_target" > "$petsc_marker" @@ -173,6 +177,43 @@ run_build() { exit 1 } + # Verify Python source files were actually installed. + # pip's wheel cache can silently serve stale .py files even with + # --no-cache-dir (version 0.0.0 problem). Compare checksums of key + # source files against the installed copies and warn if they differ. + local site_packages=$($PIXI run -e "$env" python -c \ + "import underworld3, pathlib; print(pathlib.Path(underworld3.__file__).parent)" 2>/dev/null) + + if [ -n "$site_packages" ] && [ -d "$site_packages" ]; then + local stale_files=() + # Check .py files under src/underworld3/ against installed copies + while IFS= read -r src_file; do + local rel_path="${src_file#src/underworld3/}" + local inst_file="$site_packages/$rel_path" + if [ -f "$inst_file" ]; then + local src_hash=$(shasum -a 256 "$src_file" | cut -d' ' -f1) + local inst_hash=$(shasum -a 256 "$inst_file" | cut -d' ' -f1) + if [ "$src_hash" != "$inst_hash" ]; then + stale_files+=("$rel_path") + fi + fi + done < <(find src/underworld3 -name '*.py' -not -path '*/__pycache__/*' 2>/dev/null) + + if [ ${#stale_files[@]} -gt 0 ]; then + echo "" + echo -e "${YELLOW}${BOLD}WARNING: ${#stale_files[@]} source file(s) differ from installed copies:${NC}" + for f in "${stale_files[@]}"; do + echo -e " ${YELLOW}STALE: $f${NC}" + done + echo "" + echo " Fixing: copying source files to installed location..." + for f in "${stale_files[@]}"; do + cp "src/underworld3/$f" "$site_packages/$f" + done + echo -e " ${GREEN}✓${NC} ${#stale_files[@]} file(s) updated" + fi + fi + echo "" echo -e "${GREEN}${BOLD}Build complete!${NC}" echo ""