ndiffusion

Multigroup neutron diffusion solver for 1-D and 2-D geometries. Written in C++17; exposed to Python via pybind11.

Capabilities

1-D (slab, cylinder, sphere)

• Arbitrary number of energy groups and material regions
• Vacuum, reflective, and albedo boundary conditions
• k-eigenvalue solver - matrix-free power iteration; Aφ = (1/k)Bφ
• Fixed-source solver - direct solve of Aφ = q for a user-supplied volumetric source
• Time-dependent solver - theta-weighted time stepping, unconditionally stable
• Per-group Thomas (TDMA) tridiagonal solver inside a Gauss-Seidel group sweep
• Harmonic-mean diffusion coefficients at material interfaces, over the center-to-center distance, so non-uniform edges_x is second-order accurate

Reactor kinetics (all three dimensionalities)

• Delayed neutron precursors - any number of precursor groups, per-material delayed fractions and delayed fission spectra
• Implicit fission source, so the scheme stays unconditionally stable through a supercritical transient
• Second-order time differencing - theta = 0.5 is Crank-Nicolson; the default theta = 1 is backward Euler
• Mid-transient perturbation via update_materials, for reactivity insertions

2-D structured (Cartesian XY or axisymmetric RZ)

• Finite-difference 5-point stencil on an nx x ny Cartesian grid
• Left (x=0) and bottom (y=0) boundaries hardcoded as reflective; right and top boundaries take user-specified Robin BCs per group
• k-eigenvalue solver - line-TDMA x-sweeps inside a Gauss-Seidel outer iteration
• Fixed-source solver - same spatial sweep; solves Aφ = q directly
• Time-dependent solver - theta-weighted stepping using the same line-TDMA sweep

2-D unstructured (triangles, quadrilaterals, and higher polygons)

• Cell-centered finite-volume method (FVM); cells may be any simple polygon, so hexagonal lattices work directly
• Deferred non-orthogonal correction, so skewed and triangular meshes stay second-order rather than converging to the wrong answer
• Arbitrary Robin BCs per boundary tag; harmonic-mean interface diffusion coefficients
• Mesh generation and material assignment are separate, so one geometry can drive several material layouts (assign_materials)
• Connectivity is validated: non-manifold edges, hanging nodes, unmerged coincident vertices and degenerate cells are all rejected
• k-eigenvalue solver - power iteration with point Gauss-Seidel inner solve
• Fixed-source solver - point SOR (successive over-relaxation) inner solve
• Time-dependent solver - theta-weighted stepping with point Gauss-Seidel

Installation

Requires a C++17 compiler, Python >= 3.9, and pybind11 >= 2.12.

pip install .

For development:

pip install -e ".[dev]"

The dev extra pulls in pytest, scipy, and ruff - everything CI runs. Installing without it leaves parts of the test suite silently skipped rather than failing (see Running tests).

Python edits under src/ndiffusion/ are picked up immediately; after editing C++ sources, re-run pip install -e ".[dev]" to rebuild the extension.

Quick start

1-D k-eigenvalue

import numpy as np
import ndiffusion as nd

m = nd.Materials()
m.n_mat    = 1
m.n_groups = 1
m.D        = [3.850204978408833]
m.removal  = [0.1532]
m.scatter  = [0.0]
m.chi      = [1.0]
m.nusigf   = [0.1570]

cells = 50
edges = list(np.linspace(0.0, 100.0, cells + 1))

solver = nd.KEigenSolver(
    mats       = m,
    medium_map = [0] * cells,
    edges_x    = edges,
    geom       = nd.Geometry.Sphere,
    bc         = [nd.BoundaryCondition(A=1.0, B=0.0)],
    epsilon    = 1e-8,
    max_outer  = 500,
)
result = solver.solve()
assert result.converged
print(f"keff = {result.keff:.8f}")   # -> 1.00000475

Every result carries a converged flag; always check it before trusting the answer (an unconverged run returns the last iterate without raising).

2-D structured k-eigenvalue

solver = nd.KEigenSolver2D(
    mats       = m,
    medium_map = [0] * (nx * ny),
    edges_x    = list(np.linspace(0.0, R, nx + 1)),
    edges_y    = list(np.linspace(0.0, R, ny + 1)),
    geom       = nd.Geometry2D.XY,
    bc_x       = [nd.BoundaryCondition(A=1.0, B=0.0)],   # vacuum right
    bc_y       = [nd.BoundaryCondition(A=1.0, B=0.0)],   # vacuum top
)
result = solver.solve()
flux = np.array(result.flux).reshape(nx, ny, m.n_groups)

