4.1.5. Leakage Rate Through a Boundary

This tutorial integrates the outgoing angular flux over a named boundary with ComputeLeakage.

4.1.5.1. Solve a two-dimensional fixed-source problem

A uniform source is placed in a homogeneous unit square with vacuum boundaries. Because leakage is evaluated from the outgoing angular flux, save_angular_flux must be enabled on the transport problem.

from mpi4py import MPI
from pyopensn.aquad import GLCProductQuadrature2DXY
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, nodes]).Execute()
mesh.SetUniformBlockID(0)

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=0.5, c=0.2)
quadrature = GLCProductQuadrature2DXY(n_polar=2, n_azimuthal=16, 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}],
    volumetric_sources=[VolumetricSource(block_ids=[0], group_strength=[1.0])],
    boundary_conditions=[
        {"name": name, "type": "vacuum"}
        for name in ("xmin", "xmax", "ymin", "ymax")
    ],
    options={"save_angular_flux": True},
)
solver = SteadyStateSourceSolver(problem=problem)
solver.Initialize()
solver.Execute()

4.1.5.2. Integrate the outgoing current

ComputeLeakage returns one value per energy group for each requested boundary. Here, we report the leakage through the xmax boundary and verify it using the global particle balance.

boundary_names = ["xmin", "xmax", "ymin", "ymax"]
leakage = problem.ComputeLeakage(boundary_names)
xmax_leakage = float(leakage["xmax"][0])
total_leakage = sum(float(leakage[name][0]) for name in boundary_names)

absorption = VolumePostprocessor(
    problem=problem, value_type="integral", xs_multiplier="sigma_a"
)
absorption.Execute()
absorption_rate = float(absorption.GetValue()[0][0])
balance_residual = 1.0 - absorption_rate - total_leakage

if rank == 0:
    print(f"XMAX leakage rate={xmax_leakage:.8e}")
    print(f"Total leakage rate={total_leakage:.8e}")
    print(f"Particle balance residual={balance_residual:.8e}")
assert abs(balance_residual) < 1.0e-6

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