4.1.1. Volume-Averaged and Integrated Flux

This tutorial compares two ways to integrate the same scalar-flux field: VolumePostprocessor, which operates directly on the transport problem, and FieldFunctionInterpolationVolume, which operates on an explicit scalar-flux field function.

4.1.1.1. Create a uniform reference solution

A uniform unit source in a purely absorbing material with reflecting boundaries produces a spatially uniform scalar flux of one. This makes the expected averages and integrals transparent.

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
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
length = 4.0
nodes = [length * i / 40.0 for i in range(41)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes]).Execute()
mesh.SetUniformBlockID(0)

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

4.1.1.2. Approach 1: use a volume postprocessor

A VolumePostprocessor reads the scalar flux directly from problem. The whole-domain average should be one and its integral should equal the slab length. Restricting the integral to \(1.0 \le z \le 2.5\) evaluates \(\int_{1.0}^{2.5} \phi(z)\,dz\), which should give 1.5.

whole_average = VolumePostprocessor(problem=problem, value_type="avg")
whole_integral = VolumePostprocessor(problem=problem, value_type="integral")
sample_region = RPPLogicalVolume(infx=True, infy=True, zmin=1.0, zmax=2.5)
regional_integral = VolumePostprocessor(
    problem=problem, value_type="integral", logical_volumes=[sample_region]
)
for postprocessor in (whole_average, whole_integral, regional_integral):
    postprocessor.Execute()

average = float(whole_average.GetValue()[0][0])
integral = float(whole_integral.GetValue()[0][0])
regional = float(regional_integral.GetValue()[0][0])
max_error = max(abs(average - 1.0), abs(integral - 4.0), abs(regional - 1.5))
if rank == 0:
    print(f"Whole-domain average={average:.6e}")
    print(f"Whole-domain integral={integral:.6e}")
    print(f"Regional integral={regional:.6e}")
    print(f"Volume-postprocessor max error={max_error:.6e}")
assert max_error < 1.0e-6

4.1.1.3. Approach 2: integrate a field function

The same regional integral can be evaluated from an explicit scalar-flux field function. For FieldFunctionInterpolationVolume, the "sum" operation integrates the field over the selected logical volume. Reusing sample_region therefore evaluates the same \(\int_{1.0}^{2.5} \phi(z)\,dz\) as regional_integral above.

scalar_flux = problem.GetScalarFluxFieldFunction()[0]
field_function_integral = FieldFunctionInterpolationVolume()
field_function_integral.SetOperationType("sum")
field_function_integral.SetLogicalVolume(sample_region)
field_function_integral.AddFieldFunction(scalar_flux)
field_function_integral.Execute()
field_integral = float(field_function_integral.GetValue())
field_function_error = max(abs(field_integral - 1.5), abs(field_integral - regional))
if rank == 0:
    print(f"Field-function regional integral={field_integral:.6e}")
    print(f"Field-function agreement error={field_function_error:.6e}")
assert field_function_error < 1.0e-6
if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()