2-D unstructured fixed-source

# Build an unstructured mesh (vertices + connectivity + boundary faces)
mesh = nd.UnstructuredMesh2D()
mesh.vx = vx; mesh.vy = vy
mesh.cell_vertices = cell_vertices
mesh.cell_offsets  = cell_offsets
mesh.material_id   = mat_ids
mesh.bface_v0      = bface_v0
mesh.bface_v1      = bface_v1
mesh.bface_bc_tag  = bface_bc_tag   # integer tag per face

bc = [nd.BoundaryCondition(A=1.0, B=0.0)]   # tag 0 -> vacuum

solver = nd.FixedSourceSolverUnstructured2D(
    mats      = m,
    mesh      = mesh,
    bc        = bc,
    epsilon   = 1e-10,
    max_inner = 1000,
    omega     = 1.9,    # SOR relaxation factor
)
result = solver.solve([q] * n_cells)   # volumetric source per cell

See examples/k_eigenvalue.py and examples/time_dependent.py for further examples.

Meshes and material assignment

An UnstructuredMesh2D carries a material_id per cell, but nothing requires it to be decided when the geometry is built. assign_materials paints it on as a separate step, so a single geometry serves many material layouts:

mesh = nd.load_gmsh("core.msh")          # geometry + region labels

# by Gmsh physical-group name
nd.assign_materials(mesh, {"fuel": 0, "reflector": 1})

# by region id
nd.assign_materials(mesh, {1: 0, 2: 0, 3: 1})

# by position, evaluated at cell centroids
nd.assign_materials(mesh, lambda x, y: 0 if x*x + y*y < R*R else 1)

It rewrites material_id in place and returns the mesh; pass copy=True to keep the input intact. Repainting between solves is safe, because each solver takes its own copy of the mesh at construction:

for name, painter in layouts.items():
    nd.assign_materials(mesh, painter)
    results[name] = nd.KEigenSolverUnstructured2D(mats[name], mesh, bc).solve()

A dict spec must map every region the mesh uses, and names are checked against the mesh's own region_names, so adding a physical group to the .msh is an error rather than a silent shift of every index after it.

Supporting geometry queries, useful for painting and for post-processing:

cx, cy = nd.cell_centroids(mesh)   # area-weighted; matches what the solver uses
areas  = nd.cell_areas(mesh)       # e.g. for flux-weighted averages
nd.validate_mesh(mesh)             # connectivity check the solvers run anyway

On a large mesh, precomputing the assignment with NumPy is about twice as quick as a per-cell callable:

cx, cy = nd.cell_centroids(mesh)
nd.assign_materials(mesh, np.where(np.asarray(cx) > x0, 1, 0))

Preset layouts and symmetry orientations

ndiffusion.layouts composes three independent pieces: a geometry, an orientation (the symmetry sector, and which boundaries are cuts), and a layout (a painter (x, y) -> material).

from ndiffusion import layouts

mesh = layouts.cartesian_mesh(width=100.0, h=2.0, orientation="quarter")
nd.assign_materials(mesh, layouts.core_reflector(core_radius=35.0))
bc = layouts.boundary_conditions(mesh, D=[1.4, 0.4], albedo=0.0)
keff = nd.KEigenSolverUnstructured2D(mats, mesh, bc).solve().keff

Geometry	Orientations
cartesian_mesh	full, half, quarter, eighth (45° octant), infinite
hex_mesh	full, half, sector120, sector60, sector30, infinite

Cuts are reflective by default. symmetry="rotational" joins the two cuts periodically instead, imposing rotational symmetry without a mirror — the right choice for a spiral or pinwheel loading, where reflecting solves a different problem. Only whole rotational periods qualify (Cartesian half and quarter; hex sector120 and sector60); the 45° octant and 30° hex wedge are fundamental domains only by virtue of the mirror, and are rejected.

Sector cuts slice cells into smaller polygons, which the FVM solver takes directly. Boundary faces are tagged layouts.SYMMETRY or layouts.OUTER, and boundary_conditions pairs them — reflective on the cuts, Marshak with the given albedo outside. infinite tags every boundary reflective, giving k∞.

