4.2.2. Transient Particle Balance

This tutorial checks a time-dependent calculation against both an analytic solution and the particle inventory change over one timestep.

4.2.2.1. Define an infinite-medium reference problem

A homogeneous, purely absorbing slab has reflecting boundaries, unit particle speed, \(\Sigma_a=1\), and a constant unit source. Starting from zero flux, the spatially uniform solution satisfies

\[\frac{d\phi}{dt}+\phi=1,\]

with the analytic solution \(\phi(t)=1-e^{-t}\). Time-dependent discrete ordinates problems must retain angular flux, so save_angular_flux=True.

import math
from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.post import VolumePostprocessor
from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver
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]).Execute()
mesh.SetUniformBlockID(0)

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0, velocity=1.0)
quadrature = GLProductQuadrature1DSlab(n_polar=8, scattering_order=0)
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=1,
    time_dependent=True,
    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": "zmin", "type": "reflecting"},
        {"name": "zmax", "type": "reflecting"},
    ],
    options={
        "save_angular_flux": True,
        "use_precursors": False,
        "verbose_inner_iterations": False,
    },
)

4.2.2.2. Compare the final state and inventory change

Backward Euler uses end-of-step rates, so its discrete balance for the final step is

\[N^{n+1}-N^n=\Delta t\left(P-A\right)^{n+1}.\]

The reflecting boundaries have zero net leakage. For this unit-volume, unit-speed problem, the average scalar flux equals the particle inventory and the integrated source rate is one. We save the inventory immediately before the final advance and independently integrate the final absorption rate with VolumePostprocessor.

dt = 0.025
stop_time = 1.0
num_steps = round(stop_time / dt)
solver = TransientSolver(
    problem=problem,
    initial_state="zero",
    dt=dt,
    theta=1.0,
    verbose=False,
)
solver.Initialize()

average_flux = VolumePostprocessor(problem=problem, value_type="avg")
for _ in range(num_steps - 1):
    solver.Advance()

average_flux.Execute()
previous_inventory = float(average_flux.GetValue()[0][0])
solver.Advance()
average_flux.Execute()
final_flux = float(average_flux.GetValue()[0][0])

absorption = VolumePostprocessor(
    problem=problem, value_type="integral", xs_multiplier="sigma_a"
)
absorption.Execute()
absorption_rate = float(absorption.GetValue()[0][0])
source_rate = 1.0
actual_inventory_change = final_flux - previous_inventory
predicted_inventory_change = dt * (source_rate - absorption_rate)
inventory_residual = actual_inventory_change - predicted_inventory_change
analytic_flux = 1.0 - math.exp(-stop_time)
relative_error = abs(final_flux - analytic_flux) / analytic_flux

if rank == 0:
    print(f"Transient scalar flux={final_flux:.8e}")
    print(f"Transient relative error={relative_error:.8e}")
    print(f"Transient inventory residual={inventory_residual:.8e}")
assert relative_error < 1.0e-2
assert abs(inventory_residual) < 1.0e-12

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