Repository navigation
docs/doxy rewrite - #46
Merged
Merged
Conversation
satfra
commented
Oct 12, 2026
satfra
commented
Owner
- Fix physics function bugs found in the docs audit
- docs: fix math and diagram rendering, hide internals
- docs: rewrite the physics API reference
- doc fixes and minor bugfixes
- Fix timestepping/common bugs found in the docs audit
- docs: rewrite the timestepping and common API reference
- IDA: honour implicit/minimal_dt and maximal_dt; drop implicit/max_steps
- Timesteppers refuse variables and adaptors they cannot handle
- ImplicitEuler/TRBDF2: step control from local error estimates
- Euler/TRBDF2 steppers interpolate their frames onto the output grid
- IDA steppers: carry all spatial state through mesh adaptation
- CLI overrides and ConfigTree getters: keep types, never fall back silently
- Root finders: use the seeded exponent, restartable search, bisection fix
- QuadratureProvider: fix cache lookups
- CSV: keep infinities, make header names round-trip
- timestepping/common fixes
- Fix the LLF and KT dissipation sign and the LDG constant-levels
- Fix minor discretization bugs
- HAdaptivity: don't coarsen right behind a refinement front
- CG: assemble from the constrained solution; constrained dofs are algebraic
- Tests: hide the hybrid-stepper scans, trim the zero-dim KT Gaussian
- Remove dead code, move discretization internals into internal::
- docs: rewrite the discretization and model API reference
- discretization fixes
- Mesh default axis, size-checked HDF5Input::load_map, CsvOutput memory
- KT: mesh adaptation in 1D
- IDA steppers: own IDADriver instead of deal.II's SUNDIALS::IDA
- IDA: continue warm across mesh adaptations
- HAdaptivity: slope prolongation for KT, threshold marking
- IDA: change the Jacobian to a new alpha instead of reassembling it
- HAdaptivity: tolerance marking and a default curvature indicator
- ONfiniteT batched model: non-negative refinement indicator
- HAdaptivity: tolerance marking by default
- HAdaptivity: predict, new defaults, adaptor state in snapshots
- HAdaptivity: coarsen_age 10 by default
- Tutorials, API overview: adaptivity keys of the tolerance marking
- documentation fixes
- update ONfiniteT example: split into two dirs (one batched, one standard)
- SplineInterpolator1D(Stack): the second-derivative corrections were applied to the wrong grid points (lower/upper swapped). - BosonicRegulator::dq2RB (and hence dq2RF): wrong derivative. - dddcothS: missing factor 1/T. - TFLitimSpatial::F<3>: cosh argument x/(2T) instead of x/T. - PolynomialExpRegulator: Kokkos::exp instead of std::exp (device code).
- MathJax 3 with SVG output: MathJax 2's HTML-CSS output was misaligned (inline math raised) inside the pydata theme. - Class diagrams via Graphviz as SVG, inverted in dark mode; drop the include, collaboration and directory graphs. - Hide private members and all `internal` namespaces from the reference.
Integrators, interpolators, regulators, thermodynamic and threshold
functions now document what they compute (corrected formulas: the
vacuum integrators integrate dq q^{d-1} up to sqrt(x_extent) k, not
dp^2 p^{d-2} to infinity), the kernel interface for get/map, config
keys, setters, valid dimensions, and short usage examples. All math
uses Doxygen's \f$ / \f[ delimiters.
Small code changes alongside:
- include optimize.hh where optimize_x_extent is used;
- rename constructor parameters typical_E -> k where they set k; the
Integrator_fT_p2* constructors now also start their radial cutoff at
that k (default 1, so ConfigTree-built integrators are unchanged).
Timestepping:
- IDA callbacks: deal.II 9.7 callbacks return void, so the int returned on a
recoverable failure was discarded and IDA carried on with garbage. They now
throw RecoverableUserCallbackError (internal::throw_on_ida_failure). The
FE-only residual also catches model exceptions like the other runs.
- Newton matrix of the variables was alpha - dD/dv instead of alpha + dD/dv.
- BoostRK/BoostABM skipped the frame at the final time (0.3/0.1 -> 2).
- ImplicitEuler wrote the Newton residual into the variables slot and kept a
stale Jacobian after mesh adaptation.
- Newton ran max_steps + 1 iterations; RK divided by zero for output_dt < dt
and missed the final frame; KINSOL statistics were uninitialised.
- Copy-pasted error messages, undeclared update_variables().
Common / model:
- def::AD::jacobian_variables did not compile for pure-variables models and
used the wrong loop bound for the extractor Jacobian.
- GSLMinimizer1D: unchecked bracket setup, memory leaks, dead "stuck" check;
GSLSimplexMinimizer leaked its step vector.
- ExternalDataInterpolator clamped wrongly when all x < 0.
- CubicSpline is no longer copyable (double free); undefined
Barycentric::check_consistency removed.
- compute::conj return types, GLQuadrature static_assert that never fired,
empty Polynomial underflows, make_folder(""), to_string_with_digits,
uninitialised Quadrature/QuadratureType members, powr used before its
declaration in complex_math.hh, dead concepts and commented-out code.
- *_or_warn getters now also warn when a present key can not be read.
- The generated parameter file now uses the defaults the code uses.
Timestepping: an overview in abstract_timestepper.hh with the equations solved, a table for choosing a stepper (what it integrates, step control, mesh adaptation, MPI) and every /timestepping key with its real default and the steppers that use it. Each stepper, linear solver, Newton/KINSOL and the step controllers now state what they compute; wrong statements (Boost ABM "fixed steps", inverted Newton threshold, KINSOL doc copied from Newton, ...) are fixed.
deal.II's IDA wrapper never set step-size bounds, so both keys only fed the sanity checks. The deal.II patch now adds a custom_setup(ida_mem) hook, run after every (re)initialization, through which all IDA steppers set IDASetMinStep/IDASetMaxStep. Against an unpatched deal.II this is detected at compile time and skipped with a one-time warning; the dependency bundle has to be rebuilt for the bounds to take effect. implicit/max_steps never had an effect (the patch sets IDASetMaxNumSteps to LONG_MAX) and is removed from the code, the generated parameter file, tests and examples.
ImplicitEuler, TRBDF2 and ExplicitEuler integrate only the FE/FV part and silently kept the variables at their initial values; they now throw when the flow has variables. BoostRK, BoostABM and RK accepted a mesh adaptor but never called it; they now throw when given one other than NoAdaptivity.
Both steppers chose their step from Newton's convergence measure, so
abs_tol/rel_tol never controlled the time discretization error. They now
estimate the local error, damped in stiff directions with the Newton matrix
(Shampine's filter, no extra factorization):
- ImplicitEuler: (M + dt J) e = dt/2 (R(u_n) - R(u_{n+1}));
- TRBDF2: the embedded estimate of Bank et al. (1985).
A step whose weighted RMS error exceeds 1 is rejected and retried with a
smaller step (accepted regardless at minimal_dt).
TRBDF2's BDF2 stage is divided by (2 - gamma), so both stages have the
Newton matrix M + gamma dt J / 2; before, the shared Jacobian was off by a
factor sqrt(2) in the second stage.
Newton now tests convergence on the weighted RMS of its update (<= 0.1 of
the step tolerance) instead of the mass-weighted residual, whose size
depended on the mesh and left errors far above the tolerance.
ExplicitEuler, ImplicitEuler and TRBDF2 wrote a frame at the first step past
each multiple of output_dt, so frame times drifted with the step size. They
now write frames exactly at t_start + n output_dt, interpolated within the
step that reaches them: linearly for the Euler steppers, quadratically
through u_n, u_gamma, u_{n+1} for TRBDF2 (second order, like the method).
AbstractAdaptor::adapt_at(t, sol, also_transfer) adapts and moves further vectors on the same mesh along with sol; HAdaptivity transfers them in one SolutionTransfer pass. On a refinement: - the IDA hybrids transfer their spatial bracket, pending state and checkpoints, which kept the old size and were read by the next dt_variables call; - all IDA runs transfer sol_dot instead of zeroing it, so IDA's consistent-IC solve starts from the previous derivative; - the plain IDA runs also re-initialize their solution views, and the run with variables updates its block sizes and n_FE_dofs. The hybrids' variable buffer is pruned behind each output frame (nothing earlier is queried again) and searched by bisection; it used to grow with the run and be scanned linearly.
…ently ConfigTree: - get_X(key, def) and get_X_or_warn(key, def) fall back to def only for a missing key; a present entry that cannot be read as X throws. Before, e.g. fe_order = 2.0 silently became the default. ConfigurationHelper: - -sd/-si/-sb/-ss keep the type of an existing entry (an integer entry takes integral values from -sd or -si; booleans need -sb, strings -ss), values must parse completely, -sb takes true/false/1/0 in any case, and -ss splits at the first '=' only. - A key that is not in the parameter file is created with the flag's type, with a warning naming it. - --generate-parameter-file writes TOML; --gen-toml/--gen-json choose the format. The file is named like the parameter file, written by rank 0 only, never replaces an existing file, and -p may come anywhere on the line. - append(), which never applied anything, is removed; the class is an ordinary copyable value (its copy constructor always threw).
ScalingRootFinder: - seed_exponent(beta) had no effect on the search. With a seed, the convergent model is fitted from two points when the three-point fit is not yet available or not valid, which saves an evaluation when no divergent fit stands in. - search() starts from scratch on every call (only the configuration is kept); a second call used to continue from the first one's samples. - The deep-scaling approach step uses approach_factor instead of a hard-coded 0.1. BisectionRootFinder converges on the bracket width and returns the smallest point known to succeed. It used to throw whenever the last probe failed and to return unconverged when the last probe succeeded early.
- A key slightly above every cached key was never matched: lower_bound returns end() there and the predecessor was not checked, so a new rule was built instead. Both neighbours are now checked independently. - Keys match to a relative 1e-6. is_close(a, b, 1e-6 T) used 1e-6 T as both the absolute and the relative tolerance (1e-3 relative at T = 1000), and never matched a negative key, whose entry was then re-initialized on every lookup.
- inf and -inf are read as infinities instead of NaN; ExternalDataInterpolator now drops rows with any non-finite value. - write_csv_header turns line breaks in names into spaces (the reader splits records at line breaks) and quotes names that start with the comment character, which used to turn the whole header line into a comment.
- ExplicitEuler, ImplicitEuler and TRBDF2 also write a frame at the final time when it is not on the output grid (OutputGrid::finish). - Newton decides on a Jacobian rebuild from the residual reduction of the previous iteration. - The local-error solve of ImplicitEuler/TRBDF2 asks for 1e-3 relative accuracy instead of abs_tol, which an iterative solver may not reach.
- def::LLFFlux added -alpha/2 (u_n - u_s) to every Cartesian component of the numerical flux instead of -alpha/2 (u_n - u_s) n_d. - The KT numerical flux had the same defect, H.n = avg.n - (a.n) du / 2: anti-dissipative on every lower boundary face. The dissipation is now a_d n_d du / 2 per direction, in the residual and the jacobian. - LDG: a level declared constant was built as j_gu[k] u even if a level below it was nonlinear, and j_gu[k] was never rebuilt. - ComponentDescriptor: couplings between non-adjacent levels are the closure of the adjacent ones, independent of the order of add_dependency calls; jacobians_constant(k, j) is true iff all levels in between are.
Mesh and coordinates - HAdaptivity never refines past adaptivity/level and never coarsens below the initial mesh. - LogarithmicCoordinates1D::backward_derivative: denominator b + y - c. - Linear/LogarithmicCoordinates1D reject grid_extent < 2 (and bias == 0). - to_string() of all 1D grids and Matsubara values: grids differing below 1e-6 no longer share an HDF5 coordinate key. - CoordinatePackND::operator== is false for packs of different axes instead of failing to compile. FV - TVDReconstructor::compute_gradient_at_point compares the point with the cell centre up to rounding, so a source point at the centre gets the limited slope even if FEValues and the cell disagree in the last bit. - van Albada: sign test and denominator without underflow (float and tiny slopes gave 0 or 0/0). Model / assembly / output - Readouts locate their EoM with EoM_postprocess, as extract() does. - The readout raw potential (a global solve per output frame) is skipped if no sink writes it and the model opts out with readouts_use_potential = false (and its extract() does not read it). - OutputSession(config) without /output/folder writes to ./ like OutputPath(config) instead of a temporary folder deleted at exit. - FlowingVariables::variable_data() throws if there are no variables. - ADjacobian_extractors<2> checks the hessian size. - 1D EoM root polish: moves the reconstructed minimum onto the zero of the EoM field near it, picking the candidate with the lowest potential, so it never jumps to a metastable minimum.
Coarsening re-interpolates the solution onto the parent cell. Right behind a moving front that overshoots: in ONfiniteT (CG) the flat, convexity-restoring region sits ~1e-8 above the propagator pole m2 = -k^2, the interpolated state crossed it, and IDA's restart failed with a NaN residual at t ~ 3.5. A cell is now never coarsened if it lies within /discretization/adaptivity/coarsen_buffer (default 2) layers of active neighbours of a cell flagged for refinement, or if one of the last /discretization/adaptivity/coarsen_age (default 2) adaptations created it. ONfiniteT CG now runs to t = 5 and follows a uniform mesh of the same resolution.
…braic The timesteppers never apply the hanging-node (and model affine) constraints, and the constrained rows of the residual are empty, so after a refinement in 2D/3D the hanging-node entries of IDA's vector stayed at their transferred values while their masters evolved -- and the neighbouring cells read them. - CG mass, residual and jacobians, the refinement indicator, the readouts and output, and the variables' residual read u with the constraints applied and u_dot with their homogeneous part (D(u + u_dot) - D(u), D affine). - get_differential_indices excludes constrained dofs: they are fixed by the constraints, not by an equation. - HAdaptivity applies the constraints before transferring vectors (that hunk went into e7abe1f).
- sundials_ida_boost: the six cur_dt / coupling-mode scans (46 full runs that assert nothing and only print errors) are hidden ([.][scan]); run them with "[scan]". The seven gating cases are unchanged. 125 s -> 18 s. - zero_dim_2d_kt gaussian: drop the 400x400 run and the separate EoM case; the 40x40 matrix of 16 scenarios checks the same observables (including that scenario). The model skips the readout raw potential, which only the output would read. Still [slow]. 429 s -> ~85 s.
Removed - AbstractModel::mass(m_ij, x). - ADjacobian_variables::jacobian_variables<1> (never called; <0> asserts). - SubDescriptor::names_tuple, ComponentDescriptor::get_component_block_structure. Renamed - GridAxis::get_stepwiths -> get_step_widths, ConfigurationMesh::get_step_withs_for_triangulation -> get_step_widths_for_triangulation. Moved into internal:: - la_policy.hh free functions, ParallelDoFs, serial_mirror, MeshIsRectangular. - batch.hh: PointState, PointJacobian, FaceJacobian, SeedStackWorkspace, batch_terms, contains(Term, Term) and the evaluate_* dispatchers. Models keep PointBatch, LDGPointBatch, BatchOutput, Term, batch_tie. - CG/DG i_tie, Variables fe_tie, def::BoundaryStencilDetail. - HDF5 frame/writer types and the output timing helpers. - FV::KurganovTadmor::Assembler: helpers private; public are the AbstractAssembler interface, the readout/extractor hooks the core calls and the few members the KT unit tests use.
- Constrained solution: internal::ConstrainedSolution applies the constraints once per entry point (u and C u_dot) and copies nothing when there are no constraints; replaces the constrained()/constrained_dot() pair, which made three copies and three distribute() calls per evaluation. - LDG: a level declared constant in the level below, while a lower level is not constant, reuses its own jacobian; that reuse is now checked against the assembled level on first use as well (the docs promised it). Test added. - EoM polish: locate the polished point with the mapping, as everywhere else. - Cleanups: jacobian_variables without the leftover block, one shared variables_block() helper, clearer names and comments in HAdaptivity and KT, no stale diagnostic prints in the o2 tests.
- A missing /discretization/grid/{x,y,z}_grid falls back to the axis that
ConfigurationMesh::get_defaults() writes (GridAxis::get_default_range())
instead of a different hard-coded "0:0.1:1".
- HDF5Input::load_map takes the buffer as std::span<double> and throws if its
size differs from the map's; the raw-pointer overloads remain, deprecated.
Examples and tut4 pass their vectors directly.
- CsvOutput keeps only the current frame's values instead of every value
ever written.
- Drop the declared-but-undefined FEOutput() default constructor.
- The KT assembler finds the active cell across a face also when it is refined (FV::KurganovTadmor::internal::active_face_neighbor: its child at the face, 1D). Stencils, face descriptors, traces and the jacobian sparsity all go through it; the reconstruction and diffusion gradients already use the true cell centres. - refinement_indicator: the model's cell_indicator at each cell centre, with the cell average and first/second derivatives from the averages of the cell and its neighbours, times the cell width. - Assembler::supports_mesh_adaptation (true in 1D); HAdaptivity refuses at compile time an assembler that sets it false, i.e. KT in 2D. - Test kt_adaptivity: on a 1D mesh with level jumps in both directions the KT residual is conservative, no less accurate than on the unrefined mesh, and its jacobian matches finite differences; HAdaptivity refines where the indicator is large and conserves the integral of u. - Examples/ONfiniteT: KT.cc uses HAdaptivity, model_KT.hh gets a cell_indicator; both models use |m2''| (the CG one was signed, which the fixed-fraction marking does not expect).
The IDA timesteppers drive SUNDIALS IDA through DiFfRG::IDADriver (timestepping/ida_driver.hh), which follows deal.II's interface but - bounds the step size to [minimal_dt, maximal_dt] directly (the deal.II patch for step bounds, ida_step_bounds.hh, is no longer needed); - uses its own Newton convergence test that measures the rate in every step instead of carrying it over; - exposes IDA's statistics and error test vectors. The default implicit minimal_dt is 1e-14.
After a restart (mesh adaptation, tolerance refresh) IDADriver now restores IDA's step history instead of starting over at order 1 with IDACalcIC: the divided-difference vectors phi are handed to the adaptor (history()) and transferred to the new mesh with the solution, IDA's step-control scalars are saved and put back. Setting /timestepping/implicit/warm_restart = false restores the cold restart.
KT (1D): the children of a refined cell now get the cell average plus a slope times the offset of their centre, instead of the bare average. The slope is that of the parabola through the cell's and its neighbours' averages (exact child averages for a monotone quadratic profile), scaled down where needed so that the children stay within the range of these averages. New /discretization/adaptivity/marking = "threshold" (default "fraction"): refine above refine_threshold * max indicator, coarsen below coarsen_threshold * max. A mesh that resolves the solution is then kept, and adaptations that only coarsen are carried out.
IDA asks for a new Jacobian F_u + alpha F_udot whenever a change of step size or order moves alpha out of [0.6, 1.67] times its value at the last setup. IDADriver now calls the optional update_jacobian_alpha for such requests, which adds (alpha_new - alpha) * F_udot (AbstractAssembler::jacobian_mass, at the state of the last assembly) and refactorizes. It still reassembles after a Newton failure since the last assembly, inside IDACalcIC, and 5 steps after the last assembly. JacobianAlphaUpdate implements it for TimeStepperSUNDIALS_IDA without variables and both IDA hybrids. /timestepping/implicit/reuse_jacobian = false turns it off.
marking = "tolerance" refines cells with an indicator above tolerance and coarsens those below coarsen_tolerance (default tolerance / 8). With it, a mesh that resolves the solution is kept, and refinement starts only where a feature becomes unresolved. It replaces the relative "threshold" marking, which always refined the cells near the largest indicator. indicator = "curvature" uses eta_K = max_c h_K^2 |u_c''| / range(u_c), about 8 times the relative deviation of u from its linear interpolant on the cell. It grows where a kink, shock or narrow extremum forms and drops by 4 per refinement. FV: the parabola curvature of the KT assembler (cell_curvatures; boundary cells take their neighbour's), FEM: DerivativeApproximation. With the default indicator = "auto" it is used for models that define neither cell_indicator nor face_indicator; AbstractModel's defaults now return def::DefaultIndicator so that HAdaptivity can tell.
/discretization/adaptivity/marking now defaults to "tolerance" (1e-3); refine_percent and coarsen_percent only apply with marking = "fraction". Since the tolerance is absolute, a model's own indicator needs a known scale. Tests about the transfer and prolongation, not the marking, keep "fraction". tut3's indicator was signed; now its absolute value.
predict (default on) adds each cell's indicator growth since the previous adaptation, so a forming kink is refined one adaptation earlier. Defaults now: tolerance 3e-4, coarsen_tolerance tolerance / 8, coarsen_age 5 (counted in mesh changes), adapt_dt 0.2. Snapshots now also store the adaptor's state (AbstractAdaptor::save_state / load_state; HAdaptivity: mesh-change count, cell ages, last indicators), in an optional "adaptor" group. Without it, a restarted adaptive flow made different coarsening decisions than the uninterrupted one.
At the default adapt_dt 0.2, coarsen_age 5 still coarsened the flat region of fRG potentials that is refined again shortly after: tut7 CG needed 2507 instead of 1302 residual evaluations, DG 3306 instead of 2160, ONfiniteT 10-15% more. With 10, tolerance marking needs fewer residuals than fraction marking in 4 of 5 measured cases (tut7 CG +12%) and still coarsens (ONfiniteT meshes 12-17% smaller than without coarsening).
The tutorial parameter files list the current adaptivity defaults (adapt_dt 0.2, tolerance 3e-4) instead of refine_percent/coarsen_percent, and the discretization overview describes the tolerance marking.
satfra
force-pushed
the
docs/doxy_rewrite
branch
from
October 12, 2026 16:23
8a22dc8 to
cdb063b
Compare
…ialExpRegulator Kokkos' overflow-safe complex division costs ~6 real divisions per reciprocal; conj(z)/|z|^2 needs one. The regulator series sum_m x^m/m becomes Horner with compile-time reciprocals.
satfra
force-pushed
the
docs/doxy_rewrite
branch
from
October 12, 2026 16:51
cdb063b to
6235ed2
Compare
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.