5.2.3. Inner-Solver Convergence with Scattering

This tutorial compares how OpenSn’s four inner transport methods respond as scattering becomes stronger.

5.2.3.1. Scattering-ratio study

Each calculation uses the same homogeneous reflecting slab, unit source, 400-cell mesh, and 16-direction quadrature. Only the scattering ratio and inner method change. For \(\Sigma_t=1\), the analytic scalar flux is \(1/(1-c)\).

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, SteadyStateSourceSolver
from pyopensn.source import VolumetricSource
from pyopensn.xs import MultiGroupXS

comm = MPI.COMM_WORLD
rank = comm.rank
methods = {
    "classic_richardson": "Classic Richardson",
    "petsc_richardson": "PETSc Richardson",
    "petsc_gmres": "PETSc GMRES",
    "petsc_bicgstab": "PETSc BiCGStab",
}
scattering_ratios = [0.0, 0.5, 0.8, 0.95]


def solve(method, scattering_ratio):
    mesh = OrthogonalMeshGenerator(
        node_sets=[[i / 400.0 for i in range(401)]]
    ).Execute()
    mesh.SetUniformBlockID(0)
    xs = MultiGroupXS()
    xs.CreateSimpleOneGroup(sigma_t=1.0, c=scattering_ratio)
    quadrature = GLProductQuadrature1DSlab(
        n_polar=16, scattering_order=0
    )
    groupset = {
        "groups_from_to": (0, 0),
        "angular_quadrature": quadrature,
        "inner_linear_method": method,
        "l_abs_tol": 1.0e-8,
        "l_max_its": 1500,
    }
    if method == "petsc_gmres":
        groupset["gmres_restart_interval"] = 30

    problem = DiscreteOrdinatesProblem(
        mesh=mesh,
        num_groups=1,
        groupsets=[groupset],
        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={"verbose_inner_iterations": False},
    )

    solver = SteadyStateSourceSolver(problem=problem)
    comm.Barrier()
    start = MPI.Wtime()
    solver.Initialize()
    solver.Execute()
    comm.Barrier()
    elapsed = comm.allreduce(MPI.Wtime() - start, op=MPI.MAX)

    average = VolumePostprocessor(problem=problem, value_type="avg")
    average.Execute()
    flux = sum(average.GetValue()[0])
    return flux, solver.GetNumSweeps(), elapsed

5.2.3.2. Compare accuracy and work

The analytic flux checks correctness. Sweep counts show how scattering changes transport work, while wall time provides a machine-dependent performance measure.

results = {
    (method, ratio): solve(method, ratio)
    for ratio in scattering_ratios
    for method in methods
}
relative_errors = {}
for ratio in scattering_ratios:
    analytic_flux = 1.0 / (1.0 - ratio)
    for method in methods:
        flux = results[method, ratio][0]
        relative_errors[method, ratio] = abs(flux - analytic_flux) / analytic_flux
maximum_relative_error = max(relative_errors.values())

if rank == 0:
    for ratio in scattering_ratios:
        for method, label in methods.items():
            flux, sweeps, elapsed = results[method, ratio]
            error = relative_errors[method, ratio]
            print(f"{label} c={ratio:.2f} total average flux={flux:.12e}")
            print(f"{label} c={ratio:.2f} relative flux error={error:.12e}")
            print(f"{label} c={ratio:.2f} sweeps={sweeps}")
            print(f"{label} c={ratio:.2f} wall time (s)={elapsed:.6f}")
    print(f"Maximum scattering-study relative flux error={maximum_relative_error:.12e}")

assert maximum_relative_error < 1.0e-6

A representative one-process run gives sweep count / wall time in seconds:

Method	\(c=0.00\)	\(c=0.50\)	\(c=0.80\)	\(c=0.95\)
Classic Richardson	11 / 0.0135	37 / 0.0320	107 / 0.1670	453 / 0.5483
PETSc Richardson	11 / 0.0170	38 / 0.0280	108 / 0.1384	454 / 0.5140
PETSc GMRES	8 / 0.0065	9 / 0.0092	10 / 0.0139	11 / 0.0087
PETSc BiCGStab	10 / 0.0072	12 / 0.0161	14 / 0.0101	14 / 0.0102

The maximum relative flux error is \(9.24\times10^{-9}\). Richardson slows sharply as \(c\) approaches one, while GMRES and BiCGStab remain effective. Exact timings depend on the machine.

5.2.3.3. Finalize (for Jupyter Notebook only)

In script mode, PyOpenSn handles finalization automatically. In a Jupyter kernel, finalize OpenSn before MPI.

if "opensn_console" not in globals():
    from IPython import get_ipython

    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()