Layout
homogeneous()	one material; baseline
core_reflector(core_radius=…) or (core_half_width=…)	two regions, circular or square core
checkerboard(pitch)	alternating materials; maximal heterogeneity
annular(radii, materials)	concentric rings
hex_rings(pitch, ring_materials)	material per hexagonal ring
with_rods(base, positions, radius, material)	overlays rods on another layout

Painters are plain functions of position, so they compose and work on any geometry. with_rods gives the perturbation pair the mesh/material split was built for — one geometry, two layouts, a rod worth:

unrodded = layouts.core_reflector(core_radius=35.0)
rodded   = layouts.with_rods(unrodded, [(0.0, 0.0)], 5.0, material=2)

nd.assign_materials(mesh, unrodded); k0 = solve(mesh)
nd.assign_materials(mesh, rodded);   k1 = solve(mesh)
rho = (k1 - k0) / (k1 * k0)

Periodic boundaries

Edges listed in mesh.periodic_a0/a1 are joined to periodic_b0/b1 — vertex to corresponding vertex — and become interior faces rather than boundaries, so the flux is continuous across them. That correspondence fixes the rigid transform, so both a translation (a repeating lattice) and a rotation (a symmetry sector) are expressible, and the solver derives which from the geometry.

A fully periodic homogeneous square has no leakage anywhere, so it reproduces k∞ exactly. layouts builds the rotational case for you via symmetry="rotational".

Named cross sections and benchmarks

ndiffusion.materials holds builders and published tables. Geometry is deliberately not bundled, so a table runs on any mesh:

mats = nd.materials.BIBLIS.materials()
nd.assign_materials(mesh, nd.materials.BIBLIS.layout())
assert abs(keff - nd.materials.BIBLIS.reference_keff) < 1e-3

one_group(D, sigma_a, nusigf, n_mat)	scalars or per-material sequences
two_group(rows, axial_buckling=0)	rows of (D1, D2, Sa1, Sa2, S12, nuSf1, nuSf2)
from_assembly_map(amap, pitch)	painter from a published assembly grid

Benchmarks: RINGHALS (1.0037), TWIGL (0.9133), IAEA (1.0296), BIBLIS (1.02535) — each with .materials(), .reference_keff, .source, and .layout() where the loading is an assembly map. tests/test_benchmarks.py solves these tables to their published eigenvalues.

Transport cross sections

Multigroup transport libraries tabulate a total (or absorption) cross section, a scatter matrix and fission data; the solvers want D, a removal cross section and a diagonal-free scatter matrix. make_materials_from_transport does the conversion and returns a ready-to-use Materials:

mats = nd.make_materials_from_transport([uo2, mox, moderator], G=7)

Each input dict holds SigT or Siga (the other is derived from the scatter matrix), a Scat matrix, nuSigf, and optionally chi, SigTr or Scat1:

D[g] = 1 / (3 Sigma_tr[g])	Sigma_tr from a tabulated SigTr, else the P1 outflow correction SigT - sum_g' Scat1[g->g'] (transport_correction="none" gives the uncorrected 1/(3 SigT))
Sigma_r[g] = SigT[g] - Scat[g->g]	the outflow correction cancels here, so removal uses the uncorrected total
scatter[g_to][g_from]	input is assumed Scat[g_from][g_to] (the transport convention) and is transposed; pass scatter_orientation="to_from" for data already in solver order

transport_to_diffusion(data, G) exposes the same transform for a single material and returns a plain dict, useful for inspecting the derived D, Removal and SigTr before building a Materials.

See examples/transport_cross_sections.py for a runnable end-to-end example.

Reactor kinetics

The time-dependent solvers model delayed neutron precursors:

(1/v_g) dphi_g/dt = -A_g phi_g + scatter
                  + (1-beta) chi_p,g F + sum_i chi_d,i,g lambda_i C_i
        dC_i/dt   = beta_i F - lambda_i C_i,   F = sum_g' nuSigf_g' phi_g'

The time discretization weights the right-hand side between the two time levels with a factor theta (see Time differencing below). Whatever the weight, C^{n+1} still eliminates in closed form, which folds the delayed source into a dt-dependent effective fission spectrum plus a source known from the old precursors:

