4.2.3. An Infinite-Medium Eigenvalue

This tutorial verifies a one-group \(k\)-eigenvalue calculation against an analytic infinite-medium result and checks its normalized particle balance.

4.2.3.1. Define an exactly critical material

The one-group material has \(\Sigma_t=\Sigma_a=1\), \(\Sigma_f=0.5\), and \(\nu=2\). With reflecting boundaries there is no net leakage, so

\[k_\infty=\frac{\nu\Sigma_f}{\Sigma_a}=1.\]

The analytic value is formed directly from the cross-section properties exposed by MultiGroupXS.

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
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("simple_fissile_1g.xs")
analytic_keff = float(xs.nu_sigma_f[0] / xs.sigma_a[0])
quadrature = GLProductQuadrature1DSlab(n_polar=8, scattering_order=0)
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=1,
    groupsets=[{"groups_from_to": (0, 0), "angular_quadrature": quadrature}],
    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 and balance

For a \(k\)-eigenvalue solve, ComputeBalanceTable divides fission production by the converged eigenvalue before forming the balance. The separate fission-rate and fission-production methods also make the material’s average neutron yield visible through their ratio.

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"]))
fission_rate = float(problem.ComputeFissionRate("new"))
fission_production = float(problem.ComputeFissionProduction("new"))
average_nu = fission_production / fission_rate

if rank == 0:
    print(f"Analytic infinite-medium eigenvalue={analytic_keff:.8e}")
    print(f"Computed infinite-medium eigenvalue={keff:.8e}")
    print(f"Eigenvalue balance residual={balance_residual:.8e}")
    print(f"Fission production-to-rate ratio={average_nu:.8e}")
assert abs(keff - analytic_keff) < 1.0e-10
assert balance_residual < 1.0e-10
assert abs(average_nu - 2.0) < 1.0e-12

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