4.2.3. A Multigroup Infinite-Medium Eigenvalue

This tutorial verifies a four-group \(k\)-eigenvalue calculation against the exact infinite-medium solution. It checks the eigenvalue, the converged energy spectrum, the spectrum-averaged neutron yield, and the normalized particle balance.

4.2.3.1. The infinite-medium reference

In an infinite homogeneous medium with isotropic scattering, the flux is uniform in space and isotropic, so the multigroup transport equation reduces to a balance for each group:

\[\Sigma_{t,g}\phi_g=\sum_{g'}\Sigma_s(g'\to g)\,\phi_{g'}+\frac{\chi_g}{k}\sum_{g'}\nu\Sigma_{f,g'}\phi_{g'}.\]

In matrix form this is \(A\phi=k^{-1}\chi\,(\nu\Sigma_f)^{T}\phi\) with \(A=\operatorname{diag}(\Sigma_t)-\Sigma_s^{T}\), where \(\Sigma_s[g',g]=\Sigma_s(g'\to g)\). A single fission spectrum \(\chi\) makes the production operator rank one, so the fundamental mode has the closed form

\[\phi_\infty\propto A^{-1}\chi,\qquad k_\infty=(\nu\Sigma_f)^{T}A^{-1}\chi.\]

Here \(A^{-1}\chi\) is the group flux produced by one fission neutron, and \(k_\infty\) is the number of fission neutrons that flux produces. With one group, the result reduces to \(k_\infty=\nu\Sigma_f/\Sigma_a\).

The four-group material below includes fast fission, downscatter across several groups, upscatter between the two thermal groups, and a group-dependent \(\nu\):

Group	\(\Sigma_t\)	\(\Sigma_s(g\to 0,1,2,3)\)	\(\Sigma_a\)	\(\Sigma_f\)	\(\nu\)	\(\chi\)
0 (fast)	0.20	0.150, 0.040, 0.002, 0	0.008	0.003	2.70	0.95
1 (epithermal)	0.35	0, 0.300, 0.030, 0.005	0.015	0.006	2.50	0.05
2 (thermal)	0.90	0, 0, 0.650, 0.150	0.100	0.060	2.44	0
3 (cold thermal)	1.20	0, 0, 0.060, 0.900	0.240	0.100	2.43	0

Cross sections are in \(\text{cm}^{-1}\), and \(\Sigma_a=\Sigma_t-\sum_{g'}\Sigma_s(g\to g')\). A slab with reflecting boundaries on both faces reproduces the infinite medium exactly, so neither the mesh nor the angular quadrature introduces discretization error, and OpenSn should match the closed form to the solver tolerances.

The cell below evaluates the reference directly from these values with NumPy, independently of OpenSn.

import numpy as np

sigma_t = np.array([0.20, 0.35, 0.90, 1.20])
# sigma_s[g_from, g_to]
sigma_s = np.array(
    [
        [0.150, 0.040, 0.002, 0.000],
        [0.000, 0.300, 0.030, 0.005],
        [0.000, 0.000, 0.650, 0.150],
        [0.000, 0.000, 0.060, 0.900],
    ]
)
sigma_f = np.array([0.003, 0.006, 0.060, 0.100])
nu = np.array([2.70, 2.50, 2.44, 2.43])
chi = np.array([0.95, 0.05, 0.0, 0.0])
sigma_a = sigma_t - sigma_s.sum(axis=1)

A = np.diag(sigma_t) - sigma_s.T
phi_ref = np.linalg.solve(A, chi)
analytic_keff = float((nu * sigma_f) @ phi_ref)
spectrum_ref = phi_ref / phi_ref.sum()
analytic_average_nu = float((nu * sigma_f) @ phi_ref / (sigma_f @ phi_ref))

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.post import VolumePostprocessor
from pyopensn.solver import DiscreteOrdinatesProblem, PowerIterationKEigenSolver
from pyopensn.xs import MultiGroupXS

