Taylor-Green order-reduction benchmark in FEniCS, and a remedy (supersedes #650) - #684
Open
pancetta wants to merge 7 commits into
Open
Taylor-Green order-reduction benchmark in FEniCS, and a remedy (supersedes #650)#684pancetta wants to merge 7 commits into
pancetta wants to merge 7 commits into
Conversation
Reimplements the benchmark from PR #650 with legacy FEniCS in the existing StroemungsRaum project instead of FEniCSx in a new project. Periodicity comes from a constrained_domain on the function space, so dolfinx_mpc is not needed, and the project is already covered by the CI matrix and the `fenics` marker. The manufactured solution and forcing term of #650 were verified symbolically and carried over unchanged; the implementation issues raised in review were not: - No mass matrix inversion. The right-hand side is subtracted from the Newton residual as a vector, so the Mf^-1 / Mf round trip in #650 (and in the merged NavierStokes_2D_monolithic_FEniCS.py) disappears together with its per-node solve and its accuracy floor. - The residual is only fixed on the actual Dirichlet boundary. In #650 the periodic variant zeroed the residual on the periodic boundaries too, which let it stop iterating earlier than the Dirichlet variant and biased exactly the comparison the benchmark exists for. - The form, its Jacobian, the Newton problem and the solver are built once, with the step size as a Constant, rather than rebuilt per solve_system call. - Scaling by `factor` instead of dividing by it, so there is no factor == 0 case and the residual scale, and with it snes/newton tolerances, no longer depends on dt. - Domain extents are not derived from rank-local coordinates. Adds the order study that was missing: with periodic conditions the method attains its design order 2M-1 = 5, with time-dependent Dirichlet conditions it does not. At CI-affordable resolution the order gap is modest (~4.7 vs ~4.3 in the pressure, the Radau IIA reduction from 2M-1 to M+1); the error constant gap is close to an order of magnitude, and that is what the test asserts. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
solve_system applied Mf^-1 to the right-hand side handed over by the sweeper, only for the residual form to apply Mf straight back to it via the dot(rhs, v)*dx terms. The round trip is exact, so it bought nothing while costing an LU solve per node per sweep and capping attainable accuracy at the tolerance of that solve. It was there because dolfin's solve(F == 0, ...) accepts a form and not a vector. Subtract the right-hand side vector in a NonlinearProblem instead. The same pattern is needed by the Taylor-Green benchmark, so it lives in newton_step.py and both problem classes use it. Building the Newton problem this way also means the form, its Jacobian and the solver are set up once, with the step size carried by a Constant, rather than rebuilt on every solve_system call. Previously each distinct value of `factor` was baked into the form as a literal, giving a different form signature per step size (verified) and so a separate compilation. FFC caches those across runs, so this shows up on a cold cache or when new step sizes keep appearing, not in steady-state runtime: on a warm cache the drag/lift benchmark test is unchanged at 2.2s against 2.1s before. Numerically neutral: drag and lift in that benchmark agree with the previous implementation to 13 significant digits (relative difference 9e-15 and 2e-12), which is the round-off of the LU solve that has been removed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
RADAU-RIGHT with M nodes has design order 2M-1 and falls back to the stiff
order M+1 when the boundary data is time dependent, so the gap this benchmark
can show is M-2. That is *identically zero for M = 2*, which is what the
original FEniCSx version ran with: that configuration cannot exhibit the
phenomenon it was written to demonstrate, whatever else is fixed. Measured at
nelems=24, nu=0.1, Tend=0.2, orders from consecutive step sizes:
M design stiff periodic order(p) Dirichlet order(p)
2 3 3 2.85 2.80
3 5 4 4.69 4.28
4 7 5 7.35 4.96
M = 4 is the cheapest setting that separates the two unmistakably: order 7
against 5, with the pressure error 52x larger at the finest step size. The
study script and the order test now use it, and the rule is recorded in the
problem class docstring and the project README so the M = 2 trap is not
stepped in again.
The test asserts the gap and the error ratio rather than absolute orders, and
still runs in about a minute.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Removing the time-dependent pressure condition on x = +-0.5 drops the observed pressure order from M+1 to M and grows the error by ~24x (measured at M=4: 4.96 -> 4.11). Prescribing it therefore supplies constraint information rather than being neutral, which places this benchmark at the "constraint lifting" rung of PR #641's ladder rather than the plain algebraic one, and means the 7-vs-5 gap understates a setup without prescribed boundary pressure. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Ships the remedy prototyped against the order reduction the benchmark measures.
Instead of evaluating the time-dependent boundary data pointwise at the node,
u_B(tau_m) = g(tau_m), it is imposed on the derivative and the stage value is
recovered by the collocation quadrature,
u_B(tau_m) = g(t0) + dt * sum_j Q[m,j] * gdot(tau_j).
The two differ by the quadrature error O(dt^(M+1)), but the second is consistent
with the collocation polynomial rather than pointwise exact. This is the
boundary-condition analogue of the differentiated-constraint remedy explored for
a time-dependent constraint in #641.
Measured at M=4, nelems=24, nu=0.1, pressure order and error at dt=0.1:
periodic (best possible) 6.32 7.4e-08
pointwise 5.74 9.7e-07
differentiated 6.30 1.3e-07
so it removes most of the penalty, leaving a factor of 1.8 against the periodic
case that is a constant rather than a rate.
No new sweeper was needed beyond a hook: generic_implicit_mass is already the
y-formulation, so generic_implicit_mass_diffbc only hands the problem the
collocation data of the step, which a problem class cannot see on its own. The
self-consistency trap #641 warns about does not apply here, because eval_f never
sees the boundary condition -- fix_residual zeroes those rows.
The observed orders are pre-asymptotic: the periodic reference does not reach
its design order 7 either, and on finer step sizes every variant collapses
against a solver floor near 1e-10. The tests therefore assert the error rather
than the order, and the docstrings say so; separating 2M-1 from 2M-2 needs a
better conditioned testbed than a 2D nonlinear NSE benchmark.
Verified that the differentiated variant solves the same problem: its solution
converges to the pointwise one at O(dt^5) = O(dt^(M+1)), the size of the
perturbation. Basing the quadrature on the exact g(t0) or chaining it from the
incoming numerical value gives results identical to three digits.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The ministry is now the Federal Ministry of Research, Technology and Space. The main README already carried the BMFTR logo while its text still said BMBF, so both it and the StroemungsRaum project README are updated, along with the link, which now points at bmftr.bund.de (verified to resolve). Grant numbers are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Coverage report
Click to see where and how coverage changed
This report was generated by python-coverage-comment-action |
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Nine findings from a review of this branch. Boundary conditions per node are now looked up by exact time rather than by nearest match, so a `solve_system` call at anything but a node of the prepared step fails instead of quietly using a neighbouring node's data. The arithmetic on both sides is identical, so the exact lookup costs nothing. `prepare_step` moves from `update_nodes` to `predict`, which the controller calls once per step rather than once per sweep -- what both docstrings already claimed. Measured over three steps at twelve iterations each: three calls instead of thirty-six. The Newton solver drops its `relative_tolerance = Sol_tol` and its reduced iteration budget, both of which pushed it towards raising on a marginal stall with `error_on_nonconvergence` left at its default. It now matches the merged monolithic class. The numbers in the `prepare_step` table were regenerated with this change in place and are unmoved: 7.411e-08 at dt = 0.1, order 6.32. `prepare_step` also gains the guard for the one mismatch that was missing, the `generic_implicit_mass_diffbc` sweeper on a problem with `differentiated_bc` off, which used to die on a missing attribute. `test_eval_f` takes `du/dt` from `_boundary_derivatives` instead of re-typing it, which anchors the velocity derivative the differentiated condition is built from, and checks its pressure half against a finite difference of `p_ex`, which nothing covered before. The rest is removal: the `dt_ref` branch of `order_study` that no caller used, the `results` dictionary `main` returned to nobody, and the `degree_rise` enrichment in `relative_errors`, which is pure cost for two functions in the same space -- verified to leave the result identical to four digits. `Sol_tol` becomes a keyword on `setup` defaulting to one decade below `restol`. A tenth finding, that the table in `prepare_step` mislabels its error column, was wrong: the numbers come from a four-point study whose dt = 0.1 entry is 7.411e-08 with order 6.32, so the label is correct and stands. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Supersedes #650, which has been stalled since July. That PR added a Taylor–Green
benchmark for order reduction as a new FEniCSx project; this reimplements it with
legacy FEniCS inside the existing
StroemungsRaumproject, fixes the issues raised inits review, and adds the order study it was missing along with a remedy for the effect
it measures.
The manufactured solution and forcing term from #650 were verified symbolically and
carried over unchanged — the mathematics there was correct, the implementation is what
needed work.
Why FEniCS instead of FEniCSx
Legacy FEniCS does periodic boundary conditions with
constrained_domainon thefunction space, so
dolfinx_mpcis not needed at all. Together with living in a projectthat is already in the CI matrix and already has a registered
fenicsmarker, thatdissolves three of the blockers on #650: the missing
dolfinx_mpcdependency (whichmade every test in that PR fail at collection), the project not being in CI, and the
unregistered
fenicsxmarker.Review issues from #650, addressed
solve_systemthere appliedMf^-1to the sweeper'sright-hand side only for the residual form to apply
Mfstraight back to it. The roundtrip is exact (verified:
1.7e-15), so it bought nothing while costing a solve per nodeper sweep and capping accuracy at that solve's tolerance. The right-hand side is now
subtracted as a vector in a
NonlinearProblem. The same round trip is removed from thealready-merged
NavierStokes_2D_monolithic_FEniCS.py, where drag and lift agree with theprevious implementation to 13 significant digits.
zeroed the residual on the periodic boundaries too, letting it stop iterating earlier than
the Dirichlet variant and biasing exactly the comparison the benchmark exists for.
Constant,rather than rebuilt on every
solve_systemcall.factorinstead of dividing by it, so there is nofactor == 0case and theresidual scale no longer depends on
dt.print, no misleading comments, nounused imports, correct
.. math::directives, no prints in tests.The order study that was missing
This is what #650 was written for, and what neither review comment on it got an answer to.
The number of collocation nodes decides whether anything is visible at all: RADAU-RIGHT with
Mnodes has design order2M-1and falls back to the stiff orderM+1, so the gap isM-2. Measured atnelems=24,nu=0.1:#650ran withnum_nodes = [2], where the gap is identically zero — that configurationcould not have demonstrated the phenomenon it was built for, whatever else was fixed. The
study and the test here use
M = 4.A remedy
Following the differentiated-constraint idea explored for a time-dependent constraint in
#641,
differentiated_bc=Trueimposes the boundary data on its time derivative and recoversthe stage values by collocation quadrature,
u_B(tau_m) = g(t0) + dt * sum_j Q[m,j] gdot(tau_j),instead of evaluating it pointwise. Pressure order and error at
dt = 0.1,M = 4:It needs no new sweeper beyond a hook, since
generic_implicit_massis already they-formulation;
generic_implicit_mass_diffbconly hands the problem the collocation data ofthe step, which a problem class cannot see on its own.
On what this does and does not show: the observed orders are pre-asymptotic — the periodic
reference does not reach its design order 7 either, and on finer step sizes every variant
collapses against a solver floor near
1e-10. The tests therefore assert the error rather thanthe order, and the docstrings say so. Separating
2M-1from2M-2needs a better conditionedtestbed than a 2D nonlinear NSE benchmark. The differentiated variant was checked to be solving
the same problem: its solution converges to the pointwise one at
O(dt^(M+1)), the size of theperturbation.
Also here
Prescribing the exact pressure on
x = ±0.5turns out not to be neutral — it acts as apartial lifting, worth
MtoM+1in the pressure order. Documented, because it means the7-vs-5 gap understates a setup without prescribed boundary pressure.
The funding ministry is renamed BMBF → BMFTR in the main README (whose logo was already
BMFTR while the text was not) and in the project README.
Tests
Nine tests, all under the
fenicsmarker in a project already in CI. FullStroemungsRaumsuite: 22 passed. Beyond the benchmark itself they cover
eval_fagainst the analyticaldu/dt(anchored independently, so a sign error shared withsolve_systemcannot hide), theperiodicity of the manufactured solution that the whole comparison rests on, and that
differentiated_bcfails loudly rather than silently doing nothing without its sweeper.🤖 Generated with Claude Code