chi_eff,g = (1-beta) chi_p,g
          + sum_i chi_d,i,g beta_i lambda_i (theta dt) / (1 + lambda_i theta dt)
Q_d,g     = sum_i chi_d,i,g lambda_i
            [ C_i^n (1 - (1-theta) lambda_i dt) / (1 + theta lambda_i dt)
            + (1-theta) dt beta_i F^n           / (1 + theta lambda_i dt)
            + ((1-theta)/theta) C_i^n ]

At the default theta = 1 those collapse to the backward-Euler forms chi_eff,g = (1-beta) chi_p,g + sum_i chi_d,i,g beta_i lambda_i dt / (1 + lambda_i dt) and Q_d,g = sum_i chi_d,i,g lambda_i C_i^n / (1 + lambda_i dt).

As theta*dt -> 0 the effective spectrum tends to (1-beta) chi_p (prompt only); as theta*dt -> inf it tends to the total fission spectrum, so a critical system with equilibrium precursors is an exact fixed point at any step size and any theta. Fission is evaluated at the new time level inside the Gauss-Seidel sweep, so the scheme stays unconditionally stable.

delayed = nd.make_delayed_data(nd.DELAYED_U235_6GROUP, G=2, n_mat=3, chi=mats.chi)

# Start from a genuine steady state: a k-eigenvalue flux is only stationary
# once nusigf is divided by keff.
res = nd.KEigenSolver2D(mats, medium_map, edges_x, edges_y, geom, bc_x, bc_y).solve()
critical = nd.scale_to_critical(mats, res.keff)

solver = nd.TimeDependentSolver2D(
    mats=critical, medium_map=medium_map,
    edges_x=edges_x, edges_y=edges_y, geom=nd.Geometry2D.XY,
    bc_x=bc_x, bc_y=bc_y, initial_flux=res.flux,
    delayed=delayed,          # omit for prompt-only kinetics
)
solver.update_materials(perturbed)   # step insertion at t = 0
out = solver.run(dt=2e-3, n_steps=250)
out.precursors                        # [cells * n_precursor], per unit volume

Precursors default to equilibrium with the initial flux, which is what a transient starting from steady state needs; pass initial_precursors to override. Drive a ramp by calling update_materials once per step with interpolated cross sections.

DelayedNeutronData arrays are flat: lambda_ is [n_precursor], beta is [n_mat * n_precursor], chi_delayed is [n_mat * n_precursor * n_groups] (each spectrum must sum to 1), and chi_prompt is [n_mat * n_groups] or empty.

When chi_prompt is omitted it is derived rather than defaulted to Materials.chi:

chi_p = (chi - sum_i beta_i chi_d,i) / (1 - beta)

so the prompt and delayed parts always add back up to the total spectrum, and chi_eff still tends to chi as dt -> inf for any delayed spectrum. With the usual chi_delayed = chi this reduces to chi_p = chi. Pass chi_prompt explicitly only if Materials.chi is itself the prompt spectrum.

Fission-matrix mode is supported. There is no separable spectrum, so the split is applied to the matrix: the production cross section is the column sum P[g'] = sum_g F[g][g'] (total neutrons emitted per fission caused by a group-g' neutron), the delayed yield beta_i chi_d,i[g] P[g'] is subtracted from the tabulated matrix, and the part emitted within the step is added back. The two representations agree exactly when the matrix is separable. Note that Materials.chi is all zeros in this mode, so it cannot serve as the ChiDelayed fallback - supply one.

Time differencing. All three time-dependent solvers take a theta weight, as a constructor argument and as a settable property:

(1/(v dt)) (phi^{n+1} - phi^n) = theta * R(phi^{n+1}, C^{n+1})
                               + (1-theta) * R(phi^n,   C^n)

theta = 1 (the default) is backward Euler, first order in dt; theta = 0.5 is Crank-Nicolson, second order. The weighting is applied consistently to the flux equation and the precursor balance, so the whole transient is second order at theta = 0.5 - on the infinite-medium point-kinetics problem in tests/test_kinetics.py the observed order is 2.00, and at a fixed dt the error is ~140x smaller than backward Euler's:

solver = nd.TimeDependentSolver(..., delayed=delayed, theta=0.5)
solver.theta = 1.0    # or change it mid-transient

