From 1a17f2e7cf200107ddae05a693b1ad37f65e0472 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 18 Mar 2026 20:38:12 +1100 Subject: [PATCH 01/18] Always clean build cache in ./uw build The setuptools build/ directory (lib.*, temp.*, bdist.*) persists across runs and silently reuses stale .py artefacts when only Python sources change. --no-cache-dir only bypasses pip's wheel cache, not setuptools. Underworld development team with AI support from Claude Code --- uw | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/uw b/uw index 77929c6b..a3de4084 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" From b9ed2839c5983d297fae1d22897b480f0d47f1ed Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 18 Mar 2026 20:38:31 +1100 Subject: [PATCH 02/18] Explicit stress history for VE_Stokes solver MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replaces the old approach where SemiLagrangian re-evaluated the stress formula (psi_fn) at upstream points — which broke on order transitions because UWexpression η_eff evaluates live and the formula structure changed between order-1 and order-2. New architecture: store the actual stress tensor after each solve, advect only stored values. The constitutive formula is used only during the PETSc solve (JIT) and the post-solve stress projection. VE_Stokes.solve() flow per timestep: 1. advect_history() — pure advection of stored stress to upstream 2. PETSc solve using advected σ*, σ** 3. Save advected σ*, project constitutive_model.flux → psi_star[0] 4. Shift: psi_star[1] ← saved σ* (chained characteristic tracing) Additional fixes: - effective_order property on constitutive model (checks DDt startup) - effective_order formula: max(1, n_solves) not max(1, n_solves+1) - VE_Stokes re-setups JIT when effective_order transitions - elastic_dt → dt_elastic typo fix in VE_Stokes.solve() - Removed sympy.simplify() from 5 hot-path methods (caused hangs) - stress_deviator_stored property for post-processing access - TODO annotation for tensor vtype inference in _function.pyx Validated: Maxwell shear box test - Order-1: ~3% RMS error with dt/t_r=0.1 - Order-2: ~1% RMS, 0.09% final error (16× better than order-1) - Performance: constant 0.4-0.8s/step (was 5-10s with linear growth) Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 31 +++--- src/underworld3/function/_function.pyx | 6 ++ src/underworld3/systems/ddt.py | 136 ++++++++++++++++++++++++- src/underworld3/systems/solvers.py | 83 +++++++++++++-- 4 files changed, 233 insertions(+), 23 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 8d02acde..205f6ec7 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -971,12 +971,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 @@ -1201,7 +1201,7 @@ def ve_effective_viscosity(inner_self): # 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: + if inner_self._owning_model.effective_order != 2: el_eff_visc = ( inner_self.shear_viscosity_0 * inner_self.shear_modulus @@ -1252,6 +1252,19 @@ 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 + # The following should have no setters @property def stress_star(self): @@ -1287,7 +1300,7 @@ def E_eff(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.order != 2: + if self.effective_order != 2: stress_star = self.Unknowns.DFDt.psi_star[0].sym E += stress_star / ( 2 * self.Parameters.dt_elastic * self.Parameters.shear_modulus @@ -1468,8 +1481,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 +1503,6 @@ def stress_projection(self): / (self.Parameters.dt_elastic * self.Parameters.shear_modulus) ) - stress = sympy.simplify(stress) - return stress def stress(self): @@ -1508,7 +1517,7 @@ def stress(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.order != 2: + if self.effective_order != 2: stress_star = self.Unknowns.DFDt.psi_star[0].sym stress += ( 2 @@ -1534,8 +1543,6 @@ def stress(self): ) ) - stress = sympy.simplify(stress) - return stress # def eff_edot(self): 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..6edf74c1 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]``.""" @@ -679,7 +681,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.""" @@ -1113,7 +1117,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 +1190,7 @@ def update_post_solve( if self._n_solves_completed < self.order: self._n_solves_completed += 1 + return def update_pre_solve( @@ -1536,6 +1543,121 @@ def update_pre_solve( return + def advect_history(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) + 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 + + # 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 + def bdf(self, order=None): r"""Backwards differentiation form for calculating DuDt. @@ -1748,7 +1870,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`. @@ -2067,7 +2191,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`. diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 25c85b16..61e02328 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1246,6 +1246,20 @@ def delta_t(self): """Elastic timestep from the constitutive model.""" return self.constitutive_model.Parameters.dt_elastic + @property + def stress_deviator_stored(self): + r"""Deviatoric stress from the most recent solve (stored in history). + + Returns the actual projected stress from psi_star[0], which is + stored after each solve. Use this for post-processing, + visualization, and quantitative measurement. + + Note: the base ``stress_deviator`` property returns the constitutive + formula (needed for JIT compilation). For VE problems, always use + ``stress_deviator_stored`` to read the solved stress. + """ + return self.DFDt.psi_star[0].sym + ## Solver needs to update the stress history terms as well as call the SNES solve: @timing.routine_timer_decorator @@ -1274,12 +1288,23 @@ 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 + if timestep != self.delta_t.sym: + self._constitutive_model.Parameters.dt_elastic = timestep 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 +1314,20 @@ def solve( self._setup_discretisation(verbose) self._setup_solver(verbose) + # --- Explicit stress history management --- + # + # 1. ADVECT: trace psi_star values to upstream positions along + # characteristics. psi_star[0] contains the actual stress from + # the previous solve (stored by step 4 below). After advection, + # psi_star[0] = σ* (previous stress at upstream) and + # 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.advect_history(timestep, verbose=verbose, evalf=evalf) + + # 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 +1339,45 @@ 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) + + # Save advected σ* before anything modifies psi_star[0] + import numpy as np + _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) + + # Project actual stress into psi_star[0] + # Uses the constitutive formula (not the stored values property) 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: psi_star[1] ← saved advected σ* (for chained tracing next step) + for i in range(self.DFDt.order - 1, 0, -1): + if i == 1: + self.DFDt.psi_star[i].array[...] = _advected_sigma_star + else: + # For order > 2, shift higher levels down + 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) From b03b79923768da62ba2c7bdda51a879a88a26c6e Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 18 Mar 2026 20:38:41 +1100 Subject: [PATCH 03/18] Add VE shear box validation tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Pytest test (test_1051_VE_shear_box.py) validates Maxwell viscoelastic shear against analytical solution σ_xy = η·γ̇·(1 - exp(-t/t_r)). Tests order-1 convergence, order-2 convergence, order-2 > order-1, and monotonic approach to steady state. Quick validation scripts for development use: - run_ve_shear_quick.py: 10-step order-1 smoke test - run_ve_shear_order2_quick.py: 20-step order-2 with timing - run_ve_shear_validation.py: full comparison of both orders - run_ve_order2_debug.py: traces psi_star values per step Underworld development team with AI support from Claude Code --- tests/run_ve_order2_debug.py | 81 +++++++++++++ tests/run_ve_shear_order2_quick.py | 48 ++++++++ tests/run_ve_shear_quick.py | 48 ++++++++ tests/run_ve_shear_validation.py | 112 ++++++++++++++++++ tests/test_1051_VE_shear_box.py | 181 +++++++++++++++++++++++++++++ 5 files changed, 470 insertions(+) create mode 100644 tests/run_ve_order2_debug.py create mode 100644 tests/run_ve_shear_order2_quick.py create mode 100644 tests/run_ve_shear_quick.py create mode 100644 tests/run_ve_shear_validation.py create mode 100644 tests/test_1051_VE_shear_box.py 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_shear_order2_quick.py b/tests/run_ve_shear_order2_quick.py new file mode 100644 index 00000000..9e7d6d8c --- /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(ddt.psi_star[0].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..6be23858 --- /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 directly from psi_star[0] (the projected actual stress) + val = uw.function.evaluate(ddt.psi_star[0].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/test_1051_VE_shear_box.py b/tests/test_1051_VE_shear_box.py new file mode 100644 index 00000000..a8d2709a --- /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 psi_star[0] — the actual projected stress + sigma_xy_val = uw.function.evaluate( + stokes.DFDt.psi_star[0].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 > 0), "Stress should be monotonically increasing" + assert final_err < 0.02 From 660e6fec4e3507965f44a969acb98dd27a56a86f Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 18 Mar 2026 21:55:34 +1100 Subject: [PATCH 04/18] Add solver.tau property for computed flux access MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds a `tau` property to SolverBaseClass that returns a MeshVariable containing the projected constitutive flux. Created lazily on first access — no overhead if never used. Works for: - Stokes: deviatoric stress tensor (SYM_TENSOR projection) - Poisson: diffusive flux vector (Vector projection) - VE_Stokes: override returns psi_star[0] directly (no projection) This gives users a reliable way to access "what the solver actually computed" without manually reconstructing fluxes from the constitutive formula (which may have history terms, nonlinear corrections, etc). Usage: stokes.solve() stress_data = stokes.tau.data # numerical values stress_sym = stokes.tau.sym # symbolic (for further expressions) The base `stress_deviator` property is unchanged (returns the symbolic formula, needed for JIT compilation via F1 -> stress -> stress_deviator). Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 116 ++++++++++++++++++ src/underworld3/systems/solvers.py | 18 +-- tests/run_ve_shear_order2_quick.py | 2 +- tests/run_ve_shear_quick.py | 4 +- tests/test_1051_VE_shear_box.py | 4 +- 5 files changed, 131 insertions(+), 13 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 1f9cebea..8a003fa3 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. diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 61e02328..64acc2fc 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1247,18 +1247,20 @@ def delta_t(self): return self.constitutive_model.Parameters.dt_elastic @property - def stress_deviator_stored(self): + def tau(self): r"""Deviatoric stress from the most recent solve (stored in history). - Returns the actual projected stress from psi_star[0], which is - stored after each solve. Use this for post-processing, - visualization, and quantitative measurement. + For VE_Stokes, the stress is projected into ``psi_star[0]`` after + each solve. This override returns that variable directly — no + additional projection needed. - Note: the base ``stress_deviator`` property returns the constitutive - formula (needed for JIT compilation). For VE problems, always use - ``stress_deviator_stored`` to read the solved stress. + Returns + ------- + MeshVariable + The stress history variable containing the actual deviatoric + stress from the most recent solve. """ - return self.DFDt.psi_star[0].sym + return self.DFDt.psi_star[0] ## Solver needs to update the stress history terms as well as call the SNES solve: diff --git a/tests/run_ve_shear_order2_quick.py b/tests/run_ve_shear_order2_quick.py index 9e7d6d8c..fb40d71f 100644 --- a/tests/run_ve_shear_order2_quick.py +++ b/tests/run_ve_shear_order2_quick.py @@ -38,7 +38,7 @@ solve_t = timer.time() - t0 time_phys += dt - val = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre) + 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 " diff --git a/tests/run_ve_shear_quick.py b/tests/run_ve_shear_quick.py index 6be23858..debced41 100644 --- a/tests/run_ve_shear_quick.py +++ b/tests/run_ve_shear_quick.py @@ -37,8 +37,8 @@ solve_t = timer.time() - t0 time_phys += dt - # Read stress directly from psi_star[0] (the projected actual stress) - val = uw.function.evaluate(ddt.psi_star[0].sym[0, 1], centre) + # 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 " diff --git a/tests/test_1051_VE_shear_box.py b/tests/test_1051_VE_shear_box.py index a8d2709a..1691ee31 100644 --- a/tests/test_1051_VE_shear_box.py +++ b/tests/test_1051_VE_shear_box.py @@ -97,9 +97,9 @@ def _run_ve_shear(order, n_steps, dt_over_tr): stokes.solve(zero_init_guess=False, evalf=False) time += dt - # Read stress from psi_star[0] — the actual projected stress + # Read stress from solver.tau — the actual projected stress sigma_xy_val = uw.function.evaluate( - stokes.DFDt.psi_star[0].sym[0, 1], centre + stokes.tau.sym[0, 1], centre ) sigma_xy = float(sigma_xy_val.flatten()[0]) From 878ab873b69088329a9b5c979696822f5c243451 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Fri, 20 Mar 2026 13:35:59 +1100 Subject: [PATCH 05/18] Fix VEP order-2 plasticity and add VP/VEP shear box tests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The plastic effective viscosity in ViscoElasticPlasticFlowModel was using order-1 only for the elastic strain rate contribution, regardless of effective_order. This caused order-2 VEP to overshoot the yield stress (stress reached 0.73 instead of capping at τ_y=0.5). Fix: _plastic_effective_viscosity now uses self.E_eff.sym which already handles effective_order correctly via the BDF coefficients. New test scripts: - run_vp_shear_box.py: VP yield stress sweep (machine-precision match) - run_vep_shear_box.py: VEP elastic buildup + yield cap Order-1: caps at τ_y with ~1e-7 error Order-2: caps at τ_y with ~1e-7 error (was 22% error before fix) - run_ve_oscillatory.py: time-harmonic VE validation Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 14 +-- tests/run_ve_oscillatory.py | 123 +++++++++++++++++++++++++ tests/run_vep_shear_box.py | 121 ++++++++++++++++++++++++ tests/run_vp_shear_box.py | 90 ++++++++++++++++++ 4 files changed, 338 insertions(+), 10 deletions(-) create mode 100644 tests/run_ve_oscillatory.py create mode 100644 tests/run_vep_shear_box.py create mode 100644 tests/run_vp_shear_box.py diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 205f6ec7..6ff31729 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1378,16 +1378,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}'}", diff --git a/tests/run_ve_oscillatory.py b/tests/run_ve_oscillatory.py new file mode 100644 index 00000000..0c70ba0c --- /dev/null +++ b/tests/run_ve_oscillatory.py @@ -0,0 +1,123 @@ +"""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) = η γ̇₀ De/(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 η γ̇₀ De/√(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 * De / (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.is_setup = False # force BC re-evaluation + + 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_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") From 0a603d52e05890732f30967ecdb228bad993daef Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Fri, 20 Mar 2026 17:14:40 +1100 Subject: [PATCH 06/18] Add mesh.t time symbol and solve(time=) parameter (prototype) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Introduces mesh.t — a symbolic time coordinate that can be used in expressions for time-dependent boundary conditions and forcing terms. The symbol compiles to a C variable accessible in PETSc pointwise functions, avoiding JIT recompilation when time changes. Also adds optional time= parameter to Scalar and Stokes solve() methods, and UW_DMSetTime wrapper in petsc_compat.h. NOTE: petsc_t is hardcoded to PETSC_MIN_REAL in PETSc's SNES path (DMPlexSNESComputeResidualFEM, DMPlexSNESComputeBoundaryFEM). The planned solution is to use PetscDS constants[] array instead, which generalises to any mutable constant (viscosity, moduli, dt_elastic). This commit is the prototype; the constants wiring is the next step. New files: - TimeSymbol class in unit_aware_coordinates.py - _patch_time_units for dimensional analysis - UW_DMSetTime C wrapper (for future TS interface use) - VEP oscillatory test scripts and plot Underworld development team with AI support from Claude Code --- src/underworld3/cython/petsc_compat.h | 15 ++ src/underworld3/cython/petsc_extras.pxi | 1 + .../cython/petsc_generic_snes_solvers.pyx | 37 +++- .../discretisation/discretisation_mesh.py | 31 ++++ .../utilities/unit_aware_coordinates.py | 65 +++++++ tests/plot_ve_vep_oscillatory.py | 41 +++++ tests/run_ve_vep_oscillatory_checkpoint.py | 125 ++++++++++++++ tests/run_ve_vep_oscillatory_plot.py | 162 ++++++++++++++++++ tests/run_vep_oscillatory.py | 135 +++++++++++++++ 9 files changed, 610 insertions(+), 2 deletions(-) create mode 100644 tests/plot_ve_vep_oscillatory.py create mode 100644 tests/run_ve_vep_oscillatory_checkpoint.py create mode 100644 tests/run_ve_vep_oscillatory_plot.py create mode 100644 tests/run_vep_oscillatory.py 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 8a003fa3..4c0da8a8 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1775,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. @@ -1798,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 ------- @@ -1835,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: @@ -3871,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. @@ -3898,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 ------- @@ -3939,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..1db90856 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,29 @@ 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 solver sets the time value via ``solve(time=t)``. If ``time`` + is not provided, ``petsc_t`` defaults to 0. + + 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/utilities/unit_aware_coordinates.py b/src/underworld3/utilities/unit_aware_coordinates.py index 304aadeb..26471db6 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, positive=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_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_vep_oscillatory_checkpoint.py b/tests/run_ve_vep_oscillatory_checkpoint.py new file mode 100644 index 00000000..706cc2ce --- /dev/null +++ b/tests/run_ve_vep_oscillatory_checkpoint.py @@ -0,0 +1,125 @@ +"""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 * De / (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.is_setup = False + + 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 * De / 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..b5b50c0d --- /dev/null +++ b/tests/run_ve_vep_oscillatory_plot.py @@ -0,0 +1,162 @@ +"""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 * De / (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.is_setup = False + + 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 * De / 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..3f767579 --- /dev/null +++ b/tests/run_vep_oscillatory.py @@ -0,0 +1,135 @@ +"""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 η γ̇₀ De/√(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 * De / (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 * De / 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 + stokes.is_setup = False + + 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 * De / 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() From 0ccdddf58add28e06b533da06ccc1c89673a55c4 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sat, 21 Mar 2026 21:56:20 +1100 Subject: [PATCH 07/18] Add BDF-3 VE constitutive formulae, fix dt_elastic handling Constitutive model changes: - Add effective_order property that delegates to DDt's ramp-up counter - ve_effective_viscosity: three-way branch for BDF-1/2/3 effective viscosity (eta_eff = eta*mu*dt / (c0*eta + mu*dt) where c0 = 1, 3/2, 11/6) - E_eff: three-way branch with correct BDF coefficients for history terms - stress(): three-way branch matching the BDF discretisation VE_Stokes.solve() changes: - Fix elastic_dt -> dt_elastic parameter name bug - Remove dt_elastic overwrite: solve(timestep=) no longer clobbers the constitutive relaxation parameter with the advection timestep - Improved docstring clarifying timestep vs dt_elastic semantics Convergence verified against Maxwell analytical solution: dt=0.2 O1: 3.0e-2 O2: 4.0e-3 O3: 6.4e-3 dt=0.1 O1: 1.5e-2 O2: 9.3e-4 O3: 1.7e-3 dt=0.05 O1: 7.8e-3 O2: 2.3e-4 O3: 4.3e-4 Order 2 converges at 2nd order, Order 3 at ~2nd order (BDF-3 error constant larger for this stiffness ratio). Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- src/underworld3/constitutive_models.py | 100 +++++++++++-------------- src/underworld3/systems/solvers.py | 49 ++++++++---- 2 files changed, 79 insertions(+), 70 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 6ff31729..df7c26b3 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1198,32 +1198,19 @@ 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.effective_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 - ) - ) - - # 2nd Order version (need to ask for this one) + # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt) + # c0: 1 (BDF-1), 3/2 (BDF-2), 11/6 (BDF-3) + eta = inner_self.shear_viscosity_0 + mu = inner_self.shear_modulus + dt_e = inner_self.dt_elastic + eff_order = inner_self._owning_model.effective_order + + if eff_order == 2: + el_eff_visc = 2 * eta * mu * dt_e / (3 * eta + 2 * mu * dt_e) + elif eff_order >= 3: + el_eff_visc = 6 * eta * mu * dt_e / (11 * eta + 6 * mu * dt_e) 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 / (eta + mu * dt_e) inner_self._ve_effective_viscosity.sym = el_eff_visc @@ -1300,20 +1287,23 @@ def E_eff(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.effective_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 + eff_order = self.effective_order - else: + if eff_order == 2: 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 - ) + E += stress_star / mu_dt - stress_2star / (4 * mu_dt) + elif eff_order >= 3: + stress_star = self.Unknowns.DFDt.psi_star[0].sym + stress_2star = self.Unknowns.DFDt.psi_star[1].sym + stress_3star = self.Unknowns.DFDt.psi_star[2].sym + E += (3 * stress_star / (2 * mu_dt) + - 3 * stress_2star / (4 * mu_dt) + + stress_3star / (6 * mu_dt)) + else: + stress_star = self.Unknowns.DFDt.psi_star[0].sym + E += stress_star / (2 * mu_dt) self._E_eff.sym = E @@ -1511,30 +1501,28 @@ def stress(self): if self.Unknowns.DFDt is not None: if self.is_elastic: - if self.effective_order != 2: + mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus + eff_order = self.effective_order + + if eff_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) - ) + stress_2star = self.Unknowns.DFDt.psi_star[1].sym + stress += 2 * self.viscosity * ( + stress_star / mu_dt - stress_2star / (4 * mu_dt) ) - - else: + elif eff_order >= 3: 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) - ) + stress_3star = self.Unknowns.DFDt.psi_star[2].sym + stress += 2 * self.viscosity * ( + 3 * stress_star / (2 * mu_dt) + - 3 * stress_2star / (4 * mu_dt) + + stress_3star / (6 * mu_dt) + ) + else: + stress_star = self.Unknowns.DFDt.psi_star[0].sym + stress += 2 * self.viscosity * ( + stress_star / (2 * mu_dt) ) return stress diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 64acc2fc..3437baf8 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1290,8 +1290,9 @@ def solve( if timestep is None: timestep = self.delta_t.sym - if timestep != self.delta_t.sym: - self._constitutive_model.Parameters.dt_elastic = timestep + # 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 @@ -1359,25 +1360,45 @@ def solve( if uw.mpi.rank == 0 and verbose: print(f"VE Stokes solver - store stress and shift history", flush=True) - # Save advected σ* before anything modifies psi_star[0] import numpy as np - _advected_sigma_star = np.copy(self.DFDt.psi_star[0].array[...]) + import sympy as _sympy - # Project actual stress into psi_star[0] - # Uses the constitutive formula (not the stored values property) so - # that σ* and σ** from psi_star are read correctly during projection. + # Running average blending factor. + # When timestep == dt_elastic: phi = 1 → hard copy (standard behaviour). + # When timestep < dt_elastic: phi < 1 → gradual accumulation so each + # history slot represents approximately one dt_elastic of elapsed time. + dt_elastic_value = self.delta_t.sym + try: + _phi = float(_sympy.Min(1, timestep / dt_elastic_value)) + except (TypeError, ValueError): + _phi = 1.0 + + # Save advected psi_star values before projection modifies psi_star[0]. + # We need psi_star[i-1] for the running average shift. + _saved = [] + for i in range(self.DFDt.order): + _saved.append(np.copy(self.DFDt.psi_star[i].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: psi_star[1] ← saved advected σ* (for chained tracing next step) + # Now psi_star[0] = projected τ (actual stress from this solve). + # + # Shift history with running average: + # psi_star[i] ← phi * saved[i-1] + (1-phi) * saved[i] + # + # When phi=1 (dt==dt_elastic): straight copy (classical behaviour). + # When phi<1 (dt 2, shift higher levels down - self.DFDt.psi_star[i].array[...] = self.DFDt.psi_star[i - 1].array[...] + self.DFDt.psi_star[i].array[...] = ( + _phi * _saved[i - 1] + (1 - _phi) * _saved[i] + ) # 5. BOOKKEEPING From 30133e91afc633bc4feaa2495993735c36ff4934 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sat, 21 Mar 2026 21:57:21 +1100 Subject: [PATCH 08/18] Add post-build verification to ./uw build After pip install, checksums all .py source files against the installed copies. If any differ (stale wheel cache), copies the source files directly to the installed location and warns. This catches the recurring issue where pip serves stale Python files despite --no-cache-dir because the setuptools build directory caches them. Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- uw | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/uw b/uw index a3de4084..efd4db2a 100755 --- a/uw +++ b/uw @@ -177,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 "" From c967b1802d5015e3f6d1094e489e25dc17262469 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 22 Mar 2026 14:12:55 +1100 Subject: [PATCH 09/18] Fix analytical formula, add oscillatory benchmark docs, remove running average MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Fix Maxwell analytical solution: remove spurious De factor in prefactor. Was η γ̇₀ De/(1+De²), correct is η γ̇₀/(1+De²). The bug was invisible at De=1 (factor cancels) but caused apparent "amplitude error" at other De. - Add oscillatory shear benchmark documentation (docs/advanced/benchmarks/) with convergence tables and resolution study description. - Add plot_ve_oscillatory_validation.py: generates validation figures for De=1.5 (phase lag, amplitude) with --replot support for saved .npz data. - Remove running-average history smoothing from VE_Stokes.solve(): the BDF constitutive formula couples dt_elastic to the time between evaluations, so the exponential moving average cannot compensate for dt << dt_elastic. Reverted to simple history shift (direct copy-down). Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- docs/advanced/benchmarks/index.md | 13 + .../benchmarks/ve-oscillatory-shear.md | 109 ++++++++ src/underworld3/systems/solvers.py | 36 +-- tests/plot_ve_oscillatory_validation.py | 239 ++++++++++++++++++ tests/run_ve_oscillatory.py | 2 +- 5 files changed, 371 insertions(+), 28 deletions(-) create mode 100644 docs/advanced/benchmarks/index.md create mode 100644 docs/advanced/benchmarks/ve-oscillatory-shear.md create mode 100644 tests/plot_ve_oscillatory_validation.py 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..a7663bcd --- /dev/null +++ b/docs/advanced/benchmarks/ve-oscillatory-shear.md @@ -0,0 +1,109 @@ +--- +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 inaccurate: the BDF constitutive formula +inherently couples `dt_elastic` to the time between stress evaluations, so +the history smoothing cannot compensate for the mismatch. For problems requiring +small advection steps with a longer relaxation scale, viscosity clamping +(flooring $\eta_{\text{eff}}$ at the `dt_elastic`-derived value) is the +recommended approach. diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 3437baf8..962aa874 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1361,23 +1361,9 @@ def solve( print(f"VE Stokes solver - store stress and shift history", flush=True) import numpy as np - import sympy as _sympy - # Running average blending factor. - # When timestep == dt_elastic: phi = 1 → hard copy (standard behaviour). - # When timestep < dt_elastic: phi < 1 → gradual accumulation so each - # history slot represents approximately one dt_elastic of elapsed time. - dt_elastic_value = self.delta_t.sym - try: - _phi = float(_sympy.Min(1, timestep / dt_elastic_value)) - except (TypeError, ValueError): - _phi = 1.0 - - # Save advected psi_star values before projection modifies psi_star[0]. - # We need psi_star[i-1] for the running average shift. - _saved = [] - for i in range(self.DFDt.order): - _saved.append(np.copy(self.DFDt.psi_star[i].array[...])) + # 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 @@ -1387,18 +1373,14 @@ def solve( self.DFDt._psi_star_projection_solver.solve(verbose=verbose) # Now psi_star[0] = projected τ (actual stress from this solve). - # - # Shift history with running average: - # psi_star[i] ← phi * saved[i-1] + (1-phi) * saved[i] - # - # When phi=1 (dt==dt_elastic): straight copy (classical behaviour). - # When phi<1 (dt Date: Sun, 22 Mar 2026 14:19:37 +1100 Subject: [PATCH 10/18] Update dt_elastic guidance: running average is diffusive, advise viscosity limiting Underworld development team with AI support from Claude Code (https://claude.com/claude-code) --- docs/advanced/benchmarks/ve-oscillatory-shear.md | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/advanced/benchmarks/ve-oscillatory-shear.md b/docs/advanced/benchmarks/ve-oscillatory-shear.md index a7663bcd..03e95c43 100644 --- a/docs/advanced/benchmarks/ve-oscillatory-shear.md +++ b/docs/advanced/benchmarks/ve-oscillatory-shear.md @@ -101,9 +101,8 @@ 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 inaccurate: the BDF constitutive formula -inherently couples `dt_elastic` to the time between stress evaluations, so -the history smoothing cannot compensate for the mismatch. For problems requiring -small advection steps with a longer relaxation scale, viscosity clamping -(flooring $\eta_{\text{eff}}$ at the `dt_elastic`-derived value) is the -recommended approach. +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. From b7bf3136d52b74b10a35cc20d1dd817b5cac5636 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Mon, 23 Mar 2026 13:23:45 +1100 Subject: [PATCH 11/18] Fix stale analytical formulas and minor issues from PR review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Remove spurious De factor from Maxwell oscillatory prefactor and steady-state amplitude in 5 test/validation scripts (code and docstrings). The corrected formula is η γ̇₀/(1+De²), not η γ̇₀ De/(1+De²). - TimeSymbol: change positive=True to nonnegative=True (t=0 is valid) - test_1051: relax strict monotonicity check (diffs > 0 → >= -1e-10) to avoid flaky failures from solver roundoff - mesh.t docstring: clarify that high-level Python solve() wrappers do not yet pass time= through to PETSc Underworld development team with AI support from Claude Code --- src/underworld3/discretisation/discretisation_mesh.py | 7 +++++-- src/underworld3/utilities/unit_aware_coordinates.py | 2 +- tests/plot_ve_oscillatory_validation.py | 4 ++-- tests/run_ve_oscillatory.py | 4 ++-- tests/run_ve_vep_oscillatory_checkpoint.py | 4 ++-- tests/run_ve_vep_oscillatory_plot.py | 4 ++-- tests/run_vep_oscillatory.py | 8 ++++---- tests/test_1051_VE_shear_box.py | 2 +- 8 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 1db90856..4b11dbd6 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1513,8 +1513,11 @@ def t(self): and Jacobian functions. Use ``mesh.t`` in expressions to reference this time without forcing JIT recompilation each timestep. - The solver sets the time value via ``solve(time=t)``. If ``time`` - is not provided, ``petsc_t`` defaults to 0. + 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 diff --git a/src/underworld3/utilities/unit_aware_coordinates.py b/src/underworld3/utilities/unit_aware_coordinates.py index 26471db6..0073f283 100644 --- a/src/underworld3/utilities/unit_aware_coordinates.py +++ b/src/underworld3/utilities/unit_aware_coordinates.py @@ -96,7 +96,7 @@ class TimeSymbol(sympy.Symbol): _units = None def __new__(cls, name="t", **kwargs): - obj = super().__new__(cls, name, real=True, positive=True, **kwargs) + obj = super().__new__(cls, name, real=True, nonnegative=True, **kwargs) return obj @property diff --git a/tests/plot_ve_oscillatory_validation.py b/tests/plot_ve_oscillatory_validation.py index d4dbf525..f3ac103f 100644 --- a/tests/plot_ve_oscillatory_validation.py +++ b/tests/plot_ve_oscillatory_validation.py @@ -5,7 +5,7 @@ - Bottom panel: stress σ_xy(t) for order 1 and order 2 vs analytical Maxwell analytical solution including startup transient: - σ_xy(t) = η γ̇₀ De/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)] + σ_xy(t) = η γ̇₀/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)] Results are saved as .npz checkpoint files for re-analysis. @@ -137,7 +137,7 @@ def make_plot(results, params, output_path): shear_rate = gamma_dot_0 * np.sin(omega * t_fine) phase_lag_deg = np.degrees(np.arctan(De)) - steady_amp = ETA * gamma_dot_0 * De / np.sqrt(1 + De**2) + steady_amp = ETA * gamma_dot_0 / np.sqrt(1 + De**2) period = 2.0 * np.pi / omega fig, (ax1, ax2) = plt.subplots( diff --git a/tests/run_ve_oscillatory.py b/tests/run_ve_oscillatory.py index be5c996b..98242026 100644 --- a/tests/run_ve_oscillatory.py +++ b/tests/run_ve_oscillatory.py @@ -3,12 +3,12 @@ Maxwell material under sinusoidal boundary velocity V(t) = V0 sin(ωt). Full solution including startup transient: - σ_xy(t) = η γ̇₀ De/(1+De²) [sin(ωt) - De cos(ωt) + De exp(-t/t_r)] + σ_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 η γ̇₀ De/√(1+De²) and phase lag arctan(De). +oscillates with amplitude η γ̇₀/√(1+De²) and phase lag arctan(De). """ import time as timer diff --git a/tests/run_ve_vep_oscillatory_checkpoint.py b/tests/run_ve_vep_oscillatory_checkpoint.py index 706cc2ce..b14529d5 100644 --- a/tests/run_ve_vep_oscillatory_checkpoint.py +++ b/tests/run_ve_vep_oscillatory_checkpoint.py @@ -11,7 +11,7 @@ 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 * De / (1.0 + De**2) + 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)) @@ -88,7 +88,7 @@ def run_oscillatory(order, n_steps, dt, omega, V0, tau_y=None): n_periods = 2 n_steps = int(n_periods * period / dt) - ve_amp = ETA * gamma_dot_0 * De / np.sqrt(1.0 + De**2) + 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}") diff --git a/tests/run_ve_vep_oscillatory_plot.py b/tests/run_ve_vep_oscillatory_plot.py index b5b50c0d..79ba8dcf 100644 --- a/tests/run_ve_vep_oscillatory_plot.py +++ b/tests/run_ve_vep_oscillatory_plot.py @@ -18,7 +18,7 @@ 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 * De / (1.0 + De**2) + 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)) @@ -108,7 +108,7 @@ def run_oscillatory(order, n_steps, dt, omega, V0, t0, tau_y=None): n_periods = 2 n_steps = int(n_periods * period / dt) - ve_amp = ETA * gamma_dot_0 * De / np.sqrt(1.0 + De**2) + 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}") diff --git a/tests/run_vep_oscillatory.py b/tests/run_vep_oscillatory.py index 3f767579..9831b3e8 100644 --- a/tests/run_vep_oscillatory.py +++ b/tests/run_vep_oscillatory.py @@ -1,7 +1,7 @@ """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 η γ̇₀ De/√(1+De²). +The VE amplitude at steady state is η γ̇₀/√(1+De²). If τ_y is below this amplitude, the stress should be clipped. """ @@ -16,7 +16,7 @@ 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 * De / (1.0 + De**2) + 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)) @@ -31,7 +31,7 @@ def run_vep_oscillatory(order, n_steps, dt_over_tr, De, tau_y): dt = dt_over_tr * t_r # VE steady-state amplitude - ve_amplitude = ETA * gamma_dot_0 * De / np.sqrt(1 + De**2) + 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), @@ -100,7 +100,7 @@ def run_vep_oscillatory(order, n_steps, dt_over_tr, De, tau_y): n_steps = int(n_periods * period / (dt_ratio * t_r)) ETA, gamma_dot_0 = 1.0, 1.0 - ve_amp = ETA * gamma_dot_0 * De / np.sqrt(1 + De**2) + 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 diff --git a/tests/test_1051_VE_shear_box.py b/tests/test_1051_VE_shear_box.py index 1691ee31..1608a729 100644 --- a/tests/test_1051_VE_shear_box.py +++ b/tests/test_1051_VE_shear_box.py @@ -177,5 +177,5 @@ def test_steady_state_approach(self): f"analytical = {ana[-1]:.6f}") diffs = np.diff(num) - assert np.all(diffs > 0), "Stress should be monotonically increasing" + assert np.all(diffs >= -1.0e-10), "Stress should be monotonically increasing" assert final_err < 0.02 From 267bb7026a47b41385ed248bc20c6a303431b3a8 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 09:14:51 +1100 Subject: [PATCH 12/18] Variable-dt BDF coefficients for viscoelastic constitutive model The DDt system already computes correct variable-timestep BDF coefficients via _bdf_coefficients(order, dt_n, dt_history), but the constitutive model hardcoded uniform-step values (c0 = 3/2 for BDF-2, 11/6 for BDF-3). When dt varies between timesteps, the hardcoded coefficients are wrong, silently degrading accuracy. This commit: - Adds bdf_coefficients property to all DDt classes (SemiLagrangian, Eulerian, Lagrangian, Lagrangian_Swarm) exposing the current coefficients - Adds _get_bdf_coefficients() helper to ViscoElasticPlasticFlowModel that pulls from DDt when available, falls back to uniform otherwise - Rewrites ve_effective_viscosity, E_eff, and stress() to use DDt coefficients instead of hardcoded order-based branching - Adds square-wave forcing benchmark (tests/run_ve_square_wave.py) for validating variable-dt accuracy via Fourier superposition of Maxwell analytical solutions For uniform dt, results are identical (same coefficient values). Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 63 ++----- src/underworld3/systems/ddt.py | 25 +++ tests/run_ve_square_wave.py | 249 +++++++++++++++++++++++++ 3 files changed, 294 insertions(+), 43 deletions(-) create mode 100644 tests/run_ve_square_wave.py diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index df7c26b3..b50687fe 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 @@ -1199,18 +1200,14 @@ def ve_effective_viscosity(inner_self): return inner_self.shear_viscosity_0 # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt) - # c0: 1 (BDF-1), 3/2 (BDF-2), 11/6 (BDF-3) + # c0 from DDt's BDF coefficients (variable-dt aware) eta = inner_self.shear_viscosity_0 mu = inner_self.shear_modulus dt_e = inner_self.dt_elastic - eff_order = inner_self._owning_model.effective_order + coeffs = inner_self._owning_model._get_bdf_coefficients() + c0 = coeffs[0] - if eff_order == 2: - el_eff_visc = 2 * eta * mu * dt_e / (3 * eta + 2 * mu * dt_e) - elif eff_order >= 3: - el_eff_visc = 6 * eta * mu * dt_e / (11 * eta + 6 * mu * dt_e) - else: - el_eff_visc = eta * mu * dt_e / (eta + mu * dt_e) + el_eff_visc = eta * mu * dt_e / (c0 * eta + mu * dt_e) inner_self._ve_effective_viscosity.sym = el_eff_visc @@ -1252,6 +1249,12 @@ def effective_order(self): return min(self._order, self.Unknowns.DFDt.effective_order) return self._order + def _get_bdf_coefficients(self): + """Get BDF coefficients from DDt (variable-dt aware) or fall back to uniform.""" + if self.Unknowns is not None and self.Unknowns.DFDt is not None: + return self.Unknowns.DFDt.bdf_coefficients + return _bdf_coefficients(self.effective_order, None, []) + # The following should have no setters @property def stress_star(self): @@ -1288,22 +1291,11 @@ def E_eff(self): if self.is_elastic: mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus - eff_order = self.effective_order - - if eff_order == 2: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - E += stress_star / mu_dt - stress_2star / (4 * mu_dt) - elif eff_order >= 3: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - stress_3star = self.Unknowns.DFDt.psi_star[2].sym - E += (3 * stress_star / (2 * mu_dt) - - 3 * stress_2star / (4 * mu_dt) - + stress_3star / (6 * mu_dt)) - else: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - E += stress_star / (2 * mu_dt) + coeffs = self._get_bdf_coefficients() + + # History contribution: -Σ cᵢ·σ_star[i-1] / (2·μ·dt) + for i in range(1, len(coeffs)): + E += -coeffs[i] * self.Unknowns.DFDt.psi_star[i - 1].sym / (2 * mu_dt) self._E_eff.sym = E @@ -1502,27 +1494,12 @@ def stress(self): if self.is_elastic: mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus - eff_order = self.effective_order + coeffs = self._get_bdf_coefficients() - if eff_order == 2: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - stress += 2 * self.viscosity * ( - stress_star / mu_dt - stress_2star / (4 * mu_dt) - ) - elif eff_order >= 3: - stress_star = self.Unknowns.DFDt.psi_star[0].sym - stress_2star = self.Unknowns.DFDt.psi_star[1].sym - stress_3star = self.Unknowns.DFDt.psi_star[2].sym - stress += 2 * self.viscosity * ( - 3 * stress_star / (2 * mu_dt) - - 3 * stress_2star / (4 * mu_dt) - + stress_3star / (6 * mu_dt) - ) - else: - stress_star = self.Unknowns.DFDt.psi_star[0].sym + # History contribution: 2·η_eff · (-Σ cᵢ·σ_star[i-1]) / (2·μ·dt) + for i in range(1, len(coeffs)): stress += 2 * self.viscosity * ( - stress_star / (2 * mu_dt) + -coeffs[i] * self.Unknowns.DFDt.psi_star[i - 1].sym / (2 * mu_dt) ) return stress diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 6edf74c1..292de866 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -408,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] @@ -806,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. @@ -1658,6 +1668,11 @@ def advect_history(self, dt, evalf=False, verbose=False): 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. @@ -1978,6 +1993,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. @@ -2300,6 +2320,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/tests/run_ve_square_wave.py b/tests/run_ve_square_wave.py new file mode 100644 index 00000000..788b7551 --- /dev/null +++ b/tests/run_ve_square_wave.py @@ -0,0 +1,249 @@ +"""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: oscillatory shear with square-wave envelope + x, y = mesh.X + t_sym = sympy.Symbol("t_sim") + t_expr = expression(R"t_{\mathrm{sim}}", t_sym, "Simulation time") + + # Build the Fourier-series velocity symbolically (truncated) + v_top_sym = sympy.Integer(0) + for k in range(1, n_harmonics + 1): + n = 2 * k - 1 + a_k = sympy.Rational(4, 1) / (sympy.pi * n) + v_top_sym += a_k * sympy.sin(n * omega * t_sym) + v_top_sym *= V0 + + v_bc_sym = v_top_sym * (y / (H / 2)) + v_bc = expression(R"v_{\mathrm{bc}}", v_bc_sym, "Square-wave shear velocity") + + stokes.add_dirichlet_bc((v_bc.sym, 0.0), "Top") + stokes.add_dirichlet_bc((-v_bc.sym, 0.0), "Bottom") + stokes.add_dirichlet_bc((None, 0.0), "Left") + stokes.add_dirichlet_bc((None, 0.0), "Right") + + # 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) + + # Update simulation time and dt_elastic + t_expr.sym = t_current + dt + stokes.constitutive_model.Parameters.dt_elastic = dt + + stokes.solve(timestep=dt) + + t_current += dt + step += 1 + + # Extract stress at mesh centre + tau = stokes.tau + centre = np.array([[0.0, 0.0]]) + tau_xy = uw.function.evaluate(tau.sym[0, 1], centre, mesh) + 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 = 15 + + 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") From e3412c3b92d6af44e57b4705195e36b5e51dff1b Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 14:00:12 +1100 Subject: [PATCH 13/18] Fix time-varying BCs: pass UWexpression directly, not .sym MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Passing expr.sym to add_dirichlet_bc strips the UWexpression wrapper, causing the value to be baked as a C literal in the JIT. The PetscDS constants mechanism only works when the UWexpression atom is visible to _extract_constants during compilation. - Square-wave benchmark: pass V_top (not V_top.sym) to BCs - Remove is_setup=False workaround from 4 oscillatory test scripts (run_ve_oscillatory, run_vep_oscillatory, run_ve_vep_oscillatory_plot, run_ve_vep_oscillatory_checkpoint) — unnecessary when UWexpression is passed directly Underworld development team with AI support from Claude Code --- tests/run_ve_oscillatory.py | 1 - tests/run_ve_square_wave.py | 43 ++++++++++------------ tests/run_ve_vep_oscillatory_checkpoint.py | 1 - tests/run_ve_vep_oscillatory_plot.py | 1 - tests/run_vep_oscillatory.py | 1 - 5 files changed, 19 insertions(+), 28 deletions(-) diff --git a/tests/run_ve_oscillatory.py b/tests/run_ve_oscillatory.py index 98242026..2cf7e21d 100644 --- a/tests/run_ve_oscillatory.py +++ b/tests/run_ve_oscillatory.py @@ -67,7 +67,6 @@ def run_oscillatory(order, n_steps, dt_over_tr, De): # Update boundary velocity for this timestep V_t = V0 * np.sin(omega * time_phys) V_bc.sym = V_t - stokes.is_setup = False # force BC re-evaluation stokes.solve(zero_init_guess=False, evalf=False) diff --git a/tests/run_ve_square_wave.py b/tests/run_ve_square_wave.py index 788b7551..36b2a833 100644 --- a/tests/run_ve_square_wave.py +++ b/tests/run_ve_square_wave.py @@ -96,26 +96,15 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA stokes.constitutive_model.Parameters.shear_modulus = MU - # Boundary conditions: oscillatory shear with square-wave envelope - x, y = mesh.X - t_sym = sympy.Symbol("t_sim") - t_expr = expression(R"t_{\mathrm{sim}}", t_sym, "Simulation time") + # 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") - # Build the Fourier-series velocity symbolically (truncated) - v_top_sym = sympy.Integer(0) - for k in range(1, n_harmonics + 1): - n = 2 * k - 1 - a_k = sympy.Rational(4, 1) / (sympy.pi * n) - v_top_sym += a_k * sympy.sin(n * omega * t_sym) - v_top_sym *= V0 - - v_bc_sym = v_top_sym * (y / (H / 2)) - v_bc = expression(R"v_{\mathrm{bc}}", v_bc_sym, "Square-wave shear velocity") - - stokes.add_dirichlet_bc((v_bc.sym, 0.0), "Top") - stokes.add_dirichlet_bc((-v_bc.sym, 0.0), "Bottom") - stokes.add_dirichlet_bc((None, 0.0), "Left") - stokes.add_dirichlet_bc((None, 0.0), "Right") + 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 = [] @@ -131,19 +120,25 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, else: dt = adaptive_dt(t_current, omega, dt_min, dt_max) - # Update simulation time and dt_elastic - t_expr.sym = t_current + dt + 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(timestep=dt) + stokes.solve(zero_init_guess=False, timestep=dt) t_current += dt step += 1 # Extract stress at mesh centre - tau = stokes.tau centre = np.array([[0.0, 0.0]]) - tau_xy = uw.function.evaluate(tau.sym[0, 1], centre, mesh) + tau_xy = uw.function.evaluate(stokes.tau.sym[0, 1], centre) sigma_xy = float(tau_xy.flatten()[0]) times.append(t_current) diff --git a/tests/run_ve_vep_oscillatory_checkpoint.py b/tests/run_ve_vep_oscillatory_checkpoint.py index b14529d5..8ae07f01 100644 --- a/tests/run_ve_vep_oscillatory_checkpoint.py +++ b/tests/run_ve_vep_oscillatory_checkpoint.py @@ -57,7 +57,6 @@ def run_oscillatory(order, n_steps, dt, omega, V0, tau_y=None): for step in range(n_steps): time_phys += dt V_bc.sym = V0 * np.sin(omega * time_phys) - stokes.is_setup = False stokes.solve(zero_init_guess=False, evalf=False) diff --git a/tests/run_ve_vep_oscillatory_plot.py b/tests/run_ve_vep_oscillatory_plot.py index 79ba8dcf..b5ee4907 100644 --- a/tests/run_ve_vep_oscillatory_plot.py +++ b/tests/run_ve_vep_oscillatory_plot.py @@ -77,7 +77,6 @@ def run_oscillatory(order, n_steps, dt, omega, V0, t0, tau_y=None): for step in range(n_steps): time_phys += dt V_bc.sym = V0 * np.sin(omega * time_phys) - stokes.is_setup = False stokes.solve(zero_init_guess=False, evalf=False) diff --git a/tests/run_vep_oscillatory.py b/tests/run_vep_oscillatory.py index 9831b3e8..e76e212e 100644 --- a/tests/run_vep_oscillatory.py +++ b/tests/run_vep_oscillatory.py @@ -70,7 +70,6 @@ def run_vep_oscillatory(order, n_steps, dt_over_tr, De, tau_y): time_phys += dt V_t = V0 * np.sin(omega * time_phys) V_bc.sym = V_t - stokes.is_setup = False t0s = __import__('time').time() stokes.solve(zero_init_guess=False, evalf=False) From db3caaa3a3ffa116540a347dd1b6a7a0cd5df971 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 15:22:59 +1100 Subject: [PATCH 14/18] Replace advect_history() with standard DDt update_pre_solve pathway The VE_Stokes solver now uses update_pre_solve(advect_only=True) instead of the custom advect_history() method. This reuses the battle-tested advection code (coordinate handling, global_evaluate, return_coords_to_bounds) rather than reimplementing it. advect_only=True skips the history shift and psi_fn evaluation steps, since the VE solver manages those in its post-solve sequence. Only the upstream advection of existing psi_star levels is performed. The old advect_history() is parked as _advect_history_PARKED for reference. Also fixes: return_coords_to_bounds calls in advect path, benchmark BC pattern (pass UWexpression directly, 3 harmonics). Regression: 4/4 VE shear box tests pass (order 1 + order 2). Underworld development team with AI support from Claude Code --- src/underworld3/systems/ddt.py | 67 +++++++++++++++++++----------- src/underworld3/systems/solvers.py | 17 ++++---- tests/run_ve_square_wave.py | 2 +- 3 files changed, 53 insertions(+), 33 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 292de866..8d4e89d4 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1209,11 +1209,21 @@ def update_pre_solve( evalf: Optional[bool] = False, verbose: Optional[bool] = False, dt_physical: Optional[float] = None, + advect_only: Optional[bool] = False, ): """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 + ---------- + advect_only : bool, optional + If True, skip the psi_fn evaluation into psi_star[0] (step 2) + and only perform history shift (step 1) and upstream advection + (step 3). Used by VE_Stokes where psi_star[0] already contains + the projected actual stress from the previous solve — we want to + advect that stored stress, not overwrite it with the flux expression. """ if not self._history_initialised: @@ -1241,15 +1251,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 not advect_only: + 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 advect_only=True (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: @@ -1288,28 +1302,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 not advect_only: + 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 @@ -1553,7 +1568,7 @@ def update_pre_solve( return - def advect_history(self, dt, evalf=False, verbose=False): + 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. @@ -1636,6 +1651,7 @@ def advect_history(self, dt, evalf=False, verbose=False): # 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, :] @@ -1650,6 +1666,7 @@ def advect_history(self, dt, evalf=False, verbose=False): # 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 diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 962aa874..a24f3b1c 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1317,18 +1317,21 @@ def solve( self._setup_discretisation(verbose) self._setup_solver(verbose) - # --- Explicit stress history management --- + # --- Stress history management via standard DDt pathway --- # - # 1. ADVECT: trace psi_star values to upstream positions along - # characteristics. psi_star[0] contains the actual stress from - # the previous solve (stored by step 4 below). After advection, - # psi_star[0] = σ* (previous stress at upstream) and - # psi_star[1] = σ** (stress from 2 steps ago, double-traced). + # update_pre_solve(advect_only=True) performs: + # 1. History shift: psi_star[i] ← psi_star[i-1] + # 2. Skip psi_fn evaluation (psi_star[0] already has projected stress) + # 3. Advect all history levels to upstream positions along characteristics + # + # After this: psi_star[0] = σ* (previous stress at upstream), + # psi_star[1] = σ** (stress from 2 steps ago, double-traced). if uw.mpi.rank == 0 and verbose: print(f"VE Stokes solver - advect stress history", flush=True) - self.DFDt.advect_history(timestep, verbose=verbose, evalf=evalf) + self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, + advect_only=True) # 2. SOLVE: PETSc uses the advected σ*, σ** via the constitutive model diff --git a/tests/run_ve_square_wave.py b/tests/run_ve_square_wave.py index 36b2a833..d56ffc73 100644 --- a/tests/run_ve_square_wave.py +++ b/tests/run_ve_square_wave.py @@ -183,7 +183,7 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, De = 1.5 order = 2 n_periods = 3 - n_harmonics = 15 + n_harmonics = 3 print("=" * 60) print(f"Square-wave VE benchmark: De={De}, order={order}") From 55962d10ed3d7ebc4f03e8a0dbbf8bde09f717c5 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 15:36:32 +1100 Subject: [PATCH 15/18] Rename advect_only parameter to store_result (inverted sense) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit store_result=True (default): normal DDt behaviour, evaluate psi_fn store_result=False: skip evaluation and history shift, advect only Clearer name — says what it does rather than describing a mode. Underworld development team with AI support from Claude Code --- src/underworld3/systems/ddt.py | 21 +++++++++++---------- src/underworld3/systems/solvers.py | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 8d4e89d4..400c33b8 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1209,7 +1209,7 @@ def update_pre_solve( evalf: Optional[bool] = False, verbose: Optional[bool] = False, dt_physical: Optional[float] = None, - advect_only: Optional[bool] = False, + store_result: Optional[bool] = True, ): """Sample upstream values along characteristics before solve. @@ -1218,12 +1218,13 @@ def update_pre_solve( Parameters ---------- - advect_only : bool, optional - If True, skip the psi_fn evaluation into psi_star[0] (step 2) - and only perform history shift (step 1) and upstream advection - (step 3). Used by VE_Stokes where psi_star[0] already contains - the projected actual stress from the previous solve — we want to - advect that stored stress, not overwrite it with the flux expression. + 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: @@ -1251,7 +1252,7 @@ def update_pre_solve( else: phi = sympy.sympify(1) - if not advect_only: + 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[...] @@ -1261,7 +1262,7 @@ def update_pre_solve( # Note the need to do a try/except to handle unsupported evaluations # (e.g. of derivatives) # - # When advect_only=True (VE stress history), skip this step — + # 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. @@ -1302,7 +1303,7 @@ def update_pre_solve( :, : ] - if not advect_only: + if store_result: try: # Use shifted ND coords to avoid quad mesh boundary issues # node_coords_nd is slightly shifted toward cell centroids diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index a24f3b1c..c85568f6 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1331,7 +1331,7 @@ def solve( print(f"VE Stokes solver - advect stress history", flush=True) self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, - advect_only=True) + store_result=False) # 2. SOLVE: PETSc uses the advected σ*, σ** via the constitutive model From 180986ad94e85098fd88a16d885f5f198d388ca6 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 19:17:46 +1100 Subject: [PATCH 16/18] Route BDF coefficients through PetscDS constants as UWexpressions The BDF coefficients (c0, c1, c2, c3) are now UWexpressions on the constitutive model, updated each step by _update_bdf_coefficients(). They flow through PetscDS constants[] to the compiled pointwise functions without JIT recompilation. Previously, _bdf_coefficients() returned numeric values that were baked into the C code at JIT time. Changing dt between steps had no effect on the compiled coefficients. Now: ve_effective_viscosity, E_eff, and stress() reference the UWexpression _bdf_c0.._bdf_c3 symbolically. The solver calls _update_bdf_coefficients() before each solve, which computes the variable-dt BDF coefficients and sets the .sym values. PetscDSSetConstants pushes these to PETSc at solve time. Verified: alternating dt test shows Order 2 (err=0.020) outperforms Order 1 (err=0.034), with variable/uniform ratio ~1.06 (near-ideal). Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 56 ++++++++++++++++++++------ src/underworld3/systems/solvers.py | 6 +++ 2 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index b50687fe..47293a68 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1087,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) @@ -1200,12 +1208,12 @@ def ve_effective_viscosity(inner_self): return inner_self.shear_viscosity_0 # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt) - # c0 from DDt's BDF coefficients (variable-dt aware) + # 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 - coeffs = inner_self._owning_model._get_bdf_coefficients() - c0 = coeffs[0] + c0 = inner_self._owning_model._bdf_c0 el_eff_visc = eta * mu * dt_e / (c0 * eta + mu * dt_e) @@ -1249,11 +1257,32 @@ def effective_order(self): return min(self._order, self.Unknowns.DFDt.effective_order) return self._order - def _get_bdf_coefficients(self): - """Get BDF coefficients from DDt (variable-dt aware) or fall back to uniform.""" + 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: - return self.Unknowns.DFDt.bdf_coefficients - return _bdf_coefficients(self.effective_order, 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 @@ -1291,11 +1320,12 @@ def E_eff(self): if self.is_elastic: mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus - coeffs = self._get_bdf_coefficients() + # BDF history coefficients as UWexpressions (route through constants[]) + bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] # History contribution: -Σ cᵢ·σ_star[i-1] / (2·μ·dt) - for i in range(1, len(coeffs)): - E += -coeffs[i] * self.Unknowns.DFDt.psi_star[i - 1].sym / (2 * mu_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 @@ -1494,12 +1524,12 @@ def stress(self): if self.is_elastic: mu_dt = self.Parameters.dt_elastic * self.Parameters.shear_modulus - coeffs = self._get_bdf_coefficients() + bdf_cs = [self._bdf_c1, self._bdf_c2, self._bdf_c3] # History contribution: 2·η_eff · (-Σ cᵢ·σ_star[i-1]) / (2·μ·dt) - for i in range(1, len(coeffs)): + for i in range(self.Unknowns.DFDt.order): stress += 2 * self.viscosity * ( - -coeffs[i] * self.Unknowns.DFDt.psi_star[i - 1].sym / (2 * mu_dt) + -bdf_cs[i] * self.Unknowns.DFDt.psi_star[i].sym / (2 * mu_dt) ) return stress diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index c85568f6..0bf156ed 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1333,6 +1333,12 @@ def solve( 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: From c94e02c33662453364829632c9278be66b3a5604 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 20:32:58 +1100 Subject: [PATCH 17/18] Remove deprecated constitutive_models_new.py Failed rebuild attempt, scheduled for removal, not imported anywhere. Contains stale hardcoded BDF patterns that would confuse future work. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models_new.py | 1702 -------------------- 1 file changed, 1702 deletions(-) delete mode 100644 src/underworld3/constitutive_models_new.py diff --git a/src/underworld3/constitutive_models_new.py b/src/underworld3/constitutive_models_new.py deleted file mode 100644 index d58be059..00000000 --- a/src/underworld3/constitutive_models_new.py +++ /dev/null @@ -1,1702 +0,0 @@ -r""" -Alternative constitutive model implementations. - -.. deprecated:: - This module contains a failed attempt to rebuild the constitutive models - and is scheduled for removal. Use :mod:`underworld3.constitutive_models` - instead for all constitutive model needs. - -See Also --------- -underworld3.constitutive_models : Primary constitutive model module (use this). -""" -from typing_extensions import Self -import sympy -from sympy import sympify -from sympy.vector import gradient, divergence -import numpy as np - -from typing import Optional, Callable -from typing import NamedTuple, Union - -from petsc4py import PETSc - -import underworld3 as uw -from underworld3.utilities._api_tools import uw_object -from underworld3.systems import SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt -from underworld3.swarm import IndexSwarmVariable -import underworld3.timing as timing - - -class Constitutive_Model(uw_object): - r""" - Base class for constitutive models. - - Constitutive laws relate gradients in the unknowns to fluxes of quantities - (e.g., heat fluxes related to temperature gradients via thermal conductivity). - - For scalar problems: - - .. math:: - - q_i = k_{ij} \frac{\partial T}{\partial x_j} - - where :math:`k_{ij}` are the constitutive parameters. The template assumes - :math:`k_{ij} = \delta_{ij}` (identity). - - For vector problems (e.g., Stokes): - - .. math:: - - t_{ij} = c_{ijkl} \frac{\partial u_k}{\partial x_l} - - Usually written with symmetrized gradients: - - .. math:: - - t_{ij} = c_{ijkl} \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} - + \frac{\partial u_l}{\partial x_k} \right] - - where :math:`c_{ijkl}` are the constitutive parameters. The template assumes - :math:`c_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})`, - the 4th-rank identity tensor with flux and gradient symmetry. - """ - - @timing.routine_timer_decorator - def __init__( - self, - dim: int, - u_dim: int, - ): - # Define / identify the various properties in the class but leave - # the implementation to child classes. The constitutive tensor is - # defined as a template here, but should be instantiated via class - # properties as required. - - # We provide a function that converts gradients / gradient history terms - # into the relevant flux term. - - self.dim = dim - self.u_dim = u_dim - self._solver = None - - self.Parameters = self._Parameters() - self.Parameters._solver = None - self.Parameters._reset = self._reset - - self._material_properties = None - - ## Default consitutive tensor is the identity - - if self.u_dim == 1: - self._c = sympy.Matrix.eye(self.dim) - else: # vector problem - self._c = uw.maths.tensor.rank4_identity(self.dim) - - self._C = None - - super().__init__() - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - """ - - def __init__(inner_self, k=1): - inner_self._k = k - inner_self._solver = None - - @property - def k(inner_self): - return inner_self._k - - @k.setter - def k(inner_self, value): - inner_self._k = value - inner_self._reset() - - """ - @property - def material_properties(self): - - return self._material_properties - - @material_properties.setter - def material_properties(self, properties): - - if isinstance(properties, self.Parameters): - self._material_properties = properties - else: - name = self.__class__.__name__ - raise RuntimeError(f"Use {name}.material_properties = {name}.Parameters(...) ") - - d = self.dim - self._build_c_tensor() - - if isinstance(self._solver, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)): - self._solver.is_setup = False - - return - """ - - @property - def solver(self): - """Each constitutive relationship can, optionally, be associated with one solver object. - and a solver object _requires_ a constitive relationship to be defined.""" - return self._solver - - @solver.setter - def solver(self, solver_object): - if isinstance(solver_object, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)): - self._solver = solver_object - self.Parameters._solver = solver_object - self._solver.is_setup = False - - ## Properties on all sub-classes - - @property - def C(self): - """The matrix form of the constitutive model (the `c` property) - that relates fluxes to gradients. - For scalar problem, this is the matrix representation of the rank 2 tensor. - For vector problems, the Mandel form of the rank 4 tensor is returned. - NOTE: this is an immutable object that is _a view_ of the underlying tensor - """ - - d = self.dim - rank = len(self.c.shape) - - if rank == 2: - return sympy.Matrix(self._c).as_immutable() - else: - return uw.maths.tensor.rank4_to_mandel(self._c, d).as_immutable() - - @property - def c(self): - """The tensor form of the constitutive model that relates fluxes to gradients. In scalar - problems, `c` and `C` are equivalent (matrices), but in vector problems, `c` is a - rank 4 tensor. NOTE: `c` is the canonical form of the constitutive relationship. - """ - - return self._c.as_immutable() - - def flux( - self, - ddu: sympy.Matrix = None, - ddu_dt: sympy.Matrix = None, - u: sympy.Matrix = None, # may be needed in the case of cylindrical / spherical - u_dt: sympy.Matrix = None, - ): - """Computes the effect of the constitutive tensor on the gradients of the unknowns. - (always uses the `c` form of the tensor). In general cases, the history of the gradients - may be required to evaluate the flux. - """ - - c = self.c - rank = len(c.shape) - - # tensor multiplication - - if rank == 2: - flux = c * ddu.T - else: # rank==4 - flux = sympy.tensorcontraction( - sympy.tensorcontraction(sympy.tensorproduct(c, ddu), (3, 5)), (0, 1) - ) - - return sympy.Matrix(flux) - - def flux_1d( - self, - ddu: sympy.Matrix = None, - ddu_dt: sympy.Matrix = None, - u: sympy.Matrix = None, # may be needed in the case of cylindrical / spherical - u_dt: sympy.Matrix = None, - ): - """Computes the effect of the constitutive tensor on the gradients of the unknowns. - (always uses the `c` form of the tensor). In general cases, the history of the gradients - may be required to evaluate the flux. Returns the Voigt form that is flattened so as to - match the PETSc field storage pattern for symmetric tensors. - """ - - flux = self.flux(ddu, ddu_dt, u, u_dt) - - assert ( - flux.is_symmetric() - ), "The conversion to Voigt form is only defined for symmetric tensors in underworld" - - return uw.maths.tensor.rank2_to_voigt(flux, dim=self.dim) - - def _reset(self): - d = self.dim - self._build_c_tensor() - - if isinstance(self._solver, (SNES_Scalar, SNES_Vector, SNES_Stokes_SaddlePt)): - self._solver.is_setup = False - - return - - def _build_c_tensor(self): - """Return the identity tensor of appropriate rank (e.g. for projections)""" - - d = self.dim - self._c = self.Parameters.k * uw.maths.tensor.rank4_identity(d) - - return - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - from textwrap import dedent - - display( - Markdown( - rf"This consititutive model is formulated for {self.dim} dimensional equations" - ) - ) - - -class ViscousFlowModel(Constitutive_Model): - r""" - Viscous flow constitutive model. - - .. math:: - - \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} - + \frac{\partial u_l}{\partial x_k} \right] - - where :math:`\eta` is the viscosity—a scalar constant, sympy function, - or mesh variable. This gives an isotropic (but not necessarily homogeneous - or linear) relationship between :math:`\tau` and velocity gradients. You can - also supply :math:`\eta_{IJ}` (Mandel form) or :math:`\eta_{ijkl}` (rank-4 tensor). - - The Mandel constitutive matrix is in ``viscous_model.C`` and the rank-4 - tensor form is in ``viscous_model.c``. - - Examples - -------- - >>> viscous_model = ViscousFlowModel(dim) - >>> viscous_model.material_properties = viscous_model.Parameters(viscosity=viscosity_fn) - >>> solver.constitutive_model = viscous_model - >>> tau = viscous_model.flux(gradient_matrix) - """ - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - """ - - def __init__( - inner_self, - viscosity: Union[float, sympy.Function] = None, - ): - if viscosity is None: - viscosity = sympy.sympify(1) - - inner_self._viscosity = sympy.sympify(viscosity) - - @property - def viscosity(inner_self): - return inner_self._viscosity - - @viscosity.setter - def viscosity(inner_self, value: Union[float, sympy.Function]): - inner_self._viscosity = value - inner_self._reset() - - def __init__(self, dim): - u_dim = dim - super().__init__(dim, u_dim) - - return - - def _build_c_tensor(self): - """For this constitutive law, we expect just a viscosity function""" - - d = self.dim - viscosity = self.Parameters.viscosity - - try: - self._c = 2 * uw.maths.tensor.rank4_identity(d) * viscosity - except: - d = self.dim - dv = uw.maths.tensor.idxmap[d][0] - if isinstance(viscosity, sympy.Matrix) and viscosity.shape == (dv, dv): - self._c = 2 * uw.maths.tensor.mandel_to_rank4(viscosity, d) - elif isinstance(viscosity, sympy.Array) and viscosity.shape == (d, d, d, d): - self._c = 2 * viscosity - else: - raise RuntimeError( - "Viscosity is not a known type (scalar, Mandel matrix, or rank 4 tensor" - ) - return - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - - super()._object_viewer() - - ## feedback on this instance - display(Latex(r"$\quad\eta = $ " + sympy.sympify(self.Parameters.viscosity)._repr_latex_())) - - -class ViscoPlasticFlowModel(ViscousFlowModel): - r""" - Viscoplastic flow constitutive model with yield stress. - - .. math:: - - \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} - + \frac{\partial u_l}{\partial x_k} \right] - - where :math:`\eta` is the viscosity—a scalar constant, sympy function, - or mesh variable. This gives an isotropic relationship between :math:`\tau` - and velocity gradients. You can also supply :math:`\eta_{IJ}` (Mandel form) - or :math:`\eta_{ijkl}` (rank-4 tensor). - - In a viscoplastic model, the viscosity is defined to cap the overall stress - at the *yield stress*. This assumes the yield stress is a scalar limit on - the 2nd invariant of the stress. Anisotropic models require careful yield - surface definition—only a subset of cases is available. - - The Mandel matrix is in ``viscoplastic_model.C`` and the rank-4 tensor in - ``viscoplastic_model.c``. - - Notes - ----- - If ``not~yet~defined`` appears in the effective viscosity, not all required - functions have been set. The model defaults to standard viscous behavior - if yield terms are not specified. - - Examples - -------- - >>> viscoplastic_model = ViscoPlasticFlowModel(dim) - >>> viscoplastic_model.material_properties = viscoplastic_model.Parameters( - ... viscosity=viscosity_fn, - ... yield_stress=yieldstress_fn, - ... min_viscosity=min_viscosity_fn, - ... max_viscosity=max_viscosity_fn, - ... strain_rate_II=strain_rate_inv_fn - ... ) - >>> solver.constitutive_model = viscoplastic_model - >>> tau = viscoplastic_model.flux(gradient_matrix) - """ - - # Init for VP class (not needed ??) - def __init__(self, dim): - super().__init__(dim) - - return - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - - `sympy.oo` (infinity) for default values ensures that sympy.Min simplifies away - the conditionals when they are not required - """ - - def __init__( - inner_self, - materialIndex: Union[uw.swarm.SwarmVariable, uw.discretisation.MeshVariable] = None, - shear_viscosity_0: Union[list, sympy.Function] = [1], - shear_viscosity_min: Union[list, sympy.Function] = [-sympy.oo], - shear_viscosity_max: Union[list, sympy.Function] = [sympy.oo], - yield_stress: Union[list, sympy.Function] = [sympy.oo], - yield_stress_min: Union[list, sympy.Function] = [-sympy.oo], - strainrate_inv_II: sympy.Function = sympy.oo, - strainrate_inv_II_min: float = 0.0, - averaging_method: str = "HA", - ): - if strainrate_inv_II is sympy.oo: - strainrate_inv_II = sympy.symbols( - r"\left|\dot\epsilon\right|\rightarrow\textrm{not\ defined}" - ) - - inner_self._shear_viscosity_0 = sympy.sympify(shear_viscosity_0) - inner_self._shear_viscosity_min = sympy.sympify(shear_viscosity_min) - inner_self._shear_viscosity_max = sympy.sympify(shear_viscosity_max) - - inner_self._yield_stress = sympy.sympify(yield_stress) - inner_self._yield_stress_min = sympy.sympify(yield_stress_min) - - inner_self._strainrate_inv_II = sympy.sympify(strainrate_inv_II) - inner_self._strainrate_inv_II_min = sympy.sympify(strainrate_inv_II_min) - - inner_self._averaging_method = averaging_method - inner_self._materialIndex = materialIndex - - return - - @property - def shear_viscosity_0(inner_self): - return inner_self._shear_viscosity_0 - - @shear_viscosity_0.setter - def shear_viscosity_0(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_0 = value - inner_self._reset() - - @property - def shear_viscosity_min(inner_self): - return inner_self._shear_viscosity_min - - @shear_viscosity_min.setter - def shear_viscosity_min(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_min = value - inner_self._reset() - - @property - def shear_viscosity_max(inner_self): - return inner_self._shear_viscosity_max - - @shear_viscosity_max.setter - def shear_viscosity_max(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_max = value - inner_self._reset() - - @property - def yield_stress(inner_self): - return inner_self._yield_stress - - @yield_stress.setter - def yield_stress(inner_self, value: Union[list, sympy.Function]): - inner_self._yield_stress = value - inner_self._reset() - - @property - def yield_stress_min(inner_self): - return inner_self._yield_stress_min - - @yield_stress_min.setter - def yield_stress_min(inner_self, value: Union[list, sympy.Function]): - inner_self._yield_stress_min = value - inner_self._reset() - - @property - def strainrate_inv_II(inner_self): - return inner_self._strainrate_inv_II - - @strainrate_inv_II.setter - def strainrate_inv_II(inner_self, value: sympy.Function): - inner_self._strainrate_inv_II = value - inner_self._reset() - - @property - def strainrate_inv_II_min(inner_self): - return inner_self._strainrate_inv_II_min - - @strainrate_inv_II_min.setter - def strainrate_inv_II_min(inner_self, value: float): - inner_self._strainrate_inv_II_min = sympy.sympify(value) - inner_self._reset() - - @property - def averaging_method(inner_self): - return inner_self._averaging_method - - @averaging_method.setter - def averaging_method(inner_self, value: str): - inner_self._averaging_method = value - inner_self._reset() - - ### Getter and setter for internel mask Variable - @property - def materialIndex(inner_self): - return inner_self._materialIndex - - @materialIndex.setter - def materialIndex(inner_self, indexVar): - # error checking, only support IndexSwarmVariables for now - if isinstance(indexVar, uw.swarm.IndexSwarmVariable): - inner_self._materialIndex = indexVar - inner_self._reset() - - # This has no setter !! - @property - def plastic_eff_viscosity(inner_self): - import warnings - - if ( - type(inner_self.yield_stress) != np.ndarray - and type(inner_self.yield_stress) != list - ): - inner_self.yield_stress = list([inner_self.yield_stress]) - if ( - type(inner_self.yield_stress_min) != np.ndarray - and type(inner_self.yield_stress_min) != list - ): - inner_self.yield_stress_min = list([inner_self.yield_stress_min]) - - if inner_self.materialIndex == None: - warnings.warn( - "materialIndex not specified, using the first value for each parameter", - stacklevel=2, - ) - - yield_stress = inner_self.yield_stress[0] - yield_stress_min = inner_self.yield_stress_min[0] - - if yield_stress_min == -sympy.oo: - yield_stress_fn = yield_stress - else: - yield_stress_fn = sympy.Max(yield_stress, yield_stress_min) - - if yield_stress_fn == sympy.oo or inner_self.strainrate_inv_II == sympy.oo: - pl_effective_viscosity = sympy.oo - - else: - pl_effective_viscosity = yield_stress_fn / ( - (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min - ) - - else: - ### creates list of values that has the same length as the material index - if len(inner_self.yield_stress) != inner_self.materialIndex.indices: - if len(inner_self.yield_stress) > 1: - warnings.warn( - f"Number of values in yield_stress ({len(inner_self.yield_stress)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.yield_stress = list( - np.repeat(inner_self.yield_stress[0], inner_self.materialIndex.indices) - ) - - if len(inner_self.yield_stress_min) != inner_self.materialIndex.indices: - if len(inner_self.yield_stress_min) > 1: - warnings.warn( - f"Number of values in yield_stress ({len(inner_self.yield_stress_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.yield_stress_min = list( - np.repeat( - inner_self.yield_stress_min[0], - inner_self.materialIndex.indices, - ) - ) - - yield_stress = sympy.Matrix(inner_self.yield_stress) - yield_stress_min = sympy.Matrix(inner_self.yield_stress_min) - - if yield_stress_min[0] == -sympy.oo: - yield_stress_fn = yield_stress - else: - yield_stress_list = [] - - for i in range(len(yield_stress)): - yield_stress_list.append(sympy.Max(yield_stress[i], yield_stress_min[i])) - - yield_stress_fn = sympy.Matrix(yield_stress_list) - - if yield_stress_fn[0] == sympy.oo or inner_self.strainrate_inv_II == sympy.oo: - pl_effective_viscosity = list( - np.repeat(sympy.oo, inner_self.materialIndex.indices) - ) - - else: - pl_effective_viscosity = yield_stress_fn / ( - (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min - ) - - return pl_effective_viscosity - - # This has no setter !! - @property - def viscosity(inner_self): - import warnings - - if ( - type(inner_self.shear_viscosity_0) != np.ndarray - and type(inner_self.shear_viscosity_0) != list - ): - inner_self.shear_viscosity_0 = list([inner_self.shear_viscosity_0]) - - if ( - type(inner_self.shear_viscosity_min) != np.ndarray - and type(inner_self.shear_viscosity_min) != list - ): - inner_self.shear_viscosity_min = list([inner_self.shear_viscosity_min]) - if ( - type(inner_self.shear_viscosity_max) != np.ndarray - and type(inner_self.shear_viscosity_max) != list - ): - inner_self.shear_viscosity_max = list([inner_self.shear_viscosity_max]) - - if inner_self.materialIndex == None: - warnings.warn( - "materialIndex not specified, using the first value for each parameter", - stacklevel=2, - ) - - shear_viscosity_min = inner_self.shear_viscosity_min[0] - shear_viscosity_max = inner_self.shear_viscosity_max[0] - - shear_viscosity_0 = inner_self.shear_viscosity_0[0] - - yield_visc = inner_self.plastic_eff_viscosity - - if inner_self.averaging_method.casefold() == "min": - if yield_visc != 0: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max( - shear_viscosity_min, - sympy.Min(yield_visc, shear_viscosity_0), - ), - ) - else: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max(shear_viscosity_min, shear_viscosity_0), - ) - else: - if yield_visc != 0: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max( - shear_viscosity_min, - 1.0 / ((1.0 / shear_viscosity_0) + (1.0 / yield_visc)), - ), - ) - else: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max(shear_viscosity_min, shear_viscosity_0), - ) - - else: - if len(inner_self.shear_viscosity_0) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_0) > 1: - warnings.warn( - f"Number of values in shear_viscosity_0 ({len(inner_self.shear_viscosity_0)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_0 = list( - np.repeat( - inner_self.shear_viscosity_0[0], - inner_self.materialIndex.indices, - ) - ) - - if len(inner_self.shear_viscosity_min) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_min) > 1: - warnings.warn( - f"Number of values in shear_viscosity_min ({len(inner_self.shear_viscosity_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_min = list( - np.repeat( - inner_self.shear_viscosity_min[0], - inner_self.materialIndex.indices, - ) - ) - - if len(inner_self.shear_viscosity_max) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_max) > 1: - warnings.warn( - f"Number of values in shear_viscosity_max ({len(inner_self.shear_viscosity_max)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_max = list( - np.repeat( - inner_self.shear_viscosity_max[0], - inner_self.materialIndex.indices, - ) - ) - - shear_viscosity_min = sympy.Matrix(inner_self.shear_viscosity_min) - shear_viscosity_max = sympy.Matrix(inner_self.shear_viscosity_max) - - shear_viscosity_0 = inner_self.shear_viscosity_0 - - yield_visc = inner_self.plastic_eff_viscosity - - viscosity_list = [] - for i in range(inner_self.materialIndex.indices): - if inner_self.averaging_method.casefold() == "min": - if yield_visc[i] != 0: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max( - shear_viscosity_min[i], - sympy.Min(yield_visc[i], shear_viscosity_0[i]), - ), - ) - ) - else: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max(shear_viscosity_min[i], shear_viscosity_0[i]), - ) - ) - else: - if yield_visc[i] != 0: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max( - shear_viscosity_min[i], - 1 / ((1 / shear_viscosity_0[i]) + (1 / yield_visc[i])), - ), - ) - ) - else: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max(shear_viscosity_min[i], shear_viscosity_0[i]), - ) - ) - - viscosity_fn = sympy.Matrix(viscosity_list) - - effective_viscosity = inner_self.materialIndex.sym.T.dot(viscosity_fn) - - return effective_viscosity - - ## ===== End of parameters sub_class - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - - super()._object_viewer() - - ## feedback on this instance - display( - Latex( - r"$\quad\eta_\textrm{0} = $ " - + sympy.sympify(self.Parameters.shear_viscosity_0)._repr_latex_() - ), - Latex( - r"$\quad\tau_\textrm{y} = $ " - + sympy.sympify(self.Parameters.yield_stress)._repr_latex_(), - ), - Latex( - r"$\quad|\dot\epsilon| = $ " - + sympy.sympify(self.Parameters.strainrate_inv_II)._repr_latex_(), - ), - ## Todo: add all the other properties in here - ) - - -class ViscoElasticPlasticFlowModel(Constitutive_Model): - r""" - Viscoelastic-plastic flow constitutive model. - - .. math:: - - \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} - + \frac{\partial u_l}{\partial x_k} \right] - - where :math:`\eta` is the viscosity—a scalar constant, sympy function, - or mesh variable. This gives an isotropic relationship between :math:`\tau` - and velocity gradients. You can also supply :math:`\eta_{IJ}` (Mandel form) - or :math:`\eta_{ijkl}` (rank-4 tensor). - - The Mandel matrix is in ``viscoelastic_model.C`` and the rank-4 tensor in - ``viscoelastic_model.c``. - - Examples - -------- - >>> viscoelastic_model = ViscoElasticFlowModel(dim) - >>> viscoelastic_model.material_properties = viscoelastic_model.Parameters( - ... viscosity=viscosity_fn - ... ) - >>> solver.constitutive_model = viscoelastic_model - >>> tau = viscoelastic_model.flux(gradient_matrix) - """ - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - """ - - def __init__( - inner_self, - materialIndex: Union[uw.swarm.SwarmVariable, uw.discretisation.MeshVariable] = None, - shear_viscosity_0: Union[list, sympy.Function] = [1], - shear_modulus: Union[list, sympy.Function] = [sympy.oo], - shear_viscosity_min: Union[list, sympy.Function] = [-sympy.oo], - shear_viscosity_max: Union[list, sympy.Function] = [sympy.oo], - yield_stress: Union[list, sympy.Function] = [sympy.oo], - yield_stress_min: Union[list, sympy.Function] = [-sympy.oo], - strainrate_inv_II: sympy.Function = sympy.oo, - strainrate_inv_II_min: float = 0.0, - averaging_method: str = "HA", - stress_star: sympy.Function = None, - stress_star_star: sympy.Function = None, - dt_elastic: Union[float, sympy.Function] = [sympy.oo], - ): - if strainrate_inv_II is None: - strainrate_inv_II = sympy.symbols( - r"\left|\dot\epsilon\right|\rightarrow\textrm{not\ defined}" - ) - - if stress_star is None: - stress_star = sympy.symbols(r"\sigma^*\rightarrow\textrm{not\ defined}") - - inner_self._shear_viscosity_0 = sympy.sympify(shear_viscosity_0) - inner_self._shear_modulus = sympy.sympify(shear_modulus) - inner_self._dt_elastic = sympy.sympify(dt_elastic) - inner_self._yield_stress = sympy.sympify(yield_stress) - inner_self._yield_stress_min = sympy.sympify(yield_stress_min) - inner_self._shear_viscosity_min = sympy.sympify(shear_viscosity_min) - inner_self._shear_viscosity_max = sympy.sympify(shear_viscosity_max) - inner_self._strainrate_inv_II = sympy.sympify(strainrate_inv_II) - inner_self._stress_star = sympy.sympify(stress_star) - inner_self._stress_star_star = sympy.sympify(stress_star_star) - inner_self._strainrate_inv_II_min = sympy.sympify(strainrate_inv_II_min) - - inner_self._averaging_method = averaging_method - inner_self._materialIndex = materialIndex - - return - - @property - def shear_viscosity_0(inner_self): - return inner_self._shear_viscosity_0 - - @shear_viscosity_0.setter - def shear_viscosity_0(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_0 = value - inner_self._reset() - - @property - def shear_modulus(inner_self): - return inner_self._shear_modulus - - @shear_modulus.setter - def shear_modulus(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_modulus = value - inner_self._reset() - - @property - def dt_elastic(inner_self): - return inner_self._dt_elastic - - @dt_elastic.setter - def dt_elastic(inner_self, value: Union[list, sympy.Function]): - inner_self._dt_elastic = value - inner_self._reset() - - @property - def shear_viscosity_min(inner_self): - return inner_self._shear_viscosity_min - - @shear_viscosity_min.setter - def shear_viscosity_min(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_min = value - inner_self._reset() - - @property - def shear_viscosity_max(inner_self): - return inner_self._shear_viscosity_max - - @shear_viscosity_max.setter - def shear_viscosity_max(inner_self, value: Union[list, sympy.Function]): - inner_self._shear_viscosity_max = value - inner_self._reset() - - @property - def yield_stress(inner_self): - return inner_self._yield_stress - - @yield_stress.setter - def yield_stress(inner_self, value: Union[list, sympy.Function]): - inner_self._yield_stress = value - inner_self._reset() - - @property - def yield_stress_min(inner_self): - return inner_self._yield_stress_min - - @yield_stress_min.setter - def yield_stress_min(inner_self, value: Union[list, sympy.Function]): - inner_self._yield_stress_min = value - inner_self._reset() - - @property - def strainrate_inv_II(inner_self): - return inner_self._strainrate_inv_II - - @strainrate_inv_II.setter - def strainrate_inv_II(inner_self, value: sympy.Function): - inner_self._strainrate_inv_II = value - inner_self._reset() - - @property - def stress_star(inner_self): - return inner_self._stress_star - - @stress_star.setter - def stress_star(inner_self, value: sympy.Function): - inner_self._stress_star = value - inner_self._reset() - - @property - def stress_star_star(inner_self): - return inner_self._stress_star_star - - @stress_star_star.setter - def stress_star_star(inner_self, value: sympy.Function): - inner_self._stress_star_star = value - inner_self._reset() - - @property - def strainrate_inv_II_min(inner_self): - return inner_self._strainrate_inv_II_min - - @strainrate_inv_II_min.setter - def strainrate_inv_II_min(inner_self, value: float): - inner_self._strainrate_inv_II_min = sympy.sympify(value) - inner_self._reset() - - @property - def averaging_method(inner_self): - return inner_self._averaging_method - - @averaging_method.setter - def averaging_method(inner_self, value: str): - inner_self._averaging_method = value - inner_self._reset() - - ### Getter and setter for internel mask Variable - @property - def materialIndex(inner_self): - return inner_self._materialIndex - - @materialIndex.setter - def materialIndex(inner_self, indexVar): - # error checking, only support IndexSwarmVariables for now - if isinstance(indexVar, uw.swarm.IndexSwarmVariable): - inner_self._materialIndex = indexVar - inner_self._reset() - - @property - def t_relax(inner_self): - # shear modulus defaults to infinity so t_relax goes to zero - # in the viscous limit - - if inner_self.materialIndex == None: - shear_viscosity_0 = inner_self.shear_viscosity_0[0] - shear_modulus = inner_self.shear_modulus[0] - t_relax = shear_viscosity_0 / shear_modulus - else: - shear_viscosity_0 = np.array(inner_self.shear_viscosity_0) - shear_modulus = np.array(inner_self.shear_modulus) - t_relax = shear_viscosity_0 / shear_modulus - - return t_relax - - @property - def ve_effective_viscosity(inner_self): - # the dt_elastic defaults to infinity, t_relax to zero, - # so this should be well behaved in the viscous limit - - import warnings - - if ( - type(inner_self.shear_modulus) != np.ndarray - and type(inner_self.shear_modulus) != list - ): - inner_self.shear_modulus = list([inner_self.shear_modulus]) - if type(inner_self.dt_elastic) != np.ndarray and type(inner_self.dt_elastic) != list: - inner_self.dt_elastic = list([inner_self.dt_elastic]) - if ( - type(inner_self.shear_viscosity_0) != np.ndarray - and type(inner_self.shear_viscosity_0) != list - ): - inner_self.shear_viscosity_0 = list([inner_self.shear_viscosity_0]) - - if inner_self.materialIndex == None: - warnings.warn( - "materialIndex not specified, using the first value for each parameter", - stacklevel=2, - ) - - mu = inner_self.shear_modulus[0] - dt = inner_self.dt_elastic[0] - eta = inner_self.shear_viscosity_0[0] - - if mu == sympy.oo or dt == sympy.oo: - return eta - - ## The effective viscosity depends on the number of history terms - if inner_self.stress_star_star is None: - el_eff_visc = eta * mu * dt / (mu * dt + eta) - else: - el_eff_visc = 2 * eta * mu * dt / (2 * mu * dt + 3 * eta) - - return sympy.simplify(el_eff_visc) - - else: - ### creates list of values that has the same length as the material index - if len(inner_self.shear_modulus) != inner_self.materialIndex.indices: - if len(inner_self.shear_modulus) > 1: - warnings.warn( - f"Number of values in shear_modulus ({len(inner_self.shear_modulus)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_modulus = list( - np.repeat( - inner_self.shear_modulus[0], - inner_self.materialIndex.indices, - ) - ) - - if len(inner_self.dt_elastic) != inner_self.materialIndex.indices: - if len(inner_self.dt_elastic) > 1: - warnings.warn( - f"Number of values in dt_elastic ({len(inner_self.dt_elastic)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.dt_elastic = list( - np.repeat(inner_self.dt_elastic[0], inner_self.materialIndex.indices) - ) - - if len(inner_self.shear_viscosity_0) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_0) > 1: - warnings.warn( - f"Number of values in shear_viscosity_0 ({len(inner_self.shear_viscosity_0)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_0 = list( - np.repeat( - inner_self.shear_viscosity_0[0], - inner_self.materialIndex.indices, - ) - ) - - mu = np.array(inner_self.shear_modulus) - dt = np.array(inner_self.dt_elastic) - eta = np.array(inner_self.shear_viscosity_0) - - if mu[0] == sympy.oo or dt[0] == sympy.oo: - return eta - - ## The effective viscosity depends on the number of history terms - if inner_self.stress_star_star is None: - el_eff_visc = eta * mu * dt / (mu * dt + eta) - else: - el_eff_visc = 2 * eta * mu * dt / (2 * mu * dt + 3 * eta) - - return sympy.simplify(el_eff_visc) - - @property - def plastic_eff_viscosity(inner_self): - import warnings - - if ( - type(inner_self.yield_stress) != np.ndarray - and type(inner_self.yield_stress) != list - ): - inner_self.yield_stress = list([inner_self.yield_stress]) - if ( - type(inner_self.yield_stress_min) != np.ndarray - and type(inner_self.yield_stress_min) != list - ): - inner_self.yield_stress_min = list([inner_self.yield_stress_min]) - - if inner_self.materialIndex == None: - warnings.warn( - "materialIndex not specified, using the first value for each parameter", - stacklevel=2, - ) - - yield_stress = inner_self.yield_stress[0] - yield_stress_min = inner_self.yield_stress_min[0] - - if yield_stress_min == -sympy.oo: - yield_stress_fn = yield_stress - else: - yield_stress_fn = sympy.Max(yield_stress, yield_stress_min) - - if yield_stress_fn == sympy.oo or inner_self.strainrate_inv_II == sympy.oo: - pl_effective_viscosity = sympy.oo - - else: - pl_effective_viscosity = yield_stress_fn / ( - (2 * inner_self.strainrate_inv_II) - # + inner_self.strainrate_inv_II_min - ) - - else: - ### creates list of values that has the same length as the material index - if len(inner_self.yield_stress) != inner_self.materialIndex.indices: - if len(inner_self.yield_stress) > 1: - warnings.warn( - f"Number of values in yield_stress ({len(inner_self.yield_stress)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.yield_stress = list( - np.repeat(inner_self.yield_stress[0], inner_self.materialIndex.indices) - ) - - if len(inner_self.yield_stress_min) != inner_self.materialIndex.indices: - if len(inner_self.yield_stress_min) > 1: - warnings.warn( - f"Number of values in yield_stress ({len(inner_self.yield_stress_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.yield_stress_min = list( - np.repeat( - inner_self.yield_stress_min[0], - inner_self.materialIndex.indices, - ) - ) - - yield_stress = sympy.Matrix(inner_self.yield_stress) - yield_stress_min = sympy.Matrix(inner_self.yield_stress_min) - - if yield_stress_min[0] == -sympy.oo: - yield_stress_fn = yield_stress - else: - yield_stress_list = [] - - for i in range(len(yield_stress)): - yield_stress_list.append(sympy.Max(yield_stress[i], yield_stress_min[i])) - - yield_stress_fn = sympy.Matrix(yield_stress_list) - - if yield_stress_fn[0] == sympy.oo or inner_self.strainrate_inv_II == sympy.oo: - pl_effective_viscosity = list( - np.repeat(sympy.oo, inner_self.materialIndex.indices) - ) - - else: - pl_effective_viscosity = yield_stress_fn / ( - (2 * inner_self.strainrate_inv_II) + inner_self.strainrate_inv_II_min - ) - - return pl_effective_viscosity - - # This has no setter !! - @property - def viscosity(inner_self): - # detect if values we need are defined or are placeholder symbols - - ve_eff_visc = inner_self.ve_effective_viscosity - yield_visc = inner_self.plastic_eff_viscosity - - import warnings - - if inner_self.materialIndex == None: - shear_viscosity_min = inner_self.shear_viscosity_min[0] - shear_viscosity_max = inner_self.shear_viscosity_max[0] - - if inner_self.averaging_method.casefold() == "min": - if yield_visc != 0: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max(shear_viscosity_min, sympy.Min(yield_visc, ve_eff_visc)), - ) - else: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max(shear_viscosity_min, ve_eff_visc), - ) - else: - if yield_visc != 0: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max( - shear_viscosity_min, - 1.0 / ((1.0 / ve_eff_visc) + (1.0 / yield_visc)), - ), - ) - else: - effective_viscosity = sympy.Min( - shear_viscosity_max, - sympy.Max(shear_viscosity_min, ve_eff_visc), - ) - - else: - if len(inner_self.shear_viscosity_min) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_min) > 1: - warnings.warn( - f"Number of values in shear_viscosity_min ({len(inner_self.shear_viscosity_min)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_min = list( - np.repeat( - inner_self.shear_viscosity_min[0], - inner_self.materialIndex.indices, - ) - ) - - if len(inner_self.shear_viscosity_max) != inner_self.materialIndex.indices: - if len(inner_self.shear_viscosity_max) > 1: - warnings.warn( - f"Number of values in shear_viscosity_max ({len(inner_self.shear_viscosity_max)}) does not match the number of material indices ({inner_self.materialIndex.indices}). Using the first value for all materials.", - stacklevel=2, - ) - inner_self.shear_viscosity_max = list( - np.repeat( - inner_self.shear_viscosity_max[0], - inner_self.materialIndex.indices, - ) - ) - - shear_viscosity_min = sympy.Matrix(inner_self.shear_viscosity_min) - shear_viscosity_max = sympy.Matrix(inner_self.shear_viscosity_max) - - viscosity_list = [] - for i in range(len(yield_visc)): - if inner_self.averaging_method.casefold() == "min": - if yield_visc[i] != 0: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max( - shear_viscosity_min[i], - sympy.Min(yield_visc[i], ve_eff_visc[i]), - ), - ) - ) - else: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max(shear_viscosity_min[i], ve_eff_visc[i]), - ) - ) - else: - if yield_visc[i] != 0: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max( - shear_viscosity_min[i], - 1 / ((1 / ve_eff_visc[i]) + (1 / yield_visc[i])), - ), - ) - ) - else: - viscosity_list.append( - sympy.Min( - shear_viscosity_max[i], - sympy.Max(shear_viscosity_min[i], ve_eff_visc[i]), - ) - ) - - viscosity_fn = sympy.Matrix(viscosity_list) - - effective_viscosity = inner_self.materialIndex.sym.T.dot(viscosity_fn) - - return effective_viscosity - - def __init__(self, dim): - u_dim = dim - super().__init__(dim, u_dim) - - return - - ## Is this really different from the original ? - - def _build_c_tensor(self): - """For this constitutive law, we expect just a viscosity function""" - - d = self.dim - viscosity = self.Parameters.viscosity - shear_modulus = self.Parameters.shear_modulus - dt_elastic = self.Parameters.dt_elastic - - try: - self._c = 2 * uw.maths.tensor.rank4_identity(d) * viscosity - except: - d = self.dim - dv = uw.maths.tensor.idxmap[d][0] - if isinstance(viscosity, sympy.Matrix) and viscosity.shape == (dv, dv): - self._c = 2 * uw.maths.tensor.mandel_to_rank4(viscosity, d) - elif isinstance(viscosity, sympy.Array) and viscosity.shape == (d, d, d, d): - self._c = 2 * viscosity - else: - raise RuntimeError( - "Viscosity is not a known type (scalar, Mandel matrix, or rank 4 tensor" - ) - return - - # Modify flux to use the stress history term - # This may be preferable to using strain rate which can be discontinuous - # and harder to map back and forth between grid and particles without numerical smoothing - - def flux( - self, - ddu: sympy.Matrix = None, - ddu_dt: sympy.Matrix = None, - u: sympy.Matrix = None, # may be needed in the case of cylindrical / spherical - u_dt: sympy.Matrix = None, - ): - """Computes the effect of the constitutive tensor on the gradients of the unknowns. - (always uses the `c` form of the tensor). In general cases, the history of the gradients - may be required to evaluate the flux. For viscoelasticity, the - """ - - c = self.c - rank = len(c.shape) - - # tensor multiplication - - if rank == 2: - flux = c * ddu.T - else: # rank==4 - flux = sympy.tensorcontraction( - sympy.tensorcontraction(sympy.tensorproduct(c, ddu), (3, 5)), (0, 1) - ) - - # Now add in the stress history. In the - # viscous limit, this term is not well behaved - # and we need to check that - - if self.is_elastic: - if self.Parameters.materialIndex == None: - eta = self.Parameters.shear_viscosity_0[0] - mu = self.Parameters.shear_modulus[0] - dt = self.Parameters.dt_elastic[0] - s_star = self.Parameters.stress_star - s_star_star = self.Parameters.stress_star_star - - if s_star_star is None: # 1st order - flux = sympy.Matrix(flux) + eta * s_star / (dt * mu + eta) - else: # 2nd order - flux = ( - sympy.Matrix(flux) - + 4 * eta * s_star / (2 * dt * mu + 3 * eta) - - eta * s_star_star / (2 * dt * mu + 3 * eta) - ) - else: - eta = self.Parameters.materialIndex.createMask(self.Parameters.shear_viscosity_0) - mu = self.Parameters.materialIndex.createMask(self.Parameters.shear_modulus) - dt = np.array(self.Parameters.dt_elastic) - s_star = self.Parameters.stress_star - s_star_star = self.Parameters.stress_star_star - - if s_star_star is None: # 1st order - flux = sympy.Matrix(flux) + eta * s_star / (dt * mu + eta) - else: # 2nd order - flux = ( - sympy.Matrix(flux) - + 4 * eta * s_star / (2 * dt * mu + 3 * eta) - - eta * s_star_star / (2 * dt * mu + 3 * eta) - ) - - return sympy.simplify(sympy.Matrix(flux)) - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - - super()._object_viewer() - - display(Markdown(r"### Viscous deformation")) - display( - Latex( - r"$\quad\eta_\textrm{0} = $ " - + sympy.sympify(self.Parameters.shear_viscosity_0[0])._repr_latex_() - ), - ) - - ## If elasticity is active: - display(Markdown(r"#### Elastic deformation")) - display( - Latex( - r"$\quad\mu = $ " + sympy.sympify(self.Parameters.shear_modulus[0])._repr_latex_(), - ), - Latex( - r"$\quad\Delta t_e = $ " - + sympy.sympify(self.Parameters.dt_elastic[0])._repr_latex_(), - ), - Latex( - r"$\quad \sigma^* = $ " + sympy.sympify(self.Parameters.stress_star)._repr_latex_(), - ), - ) - if self.Parameters.stress_star_star is not None: - display( - Latex( - r"$\quad \sigma^{**} = $ " - + sympy.sympify(self.Parameters.stress_star_star)._repr_latex_(), - ), - ) - - # If plasticity is active - display(Markdown(r"#### Plastic deformation")) - display( - Latex( - r"$\quad\tau_\textrm{y} = $ " - + sympy.sympify(self.Parameters.yield_stress[0])._repr_latex_(), - ), - Latex( - r"$\quad|\dot\epsilon| = $ " - + sympy.sympify(self.Parameters.strainrate_inv_II)._repr_latex_(), - ), - ## Todo: add all the other properties in here - ) - - @property - def is_elastic(self): - # If any of these is not defined, elasticity is switched off - - if self.Parameters.dt_elastic[0] is sympy.oo: - return False - - if self.Parameters.shear_modulus[0] is sympy.oo: - return False - - if self.Parameters.stress_star is None: - return False - - return True - - @property - def is_viscoplastic(self): - if self.Parameters.yield_stress[0] == sympy.oo: - return False - - if isinstance(self.Parameters.strainrate_inv_II, sympy.core.symbol.Symbol): - return False - - return True - - -### - - -class DiffusionModel(Constitutive_Model): - r""" - Diffusion constitutive model for scalar transport. - - .. math:: - - q_{i} = \kappa_{ij} \cdot \frac{\partial \phi}{\partial x_j} - - where :math:`\kappa` is a diffusivity—a scalar constant, sympy function, - or mesh variable. - - Examples - -------- - >>> diffusion_model = DiffusionModel(dim) - >>> diffusion_model.material_properties = diffusion_model.Parameters( - ... diffusivity=diffusivity_fn - ... ) - >>> scalar_solver.constitutive_model = diffusion_model - >>> flux = diffusion_model.flux(gradient_matrix) - """ - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - """ - - def __init__( - inner_self, - diffusivity: Union[float, sympy.Function] = 1, - ): - inner_self._diffusivity = diffusivity - - @property - def diffusivity(inner_self): - return inner_self._diffusivity - - @diffusivity.setter - def diffusivity(inner_self, value: Union[float, sympy.Function]): - inner_self._diffusivity = value - inner_self._reset() - - def __init__(self, dim): - self.u_dim = 1 - super().__init__(dim, self.u_dim) - - return - - def _build_c_tensor(self): - """For this constitutive law, we expect just a diffusivity function""" - - d = self.dim - kappa = self.Parameters.diffusivity - self._c = sympy.Matrix.eye(d) * kappa - - return - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - - super()._object_viewer() - - ## feedback on this instance - display( - Latex(r"$\quad\kappa = $ " + sympy.sympify(self.Parameters.diffusivity)._repr_latex_()) - ) - - return - - -class TransverseIsotropicFlowModel(Constitutive_Model): - r""" - Transversely isotropic viscous flow model. - - .. math:: - - \tau_{ij} = \eta_{ijkl} \cdot \frac{1}{2} \left[ \frac{\partial u_k}{\partial x_l} - + \frac{\partial u_l}{\partial x_k} \right] - - where the viscosity tensor :math:`\eta` is defined as: - - .. math:: - - \eta_{ijkl} = \eta_0 I_{ijkl} + (\eta_0 - \eta_1) \left[ - \frac{1}{2}\left( n_i n_l \delta_{jk} + n_j n_k \delta_{il} - + n_i n_l \delta_{jk} + n_j n_l \delta_{ik} \right) - - 2 n_i n_j n_k n_l \right] - - and :math:`\hat{\mathbf{n}} \equiv \{n_i\}` is the unit vector defining the - local orientation of the weak plane (the director). - - The Mandel matrix is in ``viscous_model.C`` and the rank-4 tensor in - ``viscous_model.c``. - - Examples - -------- - >>> viscous_model = TransverseIsotropicFlowModel(dim) - >>> viscous_model.material_properties = viscous_model.Parameters( - ... eta_0=viscosity_fn, - ... eta_1=weak_viscosity_fn, - ... director=orientation_vector_fn - ... ) - >>> solver.constitutive_model = viscous_model - >>> tau = viscous_model.flux(gradient_matrix) - """ - - class _Parameters: - """Any material properties that are defined by a constitutive relationship are - collected in the parameters which can then be defined/accessed by name in - individual instances of the class. - """ - - def __init__( - inner_self, - eta_0: Union[float, sympy.Function] = 1, - eta_1: Union[float, sympy.Function] = 1, - director: Union[sympy.Matrix, sympy.Function] = sympy.Matrix([0, 0, 1]), - ): - inner_self._eta_0 = eta_0 - inner_self._eta_1 = eta_1 - inner_self._director = director - # inner_self.constitutive_model_class = const_model - - ## Note the inefficiency below if we change all these values one after the other - - @property - def eta_0(inner_self): - return inner_self._eta_0 - - @eta_0.setter - def eta_0( - inner_self, - value: Union[float, sympy.Function], - ): - inner_self._eta_0 = value - inner_self._reset() - - @property - def eta_1(inner_self): - return inner_self._eta_1 - - @eta_1.setter - def eta_1( - inner_self, - value: Union[float, sympy.Function], - ): - inner_self._eta_1 = value - inner_self._reset() - - @property - def director(inner_self): - return inner_self._director - - @director.setter - def director( - inner_self, - value: Union[sympy.Matrix, sympy.Function], - ): - inner_self._director = value - inner_self._reset() - - def __init__(self, dim): - u_dim = dim - super().__init__(dim, u_dim) - - # default values ... maybe ?? - return - - def _build_c_tensor(self): - """For this constitutive law, we expect two viscosity functions - and a sympy matrix that describes the director components n_{i}""" - - d = self.dim - dv = uw.maths.tensor.idxmap[d][0] - - eta_0 = self.Parameters.eta_0 - eta_1 = self.Parameters.eta_1 - n = self.Parameters.director - - Delta = eta_1 - eta_0 - - lambda_mat = uw.maths.tensor.rank4_identity(d) * eta_0 - - for i in range(d): - for j in range(d): - for k in range(d): - for l in range(d): - lambda_mat[i, j, k, l] += Delta * ( - ( - n[i] * n[k] * int(j == l) - + n[j] * n[k] * int(l == i) - + n[i] * n[l] * int(j == k) - + n[j] * n[l] * int(k == i) - ) - / 2 - - 2 * n[i] * n[j] * n[k] * n[l] - ) - - lambda_mat = sympy.simplify(uw.maths.tensor.rank4_to_mandel(lambda_mat, d)) - - self._c = uw.maths.tensor.mandel_to_rank4(lambda_mat, d) - - def _object_viewer(self): - from IPython.display import Latex, Markdown, display - - super()._object_viewer() - - ## feedback on this instance - display(Latex(r"$\quad\eta_0 = $ " + sympy.sympify(self.Parameters.eta_0)._repr_latex_())) - display(Latex(r"$\quad\eta_1 = $ " + sympy.sympify(self.Parameters.eta_1)._repr_latex_())) - display( - Latex( - r"$\quad\hat{\mathbf{n}} = $ " - + sympy.sympify(self.Parameters.director.T)._repr_latex_() - ) - ) - - -class MultiMaterial(Constitutive_Model): - r""" - Manage multiple materials in a constitutive framework. - - Bundles multiple materials into a single consitutive law. The expectation - is that these all have compatible flux terms. - """ - - def __init__( - self, - material_swarmVariable: Optional[IndexSwarmVariable] = None, - constitutive_models: Optional[list] = [], - ): - self._constitutive_models = constitutive_models - self._material_var = material_swarmVariable - - return From 7199d4a99b80633b2724358195ae68e868d9f713 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 25 Mar 2026 07:20:37 +1100 Subject: [PATCH 18/18] Move square-wave benchmark to docs, remove parked advect_history - Move run_ve_square_wave.py to docs/advanced/benchmarks/ - Add ve-square-wave-shear.md benchmark documentation (10 harmonics) - Add to benchmarks toctree index - Remove _advect_history_PARKED (replaced by update_pre_solve) Underworld development team with AI support from Claude Code --- docs/advanced/benchmarks/index.md | 1 + .../benchmarks}/run_ve_square_wave.py | 6 +- .../benchmarks/ve-square-wave-shear.md | 79 ++++++++++++ src/underworld3/systems/ddt.py | 117 ------------------ 4 files changed, 83 insertions(+), 120 deletions(-) rename {tests => docs/advanced/benchmarks}/run_ve_square_wave.py (98%) create mode 100644 docs/advanced/benchmarks/ve-square-wave-shear.md diff --git a/docs/advanced/benchmarks/index.md b/docs/advanced/benchmarks/index.md index 774102ea..56ee9ae0 100644 --- a/docs/advanced/benchmarks/index.md +++ b/docs/advanced/benchmarks/index.md @@ -10,4 +10,5 @@ Validation benchmarks comparing Underworld3 solvers against analytical solutions :maxdepth: 1 ve-oscillatory-shear +ve-square-wave-shear ``` diff --git a/tests/run_ve_square_wave.py b/docs/advanced/benchmarks/run_ve_square_wave.py similarity index 98% rename from tests/run_ve_square_wave.py rename to docs/advanced/benchmarks/run_ve_square_wave.py index d56ffc73..56a988b8 100644 --- a/tests/run_ve_square_wave.py +++ b/docs/advanced/benchmarks/run_ve_square_wave.py @@ -183,7 +183,7 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, De = 1.5 order = 2 n_periods = 3 - n_harmonics = 3 + n_harmonics = 10 print("=" * 60) print(f"Square-wave VE benchmark: De={De}, order={order}") @@ -231,7 +231,7 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, # Save results np.savez( - "tests/ve_square_wave_benchmark.npz", + f"tests/ve_square_wave_{n_harmonics}h.npz", adaptive_times=result_adaptive["times"], adaptive_numerical=result_adaptive["numerical"], adaptive_analytical=result_adaptive["analytical"], @@ -241,4 +241,4 @@ def run_square_wave(order, De, n_periods, dt_min_over_tr, dt_max_over_tr, uniform_analytical=result_uniform["analytical"], uniform_timesteps=result_uniform["timesteps"], ) - print("\nResults saved to tests/ve_square_wave_benchmark.npz") + print(f"\nResults saved to tests/ve_square_wave_{n_harmonics}h.npz") diff --git a/docs/advanced/benchmarks/ve-square-wave-shear.md b/docs/advanced/benchmarks/ve-square-wave-shear.md new file mode 100644 index 00000000..6d80998a --- /dev/null +++ b/docs/advanced/benchmarks/ve-square-wave-shear.md @@ -0,0 +1,79 @@ +--- +title: "Viscoelastic Square-Wave Shear Benchmark" +--- + +# Maxwell Square-Wave Shear + +This benchmark validates the viscoelastic Stokes solver with **variable timesteps** +against an analytical solution for a Maxwell material under square-wave shear forcing. + +It tests both the variable-dt BDF-2 coefficients and the PetscDS constants mechanism +that routes these coefficients to the compiled pointwise functions at runtime. + +## Problem Setup + +Same geometry as the oscillatory shear benchmark: a box with height $H$ and width $2H$, +sheared by top/bottom boundary velocities. The shear rate is a truncated Fourier series +approximation of a square wave: + +$$\dot\gamma(t) = \frac{4\dot\gamma_0}{\pi} +\sum_{k=1}^{N} \frac{\sin\bigl((2k-1)\omega t\bigr)}{2k-1}$$ + +The sharp transitions between positive and negative shear demand small timesteps +near the transition points, while the plateaux can use much larger steps. This +makes it a natural test for adaptive (variable) timestepping. + +## Analytical Solution + +Since the Maxwell equation is linear, the stress is the superposition of +single-frequency Maxwell solutions at each Fourier harmonic: + +$$\sigma_{xy}(t) = \sum_{k=1}^{N} \sigma_k(t)$$ + +where each $\sigma_k$ is the oscillatory Maxwell solution with amplitude +$a_k = 4\dot\gamma_0 / (\pi(2k-1))$ and frequency $\omega_k = (2k-1)\omega$: + +$$\sigma_k(t) = \frac{\eta\, a_k}{1 + \text{De}_k^2} +\left[\sin(\omega_k t) - \text{De}_k\cos(\omega_k t) ++ \text{De}_k\,e^{-t/t_r}\right]$$ + +with $\text{De}_k = \omega_k t_r$. + +## Adaptive Timestep Strategy + +The timestep varies between `dt_min` near transitions and `dt_max` on plateaux, +based on distance to the nearest square-wave transition point: + +$$\Delta t = \Delta t_{\min} + (\Delta t_{\max} - \Delta t_{\min})\, f^2$$ + +where $f \in [0,1]$ measures the normalised distance from the nearest transition. + +## Variable-dt BDF Coefficients + +With uniform timesteps, BDF-2 uses constant coefficients $[3/2, -2, 1/2]$. +With variable timesteps (ratio $r = \Delta t_n / \Delta t_{n-1}$), the +coefficients become: + +$$c_0 = \frac{1+2r}{1+r}, \quad c_1 = -(1+r), \quad c_2 = \frac{r^2}{1+r}$$ + +These coefficients are stored as UWexpressions and updated each step via +`_update_bdf_coefficients()`, flowing through PetscDS `constants[]` to the +compiled pointwise functions without JIT recompilation. + +## Results (De=1.5, BDF-2, 10 harmonics) + +| Run | Steps | L2 Error | Ratio | +|-----|-------|----------|-------| +| Adaptive dt | 295 | 9.55e-04 | 1.78x | +| Uniform dt | 629 | 5.35e-04 | 1.0x | + +The adaptive run uses 53% fewer steps at only 1.78x the error. + +## Running the Benchmark + +```bash +pixi run -e default python docs/advanced/benchmarks/run_ve_square_wave.py +``` + +The script runs both adaptive and uniform timestep cases, prints convergence +data, and saves results to `.npz` files. diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 400c33b8..6e124d58 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1569,123 +1569,6 @@ 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."""