rank = MPI.COMM_WORLD.rank
nodes = [i / 10.0 for i in range(11)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes]).Execute()
mesh.SetUniformBlockID(0)

xs = MultiGroupXS()
xs.LoadFromOpenSn("infinite_medium_4g.xs")

# The file and the reference arrays must describe the same material.
for name, expected in [
    ("sigma_t", sigma_t),
    ("sigma_a", sigma_a),
    ("sigma_f", sigma_f),
    ("nu_sigma_f", nu * sigma_f),
    ("chi", chi),
]:
    assert np.allclose(getattr(xs, name), expected, rtol=0.0, atol=1.0e-12), name

quadrature = GLProductQuadrature1DSlab(n_polar=8, scattering_order=0)
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=4,
    groupsets=[
        {
            "groups_from_to": (0, 3),
            "angular_quadrature": quadrature,
            "inner_linear_method": "petsc_gmres",
            "l_abs_tol": 1.0e-12,
            "l_max_its": 500,
        }
    ],
    xs_map=[{"block_ids": [0], "xs": xs}],
    boundary_conditions=[
        {"name": "zmin", "type": "reflecting"},
        {"name": "zmax", "type": "reflecting"},
    ],
    options={
        "use_precursors": False,
        "verbose_inner_iterations": False,
        "verbose_outer_iterations": False,
    },
)

4.2.3.2. Verify the eigenvalue, spectrum, and balance

The solve converges tightly (k_tol=1e-12 and an inner GMRES tolerance of \(10^{-12}\)). Each check below requires agreement to \(10^{-10}\), a margin above the solver tolerances; by comparison, changing a single transfer entry by \(0.001\) \(\text{cm}^{-1}\) shifts \(k\) by about \(2\times10^{-3}\). The checks are:

• the computed eigenvalue against \(k_\infty\);

• the group spectrum from a VolumePostprocessor, normalized to unit sum, against the normalized \(A^{-1}\chi\);

• the spectrum-averaged neutron yield \(\bar\nu=\sum_g\nu\Sigma_{f,g}\phi_g/\sum_g\Sigma_{f,g}\phi_g\), formed from ComputeFissionProduction and ComputeFissionRate, which tests the energy weighting of the fission source;

• the normalized balance from ComputeBalanceTable, which divides fission production by the converged eigenvalue before forming the balance.

solver = PowerIterationKEigenSolver(
    problem=problem, k_tol=1.0e-12, compute_balance=True
)
solver.Initialize()
solver.Execute()

keff = float(solver.GetEigenvalue())
balance = solver.ComputeBalanceTable()
balance_residual = abs(float(balance["balance"]))

group_average = VolumePostprocessor(problem=problem, value_type="avg")
group_average.Execute()
group_flux = np.asarray(group_average.GetValue()[0])
spectrum = group_flux / group_flux.sum()
spectrum_error = float(np.max(np.abs(spectrum - spectrum_ref) / spectrum_ref))

fission_rate = float(problem.ComputeFissionRate("new"))
fission_production = float(problem.ComputeFissionProduction("new"))
average_nu = fission_production / fission_rate
average_nu_error = abs(average_nu - analytic_average_nu) / analytic_average_nu

if rank == 0:
    print("Group  Computed spectrum  Analytic spectrum")
    for g in range(len(spectrum)):
        print(f"{g:5d}  {spectrum[g]:.12f}     {spectrum_ref[g]:.12f}")
    print(f"Analytic infinite-medium eigenvalue={analytic_keff:.12e}")
    print(f"Computed infinite-medium eigenvalue={keff:.12e}")
    print(f"Spectrum max relative error={spectrum_error:.8e}")
    print(f"Average nu relative error={average_nu_error:.8e}")
    print(f"Eigenvalue balance residual={balance_residual:.8e}")
assert abs(keff - analytic_keff) < 1.0e-10
assert spectrum_error < 1.0e-10
assert average_nu_error < 1.0e-10
assert balance_residual < 1.0e-10

if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()