Values outside [0.5, 1] raise ValueError: that is exactly the A-stable range, and below it the fast spatial modes - the ones quantified below, decaying at ~1e5 to 1e6 per second - would diverge at any useful step size.

Both ends of the range are unconditionally stable, but only backward Euler damps the stiff modes. theta = 0.5 is A-stable and not L-stable: a mode with zeta = |lambda| dt >> 1 has amplification (1 - (1-theta) zeta)/(1 + theta zeta), which is ~1/zeta at theta = 1 but tends to -1 at theta = 0.5, so it decays slowly with an alternating sign instead of being killed.

Whether that ever shows up depends on how much stiff content the state actually contains. The stiff modes here are the mesh-scale spatial harmonics: on the 8-cell slab in tests/test_kinetics.py the checkerboard mode decays at ~6e5 / s, so dt = 1e-3 puts it ~560x beyond the resolved range and a localized flux bump still retains ~80% of its size after 20 steps, reversing sign on every one, while backward Euler annihilates it in a single step. A smooth, mode-shaped perturbation excites almost none of that - the TWIGL step insertion in tests/test_benchmarks.py moves the flux shape by only ~1e-5 and shows no ringing at all. Nothing diverges either way. When you do have a stiff perturbation, damp it first and then switch:

solver.update_materials(perturbed)   # step insertion at t = 0
solver.theta = 1.0
solver.run(dt, 2)                    # two damped steps
solver.theta = 0.5
solver.run(dt, n_steps)              # second order from here on

update_materials takes effect at the start of the next step, which is what a step insertion at t_n means: the new cross sections hold across the whole of [t_n, t_n + dt], the explicitly weighted term included.

Choosing dt and max_inner. An implicit fission source means the inner Gauss-Seidel sweep resolves the multiplication as well as the scatter coupling, and only the 1/(v_g*dt) diagonal term keeps that iteration contracting. When 1/(v*dt) << Sigma_r a near-critical problem converges at roughly k per sweep. The solvers apply Aitken extrapolation to that fixed point: the convergence ratio is estimated from successive iterate changes and, once it has held steady, the iterate jumps to the limit of the geometric series. In the worst case tested - exactly critical, zero leakage, 1/(v*dt) at 3% of Sigma_a - this cuts the iterations per step from thousands to a few dozen. The extrapolation is safeguarded (ratio stability judged against 1 - sigma, and the jump capped relative to ||phi||) and is self-correcting, since convergence is still measured across the sweep. Any step that nonetheless hits max_inner prints a warning to stderr naming the solver and the residual - never trust a transient that warned.

See examples/kinetics.py for a runnable end-to-end transient.

Boundary conditions

Type	A	B
Zero-flux (approx. vacuum)	1.0	0.0
Marshak vacuum	(1-α)/(4(1+α))	D/2
Reflective	0.0	1.0

The ndiffusion.boundary_conditions(Dg, alpha) helper constructs the coefficient array from an albedo value alpha (0 = vacuum, 1 = reflective).

The 1-D and 2-D structured solvers take one BoundaryCondition per energy group per boundary. The unstructured solvers index bc by boundary tag:

bc[tag * n_groups + g]      # length = n_bc_types * n_groups

so a two-group problem needs two entries even for a single tag. The constructor rejects a length that is not a multiple of n_groups, or that fails to cover every tag the mesh uses - otherwise those boundaries would silently behave as reflective, and an all-reflective system just reads as k∞.

Adjoint & solution verification

Two Python helpers layer on top of the compiled solvers (they reuse the existing solver classes, so no rebuild is involved):

• Adjoint - ndiffusion.make_adjoint_materials(mats) returns the adjoint cross sections (group scatter transposed, chi/nusigf swapped). Running any solver on them solves the adjoint (importance) problem; the k-eigenvalue is identical to the forward one, and the flux is the neutron importance function.

• Method of nearby problems (MNP) - a discretization-error estimator (ndiffusion.nearby_fixed_source, ndiffusion.nearby_k_eigenvalue). It fits a smooth curve through the numerical flux, substitutes it into the continuous diffusion operator to form a residual source, and re-solves the resulting "nearby problem" whose exact solution is the fit - so nearby - fit estimates the true error. Works in 1-D, 2-D structured, and 2-D unstructured (a high-order least-squares reconstruction supplies the Laplacian on the FVM mesh). Requires SciPy: pip install ndiffusion[nearby].

