4.1.5. Leakage Rate Through a Boundary

This tutorial computes the net leakage through 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. Leakage is evaluated from the boundary currents that every sweep records, so save_angular_flux is not required.

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")
    ],
)
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()