2.5.1. Multiple Groupsets in an Infinite-Medium

This tutorial divides a realistic 30-group concrete calculation into three groupsets and shows how OpenSn resolves scattering between energy ranges.

2.5.1.1. Define the groupsets

The problem uses the LANL concrete mixture represented with the 30-group LANL energy structure (LANL30). Its cross sections, including the scattering transfer matrix, are imported from OpenMC.

The groupsets cover groups 0–9, 10–19, and 20–29. Each groupset has its own inner solve, while across-groupset iteration (AGS) reconciles scattering contributions that cross groupset boundaries.

A unit source is applied to group 0. Reflecting boundaries make the slab an infinite-medium model, keeping the example focused on energy coupling rather than spatial leakage.

Data note: When importing OpenMC data, OpenSn reconstructs the absorption cross section from the total cross section and transfer matrix. Small statistical differences in this concrete dataset produce a negative reconstructed absorption cross section for group 13, causing OpenSn to issue a warning. This tutorial uses the supplied data without modification.

from pathlib import Path

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, SteadyStateSourceSolver
from pyopensn.source import VolumetricSource
from pyopensn.xs import MultiGroupXS

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

xs = MultiGroupXS()
xs.LoadFromOpenMC(str(Path("LANL30/OpenMC/Concrete.h5")), "set1", 294.0)
num_groups = xs.num_groups

source_strength = [0.0] * num_groups
source_strength[0] = 1.0
source = VolumetricSource(block_ids=[0], group_strength=source_strength)
quadrature = GLProductQuadrature1DSlab(n_polar=16, scattering_order=3)

groupset_ranges = [(0, 9), (10, 19), (20, 29)]
groupsets = [
    {
        "groups_from_to": group_range,
        "angular_quadrature": quadrature,
        "inner_linear_method": "petsc_gmres",
        "l_abs_tol": 1.0e-10,
        "l_max_its": 100,
    }
    for group_range in groupset_ranges
]

problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=num_groups,
    groupsets=groupsets,
    xs_map=[{"block_ids": [0], "xs": xs}],
    volumetric_sources=[source],
    boundary_conditions=[
        {"name": "zmin", "type": "reflecting"},
        {"name": "zmax", "type": "reflecting"},
    ],
    options={
        "max_ags_iterations": 100,
        "ags_tolerance": 1.0e-8,
        "verbose_inner_iterations": False,
    },
)
solver = SteadyStateSourceSolver(problem=problem)
solver.Initialize()
solver.Execute()

2.5.1.2. Verify the coupled solution

Scattering redistributes the group-0 source across the concrete spectrum. To summarize that spectrum, the code reports the flux integrated over each groupset. Because the boundaries are reflecting, no particles leak from the slab; at convergence, the absorption rate must therefore equal the unit source rate.

group_averages = []
for group in range(num_groups):
    postprocessor = VolumePostprocessor(problem=problem, value_type="avg", group=group)
    postprocessor.Execute()
    group_averages.append(float(postprocessor.GetValue()[0][0]))

groupset_fluxes = [
    sum(group_averages[first : last + 1])
    for first, last in groupset_ranges
]
absorption_rate = sum(
    xs.sigma_a[group] * group_averages[group]
    for group in range(num_groups)
)
balance_error = abs(absorption_rate - 1.0)
if rank == 0:
    for (first, last), flux in zip(groupset_ranges, groupset_fluxes):
        print(f"Groupset {first}-{last} integrated flux={flux:.8e}")
    print(f"Concrete absorption rate={absorption_rate:.6e}")
    print(f"Groupset balance error={balance_error:.6e}")
assert min(group_averages) > -1.0e-10
assert balance_error < 1.0e-6
if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()