Route BDF/AM coefficients through PetscDS constants array#88
Route BDF/AM coefficients through PetscDS constants array#88lmoresi merged 3 commits intodevelopmentfrom
Conversation
BDF and Adams-Moulton coefficients in all 5 DDt classes are now UWexpression objects that the JIT routes through PetscDSSetConstants(). This fixes a bug where the NS solver's F0/F1 expressions were frozen at order-1 coefficients from the first solve — effective_order ramping had no effect on the compiled code. Changes: - Add module-level helpers: _create_coefficients(), _update_bdf_values(), _update_am_values(), _build_weighted_sum() - All DDt classes create coefficient UWexpressions in __init__ - bdf() and adams_moulton_flux() return fixed-structure expressions - Coefficient values update each step in update_pre_solve() - NS solver F0: bdf(1) -> bdf() (was hardcoded to order 1) The JIT compiles the full expression once (all terms present, inactive ones have coefficient=0). When effective_order ramps or dt changes, only PetscDSSetConstants() is called — no recompilation. Underworld development team with AI support from Claude Code
Clamp mid_pt_coords and end_pt_coords to domain bounds in SemiLagrangian.update_pre_solve() to prevent extrapolation blowup near boundary discontinuities. Underworld development team with AI support from Claude Code
Allows independent control of Adams-Moulton flux order via flux_order= parameter. When set, AM coefficients are capped at the specified order while BDF follows effective_order normally. Default (flux_order=None) keeps BDF and AM matched. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
This PR fixes a Navier–Stokes time-integration bug where BDF/Adams–Moulton coefficients were effectively frozen at order 1 after the first JIT compile, by routing these coefficients through PETSc constants[] so they can be updated per-solve without recompilation.
Changes:
- Refactors DDt implementations to build fixed-structure BDF/AM symbolic expressions with coefficient values carried via
UWexpressionconstants updated each timestep. - Updates Navier–Stokes momentum source to use
DuDt.bdf()(effective order) and introduces an optionalflux_orderoverride for AM flux integration. - Adds SemiLagrangian departure-point clamping to keep upstream sampling inside the domain.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
src/underworld3/systems/ddt.py |
Introduces constant-routed BDF/AM coefficient expressions and updates all DDt classes to use fixed-structure sums; adds SemiLagrangian coordinate clamping. |
src/underworld3/systems/solvers.py |
Switches NS inertia term to bdf() and adds flux_order logic to override AM coefficients during solve. |
Comments suppressed due to low confidence (1)
src/underworld3/systems/solvers.py:2707
flux_orderis added to the constructor signature, but the class docstring’s Parameters section doesn’t mention it and there’s no input validation. As written, negative values or values >3 will silently map to order0/3 behavior in_update_am_values. Please document the new parameter and validate/clamp it to the supported Adams–Moulton orders (0–3) (or raise a clear error).
rho: Optional[float] = 0.0,
restore_points_func: Callable = None,
order: Optional[int] = 2,
flux_order: Optional[int] = None,
p_continuous: Optional[bool] = False,
verbose: Optional[bool] = False,
DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None,
DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None,
):
## Parent class will set up default values and load u_Field into the solver
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| from underworld3.systems.ddt import _update_am_values | ||
| fo = min(self._flux_order, self.DFDt.effective_order) | ||
| _update_am_values(self.DFDt._am_coeffs, fo, 0.5) |
There was a problem hiding this comment.
The flux_order override reaches into DFDt._am_coeffs and imports the private helper _update_am_values, which tightly couples the solver to ddt.py internals. It also hard-codes theta=0.5, ignoring the solver’s configured theta (passed when constructing DFDt). Consider adding a public method on the DDt object (e.g. set_flux_order(order) or update_am_coeffs(order, theta=self.theta)) and use self.DFDt.theta (or equivalent) instead of 0.5.
| from underworld3.systems.ddt import _update_am_values | |
| fo = min(self._flux_order, self.DFDt.effective_order) | |
| _update_am_values(self.DFDt._am_coeffs, fo, 0.5) | |
| fo = min(self._flux_order, self.DFDt.effective_order) | |
| theta = getattr(self.DFDt, "theta", 0.5) | |
| update_am = getattr(self.DFDt, "update_am_coeffs", None) | |
| if callable(update_am): | |
| # Prefer public API on DFDt when available | |
| update_am(order=fo, theta=theta) | |
| else: | |
| # Fallback to legacy private helper and internal coefficients | |
| from underworld3.systems.ddt import _update_am_values | |
| _update_am_values(self.DFDt._am_coeffs, fo, theta) |
|
|
||
| # This is actually calculated over several steps so scaling is required | ||
| return bdf0 | ||
| return _build_weighted_sum(self._bdf_coeffs, self.psi_fn, [ps.sym for ps in self.psi_star]) |
There was a problem hiding this comment.
Lagrangian_Swarm.bdf() previously scaled the BDF expression by self.step_averaging (see the class docstring and update_post_solve’s use of phi = 1/self.step_averaging). That scaling is now missing, so the returned Δψ will be too large when step_averaging != 1. Reintroduce the 1/self.step_averaging scaling (either inside _build_weighted_sum for this class or by dividing the returned expression here).
| return _build_weighted_sum(self._bdf_coeffs, self.psi_fn, [ps.sym for ps in self.psi_star]) | |
| return _build_weighted_sum( | |
| self._bdf_coeffs, | |
| self.psi_fn, | |
| [ps.sym for ps in self.psi_star], | |
| ) / self.step_averaging |
| def adams_moulton_flux(self, order: Optional[int] = None): | ||
| r"""Adams-Moulton flux approximation for implicit time integration. | ||
|
|
||
| When ``order`` is not specified, uses ``effective_order`` which | ||
| ramps up from 1 during startup. | ||
| Returns a fixed-structure symbolic expression using UWexpression | ||
| coefficients. Values are updated each step in ``update_pre_solve``. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| order : int, optional | ||
| Order of the approximation (1-3). Defaults to ``effective_order``. | ||
|
|
||
| Returns | ||
| ------- | ||
| sympy.Matrix | ||
| Weighted average of :math:`\psi` and history terms. | ||
|
|
||
| Notes | ||
| ----- | ||
| The Adams-Moulton formulas for order 1-3 are: | ||
|
|
||
| - Order 1: :math:`\theta \psi + (1-\theta) \psi^*` | ||
| - Order 2: :math:`\frac{5\psi + 8\psi^* - \psi^{**}}{12}` | ||
| - Order 3: :math:`\frac{9\psi + 19\psi^* - 5\psi^{**} + \psi^{***}}{24}` | ||
| Ignored (kept for API compatibility). | ||
| """ | ||
| if order is None: | ||
| order = self.effective_order | ||
| else: | ||
| order = max(1, min(self.effective_order, order)) | ||
|
|
||
| with sympy.core.evaluate(False): | ||
| if order == 1: | ||
| am = self.theta * self.psi_fn + (1.0 - self.theta) * self.psi_star[0] | ||
| elif order == 2: | ||
| am = (5 * self.psi_fn + 8 * self.psi_star[0] - self.psi_star[1]) / 12 | ||
| elif order == 3: | ||
| am = ( | ||
| 9 * self.psi_fn | ||
| + 19 * self.psi_star[0] | ||
| - 5 * self.psi_star[1] | ||
| + self.psi_star[2] | ||
| ) / 24 | ||
| return am | ||
| return _build_weighted_sum(self._am_coeffs, self.psi_fn, self.psi_star) |
There was a problem hiding this comment.
The updated adams_moulton_flux() docstrings no longer mention that the implemented Adams–Moulton coefficients are constant-Δt formulas for order>1 (and _update_am_values() still does not account for variable timesteps). Please add an explicit note/warning here (or reference _update_am_values) so users don’t assume variable-Δt AM2/AM3 is handled just because coefficients are now updated via constants[].
|
Bug fix - see the response to PR #85. 2nd order solver now works. |
Summary
Fixes a bug where the Navier-Stokes solver's BDF and Adams-Moulton coefficients were frozen at order 1 from the first solve. The
effective_orderramp (1 to 2 over successive solves) had no effect on the compiled pointwise functions because the coefficients were baked as numerical literals during JIT compilation.Root cause: The NS solver's
F0andF1properties callDuDt.bdf()andDFDt.adams_moulton_flux()during_setup_pointwise_functions, which runs once (guarded byself.is_setup). The resulting SymPy expressions are JIT-compiled to C code with the coefficient values as literals. Wheneffective_orderramps from 1 to 2 on subsequent solves, the compiled code is never rebuilt.Fix: Make BDF/AM coefficients
UWexpressionobjects that the JIT routes through PETSc'sPetscDSSetConstants()mechanism (from PR #87). The symbolic expression structure is built once at max order (all terms present, inactive ones have coefficient = 0). Wheneffective_orderramps or the timestep changes, onlyPetscDSSetConstants()is called -- no recompilation.Changes
ddt.py: Add module-level helpers (_create_coefficients,_update_bdf_values,_update_am_values,_build_weighted_sum). All 5 DDt classes (Symbolic, Eulerian, SemiLagrangian, Lagrangian, Lagrangian_Swarm) create coefficientUWexpressionobjects in__init__, update them inupdate_pre_solve, and return fixed-structure expressions frombdf()/adams_moulton_flux().solvers.py: NS solverF0usesbdf()instead ofbdf(1). Newflux_orderparameter allows independent AM order control.Verification
Before fix: All order-2 BDF/AM combinations (BDF1+AM2, BDF2+AM2, BDF2+AM1) produced identical results to 6 decimal places (RMS = 0.041469 vs Ghia) -- the coefficients were frozen at order 1 with extra history noise.
After fix: Lid-driven cavity Re=100 (cellSize=0.04, 200 steps, dt=0.05):
All solutions are clean with no instabilities.
Test plan
flux_orderparameter works (BDF2+AM1 differs from BDF2+AM2)Underworld development team with AI support from Claude Code