3.1.7. A Striped Fixed-Source Problem

This two-dimensional problem places three source-bearing material stripes in a scattering background. Two reflecting boundaries reproduce neighboring copies of the modeled quadrant, while particles can escape through the top and right vacuum boundaries.

3.1.7.1. Import the OpenSn objects

if "opensn_console" not in globals():
    from mpi4py import MPI
    from pyopensn.aquad import GLCProductQuadrature2DXY
    from pyopensn.context import Finalize, UseColor
    from pyopensn.logvol import RPPLogicalVolume
    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
    UseColor(False)

3.1.7.2. Create the mesh

The domain is \([0,10]\times[0,10]\) cm. A uniform spacing of 0.1 cm places mesh nodes on every stripe edge and on the horizontal interface at \(y=9\) cm. The two-dimensional calculation represents a unit thickness in the \(z\) direction.

domain_length = 10.0
cells_per_cm = 10
num_cells = int(domain_length * cells_per_cm)
nodes = [i / cells_per_cm for i in range(num_cells + 1)]

mesh = OrthogonalMeshGenerator(node_sets=[nodes, nodes]).Execute()
mesh.SetOrthogonalBoundaries()

3.1.7.3. Assign the material stripes

Material M fills the domain. Three F regions occupy \(1\leq x\leq2\), \(4\leq x\leq5\), and \(7\leq x\leq8\) cm, each extending from \(y=0\) to \(y=9\) cm.

The F regions share one block ID, allowing a single material definition and source to be applied to all three disconnected stripes.

MATERIAL_M = 0
MATERIAL_F = 1
stripe_intervals = [(1.0, 2.0), (4.0, 5.0), (7.0, 8.0)]

mesh.SetUniformBlockID(MATERIAL_M)

for xmin, xmax in stripe_intervals:
    stripe = RPPLogicalVolume(
        xmin=xmin, xmax=xmax,
        ymin=0.0, ymax=9.0,
        infz=True,
    )
    mesh.SetBlockIDFromLogicalVolume(stripe, MATERIAL_F, True)

3.1.7.4. Define the materials and source

This tutorial uses the total and isotropic-scattering cross sections from the reference material table. The tabulated fission terms for material F are deliberately set to zero, and a unit external source is applied to the F stripes instead.

Material	\(\sigma_t\) (cm\(^{-1}\))	\(\sigma_{s0}\) (cm\(^{-1}\))	\(c=\sigma_{s0}/\sigma_t\)	External source (cm\(^{-3}\) s\(^{-1}\))
M	1.0	0.93	0.93	0.0
F	1.5	1.35	0.90	1.0

Because CreateSimpleOneGroup creates non-fissile one-group data, no fission source is present in either material.

xs_m = MultiGroupXS()
xs_m.CreateSimpleOneGroup(sigma_t=1.0, c=0.93)

xs_f = MultiGroupXS()
xs_f.CreateSimpleOneGroup(sigma_t=1.5, c=1.35 / 1.5)

source = VolumetricSource(
    block_ids=[MATERIAL_F], group_strength=[1.0]
)

3.1.7.5. Assemble and solve the transport problem

The xmin and ymin faces are reflecting, while xmax and ymax are vacuum. A Gauss–Legendre–Chebyshev quadrature represents directions in the \(xy\) plane.

quadrature = GLCProductQuadrature2DXY(
    n_polar=2, n_azimuthal=32, scattering_order=0
)

problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=1,
    groupsets=[
        {
            "groups_from_to": (0, 0),
            "angular_quadrature": quadrature,
            "inner_linear_method": "petsc_gmres",
            "l_abs_tol": 1.0e-9,
            "l_max_its": 300,
            "gmres_restart_interval": 30,
        }
    ],
    xs_map=[
        {"block_ids": [MATERIAL_M], "xs": xs_m},
        {"block_ids": [MATERIAL_F], "xs": xs_f},
    ],
    volumetric_sources=[source],
    boundary_conditions=[
        {"name": "xmin", "type": "reflecting"},
        {"name": "xmax", "type": "vacuum"},
        {"name": "ymin", "type": "reflecting"},
        {"name": "ymax", "type": "vacuum"},
    ],
)

solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

3.1.7.6. Check the source and particle balance

Each F stripe has area \(1\times9=9\) cm\(^2\). The three unit-strength source regions therefore produce 27 particles per second per unit depth. We also compute the average scalar flux and absorption rate over the combined F-region block ID.

production_rate measures the external volumetric source. The balance table also reports particles leaving and reentering through reflecting faces as outflow and inflow, respectively, so the complete balance includes both boundary terms.

f_flux_postprocessor = VolumePostprocessor(
    problem=problem, value_type="avg", block_ids=[MATERIAL_F]
)
f_flux_postprocessor.Execute()
f_region_average_flux = f_flux_postprocessor.GetValue()[0][0]

f_absorption_postprocessor = VolumePostprocessor(
    problem=problem,
    value_type="integral",
    block_ids=[MATERIAL_F],
    xs_multiplier="sigma_a",
)
f_absorption_postprocessor.Execute()
f_region_absorption_rate = f_absorption_postprocessor.GetValue()[0][0]

balance = solver.ComputeBalanceTable()
external_source_rate = balance["production_rate"]
total_gain_rate = external_source_rate + balance["inflow_rate"]
loss_rate = balance["absorption_rate"] + balance["outflow_rate"]
relative_balance_error = abs(total_gain_rate - loss_rate) / total_gain_rate

if rank == 0:
    print(f"STRIPED_SOURCE_RATE={external_source_rate:.8e}")
    print(f"STRIPED_SOURCE_AVG_FLUX={f_region_average_flux:.8e}")
    print(f"STRIPED_SOURCE_ABSORPTION_RATE={f_region_absorption_rate:.8e}")
    print(f"STRIPED_SOURCE_BALANCE_ERROR={relative_balance_error:.8e}")

assert abs(external_source_rate - 27.0) < 1.0e-8
assert relative_balance_error < 1.0e-8

3.1.7.7. Export and visualize the scalar flux

After running the generated Python input, use the following code to export the scalar flux for visualization. It is kept in Markdown so the regression test does not create output files.

from pyopensn.fieldfunc import FieldFunctionGridBased

scalar_flux = problem.GetScalarFluxFieldFunction()[0]
FieldFunctionGridBased.ExportMultipleToPVTU(
    [scalar_flux], "Flux/StripedSource_Phi"
)

3.1.7.8. Next steps

Compare the flux in the three source stripes, replace the external source with the original fission data to recover an eigenvalue problem, or change the reflecting faces to vacuum and examine the resulting leakage.

if "opensn_console" not in globals():
    from IPython import get_ipython

    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()