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 - backward-Euler time stepping, unconditionally stable
• Per-group Thomas (TDMA) tridiagonal solver inside a Gauss-Seidel group sweep
• Harmonic-mean diffusion coefficients at material interfaces

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 - backward-Euler stepping using the same line-TDMA sweep

2-D unstructured (triangles and/or quadrilaterals)

• Cell-centred finite-volume method (FVM)
• Arbitrary Robin BCs per boundary tag; harmonic-mean interface diffusion coefficients
• k-eigenvalue solver - power iteration with point Gauss-Seidel inner solve
• Fixed-source solver - point SOR (successive over-relaxation) inner solve
• Time-dependent solver - backward-Euler 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 .

Python edits under src/ndiffusion/ are picked up immediately; after editing C++ sources, re-run pip install -e . 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.

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).

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 - Gmsh .msh import for unstructured meshes

tests/
  test_1d_k_eigenvalue.py
  test_1d_time_dependent.py
  test_1d_fixed_source.py
  test_2d_k_eigenvalue.py
  test_2d_time_dependent.py
  test_2d_fixed_source.py

examples/
  k_eigenvalue.py
  time_dependent.py

Running tests

pytest

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

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

Geometry

• 3-D structured geometry (x-y-z) and 3-D unstructured (tetrahedra/hexahedra)
• General boundary conditions on all edges (1-D currently hardcodes symmetry at the left/inner edge; 2-D structured hardcodes left and bottom as reflective)
• Non-orthogonal correction for the unstructured FVM two-point flux approximation (accuracy degrades on skewed meshes)

Physics

• Delayed neutron precursor groups in the time-dependent solver (currently prompt-only), and implicit treatment of the fission source (currently explicit)
• Adjoint flux solver for sensitivity and perturbation analysis
• Depletion coupling - Bateman equations for nuclide inventory evolution

Solvers and performance

• Flip the default inner solver for the 2-D k-eigenvalue solvers to the within-group CG (now a use_cg constructor option; default remains Gauss-Seidel, overridable via NDIFFUSION_KEIG_CG=1); extend CG to the fixed-source and time-dependent solvers, replacing hand-tuned SOR
• Power-iteration acceleration (Wielandt shift or Chebyshev extrapolation); CMFD (Coarse Mesh Finite Difference) for unstructured k-eigenvalue convergence
• Zero-copy numpy arrays across the pybind11 boundary (fluxes and sources currently cross as Python lists)
• OpenMP parallelism for the spatial sweep loops

Testing

• Published two-group benchmark regressions are in tests/test_benchmarks.py (1-D Ringhals-4 slab, 2-D TWIGL, 2-D IAEA PWR on the stepped quarter core); still to add: BIBLIS-2D and C5G7 diffusion via tools/c5g7_fuel_mesh.py