4.2.1. Global Particle Balance

This tutorial checks that volumetric production and boundary inflow balance absorption and outward leakage in a fixed-source calculation.

4.2.1.1. Solve with balance accounting enabled

The slab contains a uniform source and an absorbing, scattering material with vacuum boundaries. Setting compute_balance=True instructs the steady-state solver to retain the terms needed by ComputeBalanceTable.

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
from pyopensn.source import VolumetricSource
from pyopensn.xs import MultiGroupXS

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

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.6)
source = VolumetricSource(block_ids=[0], group_strength=[1.0])
quadrature = GLProductQuadrature1DSlab(n_polar=64, 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-10,
        }
    ],
    xs_map=[{"block_ids": [0], "xs": xs}],
    volumetric_sources=[source],
    boundary_conditions=[
        {"name": "zmin", "type": "vacuum"},
        {"name": "zmax", "type": "vacuum"},
    ],
)
solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

4.2.1.2. Form the balance residual

For a converged steady-state solve, production plus inflow must equal absorption plus outflow. A normalized residual makes the check useful across problems with different source magnitudes.

balance = solver.ComputeBalanceTable()
source_rate = balance["production_rate"] + balance["inflow_rate"]
loss_rate = balance["absorption_rate"] + balance["outflow_rate"]
residual = abs(source_rate - loss_rate) / max(abs(source_rate), 1.0e-16)
if rank == 0:
    print(f"Production={balance['production_rate']:.6e}")
    print(f"Absorption={balance['absorption_rate']:.6e}")
    print(f"Outflow={balance['outflow_rate']:.6e}")
    print(f"Balance residual={residual:.6e}")
assert residual < 1.0e-8
if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()