Discrete Ordinates Method (SN)

Key Facts

Read this before modifying the SN solver.

  • Transport equation: \(\mu_m \frac{\partial\psi_m}{\partial x} + \Sigma_t \psi_m = Q/2\) (1D slab)

  • Curvilinear adds angular redistribution: \(\alpha\) coefficients couple ordinates

  • Diamond difference: \(\psi^a = (1+\beta)\psi_{\text{out}} - \beta\psi_{\text{in}}\), Morel-Montry sets \(\beta = 0\)

  • Scattering convention: SigS[l][g_from, g_to] — source uses transpose: Q = SigS^T @ phi

  • GL weights sum to 2; Lebedev/LS/Product sum to \(4\pi\)

  • Gotcha: 1-group tests are degenerate (k = νΣ_f/Σ_a regardless of flux shape)

  • Gotcha: homogeneous tests are blind to curvilinear redistribution bugs (flat flux → α terms vanish)

  • Gotcha: conservation holds even with wrong per-ordinate balance (telescoping sum identity)

  • The \(\alpha\) dome must be non-negative; negative → NaN/overflow

  • Fixed-source flat-flux diagnostic (Q/Σ_t) is the most powerful curvilinear bug detector

  • Key reference: Bailey, Morel & Chang (2009) — Eq. 50 (α recursion), Eq. 74 (M-M weights)

  • Verification uses synthetic cross sections, not real nuclear data

Conventions

Overview

The discrete ordinates (SN) method solves the multi-group eigenvalue problem in integro-differential form by discretising the angular variable \(\hat{\Omega}\) into a finite set of directions (ordinates). Unlike the collision probability method (which works with the integral form and the scalar flux), SN retains the angular flux \(\psi(\mathbf{r}, \hat{\Omega}, E)\) and resolves directional effects such as streaming, anisotropic scattering, and angular current at interfaces.

Three coordinate systems are supported:

  • Cartesian (slab / 2D) — the simplest case; no angular coupling between ordinates.

  • Spherical (1D radial) — angular redistribution in \(\mu\); a single dome of \(\alpha\) coefficients couples all ordinates.

  • Cylindrical (1D radial) — azimuthal redistribution per \(\mu\)-level; independent \(\alpha\) domes on each level.

All three share a single balance-equation framework with a geometry factor \(\Delta A / w\) that guarantees per-ordinate flat-flux consistency. The treatment follows [Bailey2009] for the curvilinear formulation, [LewisMiller1984] for the general framework, and [CaseZweifel1967] for the angular discretisation.

The solver is implemented in SNSolver, which satisfies the EigenvalueSolver protocol. The convenience function solve_sn() runs the full calculation and returns an SNResult.

Architecture

Two-Layer Mesh Pattern

The SN solver follows the same two-layer pattern as the CP solver. This pattern (base Mesh1D + augmented mesh) is shared with Collision Probability Method and Method of Characteristics (MOC).

  1. Base geometry — Mesh1D or Mesh2D stores cell edges, material IDs, and coordinate system.

  2. Augmented geometry — SNMesh wraps the base mesh and an angular quadrature, precomputing the coordinate-specific streaming stencil:

    • Cartesian: streaming_x[n,i] = 2|mu_x|/dx[i] and streaming_y[n,j] = 2|mu_y|/dy[j] — the diamond-difference denominator terms, precomputed to avoid per-cell division in the sweep hot loop.

    • Spherical: face_areas (\(4\pi r^2\)), delta_A, alpha_half (angular redistribution dome), redist_dAw (\(\Delta A_i / w_n\), shape (nx, N)), and tau_mm (Morel–Montry closure weights).

    • Cylindrical: face_areas (\(2\pi r\)), delta_A, alpha_per_level (per-level redistribution domes), redist_dAw_per_level (list of (nx, M) arrays), and tau_mm_per_level (per-level Morel–Montry weights).

  3. Solver — solve_sn() creates an SNMesh, builds the SNSolver, and runs power iteration.

Mesh1D / Mesh2D (base geometry)
    |
    v
SNMesh (stencil + quadrature + alpha coefficients)
    |
    v
solve_sn() --> SNResult

Quadrature Dispatch

The sweep dispatcher in transport_sweep() routes based on the SNMesh.curvature attribute and the quadrature type:

if sn_mesh.curvature == "spherical":
    return _sweep_1d_spherical(...)
elif sn_mesh.curvature == "cylindrical":
    return _sweep_1d_cylindrical(...)
elif is_gl_1d:            # ny=1, mu_y=0, no aniso source
    return _sweep_1d_cumprod(...)
else:
    return _sweep_2d_wavefront(...)

For 1D meshes (ny=1):

  • Gauss–Legendre quadrature takes the fast cumprod path (all \(\mu_y = 0\), so no y-streaming).

  • Lebedev quadrature falls through to the 2D wavefront sweep. Ordinates with \(\mu_x \neq 0\) stream along x; the y-streaming terms cancel via reflective BCs on the single-cell y-dimension. Ordinates with \(\mu_x = \mu_y = 0\) (z-directed) reduce to pure collision: \(\psi = Q \cdot w_{\text{norm}} / \Sigt{}\).

Both quadratures recover the analytical eigenvalue exactly on homogeneous problems (verified to machine precision for 1G/2G/4G).

The Transport Equation

Cartesian 1D (Slab)

The steady-state transport equation in a 1D slab:

(1)\[\mu \frac{\partial \psi(x, \mu)}{\partial x} + \Sigt{} \, \psi(x, \mu) = \frac{Q}{W}\]

where \(\mu = \cos\theta\) is the direction cosine, \(Q\) is the total isotropic source (fission + scattering), and \(W = \sum_n w_n\) is the quadrature weight sum.

Cartesian 2D

In two Cartesian dimensions the angular flux depends on two direction cosines \(\mu_x\) and \(\mu_y\):

(2)\[\mu_x \frac{\partial \psi}{\partial x} + \mu_y \frac{\partial \psi}{\partial y} + \Sigt{} \, \psi = \frac{Q}{W}\]

There is no angular coupling between ordinates — each direction is solved independently. The two streaming terms are the only difference from the 1D case.

Spherical 1D

In spherical coordinates the transport equation acquires an angular redistribution term that couples ordinates:

(3)\[\mu \frac{\partial \psi}{\partial r} + \frac{1 - \mu^2}{r} \frac{\partial \psi}{\partial \mu} + \Sigt{} \psi = \frac{Q}{W}\]

The curvature term \((1 - \mu^2)/r \cdot \partial\psi/\partial\mu\) arises because a neutron streaming radially at angle \(\mu\) rotates its direction cosine as it moves to a different radius. Discretising this term requires diamond difference in both space and angle.

Cylindrical 1D

For an infinitely long cylinder with azimuthal symmetry, the transport equation in the radial variable \(r\) is:

(4)\[\frac{\eta}{r} \frac{\partial(r\psi)}{\partial r} - \frac{1}{r} \frac{\partial(\xi\psi)}{\partial\varphi} + \Sigt{} \psi = \frac{Q}{W}\]

where the direction cosines are:

  • \(\eta = \sin\theta\cos\varphi\) — radial projection (streaming)

  • \(\xi = \sin\theta\sin\varphi\) — azimuthal component

  • \(\mu = \cos\theta\) — axial component

The constraint \(\eta^2 + \xi^2 + \mu^2 = 1\) holds. The azimuthal redistribution \(-\partial(\xi\psi)/\partial\varphi\) couples ordinates on each \(\mu\)-level.