solver = nd.FixedSourceSolver(mats, medium_map, edges_x, nd.Geometry.Slab, bc)
result = nd.nearby_fixed_source(solver, mats, source,
                                medium_map=medium_map, edges_x=edges_x,
                                geometry=nd.Geometry.Slab)
# result.error_estimate estimates (numerical - exact) flux

Standalone C++ driver

To build and run the 1-D reference problems without Python:

cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build
./build/cpp/ndiffusion_driver

Project structure

CMakeLists.txt              top-level CMake
pyproject.toml              build config (scikit-build-core)

cpp/
  CMakeLists.txt
  include/ndiffusion/
    types.hpp               shared types: Geometry, Materials, BoundaryCondition, results
    solver_1d.hpp           1-D solver class declarations
    solver_2d.hpp           2-D structured and unstructured solver declarations
    solver_3d.hpp           3-D solver declarations (placeholder)
  src/
    solver_1d.cpp               1-D solver implementation
    solver_2d_structured.cpp    structured 2-D implementation
    solver_2d_unstructured.cpp  unstructured 2-D implementation
    main.cpp                    standalone driver (1-D reference problems)
  python/
    bindings.cpp            pybind11 bindings -> ndiffusion._core

src/ndiffusion/
  __init__.py               re-exports from _core + create/mesh utilities
  create.py                 make_materials / make_medium_map / boundary_conditions
  mesh.py                   load_gmsh + assign_materials / copy_mesh
  layouts.py                preset layouts, symmetry orientations, matching BCs
  materials.py              named cross-section sets and published benchmarks
  transport.py              transport -> diffusion cross-section transform
  adjoint.py                make_adjoint_materials - forward -> adjoint transform
  kinetics.py               delayed neutron data + critical scaling helpers
  nearby.py                 method of nearby problems (fixed-source & k-eigenvalue)

tests/
  test_1d_k_eigenvalue.py       1-D k-eigenvalue, incl. non-uniform meshes
  test_1d_fixed_source.py       1-D fixed source, incl. iteration accounting
  test_1d_time_dependent.py
  test_2d_k_eigenvalue.py
  test_2d_fixed_source.py
  test_2d_time_dependent.py
  test_kinetics.py              delayed neutrons, point-kinetics limits
  test_benchmarks.py            published two-group benchmark regressions
  test_materials_validation.py  input validation: Materials, bc, mesh, dt
  test_transport.py             transport -> diffusion cross sections
  test_adjoint.py               adjoint materials transform
  test_nearby_1d.py             method of nearby problems (1-D)
  test_nearby_2d.py             method of nearby problems (2-D structured)
  test_nearby_unstructured.py   method of nearby problems (FVM)
  test_cg_ab.py                 within-group CG vs Gauss-Seidel A/B
  test_mesh_gmsh.py             Gmsh import (skipped without the gmsh package)
  test_mesh_materials.py        mesh geometry queries and material assignment
  test_layouts.py               preset layouts, orientations, periodic boundaries
  test_materials_sets.py        cross-section builders and benchmark bundles

examples/
  k_eigenvalue.py
  time_dependent.py
  kinetics.py
  transport_cross_sections.py
  c5g7_quarter_core.py

tools/
  c5g7_fuel_mesh.py         regenerates the gitignored C5G7 .msh meshes

Running tests

pip install -e ".[dev]"
pytest

The test suite lives in tests/ and is configured via pyproject.toml.

Install the extra. Optional dependencies are handled with pytest.importorskip, so a bare pip install -e . gives a run that looks clean while quietly skipping tests: without scipy that is the point-kinetics reference comparisons in tests/test_kinetics.py - the primary kinetics validation - plus all of tests/test_nearby_*.py. Run pytest -rs to see what skipped and why. (tests/test_mesh_gmsh.py needs pip install -e ".[mesh]", which is a heavier dependency and is not part of dev.)

The Testing/ directory is a generated CMake/CTest artifact and is not part of the source test suite.

API docs

Doxyfile configures Doxygen for the C++ sources under cpp/include/ndiffusion, cpp/src, and cpp/python. To generate the HTML documentation:

doxygen Doxyfile

Or, after configuring CMake:

cmake --build build --target docs

The output is written to docs/doxygen/html/.

Future work

See TODO.md for what is planned for 1.0 and after.