3.8. Delta-Forward Scattering and Galerkin Quadrature Methods
This tutorial introduces the delta-forward scattering kernel — the limiting case of forward-peaked scattering where all particles scatter exactly in their incident direction — and uses it to compare the standard, Galerkin One, and Galerkin Three operator construction methods in 3D.
The delta-forward kernel is defined by a Dirac delta in the scattering angle cosine:
Its Legendre moments are all identically one:
so the scattering cross-section moments are \(\sigma_{s,\ell} = \sigma_s\) for every \(\ell\). No finite Legendre truncation can exactly represent this kernel.
Key analytical result. Substituting the delta kernel into the transport equation yields
Delta-forward scattering is therefore exactly equivalent to pure absorption with \(\sigma_a = \sigma_t - \sigma_s\), which provides an exact reference solution.
Why 3D matters. In 1D slab geometry with Gauss-Legendre quadrature, the GL rule exactly integrates products \(P_\ell(\mu)\,P_m(\mu)\), so the standard operator already satisfies \(\mathbf{D}\mathbf{M}=\mathbf{I}\) and all three methods are equivalent. In 3D, the GL–Chebyshev product quadrature uses equally-spaced azimuthal angles that do not satisfy exact spherical harmonic integration, so the standard \(\mathbf{DM}\) product deviates measurably from the identity. Galerkin methods correct this and produce a genuinely different — and more accurate — transport solution.
To run the code: jupyter nbconvert --to python --execute delta_forward.ipynb.
To convert to a Python script: jupyter nbconvert --to python delta_forward.ipynb
[ ]:
import os
import sys
import numpy as np
from mpi4py import MPI
sys.path.append("../../../..")
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.xs import MultiGroupXS
from pyopensn.source import VolumetricSource
from pyopensn.aquad import GLCProductQuadrature3DXYZ
from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
from pyopensn.logvol import RPPLogicalVolume
from pyopensn.context import UseColor, Finalize
UseColor(False)
rank = MPI.COMM_WORLD.rank
3.8.1. Problem Setup
Parameter |
Value |
|---|---|
Geometry |
3D cube, \([0,\,10]^3\) mfp |
Mesh |
\(20 \times 20 \times 20 = 8000\) cells |
\(\sigma_t\) |
1.5 |
\(\sigma_s\) |
1.0 |
Scattering |
Delta-forward: \(\sigma_{s,\ell} = 1.0\) for all \(\ell\) |
\(\sigma_a = \sigma_t - \sigma_s\) |
0.5 |
Volumetric source |
Isotropic, \(Q = 1.0\) |
Boundary conditions |
Isotropic incoming flux \(= Q/\sigma_a = 2.0\) on all six faces |
Angular quadrature |
GL-Chebyshev, \(N_p = 4\) polar, \(N_\phi = 8\) azimuthal (64 directions) |
Analytical solution |
Uniform scalar flux \(\phi^* = Q/\sigma_a = 2.0\) |
Analytical solution. With delta-forward scattering reducing to pure absorption with cross-section \(\sigma_a = \sigma_t - \sigma_s\), the steady-state balance in an infinite medium is simply
By setting every boundary face to an isotropic incoming flux equal to \(\phi^*\), the spatially uniform solution \(\phi(\mathbf{r}) = \phi^*\) satisfies the transport equation everywhere in the domain. Any deviation of the numerical solution from this uniform value directly quantifies operator error in the scattering term.
We compare three operator methods: standard, galerkin_one, and galerkin_three at \(L = 7\) (64 moments).
[ ]:
sigma_t = 1.5
sigma_s = 1.0
sigma_a = sigma_t - sigma_s # = 0.5
Q = 1.0
analytical_phi = Q / sigma_a # = 2.0
bsrc = analytical_phi # incoming isotropic flux on every face
L_compare = 7 # scattering order for standard and Galerkin Three
L_max = 7 # highest moment stored in the XS file
L_domain = 10.0
N_mesh = 20 # cells per dimension (N_mesh^3 cells total)
N_polar = 4
N_azimuthal = 8
n_dirs = 2 * N_polar * N_azimuthal
L_auto = int(np.sqrt(n_dirs) - 1)
print(f"sigma_t = {sigma_t}")
print(f"sigma_s = {sigma_s} -> sigma_a = {sigma_a}")
print(f"Q = {Q} -> phi_analytical = {analytical_phi}")
print(f"Mesh = {N_mesh}^3 = {N_mesh**3} cells")
print(f"Directions = {n_dirs} (2 x {N_polar} polar x {N_azimuthal} azimuthal)")
print(f"Galerkin One auto L = {L_auto} ((L+1)^2 = {(L_auto+1)**2})")
# ── Delta-forward XS: sigma_s,ell = sigma_s for ALL ell ──────────────────────
xs_filename = "delta_forward.xs"
with open(xs_filename, "w") as f:
f.write("NUM_GROUPS 1\n")
f.write(f"NUM_MOMENTS {L_max + 1}\n\n")
f.write("SIGMA_T_BEGIN\n")
f.write(f"0 {sigma_t:.10f}\n")
f.write("SIGMA_T_END\n\n")
f.write("TRANSFER_MOMENTS_BEGIN\n")
for ell in range(L_max + 1):
f.write(f"M_GFROM_GTO_VAL {ell} 0 0 {sigma_s:.10f}\n")
f.write("TRANSFER_MOMENTS_END\n")
# ── Pure-absorber reference XS: sigma_t = sigma_a, no scattering ─────────────
xs_ref_filename = "pure_absorber.xs"
with open(xs_ref_filename, "w") as f:
f.write("NUM_GROUPS 1\n")
f.write("NUM_MOMENTS 1\n\n")
f.write("SIGMA_T_BEGIN\n")
f.write(f"0 {sigma_a:.10f}\n")
f.write("SIGMA_T_END\n\n")
f.write("TRANSFER_MOMENTS_BEGIN\n")
f.write("TRANSFER_MOMENTS_END\n")
# ── 3D cube mesh ──────────────────────────────────────────────────────────────
dx = L_domain / N_mesh
nodes = [i * dx for i in range(N_mesh + 1)]
meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes])
grid = meshgen.Execute()
grid.SetUniformBlockID(0)
3.8.2. Operator Construction Methods
OpenSn provides three methods for building the discrete-to-moment (\(\mathbf{D}\)) and moment-to-discrete (\(\mathbf{M}\)) operators:
Method |
D construction |
\(\mathbf{DM}=\mathbf{I}\) guaranteed? |
Scattering order |
|---|---|---|---|
|
\(D_{nm} = w_n\,Y_\ell^m(\hat\Omega_n)\) |
No in 3D |
User-specified |
|
\(D = M^{-1}\) (square system) |
Yes |
Auto: \((L+1)^2 = N_{dir}\) |
|
Orthogonalize \(Y_\ell^m\) |
Yes |
User-specified |
The scattering source in the \(S_N\) equations is \(\sigma_s\,\mathbf{M}\,(\mathbf{D}\,\psi)\). If \(\mathbf{DM}\ne\mathbf{I}\), angular moments are mixed incorrectly: even an angular flux that lies exactly in the column space of \(\mathbf{M}\) is not recovered exactly. For the delta-forward kernel — where every Legendre moment contributes equally — this mixing error accumulates across all \(\ell\) and produces a noticeably different flux profile.
We first verify the \(\mathbf{DM}\) orthogonality numerically, then compare transport solutions.
[ ]:
def check_dm(label, q):
D2M = q.GetDiscreteToMomentOperator() # stored as D^T in OpenSn
M2D = q.GetMomentToDiscreteOperator()
product = D2M.T @ M2D # = D @ M (n_mom x n_mom)
n_mom = M2D.shape[1]
diag = np.diag(product)
mask = ~np.eye(product.shape[0], product.shape[1], dtype=bool)
max_off = np.abs(product[mask]).max() if product[mask].size > 0 else 0.0
max_dev = np.abs(diag - 1.0).max()
print(f" {label}")
print(f" D2M {D2M.shape}, M2D {M2D.shape} ({n_mom} moments)")
print(f" Max |off-diagonal of D@M| : {max_off:.3e}")
print(f" Max |diag(D@M) - 1| : {max_dev:.3e}")
print(f"D@M orthogonality (N_polar={N_polar}, N_azimuthal={N_azimuthal}, {n_dirs} directions)")
print("-" * 60)
q_std = GLCProductQuadrature3DXYZ(
n_polar=N_polar, n_azimuthal=N_azimuthal,
scattering_order=L_compare, operator_method='standard')
q_g1 = GLCProductQuadrature3DXYZ(
n_polar=N_polar, n_azimuthal=N_azimuthal,
operator_method='galerkin_one')
q_g3 = GLCProductQuadrature3DXYZ(
n_polar=N_polar, n_azimuthal=N_azimuthal,
scattering_order=L_compare, operator_method='galerkin_three')
check_dm(f"Standard (L={L_compare}, {(L_compare+1)**2} moments)", q_std)
check_dm(f"Galerkin One (L=auto={L_auto}, {(L_auto+1)**2} moments)", q_g1)
check_dm(f"Galerkin Three (L={L_compare}, {(L_compare+1)**2} moments)", q_g3)
3.8.3. Transport Solutions
We solve the delta-forward transport problem with each operator method and compare the domain-wide scalar flux statistics against the exact analytical solution \(\phi^* = Q/\sigma_a = 2.0\).
For each run we extract:
Max flux — maximum cell-averaged scalar flux over the entire domain
Avg flux — volume-weighted mean scalar flux
Error in avg (%) — \(|\bar\phi - \phi^*|\,/\,\phi^* \times 100\)
Non-uniformity (%) — \((\phi_\mathrm{max} - \bar\phi)\,/\,\phi^* \times 100\)
Because the analytical solution is spatially uniform, non-uniformity is itself a measure of how much operator error is polluting the scattering source.
[ ]:
def run_transport(pquad, xs_file):
"""Run one steady-state transport solve; return (maxval, avgval) on all ranks."""
xs_mat = MultiGroupXS()
xs_mat.LoadFromOpenSn(xs_file)
mg_src = VolumetricSource(block_ids=[0], group_strength=[Q])
phys = DiscreteOrdinatesProblem(
mesh=grid,
num_groups=1,
groupsets=[{
"groups_from_to": (0, 0),
"angular_quadrature": pquad,
"angle_aggregation_num_subsets": 1,
"inner_linear_method": "petsc_gmres",
"l_abs_tol": 1.0e-8,
"l_max_its": 100,
"gmres_restart_interval": 100,
}],
xs_map=[{"block_ids": [0], "xs": xs_mat}],
volumetric_sources=[mg_src],
boundary_conditions=[
{"name": "xmin", "type": "isotropic", "group_strength": [bsrc]},
{"name": "xmax", "type": "isotropic", "group_strength": [bsrc]},
{"name": "ymin", "type": "isotropic", "group_strength": [bsrc]},
{"name": "ymax", "type": "isotropic", "group_strength": [bsrc]},
{"name": "zmin", "type": "isotropic", "group_strength": [bsrc]},
{"name": "zmax", "type": "isotropic", "group_strength": [bsrc]},
],
)
ss = SteadyStateSourceSolver(problem=phys)
ss.Initialize()
ss.Execute()
fflist = phys.GetScalarFluxFieldFunction()
vol_all = RPPLogicalVolume(infx=True, infy=True, infz=True)
ffi_max = FieldFunctionInterpolationVolume()
ffi_max.SetOperationType("max")
ffi_max.SetLogicalVolume(vol_all)
ffi_max.AddFieldFunction(fflist[0])
ffi_max.Initialize()
ffi_max.Execute()
maxval = ffi_max.GetValue()
ffi_avg = FieldFunctionInterpolationVolume()
ffi_avg.SetOperationType("avg")
ffi_avg.SetLogicalVolume(vol_all)
ffi_avg.AddFieldFunction(fflist[0])
ffi_avg.Initialize()
ffi_avg.Execute()
avgval = ffi_avg.GetValue()
return maxval, avgval
# ── Reference (pure absorber = exact delta-forward solution) ──────────────────
print("=== Reference: pure absorber ===")
pquad_ref = GLCProductQuadrature3DXYZ(
n_polar=N_polar, n_azimuthal=N_azimuthal, scattering_order=0)
max_ref, avg_ref = run_transport(pquad_ref, xs_ref_filename)
if rank == 0:
print(f" Max flux = {max_ref:.5f}, Avg flux = {avg_ref:.5f}")
print(f" Analytical = {analytical_phi:.5f}")
# ── Method comparison ─────────────────────────────────────────────────────────
runs = [
(f"Standard (L={L_compare})", q_std),
(f"Galerkin One (L=auto={L_auto})", q_g1),
(f"Galerkin Three (L={L_compare})", q_g3),
]
flux_methods = {}
for label, pquad in runs:
print(f"\n=== {label} ===")
maxval, avgval = run_transport(pquad, xs_filename)
if rank == 0:
flux_methods[label] = (maxval, avgval)
print(f" Max flux = {maxval:.5f}, Avg flux = {avgval:.5f}")
print(f"\nnumber of operator methods tested = {len(runs)}")
[ ]:
if rank == 0 and flux_methods:
col = f"{'Method':<32} {'Max φ':>8} {'Avg φ':>8} {'Err avg %':>10} {'Non-unif %':>11}"
print(col)
print("-" * len(col))
def row(label, maxv, avgv):
err = abs(avgv - analytical_phi) / analytical_phi * 100
nonu = (maxv - avgv) / analytical_phi * 100
print(f"{label:<32} {maxv:>8.5f} {avgv:>8.5f} {err:>10.4f} {nonu:>11.4f}")
row("Reference (pure absorber)", max_ref, avg_ref)
for label, (maxv, avgv) in flux_methods.items():
row(label, maxv, avgv)
print(f"\n Analytical φ* = Q/σ_a = {Q}/{sigma_a} = {analytical_phi:.4f}")
[ ]:
if rank == 0:
for fname in [xs_filename, xs_ref_filename]:
if os.path.exists(fname):
os.remove(fname)
print("Temporary files cleaned up.")
3.8.4. Discussion
3.8.4.1. The delta-forward kernel as a stress test
Because \(\sigma_{s,\ell} = \sigma_s\) for every \(\ell\), the scattering source depends on every angular moment with equal weight. Any error in \(\mathbf{DM}\) is therefore amplified by the full scattering ratio \(c = \sigma_s/\sigma_t = 1.0/1.5 \approx 0.67\), making this a sensitive test of operator consistency.
The uniform-source, uniform-boundary-condition setup provides a clean analytical reference: the exact scalar flux is \(\phi^* = Q/\sigma_a = 2.0\) everywhere. Error in the scattering operator breaks this uniformity, producing a flux that is higher in the interior (where the scattering source is over-estimated) and lower near the boundaries. The non-uniformity metric directly measures this effect without any spatial averaging artefacts.
3.8.4.2. Why standard differs from Galerkin in 3D
The GL–Chebyshev quadrature places azimuthal angles at equally spaced points \(\phi_k = 2\pi k / N_\phi\). This is not a Gaussian rule for trigonometric functions, so integrals of the form \(\sum_n w_n\,Y_{\ell}^{m}(\hat\Omega_n)\,Y_{\ell'}^{m'}(\hat\Omega_n)\) differ from the exact value \(\delta_{\ell\ell'}\delta_{mm'}\). The orthogonality check shows that the standard \(\mathbf{DM}\) has off-diagonal entries of order \(10^{0}\), while both Galerkin methods reduce them to machine precision.
As a result, Standard (:math:`L=7`) systematically mis-estimates the scattering source at each iteration, producing a non-uniform flux with measurable deviation from \(\phi^*\). Additionally, it does not converge within the 100 iteration limit, whereas the two Galerkin Methods converge. Galerkin Three (:math:`L=7`) corrects the operator at the same computational cost and recovers a nearly uniform solution. Galerkin One (auto :math:`L=7`) additionally uses the maximum feasible scattering order — auto-selected so \((L+1)^2 = N_{dir}\) — for the best overall agreement.
3.8.4.3. Practical recommendations
Situation |
Recommendation |
|---|---|
1D slab, GL quadrature |
|
3D product quadrature |
|
Optimal \(L\) unknown |
|
Strongly forward-peaked scattering |
Maximise \(L\); use |
3.8.5. Finalize (for Jupyter Notebook only)
In Python script mode, PyOpenSn automatically handles environment termination. However, this automatic finalization does not occur when running in a Jupyter notebook, so explicit finalization of the environment at the end of the notebook is required. Do not call the finalization in Python script mode, or in console mode.
Note that PyOpenSn’s finalization must be called before MPI’s finalization.
[ ]:
from IPython import get_ipython
def finalize_env():
Finalize()
MPI.Finalize()
ipython_instance = get_ipython()
if ipython_instance is not None:
ipython_instance.events.register("post_execute", finalize_env)