Multi-Group Extension

For \(G\) energy groups, each transport equation becomes a coupled system with scattering transfer \(\Sigs{g' \to g}\) between groups:

(5)\[\text{streaming} + \Sigt{g} \psi_g = \frac{1}{W} \left[ \sum_{g'} \Sigs{g' \to g} \phi_{g'} + \frac{\chi_g}{k} \sum_{g'} \nSigf{g'} \phi_{g'} \right]\]

where the streaming operator depends on the coordinate system and \(\phi_g = \sum_n w_n \psi_{g,n}\) is the scalar flux.

Angular Quadratures

ORPHEUS provides four angular quadrature types for different geometries.

Gauss–Legendre (1D)

For 1D slab and spherical geometry: \(N\) points on \(\mu \in [-1, 1]\), weights sum to 2. Optimal for polynomial integrands (degree \(2N-1\) exact). Also used for spherical 1D, where the single direction cosine \(\mu\) suffices for the angular redistribution.

Implemented in GaussLegendre1D.

Lebedev (Sphere)

For 2D/3D Cartesian geometry: \(N\) points on the unit sphere with octahedral symmetry [Lebedev1999]. Weights sum to \(4\pi\). On a 1D mesh, z-directed ordinates (\(\mu_x = \mu_y = 0\)) are handled as pure collision; all others stream along x with y-terms cancelling via reflective BCs.

Implemented in LebedevSphere.

Level-Symmetric SN

Standard triangular quadrature with \(N/2\) distinct \(\mu_z\) values per hemisphere. Ordinates on each level are permutations of the direction cosine set satisfying \(\eta^2 + \xi^2 + \mu^2 = 1\). Equal spacing in \(\mu^2\) is used with \(\mu_1^2 = 4/(N(N+2))\) [CarlsonLathrop1965].

Weights sum to \(4\pi\). Provides the level_indices structure needed by the cylindrical sweep. Unlike ProductQuadrature (which has one level per \(\mu_z\) value), the Level-Symmetric quadrature groups both \(+\mu_z\) and \(-\mu_z\) hemispheres on the same level (grouped by \(|\mu_z|\)). Within each level, ordinates are sorted by increasing \(\eta\) for the azimuthal sweep.

Implemented in LevelSymmetricSN.

Product Quadrature (GL x equispaced)

Tensor product of Gauss–Legendre in \(\mu = \cos\theta\) (polar) and equispaced points in \(\varphi\) (azimuthal). Each \(\mu\) level has the same number of azimuthal points, giving a clean level structure ideal for the cylindrical sweep. Weights:

\[w_{p,m} = w_{\text{GL}}(\mu_p) \cdot \frac{2\pi}{N_\varphi}\]

Sum to \(4\pi\). Within each level, ordinates are sorted by increasing \(\eta = \sin\theta\cos\varphi\) to match the \(\alpha\) recursion convention from [Bailey2009] Eq. 50.

Implemented in ProductQuadrature.

Reflection Index

Each quadrature implements a reflection_index() method that returns an index array mapping each ordinate \(n\) to its mirror image \(n'\) obtained by negating the direction cosine along a specified axis. For example, reflection_index("x") finds the ordinate whose direction cosines match \((-\mu_x, \mu_y, \mu_z)\).

The implementation in _find_reflections() computes the Euclidean distance between the target direction (with one component negated) and all ordinate directions, then returns the closest match:

\[n' = \arg\min_j \bigl[ (\mu_{x,j} - (-\mu_{x,n}))^2 + (\mu_{y,j} - \mu_{y,n})^2 + (\mu_{z,j} - \mu_{z,n})^2 \bigr]\]

For Gauss–Legendre (1D), the reflection in x is simply \(n' = N - 1 - n\) because the GL points are symmetric about zero. Reflection in y is the identity since \(\mu_y = 0\).

For multi-dimensional quadratures (Lebedev, Level-Symmetric, Product), the reflection indices are precomputed at construction time for all three axes (x, y, z) and stored as _ref_x, _ref_y, _ref_z. These indices are used by the sweep to implement reflective boundary conditions (see Boundary Conditions).

Comparison Table

Quadrature

Geometry

\(\sum w\)

Level structure

Best for

Gauss–Legendre

Slab, Sphere

2

No

1D problems

Lebedev

Cartesian 2D

\(4\pi\)

No

2D/3D Cartesian

Level-Symmetric

Sphere, Cylinder

\(4\pi\)

Yes

Curvilinear

Product

Cylinder

\(4\pi\)

Yes

Cylindrical 1D

The Discrete Balance Equation

This is the core of the SN method. The balance equations are presented from simplest to most complex: Cartesian geometries have no angular redistribution; curvilinear geometries add \(\alpha\) coupling and a geometry factor \(\Delta A/w\).

Cartesian 1D Balance Equation

Integrating (1) over a spatial cell \([x_{i-1/2}, x_{i+1/2}]\) of width \(\Delta x_i\) and applying the divergence theorem to the streaming term:

\[\mu_n \bigl[\psi_{i+\frac12} - \psi_{i-\frac12}\bigr] + \Sigt{} \Delta x_i\, \psi_{n,i} = S_i \Delta x_i\]

where \(S_i = Q_i / W\) and face areas are unity in slab geometry. Applying the diamond-difference closure \(\psi_{n,i} = \frac{1}{2}(\psi_{\rm in} + \psi_{\rm out})\) and \(\psi_{\rm out} = 2\psi_{n,i} - \psi_{\rm in}\), we solve for the cell-average angular flux:

(6)\[\psi_{n,i} = \frac{S_i + \dfrac{2|\mu_n|}{\Delta x_i}\, \psi_{\rm in}} {\Sigt{} + \dfrac{2|\mu_n|}{\Delta x_i}}\]

This is the simplest balance equation: no \(\alpha\) redistribution and no \(\Delta A\) factor, because slab geometry has no curvature. The streaming coefficient \(2|\mu|/\Delta x\) is precomputed by SNMesh as streaming_x[n, i].

Cartesian 2D Balance Equation

Integrating (2) over a rectangular cell \(\Delta x_i \times \Delta y_j\):

\[\mu_{x,n}\bigl[\psi_{i+\frac12,j} - \psi_{i-\frac12,j}\bigr] \Delta y_j + \mu_{y,n}\bigl[\psi_{i,j+\frac12} - \psi_{i,j-\frac12}\bigr] \Delta x_i + \Sigt{} \Delta x_i \Delta y_j\, \psi_{n,i,j} = S_{i,j}\, \Delta x_i \Delta y_j\]

Dividing through by \(\Delta x_i \Delta y_j\) and applying diamond-difference closures in both directions simultaneously:

\[\begin{split}\psi_{n,i} &= \tfrac{1}{2}(\psi^x_{\rm in} + \psi^x_{\rm out}) \qquad\text{(x-closure)} \\ \psi_{n,i} &= \tfrac{1}{2}(\psi^y_{\rm in} + \psi^y_{\rm out}) \qquad\text{(y-closure)}\end{split}\]

yields the 2D DD equation:

(7)\[\psi_{n,i,j} = \frac{S_{i,j} + s_x\, \psi^x_{\rm in} + s_y\, \psi^y_{\rm in}} {\Sigt{} + s_x + s_y}\]

where the streaming coefficients are:

\[s_x = \frac{2|\mu_{x,n}|}{\Delta x_i}, \qquad s_y = \frac{2|\mu_{y,n}|}{\Delta y_j}\]

Both outgoing face fluxes are then updated from the DD closure:

\[\psi^x_{\rm out} = 2\psi_{n,i,j} - \psi^x_{\rm in}, \qquad \psi^y_{\rm out} = 2\psi_{n,i,j} - \psi^y_{\rm in}\]

These are precomputed by SNMesh as streaming_x[n, i] and streaming_y[n, j], so the inner loop in _sweep_2d_wavefront() reduces to a single vectorised division per diagonal.

Curvilinear Balance Equation (Spherical and Cylindrical)

Derivation from the Continuous PDE

Start with the general 1D curvilinear transport equation. In conservative form for a coordinate \(r\) with face area \(A(r)\) and volume element \(V\):

(8)\[\frac{\mu_n}{V_i} \bigl[A_{i+\frac12}\psi_{i+\frac12} - A_{i-\frac12}\psi_{i-\frac12}\bigr] + \frac{1}{V_i} \bigl[\alpha_{n+\frac12}\psi_{n+\frac12} - \alpha_{n-\frac12}\psi_{n-\frac12}\bigr] + \Sigt{} \psi_{n,i} = S_i\]

where the streaming cosine is \(\mu_n\) for spherical and \(\eta_m\) for cylindrical, and \(S_i = Q_i / W\) is the isotropic source density divided by the quadrature weight sum.

Step 1: Integrate the PDE over a spatial cell.

For spherical geometry, integrating (3) over the shell \([r_{i-1/2}, r_{i+1/2}]\) and using the divergence theorem on the radial streaming gives:

\[\mu_n \bigl[A_{i+\frac12}\psi_{i+\frac12} - A_{i-\frac12}\psi_{i-\frac12}\bigr] + \int_{V_i} \frac{1-\mu^2}{r} \frac{\partial\psi}{\partial\mu}\, dV + \Sigt{} V_i \psi_{n,i} = S_i V_i\]

For cylindrical geometry, integrating (4) over the annular shell gives:

\[\eta_m \bigl[A_{i+\frac12}\psi_{i+\frac12} - A_{i-\frac12}\psi_{i-\frac12}\bigr] - \int_{V_i} \frac{1}{r} \frac{\partial(\xi\psi)}{\partial\varphi}\, dV + \Sigt{} V_i \psi_{m,i} = S_i V_i\]

Step 2: Discretise the angular redistribution.

The angular integral is discretised as a finite difference in the ordinate index. For spherical:

\[\int_{V_i} \frac{1-\mu^2}{r}\frac{\partial\psi}{\partial\mu}\, dV \;\approx\; \alpha_{n+\frac12}\psi_{n+\frac12} - \alpha_{n-\frac12}\psi_{n-\frac12}\]

For cylindrical (per \(\mu\)-level):

\[-\int_{V_i} \frac{1}{r}\frac{\partial(\xi\psi)}{\partial\varphi}\, dV \;\approx\; \alpha_{m+\frac12}\psi_{m+\frac12} - \alpha_{m-\frac12}\psi_{m-\frac12}\]

Step 3: Apply the geometry factor \(\Delta A / w\).

The raw discretisation above does NOT preserve per-ordinate flat-flux consistency. The correct form from [Bailey2009] includes the geometry factor \(\Delta A_i / w_n\):

(9)\[\mu_n \bigl[A_{i+\frac12}\psi_{i+\frac12} - A_{i-\frac12}\psi_{i-\frac12}\bigr] + \frac{\Delta A_i}{w_n} \bigl[\alpha_{n+\frac12}\psi_{n+\frac12} - \alpha_{n-\frac12}\psi_{n-\frac12}\bigr] + \Sigt{} V_i \psi_{n,i} = S_i V_i\]

where \(\Delta A_i = A_{i+1/2} - A_{i-1/2}\). This is [Bailey2009] Eq. 7–10 for spherical and Eq. 50–55 for cylindrical.

Note why (6) has no \(\alpha\) or \(\Delta A\) terms: in Cartesian geometry the face area is unity (\(A = 1\)), so \(\Delta A = 0\), and there is no curvature to redistribute angular flux.

The Alpha Redistribution Coefficients

The \(\alpha\) coefficients encode how the angular flux redistributes between neighbouring ordinates due to the geometry curvature. They are defined recursively:

(10)\[\alpha_{n+\frac12} = \alpha_{n-\frac12} - w_n \mu_n\]

with the boundary condition \(\alpha_{1/2} = 0\).

For spherical geometry, all \(N\) ordinates form a single sequence sorted by \(\mu\) (most negative to most positive). The \(\alpha\) values form a non-negative dome: they rise while \(\mu < 0\), peak near \(\mu = 0\), and fall back to zero at \(\mu = 1\). The endpoint condition \(\alpha_{N+1/2} = 0\) is guaranteed by Gauss–Legendre antisymmetry: \(\sum_n w_n \mu_n = 0\).

For cylindrical geometry, each \(\mu\)-level has its own independent \(\alpha\) sequence. On level \(p\), the ordinates are sorted by increasing \(\eta\) (radial direction cosine), and the recursion uses \(\eta\) instead of \(\mu\):

(11)\[\alpha_{p,m+\frac12} = \alpha_{p,m-\frac12} - w_m \eta_m\]

This is [Bailey2009] Eq. 50. Each level’s \(\alpha\) values form an independent dome from \(\eta = -\sin\theta\) to \(\eta = +\sin\theta\).

Dome shape properties:

  • \(\alpha_{n+1/2} \geq 0\) for all \(n\) (non-negative dome).

  • The peak occurs near the ordinate where \(\mu_n\) (or \(\eta_m\)) crosses zero.

  • The dome height scales with the quadrature weight sum: higher-order quadratures have narrower but taller domes.

  • Non-negativity ensures the denominator of the DD equation is unconditionally positive, guaranteeing numerical stability.

The code stores these in SNMesh.alpha_half (spherical, shape (N+1,)) and SNMesh.alpha_per_level (cylindrical, list of (M+1,) arrays).

The Geometry Factor and Why It Is Needed

The geometry factor \(\Delta A_i / w_n\) in (9) is the key to correct curvilinear transport. Without it, the balance equation violates per-ordinate flat-flux consistency: for a spatially uniform, isotropic flux \(\psi = \text{const}\), the streaming and redistribution terms should cancel exactly for EACH ordinate individually.

Proof of consistency.

Set \(\psi_{n,i} = \psi_{n+1/2} = \psi_{n-1/2} = \psi_0\) (flat in both space and angle) and \(\psi_{i+1/2} = \psi_{i-1/2} = \psi_0\) (flat in space). The streaming term becomes:

\[\mu_n \bigl[A_{i+\frac12} - A_{i-\frac12}\bigr] \psi_0 = \mu_n \,\Delta A_i\, \psi_0\]

The redistribution term with the \(\Delta A/w\) factor becomes:

\[\frac{\Delta A_i}{w_n} \bigl[\alpha_{n+\frac12} - \alpha_{n-\frac12}\bigr] \psi_0 = \frac{\Delta A_i}{w_n} (-w_n \mu_n) \psi_0 = -\mu_n \,\Delta A_i\, \psi_0\]

where we used the recursion (10): \(\alpha_{n+1/2} - \alpha_{n-1/2} = -w_n \mu_n\). The two terms cancel exactly, giving \(\Sigt{} \psi_0 = S_0\), which is the correct homogeneous solution.

Without the \(\Delta A/w\) factor (i.e., using \([\alpha_{n+1/2}\psi_{n+1/2} - \alpha_{n-1/2}\psi_{n-1/2}]\) directly), the redistribution term for flat flux is \((-w_n \mu_n)\psi_0\), but the streaming term is \(\mu_n \Delta A_i \psi_0\). These differ by a factor of \(\Delta A_i\), so consistency only holds in the limit \(\Delta A_i \to 0\) (i.e., at the origin or on an infinitely fine mesh).

Consequence of the missing factor: The solver creates artificial angular anisotropy that worsens with mesh refinement near \(r = 0\) (where \(\Delta A_i\) is smallest but non-zero). This manifests as a flux spike at the origin in fixed-source problems and as divergent eigenvalues in heterogeneous eigenvalue problems.

The code precomputes this factor as SNMesh.redist_dAw (spherical, shape (nx, N)) and SNMesh.redist_dAw_per_level (cylindrical, list of (nx, M) arrays).

The Morel–Montry Flux Dip

Even with the correct \(\Delta A/w\) factor, the standard diamond-difference closure (equal weight \(\tau = 0.5\)) introduces a flux error near \(r = 0\) known as the Morel–Montry flux dip [MorelMontry1984].

The standard DD angular closure is:

\[\psi_{n,i} = \frac{1}{2}(\psi_{n-\frac12} + \psi_{n+\frac12})\]

This can be rewritten as:

\[\psi_{n+\frac12} = 2\psi_{n,i} - \psi_{n-\frac12}\]

The contamination factor \(\beta\) ([Bailey2009] Eq. 41) quantifies the coupling between the leading-order scalar flux and the first-order current in the asymptotic diffusion limit. For spherical geometry:

\[\beta = \frac{1}{2} \sum_{n=1}^{N} \mu_n \bigl[\alpha_{n+\frac12}\, \mu_{n+\frac12} - \alpha_{n-\frac12}\, \mu_{n-\frac12}\bigr]\]

where \(\mu_{n\pm 1/2}\) are the angular cell-edge cosines. For cylindrical, the equivalent is a per-level sum using \(\eta\) and \(\eta_{m\pm 1/2}\). When \(\beta \neq 0\), the discrete SN equations satisfy a contaminated diffusion equation near \(r = 0\), producing the artificial flux dip (or spike).

The module derivations.sn_contamination computes \(\beta\) for any quadrature and geometry. With the correct \(\Delta A/w\) factor AND Morel–Montry weights, \(\beta \sim 10^{-16}\) (machine zero) for both spherical and cylindrical.

Weighted Diamond Difference (WDD) and Morel–Montry Weights

The Morel–Montry (M-M) angular closure replaces the equal-weight DD with position-dependent weights \(\tau_n\) [Bailey2009] Eq. 74:

(12)\[\psi_{n,i} = (1 - \tau_n)\,\psi_{n-\frac12} + \tau_n\,\psi_{n+\frac12}\]

Solving for the angular face flux:

(13)\[\psi_{n+\frac12} = \frac{\psi_{n,i} - (1 - \tau_n)\,\psi_{n-\frac12}}{\tau_n}\]

The M-M weights are defined as:

(14)\[\tau_n = \frac{\mu_n - \mu_{n-\frac12}}{\mu_{n+\frac12} - \mu_{n-\frac12}}\]

where \(\mu_{n\pm 1/2}\) are the angular cell edges.

Spherical cell edges: \(\mu_{1/2} = -1\), \(\mu_{N+1/2} = +1\), and interior edges by weight-sum: \(\mu_{n+1/2} = \mu_{n-1/2} + w_n\). This is exact for Gauss–Legendre quadrature because the weights correspond to the \(\mu\)-space widths of the angular cells.

Cylindrical cell edges: \(\eta_{1/2} = -\sin\theta\), \(\eta_{M+1/2} = +\sin\theta\), and interior edges at midpoints of consecutive \(\eta\) values: \(\eta_{m+1/2} = (\eta_m + \eta_{m+1})/2\). The weight-sum approach is NOT used for cylindrical because the quadrature weights are uniform in \(\varphi\)-space (not \(\eta\)-space): the Product quadrature spaces \(\varphi\) equally, but \(\eta = \sin\theta\cos\varphi\) is cosine-distributed, so equal \(\varphi\)-widths map to unequal \(\eta\)-widths. The midpoint approach gives a proper partition of \([-\sin\theta, +\sin\theta]\).

For the Product quadrature with equally-spaced \(\varphi\), ordinates come in pairs with the same \(|\eta|\) but opposite \(\xi\) (e.g., \(\varphi = \pi/4\) and \(\varphi = 7\pi/4\) both give \(\eta = \sin\theta/\sqrt{2}\)). The midpoint between paired ordinates equals their shared \(\eta\), creating zero-width angular cells. The resulting \(\tau\) alternates between 0.5 (DD, for the left member of each pair) and 1.0 (step, for the right member). This alternating pattern is correct but could be smoothed by using a Gauss-type azimuthal quadrature with distinct \(\eta\) values (see IMPROVEMENTS.md item DO-20260405-002).

The M-M weights force the contamination factor \(\beta\) to machine zero (verified: \(\beta \sim 10^{-16}\)), completely eliminating the Morel–Montry flux dip.

Clipping: The code clips \(\tau_n\) to \([0.5, 1.0]\) to prevent negative angular face fluxes. This preserves the stability of the standard DD while gaining the accuracy of the M-M correction.

The code stores these as SNMesh.tau_mm (spherical, shape (N,)) and SNMesh.tau_mm_per_level (cylindrical, list of (M,) arrays).

Substituting the WDD Closure into the Balance Equation

Combining the balance equation (9) with the WDD angular closure (12) and the standard spatial DD (\(\psi_{n,i} = \frac{1}{2}(\psi_{\rm in}^s + \psi_{\rm out}^s)\), \(\psi_{\rm out}^s = 2\psi_{n,i} - \psi_{\rm in}^s\)), define:

\[\begin{split}c_{\rm out} &= \frac{\alpha_{n+\frac12}}{\tau_n} \\[6pt] c_{\rm in} &= \frac{(1-\tau_n)}{\tau_n}\,\alpha_{n+\frac12} + \alpha_{n-\frac12}\end{split}\]

The cell-average angular flux is then:

(15)\[\psi_{n,i} = \frac{ S_i V_i + |\mu_n|(A_{\rm in} + A_{\rm out})\,\psi_{\rm in}^s + \dfrac{\Delta A_i}{w_n}\, c_{\rm in}\, \psi_{n-\frac12} }{ 2|\mu_n|\, A_{\rm out}^s + \dfrac{\Delta A_i}{w_n}\, c_{\rm out} + \Sigt{} V_i }\]

where the superscript \(s\) denotes spatial face fluxes, and \(A_{\rm in}\), \(A_{\rm out}\) are the cell face areas in the direction of neutron travel (see Sweep Algorithm below for their definition). This is the equation solved by both _sweep_1d_spherical() and _sweep_1d_cylindrical().

Geometry Comparison

Aspect

Cartesian

Spherical

Cylindrical

Streaming cosine

\(\mu\)

\(\mu\)

\(\eta\) (radial)

Face area \(A\)

1 (per unit area)

\(4\pi r^2\)

\(2\pi r\)

Volume \(V\)

\(\Delta x\)

\(\tfrac{4}{3}\pi(r_{\rm out}^3 - r_{\rm in}^3)\)

\(\pi(r_{\rm out}^2 - r_{\rm in}^2)\)

\(\Delta A\)

0

\(4\pi(r_{\rm out}^2 - r_{\rm in}^2)\)

\(2\pi(r_{\rm out} - r_{\rm in})\)

Redistribution

None

\(+(\Delta A/w)\,[\alpha\psi]\)

\(+(\Delta A/w)\,[\alpha\psi]\)

\(\alpha\) scope

N/A

Global (all \(N\) ordinates)

Per \(\mu\)-level

\(\alpha\) recursion variable

N/A

\(\mu\)

\(\eta\)

Quadrature required

GL or Lebedev

GL

Product or Level-Sym

Sweep Algorithm

Because each cell’s outgoing flux becomes the next cell’s incoming flux, the equations must be solved in the direction of neutron travel — this is called a transport sweep.

Cartesian 1D: Cumprod Recurrence

For the 1D slab with Gauss–Legendre quadrature, the DD equation (6) defines a recurrence for the outgoing face flux:

(16)\[\psi_{\rm out} = a_i\, \psi_{\rm in} + b_i\]

where the coefficients for cell \(i\) are:

\[a_i = \frac{2|\mu_n|/\Delta x_i - \Sigt{}} {2|\mu_n|/\Delta x_i + \Sigt{}}, \qquad b_i = \frac{S_i} {2|\mu_n|/\Delta x_i + \Sigt{}}\]

This arises from substituting the DD closure \(\psi_{\rm out} = 2\psi_{\rm avg} - \psi_{\rm in}\) into (6). The coefficient \(a_i\) is the stream-to-collision ratio: it controls how much incoming flux propagates through cell \(i\).

Unrolling the recurrence \(\psi_{\rm out}^{(i)} = a_i\, \psi_{\rm out}^{(i-1)} + b_i\) gives a linear first-order relation that can be solved analytically using cumulative products. Define:

\[C_i = \prod_{k=0}^{i} a_k, \qquad R_i = \sum_{k=0}^{i} \frac{b_k}{C_k}\]

Then the incoming face flux at cell \(i+1\) is:

\[\psi_{\rm in}^{(i+1)} = C_i \bigl(\psi_{\rm in}^{(0)} + R_i\bigr)\]

and the cell-average flux is \(\psi_{\rm avg}^{(i)} = \frac{1}{2}(\psi_{\rm in}^{(i)} + \psi_{\rm out}^{(i)})\).

The implementation in _sweep_1d_cumprod() computes \(C\) and \(R\) via np.cumprod and np.cumsum, giving an \(O(N \cdot n_x)\) vectorised sweep — all spatial cells for a given ordinate are resolved simultaneously in numpy array operations, with no Python-level cell loop. This typically runs in sub-millisecond time for practical meshes.

Exploiting GL symmetry, only positive-\(\mu\) ordinates are swept forward; negative-\(\mu\) ordinates are obtained by reversing the cell array and sweeping with the same coefficients.

Cartesian 2D: Anti-Diagonal Wavefront Sweep

In 2D, the DD equation (7) creates a data dependency: cell \((i, j)\) requires incoming face fluxes from its upwind neighbours in both \(x\) and \(y\). Cells along an anti-diagonal \(i + j = k\) are mutually independent because they share no incoming faces, so they can be solved simultaneously.

The four quadrant sweeps. Each ordinate has a sign pair \((\text{sgn}(\mu_x), \text{sgn}(\mu_y))\) that determines the sweep direction. The four combinations define four sweep patterns:

\(\mu_x\)

\(\mu_y\)

x-direction

y-direction

\(+\)

\(+\)

left \(\to\) right

bottom \(\to\) top

\(-\)

\(+\)

right \(\to\) left

bottom \(\to\) top

\(+\)

\(-\)

left \(\to\) right

top \(\to\) bottom

\(-\)

\(-\)

right \(\to\) left

top \(\to\) bottom

For each direction pair, the sweep visits anti-diagonals \(k = 0, 1, \ldots, n_x + n_y - 2\). On diagonal \(k\), the cells \((i, j)\) satisfying \(i + j = k\) (in the swept index space) are gathered into a numpy batch and solved with a single vectorised evaluation of (7).

Vectorisation within each diagonal. Each diagonal contains up to \(\min(n_x, n_y)\) cells. The incoming face fluxes psi_in_x and psi_in_y are gathered by advanced indexing; the DD equation is evaluated as one numpy operation; and the outgoing face fluxes are scattered back. There is no Python-level cell loop within a diagonal.

Reflective BCs in 2D. At each boundary face, the incoming flux for ordinate \(n\) is set to the outgoing flux of its reflected partner. For the left/right boundaries (x-reflection), the partner is ref_x[n] (negating \(\mu_x\)); for the top/bottom boundaries (y-reflection), the partner is ref_y[n] (negating \(\mu_y\)). The reflection indices are precomputed by the quadrature’s reflection_index() method.

Implemented in _sweep_2d_wavefront().

Curvilinear 1D: Sequential Ordinate Sweep

For spherical and cylindrical geometries, the angular redistribution couples successive ordinates, preventing vectorisation across the ordinate dimension. The sweep proceeds cell-by-cell, ordinate-by-ordinate:

Spherical: Ordinates are processed from most negative \(\mu\) to most positive (a single global sequence). Negative-\(\mu\) ordinates sweep inward (outer boundary to centre); positive-\(\mu\) ordinates sweep outward (centre to outer boundary).

Cylindrical: For each \(\mu\)-level, azimuthal ordinates are processed from most-inward (\(\eta = -\sin\theta\)) to most-outward (\(\eta = +\sin\theta\)). Negative-\(\eta\) ordinates sweep inward; \(\eta \approx 0\) ordinates have no radial streaming (pure redistribution); positive-\(\eta\) ordinates sweep outward.

At each cell, the sweep solves (15) for \(\psi_{n,i}\), then updates:

  1. Spatial face flux: \(\psi_{\rm out}^s = 2\psi_{n,i} - \psi_{\rm in}^s\)

  2. Angular face flux: \(\psi_{n+1/2} = (\psi_{n,i} - (1-\tau_n)\psi_{n-1/2})/\tau_n\) (using M-M weights)

  3. Scalar flux accumulation: \(\phi_i \mathrel{+}= w_n \psi_{n,i}\)

The spatial face flux propagates to the next cell; the angular face flux propagates to the next ordinate on the same cell.

Implemented in _sweep_1d_spherical() and _sweep_1d_cylindrical().

Starting Direction

Both curvilinear sweeps initialise the angular face flux \(\psi_{1/2}\) to zero for each spatial cell. This is valid because \(\alpha_{1/2} = 0\) by construction, so the product \(\alpha_{1/2} \psi_{1/2}\) in the balance equation vanishes regardless of the value of \(\psi_{1/2}\). The starting-direction treatment of [LewisMiller1984] Section 4.5.4 (which tracks \(\alpha\psi\) instead of \(\psi\) alone) is therefore unnecessary when the \(\alpha\) recursion is implemented correctly.

BiCGSTAB Alternative

Instead of sweep-based source iteration, the within-group transport problem can be solved directly using a Krylov method (BiCGSTAB).

Explicit Transport Operator

The finite-difference transport operator \(T\) is formed explicitly:

\[T\psi = \mu_n \nabla\psi + \Sigt{}\psi\]

For Cartesian geometry, this is a banded matrix with finite-difference gradients. For curvilinear geometries, the operator includes the same \(\Delta A/w\) geometry factor and Morel–Montry angular closure weights used by the sweep. The system \(T\psi = b\) is solved with scipy.sparse.linalg.bicgstab.

Consistency with the Sweep

Both the sweep and BiCGSTAB path must use the same spatial discretisation to produce identical eigenvalues. In practice:

  • The sweep uses diamond-difference (DD): \(T^{-1}\) is applied implicitly.

  • BiCGSTAB forms \(T\) explicitly using finite-difference (FD) gradients.

On coarse meshes, DD and FD have different truncation error constants, so the two paths may give slightly different \(\keff\) values. They converge to the same answer as \(h \to 0\).

The curvilinear BiCGSTAB operators read redist_dAw (spherical) and redist_dAw_per_level (cylindrical) from SNMesh, along with the M-M weights tau_mm / tau_mm_per_level. This ensures both paths share exactly the same physics.

Boundary Conditions

ORPHEUS uses reflective boundary conditions on all faces, representing an infinite lattice or infinite medium. The CP solver uses white (isotropic) BCs instead; see Numerical Evidence: White-BC Approximation Quality for a comparison showing the ~1% gap between the two approaches.

Outer Boundary

At the outer boundary \(r = R\) (or \(x = L\)), the incoming flux for ordinate \(n\) is set to the outgoing flux of its reflected partner:

(17)\[\psi_n^{\rm in} = \psi_{n'}^{\rm out}\]

where \(n'\) is the reflected partner ordinate (negating the appropriate direction cosine). Reflective partner indices are precomputed by each quadrature’s reflection_index() method.

Inner Boundary (Curvilinear)

At \(r = 0\):

  • The face area \(A(0) = 0\), so no spatial flux crosses the origin. The spatial incoming flux for outward-sweeping ordinates is zero.

  • The angular redistribution provides the inward-to-outward transition: flux entering as an inward-directed ordinate (\(\mu < 0\) or \(\eta < 0\)) is redistributed to outward ordinates (\(\mu > 0\) or \(\eta > 0\)) through the \(\alpha\) coupling.

This means the curvilinear sweep does not need an explicit boundary condition at \(r = 0\) — the geometry handles it naturally.

Scattering

Matrix Convention

The Mixture.SigS[l] matrices use the convention \(\text{SigS}[g_{\rm from}, g_{\rm to}]\):

\[\begin{split}\text{SigS}[0] = \begin{pmatrix} \Sigs{0\to0} & \Sigs{0\to1} \\ \Sigs{1\to0} & \Sigs{1\to1} \end{pmatrix}\end{split}\]

For the in-scatter source (total scattering into group \(g\) from all groups \(g'\)):

\[Q_{\rm scatter}[g] = \sum_{g'} \Sigs{g'\to g}\, \phi_{g'} = (\text{SigS}^T \cdot \phi)[g]\]

The vectorised form for batched cells is phi @ SigS (equivalent to \((\text{SigS}^T \phi^T)^T\) for row-vector \(\phi\)).

The analytical eigenvalue problem uses:

\[\begin{split}A &= \text{diag}(\Sigt{}) - \text{SigS}^T \quad\text{(removal matrix)} \\ F &= \chi \otimes \nSigf{} \quad\text{(fission matrix)} \\ \kinf &= \lambda_{\max}(A^{-1} F)\end{split}\]

Note the transpose: \(\text{SigS}^T[g, g'] = \Sigs{g'\to g}\) gives the in-scatter contribution, so \(\text{diag}(\Sigt{}) - \text{SigS}^T\) is the net removal matrix. See Scattering Matrix Convention for the full derivation of this convention.

P0 Isotropic Scattering

The default mode (scattering_order=0). A direction-independent source is added to all ordinates equally:

\[Q_{\rm scatter}(\hat{\Omega}_n) = \sum_{g'} \Sigs{g'\to g}^{(0)}\, \phi_{g'} / W\]

Implemented in SNSolver._add_scattering_source(), which performs phi @ SigS[0] per material.

PN Anisotropic Scattering

When scattering_order >= 1, per-ordinate anisotropic sources are computed from the Legendre moments of the angular flux. The full anisotropic scattering source for ordinate \(n\) and group \(g\) is:

(18)\[Q_{\rm scatter}(\hat{\Omega}_n, g) = \sum_{\ell=0}^{L} (2\ell+1) \sum_{m=-\ell}^{\ell} \sum_{g'} \Sigs{g'\to g}^{(\ell)}\, f_{\ell,g'}^m \; Y_\ell^m(\hat{\Omega}_n)\]

where \(Y_\ell^m\) are real spherical harmonics and the angular flux moments are computed by quadrature:

(19)\[f_{\ell,g}^m = \sum_{n=1}^{N} w_n \, \psi_{n,g} \, Y_\ell^m(\hat{\Omega}_n)\]

The \((2\ell+1)\) factor is the addition theorem normalisation for real spherical harmonics: it ensures that the PL expansion reproduces the angular flux moments exactly when the angular flux is a polynomial of degree \(\leq L\).

Implementation in SNSolver._build_aniso_scattering():

  1. Compute spherical harmonics at construction time: \(Y[n, \ell, \ell+m]\) for all ordinates, stored as self._Y with shape (N, L+1, 2L+1). The convention is \(Y_0^0 = 1\), \(Y_1^{-1} = \mu_z\), \(Y_1^0 = \mu_x\), \(Y_1^1 = \mu_y\).

  2. Compute flux moments via an einsum contraction over the ordinate index:

    fiL[:, :, :, l, l+m] = np.einsum(
    'n,nxyg->xyg', w * Y[:, l, l+m], angular_flux,
)

    This contracts \(\sum_n w_n Y_\ell^m(\hat{\Omega}_n) \psi_n(x,y,g)\) into a spatial-energy field of shape (nx, ny, ng).

  3. Reconstruct per-ordinate source: for each Legendre order \(\ell \geq 1\) (the \(\ell = 0\) term is handled by SNSolver._add_scattering_source()) and each \(m\), the scattered moment moment @ sig_s_l[l] is multiplied by \((2\ell+1) Y_\ell^m(\hat{\Omega}_n)\) and accumulated into Q_aniso[n, :, :, :].

  4. The resulting Q_aniso array of shape (N, nx, ny, ng) is passed to transport_sweep(), which adds it to the isotropic source on a per-ordinate basis.

Equivalence of the code to the mathematical form. Equation (18) writes the sum as \(\sum_\ell \sum_m \sum_{g'} \Sigs{}^{(\ell)} f_\ell^m Y_\ell^m\). The code separates the \(\ell = 0\) term (isotropic, handled by _add_scattering_source) from the \(\ell \geq 1\) terms (anisotropic, handled by _build_aniso_scattering). For \(\ell = 0\), \(Y_0^0 = 1\) and \((2 \cdot 0 + 1) = 1\), so the sum reduces to \(\sum_{g'} \Sigs{g' \to g}^{(0)} f_{0,g'}^0 = \sum_{g'} \Sigs{g' \to g}^{(0)} \phi_{g'}\), which is exactly the P0 source. The split is therefore exact with no double-counting.

The 421-group cross-section library provides both P0 and P1 matrices.

(n,2n) Reactions

The \((n,2n)\) reaction is a threshold reaction in which a neutron is absorbed by a nucleus, which then emits two neutrons. The net effect is a gain of one neutron per reaction (the incident neutron is consumed, two are produced).

The \((n,2n)\) cross section is stored as a group-to-group transfer matrix Mixture.Sig2 with the same [g_from, g_to] convention as the scattering matrix. The source contribution is:

\[Q_{(n,2n)}(g) = 2 \sum_{g'} \Sigma_{2,g'\to g}\, \phi_{g'}\]

The factor of 2 accounts for the two neutrons produced per reaction. The implementation in SNSolver._add_n2n_source() performs:

Q[ix, iy, :] += 2.0 * (phi[ix, iy, :] @ self.sig2[mid])

This is added to the isotropic source before the transport sweep, on the same footing as the P0 scattering source. The \((n,2n)\) contribution also enters the \(\keff\) production term in SNSolver.compute_keff(), where row sums of Sig2 (total \((n,2n)\) removal rate) are used.

Normalization Chain

The normalization chain in the code ensures consistent scaling:

  1. Fission source (SNSolver.compute_fission_source()): \(Q_f = \chi \cdot (\nSigf{} \cdot \phi) / k\) — raw, un-normalised.

  2. Scattering source (SNSolver._add_scattering_source()): \(Q_s = \text{SigS}^T \cdot \phi\) — also un-normalised.

  3. Sweep (transport_sweep()): applies \(Q_{\rm scaled} = Q \cdot w_{\rm norm}\) where \(w_{\rm norm} = 1/\sum w_n\). This is the \(1/W\) division in the SN equation.

  4. Scalar flux (inside sweep): \(\phi = \sum_n w_n \psi_n\) — standard quadrature integration.

  5. keff (SNSolver.compute_keff()): \(k = (\nSigf{} \cdot \phi \cdot V) / (\Sigma_a \cdot \phi \cdot V)\) — volume-weighted ratio.

The \(1/W\) in step 3 and the \(W\) implicit in step 4 cancel: \(\phi = \sum w_n \cdot Q/(W \Sigt{}) = Q/\Sigt{}\) for uniform isotropic source.

Convention rule: Sources passed to the sweep must NOT include \(1/W\) — the sweep applies it. The BiCGSTAB path (direct operator) must divide sources by \(W\) itself, since it solves \(T\psi = b\) without the sweep.

The Eigenvalue Problem

The eigenvalue \(\keff\) is determined by power iteration: an outer loop updates \(k\) from the production/absorption ratio, with an inner loop that solves the within-group scattering problem.

Two Inner Solvers

Source iteration (sweep-based):

  • Operator: \(T^{-1}\) (diamond-difference sweep)

  • Solution variable: scalar flux \(\phi(x, y, g)\)

  • Fixed-point: \(\phi^{(k+1)} = T^{-1}(S \cdot \phi^{(k)} + Q_f)\)

  • Convergence rate: spectral radius of \(T^{-1}S\) (~0.97 for 421 groups)

  • Cost per iteration: one transport sweep

  • Works for all geometries

BiCGSTAB (direct operator):

  • Operator: \(T = \mu \nabla + \Sigt{}\) (finite-difference gradients)

  • Solution variable: angular flux \(\psi(x, y, n, g)\) (much larger)

  • System: \(T\psi = b\) where \(b\) = fission + scattering

  • Convergence: ~100 Krylov iterations at tol=1e-4 (always converges)

  • Available for all geometries (Cartesian, spherical, cylindrical)

The two architectures use different spatial discretisations (DD sweep vs FD gradient) that converge to different \(\keff\) on coarse meshes. They agree in the limit \(h \to 0\).

Verification

Homogeneous Infinite Medium

For homogeneous geometry with reflective BCs, the flux is spatially flat and \(\keff = \lambda_{\max}(A^{-1}F)\). This is geometry-independent — Cartesian, spherical, and cylindrical must all give the same \(\keff\).

Groups

\(\kinf\)

Cartesian (GL S8)

Spherical (GL S8)

Cylindrical (Prod 4x8)

Cylindrical (LS S4)

1

1.5000

exact

exact

exact

exact

2

1.8750

exact

exact

exact

exact

4

1.4878

exact

exact

exact

exact

All entries are exact to machine precision. Spherical 2G/4G results (previously showing ~1% error) are now exact thanks to the M-M angular closure weights.

Heterogeneous Convergence

For a cylindrical fuel (r < 0.5) + moderator (r < 1.0) geometry with Product(4x8) quadrature:

Cells/region

\(\keff\) (1G)

\(\Delta k\) from previous

5

0.9769

10

0.9842

+0.0073

20

0.9874

+0.0032

\(\keff\) converges monotonically toward the CP reference (0.9955). The ~1% residual gap is the white-BC (CP) vs reflective-BC (SN) approximation difference, consistent with the slab geometry findings.

For the 2G heterogeneous resolution test, Product(4x8) and Product(8x8) agree to \(< 0.01\%\) (keff = 0.7227 for both), confirming angular convergence.

Why 1-Group Verification Is Degenerate

For 1 energy group, the eigenvalue is:

\[k = \frac{\nSigf{}}{\Sigma_a}\]

This is a scalar ratio independent of the spatial or angular flux distribution. Consequences:

  • Angular weight errors scale all flux equally — cancels in \(k\).

  • Wrong scattering convention — no inter-group coupling to distort.

  • Wrong flux shape — does not matter; \(k\) is a material property.

Only multi-group problems have a flux-shape-dependent eigenvalue: \(k = (\nSigf{} \cdot \phi) / (\Sigma_a \cdot \phi)\) where the dot product weights each group differently. A wrong group ratio (from angular errors, normalization errors, or convergence failures) directly shifts \(\keff\).

Rule: Every transport solver must be verified on at least 2-group problems. 1-group success gives false confidence.

Spatial and Angular Convergence

The diamond-difference scheme converges at \(O(h^2)\) with mesh refinement. Gauss–Legendre quadrature shows spectral convergence in angle. Both are verified in test_sn_1d.py.

Property Tests

For all geometries:

  • Particle balance: production / absorption \(= \keff\)

  • Flux non-negativity: \(\phi \geq 0\) everywhere

  • Angular flux at \(r = 0\) all positive (curvilinear only)

  • Multi-group eigenvector not flat: flux spectrum differs between fuel and moderator (catches 1G-degenerate bugs)

Run the full suite:

pytest tests/test_sn_1d.py tests/test_sn_properties.py \
       tests/test_sn_solver_components.py tests/test_sn_spherical.py \
       tests/test_sn_cylindrical.py tests/test_sn_quadrature.py \
       tests/test_sn_sweep_regression.py -v -m "not slow"

Investigation History: Curvilinear Bug

This section documents the full history of the cylindrical DD bug and its resolution. It is preserved to prevent future sessions from repeating the same dead ends.

Symptoms

  1. Homogeneous eigenvalue problems: exact (1G/2G/4G).

  2. Heterogeneous eigenvalue problems: divergent \(\keff\) with mesh refinement (5 cells: 1.15, 10 cells: 0.90, 20 cells: 0.52).

  3. Fixed-source flux range: [0.59, 5.09] (should be near-flat).

  4. \(\keff\) depended strongly on angular quadrature order (4x8 vs 8x8 gave a 67% gap on heterogeneous problems).

The Root Cause

Two bugs, both breaking per-ordinate flat-flux consistency:

Bug 1: Wrong \(\alpha\) recursion. The code used \(\alpha = \text{cumsum}(+w\xi)\) with the azimuthal cosine \(\xi\) (mu_y). The correct recursion ([Bailey2009] Eq. 50) is \(\alpha = \text{cumsum}(-w\eta)\) with the radial cosine \(\eta\) (mu_x), and ordinates must be sorted by increasing \(\eta\) within each level.

Bug 2: Missing \(\Delta A/w\) geometry factor. The redistribution term in the balance equation must include \(\Delta A_i / w_m\). Without this factor, the streaming and redistribution do NOT cancel per-ordinate for a spatially flat flux, creating artificial angular anisotropy that worsens with mesh refinement near \(r = 0\).

Six Failed Approaches

Before the correct fix was found, six approaches were tested. All failed because they addressed symptoms, not the root cause:

  1. Reverse sweep: Reversed the azimuthal sweep direction. No effect — the dome shape is symmetric.

  2. Step closure: Replaced DD (\(\tau = 0.5\)) with step (\(\tau = 1.0\)) in angle. Reduced the divergence slightly but did not eliminate it.

  3. Lewis & Miller starting direction: Tracked \(\alpha\psi\) instead of \(\psi\) alone (Section 4.5.4 of [LewisMiller1984]). Unnecessary when the \(\alpha\) recursion is correct, since \(\alpha_{1/2} = 0\).

  4. Bidirectional sweep: Swept azimuthal ordinates in both directions and averaged. Masked the asymmetry but did not fix it.

  5. Scaled \(\alpha\): Empirically scaled the \(\alpha\) coefficients. Could not find a consistent scaling.

  6. Zero redistribution: Set \(\alpha = 0\) to isolate the spatial streaming. Confirmed the bug was in the redistribution, not the spatial DD.

Why the Original Sign Convention Hypothesis Was Wrong

The initial hypothesis was that the minus sign before the redistribution in (4) required special treatment (tracking \(\alpha\psi\), using \(|\alpha|\), etc.). This was wrong. The sign convention is absorbed into the \(\alpha\) definition: \(\text{cumsum}(-\eta w)\) with a \(+\) sign in the balance equation gives the same physics as \(\text{cumsum}(+\xi w)\) with a \(-\) sign, for symmetric quadratures.

The real issue was the missing \(\Delta A/w\) geometry factor, which has nothing to do with signs.

The Fix

Applied the [Bailey2009] formulation:

  1. Corrected \(\alpha\) recursion: \(\alpha_{m+1/2} = \alpha_{m-1/2} - w_m \eta_m\) with ordinates sorted by increasing \(\eta\).

  2. Added \(\Delta A/w\) geometry factor to both the DD sweep and the BiCGSTAB operator.

  3. Added Morel–Montry angular closure weights (M-M) to eliminate the flux dip at \(r = 0\).

  4. Applied the same \(\Delta A/w\) fix to the spherical sweep (which had the same missing factor). The fixed-source flux spike at \(r = 0\) dropped from 5.1x to 1.1x.

  5. Applied the fix to both the spherical and cylindrical BiCGSTAB operators (transport_operator_matvec_spherical(), transport_operator_matvec_cylindrical()). Multi-group spherical BiCGSTAB had been unstable (keff → NaN for 2G+); the root cause was the same missing \(\Delta A/w\) factor in the explicit FD operator. After the fix, 2G and 4G spherical BiCGSTAB converge to \(< 10^{-6}\) of the analytical eigenvalue.

Results After Fix

Test

Before

After

Homogeneous 1G/2G/4G

Exact

Exact

Heterogeneous 1G, 5/10/20 cells

1.15 / 0.90 / 0.52 (diverges)

0.977 / 0.984 / 0.987 (converges)

Heterogeneous 2G, 4x8 vs 8x8

0.54 vs 0.91 (67% gap)

0.723 vs 0.723 (<0.01%)

Fixed-source flux range (40 cells)

[0.59, 5.09] (spike)

[0.51, 1.12] (bounded)

Contamination \(\beta\) (cylindrical)

~2.0

~10-16 (machine zero)

Numerical Sensitivities

\(\keff\) Sensitivity Table (421-Group Heterogeneous PWR Slab)

All cases: 10 cells, \(\delta = 0.2\) cm, material layout [fuel x 5, clad x 1, cool x 4], P0 scattering, 421 energy groups.

Configuration

\(\keff\)

Notes

1D GL S16, BiCGSTAB (FD operator)

1.03882

True 1D, 16 ordinates

1D Lebedev 110, source iteration (DD sweep)

1.04294

1D mesh, 2D quadrature

2D (10x2) Lebedev 110, source iter (DD sweep)

1.04294

Pseudo-2D, full volumes

2D (10x2) Lebedev 110, BiCGSTAB (FD)

1.04007

Pseudo-2D, full volumes

2D (10x2) Lebedev 110, BiCGSTAB, half-volumes

1.04192

MATLAB convention

MATLAB reference

1.04188

2D Lebedev, FD, half-volumes

Sources of Variation

  1. Angular quadrature (GL vs Lebedev): ~0.004 difference. GL S16 integrates 1D angular flux with 16 points on \([-1,1]\). Lebedev 110 integrates over the unit sphere — more angular resolution but different effective weights per \(\mu_x\) direction. On a coarse heterogeneous mesh, these give different eigenvalues.

  2. Spatial discretisation (DD sweep vs FD gradient): ~0.003 difference. Source iteration uses the DD wavefront sweep (\(T^{-1}\)). BiCGSTAB uses the explicit FD transport operator (\(T\)). Both are \(O(h)\) on this mesh but with different truncation error constants.

  3. Boundary volume weighting: ~0.002 difference (full vs half). The MATLAB code halves boundary cell volumes. With ny=2 and materials uniform in y, only the x-direction halving (fuel edge, coolant edge) affects \(\keff\). This is an artifact of the pseudo-2D implementation: a true 1D calculation has no y-volumes.

  4. Inner convergence: source iteration with max_inner=200, inner_tol=1e-8 does not fully converge for 421 groups (spectral radius ~0.97). BiCGSTAB fully converges the inner solve in ~100 Krylov iterations.

Matching the MATLAB Result

The MATLAB code uses: 2D Lebedev 110 on a 10x2 mesh, explicit FD operator with BiCGSTAB, boundary half-volumes, P0 scattering.

The BiCGSTAB path with half-volumes reproduces 1.04192 vs MATLAB’s 1.04188 (\(4 \times 10^{-5}\) agreement). The residual difference is from floating-point details in cross-section processing.

The cleanest reference is the 1D GL BiCGSTAB result (1.03882): no pseudo-2D artifacts, well-conditioned angular quadrature, fully converged inner solve.

References

[Bailey2009] (1,2,3,4,5,6,7,8,9)

T.S. Bailey, J.E. Morel, and J.H. Chang, “A piecewise-linear finite element discretization of the diffusion equation for arbitrary polyhedral grids,” Nuclear Science and Engineering, 162:3, 2009. (Eq. 50: \(\alpha\) recursion; Eq. 53–54: WDD; Eq. 74: Morel–Montry weights.)

[MorelMontry1984]

J.E. Morel and G.R. Montry, “Analysis and elimination of the discrete-ordinates flux dip,” Transport Theory and Statistical Physics, 13:5, 1984.

[LewisMiller1984] (1,2,3)

E.E. Lewis and W.F. Miller, Jr., Computational Methods of Neutron Transport, John Wiley & Sons, 1984.

[CaseZweifel1967]

K.M. Case and P.F. Zweifel, Linear Transport Theory, Addison-Wesley, 1967.

[Lebedev1999]

V.I. Lebedev and D.N. Laikov, “A quadrature formula for the sphere of the 131st algebraic order of accuracy,” Doklady Mathematics, 59(3):477–481, 1999.

[CarlsonLathrop1965]

B.G. Carlson and K.D. Lathrop, “Transport theory – the method of discrete ordinates,” in Computing Methods in Reactor Physics, Gordon and Breach, 1968.