From f45757355a520bd21147939a89431939adf2f733 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 fd37fe077..c4d52af79 100755 --- a/uw +++ b/uw @@ -162,8 +162,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 1415fe503fe0448dd303176822bab1c7136a1397 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 | 84 +++++++++++++-- 4 files changed, 234 insertions(+), 23 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 8d02acde9..205f6ec74 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 d9a0b825a..98f218e3b 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 25e2e3e0f..9e20739b8 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -445,7 +445,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]``.""" @@ -763,7 +765,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.""" @@ -1158,7 +1162,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`. @@ -1229,6 +1235,7 @@ def update_post_solve( if self._n_solves_completed < self.order: self._n_solves_completed += 1 + return def update_pre_solve( @@ -1595,6 +1602,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"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -1773,7 +1895,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 550632306..8d86f10d7 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,24 @@ def solve( if timestep is None: timestep = self.delta_t.sym - if timestep != self.delta_t: - self._constitutive_model.Parameters.dt_elastic = timestep # this will force an initialisation because the functions need to be updated + # dt_elastic is a constitutive parameter (relaxation timescale) — + # never overwritten by solve(). The advection timestep (for departure + # point tracing) and the elastic relaxation timescale are independent. if _force_setup: self.is_setup = False + # Re-setup when effective_order changes (e.g. DDt history ramp-up + # from order 1 to order 2). The JIT-compiled pointwise functions + # depend on the order used in the constitutive model. + _current_eff_order = self.constitutive_model.effective_order + if not hasattr(self, '_prev_effective_order'): + self._prev_effective_order = None + if _current_eff_order != self._prev_effective_order: + self.is_setup = False + self.constitutive_model._solver_is_setup = False + self._prev_effective_order = _current_eff_order + if not self.constitutive_model._solver_is_setup: self.is_setup = False self.DFDt.psi_fn = self.constitutive_model.flux.T @@ -1289,11 +1315,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 +1340,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 ceb5216619d499f7e6276639bc0e098681a0490c 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 000000000..286b73595 --- /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 000000000..9e7d6d8c8 --- /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 000000000..6be238580 --- /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 000000000..759616edb --- /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 000000000..a8d2709ad --- /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 3c685844d3ab5b20952305e3a1d643e048b59862 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 e830fdc42..92379c885 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -956,6 +956,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 8d86f10d7..757fede1d 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 9e7d6d8c8..fb40d71f9 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 6be238580..debced41d 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 a8d2709ad..1691ee319 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 329826f11718c1f826965c2a99d8262c3e70ea8f 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 205f6ec74..6ff31729a 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 000000000..0c70ba0ce --- /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 000000000..7d4fabe3f --- /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 000000000..01a0c1651 --- /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 f5970b43c33831f582254ae2f40eb7bede21114f 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 d502dd929..f2e5cc103 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 c39d9c715..edd73ca15 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 92379c885..c4b3377e0 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1784,7 +1784,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. @@ -1807,6 +1808,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 ------- @@ -1844,6 +1850,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: @@ -4153,7 +4169,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. @@ -4180,6 +4197,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 ------- @@ -4221,6 +4243,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 c111d0fdd..1db908564 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 304aadebf..26471db63 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 000000000..abc234ca6 --- /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 000000000..706cc2ce8 --- /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 000000000..b5b50c0db --- /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 000000000..3f767579e --- /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 7cd6a43b8040025d3a08efa3d8b4a242aa9a1d7c 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 | 44 ++++++++--- 2 files changed, 76 insertions(+), 68 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 6ff31729a..df7c26b3d 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 757fede1d..12c60a90d 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1360,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 5cd5580e8dfa6fdf2043e57f7040e4b639699b0d Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 22 Mar 2026 14:12:55 +1100 Subject: [PATCH 08/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 000000000..774102ea4 --- /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 000000000..a7663bcde --- /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 12c60a90d..e6f51dd87 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 09/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 a7663bcde..03e95c437 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 0cab29357e9614a366aa1290dd73cf1915a3f447 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Mon, 23 Mar 2026 13:23:45 +1100 Subject: [PATCH 10/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 1db908564..4b11dbd6e 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 26471db63..0073f2830 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 d4dbf5255..f3ac103f7 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 be5c996be..98242026e 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 706cc2ce8..b14529d51 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 b5b50c0db..79ba8dcf5 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 3f767579e..9831b3e80 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 1691ee319..1608a7291 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 ae8066b2074ac78a1331cebca77b6c3c27028af9 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 09:14:51 +1100 Subject: [PATCH 11/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 df7c26b3d..b50687fe8 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 9e20739b8..27705b914 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -526,6 +526,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"""Backward differentiation approximation of the time-derivative of ψ. @@ -894,6 +899,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"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -1717,6 +1727,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"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -2010,6 +2025,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"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -2307,6 +2327,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"""Backward differentiation approximation of the time-derivative of :math:`\psi`. diff --git a/tests/run_ve_square_wave.py b/tests/run_ve_square_wave.py new file mode 100644 index 000000000..788b7551e --- /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 8bc3b6662f13868a0bccea4a00ad6665de6031cf Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 14:00:12 +1100 Subject: [PATCH 12/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 98242026e..2cf7e21da 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 788b7551e..36b2a8336 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 b14529d51..8ae07f015 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 79ba8dcf5..b5ee4907c 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 9831b3e80..e76e212e5 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 f2204537d280b75c160bac2cc457a9c43e7c89dc Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 15:22:59 +1100 Subject: [PATCH 13/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 27705b914..7ef2f4d94 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1254,11 +1254,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. """ self._dt = dt @@ -1292,15 +1302,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: @@ -1339,28 +1353,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 @@ -1612,7 +1627,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. @@ -1695,6 +1710,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, :] @@ -1709,6 +1725,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 e6f51dd87..e536fbc43 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 36b2a8336..d56ffc73e 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 a43340cf2d88687ee32ef46ece9a9c219ad5ef48 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 15:36:32 +1100 Subject: [PATCH 14/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 7ef2f4d94..24c77efe9 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1254,7 +1254,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. @@ -1263,12 +1263,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. """ self._dt = dt @@ -1302,7 +1303,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[...] @@ -1312,7 +1313,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. @@ -1353,7 +1354,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 e536fbc43..7772d0821 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 8cd7293725d206ca089fbca664818a850ff0bb70 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 19:17:46 +1100 Subject: [PATCH 15/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 b50687fe8..47293a680 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 7772d0821..62387707a 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 995ec7672c2ee2c56b966f07081ef680e53b8ccd Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 24 Mar 2026 20:32:58 +1100 Subject: [PATCH 16/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 d58be0598..000000000 --- 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 1d575caf1af8e8a56b3e392d690b49f90892aa72 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 25 Mar 2026 07:20:37 +1100 Subject: [PATCH 17/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 774102ea4..56ee9ae04 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 d56ffc73e..56a988b82 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 000000000..6d80998a4 --- /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 24c77efe9..64c10b255 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1628,123 +1628,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.""" From 069cb722e70687d36164fee82b757d820be47d2f Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 25 Mar 2026 16:02:51 +1100 Subject: [PATCH 18/18] Update uw script from development (PIXI_PROJECT_MANIFEST fix) Underworld development team with AI support from Claude Code --- uw | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/uw b/uw index c4d52af79..fd37fe077 100755 --- a/uw +++ b/uw @@ -162,12 +162,8 @@ 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"