Method of Characteristics (MOC)

Key Facts

Read this before modifying the MOC solver.

  • Characteristic ODE: \(\frac{d\psi}{ds} + \Sigma_t \psi = Q / 4\pi\)

  • Flat-source solution: \(\bar\psi = \psi_{\text{in}} \frac{1 - e^{-\tau}}{\tau} + \frac{Q}{4\pi\Sigma_t}(1 - \frac{1-e^{-\tau}}{\tau})\)

  • Scalar flux update (Boyd Eq. 45): \(\phi = (4\pi Q + \Delta\phi / A) / \Sigma_t\)

  • Weight formula: \(4\pi \cdot \omega_a \cdot \omega_p \cdot t_s \cdot \sin(\theta_p)\) — the \(4\pi\) and \(\sin\theta_p\) are invisible to homogeneous tests (ERR-019)

  • Inverse Wigner-Seitz: pitch = mesh.edges[-1] * sqrt(pi), outermost annulus = square border

  • Tabuchi-Yamamoto polar quadrature (TY-1/2/3) × uniform azimuthal

  • Reflective BC via track linking: vertical → same direction, horizontal → reversed

  • Gotcha: homogeneous flat-source is exact regardless of weight errors — only heterogeneous multi-region exposes bugs

  • Ray-circle intersection is exact (analytical); ray-box uses parametric clipping

  • Convergence: O(ray_spacing²) for spatial, spectral for angular

  • Verification uses synthetic cross sections, not real nuclear data

Conventions

Overview

The method of characteristics (MOC) solves the multi-group eigenvalue problem in integro-differential form by tracing characteristic rays across the geometry and integrating the transport equation analytically along each ray. Unlike discrete ordinates (which discretises both angle and space on a grid) or collision probabilities (which integrates out the angular variable entirely), MOC retains the angular flux \(\psi(\mathbf{r}, \hat{\Omega}, E)\) along each characteristic while using a flat-source approximation within spatial sub-regions.

For a 2-D pin cell with concentric annuli inside a square lattice cell, MOC offers several advantages over the other ORPHEUS solvers:

  • Exact circular geometry — rays intersect circles analytically, avoiding the staircase approximation of Cartesian SN or the white-BC approximation of CP.

  • Arbitrary angular resolution — configurable azimuthal angles and ray spacing (unlike the fixed 8 directions of the pedagogical solver).

  • Reflective boundary conditions — exact, without the isotropic re-entry assumption of the CP white BC.

The implementation follows [Boyd2014] (the OpenMOC formulation) for the transport sweep and scalar flux update, and [Yamamoto2007] for the Tabuchi-Yamamoto (TY) polar quadrature. The original MOC formulation for reactor physics is due to [Askew1972]; the TY quadrature tables are also documented in [KnottYamamoto2010].

Derivation sources. The analytical eigenvalues used for verification are computed by:

  • derivations/_eigenvalue.py — kinf_homogeneous() for infinite-medium eigenvalue (all solvers share this)

  • derivations/moc.py — homogeneous cases (analytical) and heterogeneous cases (Richardson extrapolation of the MOC solver with ray-spacing refinement)

Every equation in this chapter can be verified against these scripts and the cited references.

The key equations (segment-averaged angular flux, scalar flux update weight, homogeneous consistency) are independently verified by SymPy in derivations/moc_equations.py.

Architecture

Two-Layer Mesh Pattern

The MOC solver follows the same two-layer pattern as all ORPHEUS deterministic solvers:

  1. Base geometry — Mesh1D with CoordSystem.CYLINDRICAL, typically constructed by pwr_pin_equivalent(). Stores radial cell edges, material IDs, and volumes.

  2. Augmented geometry — MOCMesh wraps the Mesh1D and an MOCQuadrature, precomputing all ray-tracing data: tracks, segments through flat-source regions, effective ray spacings, and reflective boundary condition links.

  3. Solver — MOCSolver satisfies the EigenvalueSolver protocol. solve_moc() orchestrates the calculation and returns an MoCResult.

Mesh1D (cylindrical, Wigner-Seitz)
    |
    v
MOCMesh (tracks + segments + reflective BC links)
    |
    v
MOCSolver (EigenvalueSolver protocol)
    |
    v
solve_moc() --> MoCResult

Inverse Wigner-Seitz Geometry

The base Mesh1D is constructed by pwr_pin_equivalent(), which approximates the square unit cell (side = pitch) by a cylinder of equal area:

(1)\[r_{\text{cell}} = \frac{\text{pitch}}{\sqrt{\pi}}\]

All annular regions, including the outermost (coolant), are defined within this equivalent cylinder. Material IDs are assigned via mesh.mat_ids.

MOCMesh inverts this approximation: it recovers the pitch from the outer Wigner-Seitz radius and reinterprets the outermost annular region as the square border — the physical region between the last circular boundary and the square cell walls:

(2)\[\text{pitch} = r_{\text{cell}} \cdot \sqrt{\pi} = \texttt{mesh.edges[-1]} \cdot \sqrt{\pi}\]

The circular boundaries used for ray tracing are mesh.edges[1:-1] (all radial edges except the origin and the Wigner-Seitz outer radius). The outermost region is bounded by the last circle and the cell walls.

Region areas are exact:

\[\begin{split}A_k = \begin{cases} \pi(r_{k+1}^2 - r_k^2) & k < N-1 \text{ (annular)} \\ \text{pitch}^2 - \pi r_k^2 & k = N-1 \text{ (square border)} \end{cases}\end{split}\]

By construction of the Wigner-Seitz radius, the total area is \(\sum A_k = \text{pitch}^2\), and the square-border area equals the Wigner-Seitz annular area: \(\text{pitch}^2 - \pi r_{N-2}^2 = \pi(r_\text{cell}^2 - r_{N-2}^2)\).

This convention means material IDs flow directly from mesh.mat_ids with no special treatment of the outer region. The same Mesh1D used by the CP solver can be reused by MOC.

The Transport Equation Along a Characteristic

Starting Point

The steady-state multi-group Boltzmann transport equation with isotropic scattering:

(3)\[\hat{\Omega} \cdot \nabla \psi_g(\mathbf{r}, \hat{\Omega}) + \Sigt{g}(\mathbf{r}) \, \psi_g(\mathbf{r}, \hat{\Omega}) = Q_g(\mathbf{r})\]

where \(Q_g\) is the isotropic source (fission + scattering + (n,2n)), independent of \(\hat{\Omega}\).

A characteristic (ray) is parametrised by arc length \(s\) along the direction \(\hat{\Omega}\):

\[\mathbf{r}(s) = \mathbf{r}_0 + s \, \hat{\Omega}, \quad s \geq 0\]

Along this ray, the PDE reduces to a first-order ODE:

(4)\[\frac{d\psi_g}{ds} + \Sigt{g}(s) \, \psi_g(s) = Q_g(s)\]

Flat-Source Approximation

Within each flat-source region (FSR) \(i\), both \(\Sigt{i,g}\) and \(Q_{i,g}\) are assumed spatially constant. The CP method uses the same approximation; see Flat-Source Approximation. For a ray segment of 2-D length \(\ell\) crossing region \(i\) at polar angle \(\theta_p\), the 3-D path length is \(\ell / \sin\theta_p\) and the optical thickness is:

(5)\[\tau_{g,p} = \Sigt{i,g} \cdot \frac{\ell}{\sin\theta_p}\]

The ODE (4) has the exact analytical solution:

(6)\[\psi_g^{\text{out}} = \psi_g^{\text{in}} \, e^{-\tau_{g,p}} + \frac{Q_{i,g}}{\Sigt{i,g}} \bigl(1 - e^{-\tau_{g,p}}\bigr)\]

This is the MOC attenuation formula. The first term is the uncollided flux from the incoming angular flux; the second term is the contribution from the in-region source.

The angular flux change across the segment is:

(7)\[\Delta\psi_{g,p} = \psi_g^{\text{in}} - \psi_g^{\text{out}} = \Bigl(\psi_g^{\text{in}} - \frac{Q_{i,g}}{\Sigt{i,g}}\Bigr) \bigl(1 - e^{-\tau_{g,p}}\bigr)\]

This quantity is the central building block of the MOC scalar flux update.

Note

Numerical stability. For small \(\tau < 10^{-10}\), the exponential is computed via Taylor expansion \(1 - e^{-\tau} \approx \tau(1 - \tau/2)\) to avoid catastrophic cancellation. This is implemented in the inner loop of MOCSolver.solve_fixed_source().

The Isotropic Source

For isotropic scattering and isotropic fission emission:

(8)\[Q_{i,g} = \frac{1}{4\pi} \left[ \sum_{g'=1}^{G} \Sigs{g' \to g} \, \phi_{i,g'} + 2 \sum_{g'=1}^{G} \Sigma_{2,g' \to g} \, \phi_{i,g'} + F_{i,g} \right]\]

where \(F_{i,g}\) is the fission source (from the outer iteration):

\[F_{i,g} = \frac{\chi_{i,g}}{k_{\text{eff}}} \sum_{g'=1}^{G} \nu\Sigf{i,g'} \, \phi_{i,g'}\]

The \(1/(4\pi)\) normalisation converts from volumetric source density to isotropic angular source density.

Scattering convention. SigS[0] is indexed as SigS[g_from, g_to]. The in-scatter source uses the transpose: Q += SigS.T @ phi (scattering into group \(g\) from all \(g'\)). This is the same convention used by all ORPHEUS solvers (CP, SN, MC). The \((n,2n)\) matrix follows the same convention with a factor of 2.

Angular Quadrature

MOC uses a product quadrature: azimuthal angles in the 2-D plane times polar angles from the z-axis. Contrast with the SN quadrature types, which use Gauss-Legendre (1D) and Lebedev (2D).

Azimuthal Quadrature

\(N_\varphi\) uniformly spaced angles in \([0, \pi)\):

(9)\[\varphi_m = \frac{\pi}{2 N_\varphi} + (m-1) \frac{\pi}{N_\varphi}, \quad m = 1, \ldots, N_\varphi\]

with weights \(\omega_m^a = 1 / N_\varphi\) (summing to 1).

The supplementary angle \(\varphi_m + \pi\) corresponds to the same physical track traversed in the opposite direction (the backward sweep). Thus, only \([0, \pi)\) needs to be traced explicitly; the backward direction is obtained by reversing the segment order.

Typical values: \(N_\varphi = 16, 32, 64\) for increasing accuracy.

Tabuchi-Yamamoto Polar Quadrature

The polar angle \(\theta_p\) enters the MOC only through the effective 3-D optical thickness \(\tau / \sin\theta_p\). The angular integration involves the Bickley function:

(10)\[\text{Ki}_3(\tau) = \int_0^{\pi/2} e^{-\tau/\sin\theta} \sin\theta \, d\theta\]

Yamamoto et al. [Yamamoto2007] derived optimal polar angles and weights that minimise the maximum approximation error of the Bickley function. The TY quadrature with \(N_p = 3\) polar angles per half-space achieves accuracy comparable to Gauss-Legendre with 12–16 points.

The TY weights include the \(\sin\theta\) Jacobian from the angular measure \(d\Omega = \sin\theta \, d\theta \, d\varphi\). They sum to 0.5 for one hemisphere (upper or lower); the full sphere gives \(2 \times 0.5 = 1.0\).

Tabuchi-Yamamoto quadrature tables [Yamamoto2007]

\(N_p\)

\(p\)

\(\sin\theta_p\)

\(\omega_p\)

1

1

0.798184

0.500000

2

1

0.363900

0.106427

2

2

0.899900

0.393573

3

1

0.166648

0.023117

3

2

0.537707

0.141810

3

3

0.932954

0.335074

These tables are hardcoded in MOCQuadrature from [Yamamoto2007] Table 2. Verified by test_moc_quadrature.py::test_ty3_values_match_published.

The combined angular weight normalisation satisfies:

\[2 \sum_{m=1}^{N_\varphi} \omega_m^a \sum_{p=1}^{N_p} \omega_p = 2 \cdot 1 \cdot 0.5 = 1\]

The factor of 2 accounts for the two hemispheres (upper = forward sweep, lower = backward sweep). This is verified by test_moc_quadrature.py::test_combined_weight_normalisation.

Ray Tracing

Track Generation

For each azimuthal angle \(\varphi_m\), a family of parallel rays is laid across the square cell \([0, \text{pitch}]^2\) with perpendicular spacing \(t_s\).

The perpendicular coordinate of a ray is:

\[t = -x \sin\varphi_m + y \cos\varphi_m\]

The range of \(t\) values that intersect the cell is determined by projecting the four cell corners onto the perpendicular axis. Rays are placed at \(t_k = t_{\min} + (k + \tfrac{1}{2}) t_s^{\text{eff}}\), where:

(11)\[t_s^{\text{eff}} = \frac{t_{\max} - t_{\min}}{n_{\text{rays}}}, \quad n_{\text{rays}} = \left\lceil \frac{t_{\max} - t_{\min}}{t_s} \right\rceil\]

The effective spacing is chosen so that rays exactly cover the cell width at each angle, with no gaps.

Ray-Circle Intersection

A ray starting at \((x_0, y_0)\) with direction \((\cos\varphi, \sin\varphi)\) is parametrised as:

\[x(s) = x_0 + s \cos\varphi, \quad y(s) = y_0 + s \sin\varphi\]

The intersection with a circle of radius \(R_k\) centred at \((c_x, c_y)\) satisfies:

(12)\[s^2 + 2bs + c = 0\]

where:

\[b = (x_0 - c_x)\cos\varphi + (y_0 - c_y)\sin\varphi, \quad c = (x_0 - c_x)^2 + (y_0 - c_y)^2 - R_k^2\]

The discriminant \(\Delta = b^2 - c\) determines whether the ray intersects the circle:

  • \(\Delta < 0\) — ray misses the circle (no intersection)

  • \(\Delta \geq 0\) — two intersections at \(s_{1,2} = -b \mp \sqrt{\Delta}\)

Only intersections within the cell (between entry and exit parameters) are retained.

Segment Construction

For a given ray through the cell:

  1. Compute all intersection parameters: cell-wall entry/exit (from _ray_box_intersections()) and circle crossings (from _ray_circle_intersections()).

  2. Sort all parameters and remove near-duplicates.

  3. For each consecutive pair \((s_a, s_b)\), the segment length is \(\ell = s_b - s_a\). The midpoint determines which FSR the segment belongs to (by computing distance from pin centre).

  4. Build the segment list: (region_id, length) tuples.

Region assignment. A point at distance \(d\) from the pin centre belongs to region \(k\) where \(r_{k-1} < d \leq r_k\) for \(k < N-1\), or region \(N-1\) (square border) if \(d > r_{N-2}\).

Verified by test_moc_ray_tracing.py:

  • test_ray_circle_chord_length: chord length = \(2\sqrt{R^2 - d^2}\) for impact parameter \(d\)

  • test_trace_three_regions: correct segment count (5) and region ordering for 3 concentric annuli

  • test_volume_conservation: sum of (segment length × ray spacing) approximates region areas to < 5% at ray_spacing=0.02

Reflective Boundary Conditions

For a pin cell with reflective BCs, each ray that exits through a cell wall re-enters as a reflected ray. The outgoing angular flux from one track becomes the incoming flux for the reflected track.

Reflection Rules

For azimuthal angles \(\varphi \in [0, \pi)\), all forward rays have \(\sin\varphi \geq 0\) (upward or horizontal). Two types of reflection occur:

Reflection rules for angles in \([0, \pi)\)

Sweep

Exit wall

Reflected angle

Target direction

Forward

Vertical (L/R)

\(\pi - \varphi\)

Forward

Forward

Horizontal (T)

\(\pi - \varphi\)

Backward

Backward

Horizontal (B)

\(\pi - \varphi\)

Forward

Backward

Vertical (L/R)

\(\pi - \varphi\)

Backward

Key distinction: Both reflection types map the azimuthal angle to \(\pi - \varphi\). The difference is whether the reflected ray continues in the same traversal direction (forward → forward at the reflected angle) or reverses (forward → backward at the reflected angle).

  • Vertical wall (right \(x = P\), left \(x = 0\)): reverses the \(x\)-component of \(\hat{\Omega}\). Since \((-\cos\varphi, \sin\varphi) = (\cos(\pi-\varphi), \sin(\pi-\varphi))\), this is the forward direction at the reflected angle.

  • Horizontal wall (top \(y = P\), bottom \(y = 0\)): reverses the \(y\)-component. Since \((\cos\varphi, -\sin\varphi) = -(\cos(\pi-\varphi), \sin(\pi-\varphi))\), this is the backward direction at the reflected angle.

Summary rule: vertical reflections preserve forward/backward; horizontal reflections flip it. In the code, _is_vertical() tests this and MOCMesh._link_tracks() uses it to set fwd_link_fwd and bwd_link_fwd.

Track Linking

After generating all tracks, MOCMesh._link_tracks() assigns four link fields per track:

  • fwd_link, fwd_link_fwd: where the forward sweep’s outgoing flux goes (target track index, and whether it feeds into the target’s forward or backward entry)

  • bwd_link, bwd_link_fwd: same for the backward sweep

The target track is found by matching the exit point of one track to the entry point (for forward targets) or exit point (for backward targets) of a track at the reflected azimuthal angle. The closest match is used.

Verified by test_moc_verification.py::TestL0GeometricInvariants:: test_reflective_links_form_cycles: following forward links from any track must return to the starting track after a finite number of reflections (closed cycle).

Scalar Flux Update

Derivation from First Principles

The scalar flux in FSR \(i\) for energy group \(g\) is the angular integral of the segment-averaged angular flux over all directions:

(13)\[\phi_{i,g} = \frac{1}{A_i} \int_{4\pi} \int_{A_i} \bar{\psi}_g(\mathbf{r}, \hat{\Omega}) \, dA \, d\Omega\]

where \(\bar{\psi}\) is the angular flux averaged along the characteristic segment through region \(i\).

Segment-averaged angular flux. The average angular flux along a segment of 3-D path length \(L = \ell / \sin\theta_p\) is obtained by integrating the ODE solution (6) over the segment:

\[\bar{\psi}_{k,g,p} = \frac{1}{L} \int_0^{L} \psi_g(s) \, ds = \frac{1}{L} \int_0^{L} \left[ \psi_g^{\text{in}} \, e^{-\Sigt{} s} + \frac{Q}{\Sigt{}} (1 - e^{-\Sigt{} s}) \right] ds\]

Evaluating the integral term by term:

\[\int_0^L \psi^{\text{in}} e^{-\Sigt{} s} \, ds = \frac{\psi^{\text{in}}}{\Sigt{}} (1 - e^{-\tau})\]
\[\int_0^L \frac{Q}{\Sigt{}} (1 - e^{-\Sigt{} s}) \, ds = \frac{Q}{\Sigt{}} \left[ L - \frac{1-e^{-\tau}}{\Sigt{}} \right]\]

where \(\tau = \Sigt{} L = \Sigt{i,g} \ell_k / \sin\theta_p\). Combining and dividing by \(L\):

\[\bar{\psi} = \frac{Q}{\Sigt{}} + \frac{1}{\tau} \left( \psi^{\text{in}} - \frac{Q}{\Sigt{}} \right) (1 - e^{-\tau})\]

Recognising \(\Delta\psi = (\psi^{\text{in}} - Q/\Sigt{})(1 - e^{-\tau})\) from (7):

(14)\[\boxed{ \bar{\psi}_{k,g,p} = \frac{Q_{i,g}}{\Sigt{i,g}} + \frac{\Delta\psi_{k,g,p}}{\tau_{k,g,p}} }\]

This result is exact under the flat-source approximation. The first term is the asymptotic angular flux in the region; the second term is the correction from non-equilibrium boundary flux.

Physical interpretation. In an infinite homogeneous medium, \(\psi^{\text{in}} = Q/\Sigt{}\) everywhere, so \(\Delta\psi = 0\) and \(\bar\psi = Q/\Sigt{}\). The correction term carries the spatial coupling between heterogeneous regions — it is the only term that distinguishes fuel from coolant. This is why the angular integration weight (which multiplies \(\Delta\psi\)) is invisible to homogeneous tests: the quantity it multiplies is identically zero.

Substituting (14) into (13):

\[\phi_{i,g} = \frac{1}{A_i} \int_{4\pi} \sum_{k \in i} t_s \, \ell_k \left[ \frac{Q_{i,g}}{\Sigt{i,g}} + \frac{\Delta\psi_{k,g,p} \, \sin\theta_p}{\Sigt{i,g} \ell_k} \right] d\Omega\]

where the \(\sin\theta_p\) factor arises from \(\tau = \Sigt{} \ell / \sin\theta_p\) in the denominator of (14): \(\Delta\psi / \tau = \Delta\psi \sin\theta_p / (\Sigt{} \ell)\).

The first term uses \(\sum_k t_s \ell_k \approx A_i\) (the tracks at each azimuthal angle tile the region area), and \(\int_{4\pi} d\Omega = 4\pi\):

\[\phi_{i,g} = \frac{4\pi Q_{i,g}}{\Sigt{i,g}} + \frac{1}{A_i \, \Sigt{i,g}} \int_{4\pi} \sin\theta_p \sum_{k \in i} t_s \, \Delta\psi_{k,g,p} \, d\Omega\]

Discretising the angular integral with the product quadrature (azimuthal \(\omega_m^a\), polar \(\omega_p\), factor of \(4\pi\) for the full sphere, and factor of 2 for forward + backward sweeps absorbed into the summation over both sweep directions):

(15)\[\boxed{ \phi_{i,g} = \frac{1}{\Sigt{i,g}} \left[ 4\pi \, Q_{i,g} + \frac{1}{A_i} \sum_{m,p,k \in i} 4\pi \, \omega_m^a \, \omega_p \, t_s \, \sin\theta_p \, \Delta\psi_{k,g,p} \right] }\]

This is [Boyd2014] Equation 45 (with our weight normalisation convention).

Implementation. In MOCSolver.solve_fixed_source(), the accumulator delta_phi[i, g] collects the weighted \(\Delta\psi\) contributions during the sweep. The weight per segment is:

weight = 4.0 * np.pi * omega_a * omega_p * ts * sin_p

After the sweep, the scalar flux is updated as:

phi[i, g] = (4*pi*Q[i,g] + delta_phi[i,g] / area[i]) / sig_t[i,g]

ERR-019 — Missing weight factor (caught during development)

Bug: The initial implementation used weight = omega_a * omega_p * ts, missing both the \(4\pi\) full-sphere factor and the \(\sin\theta_p\) factor from (15).

Impact: Heterogeneous keff was completely wrong: MOC gave 1.344 vs CP reference of 0.902 for a 2-region pin cell. All three homogeneous tests passed to machine precision because \(\Delta\psi = 0\) when the angular flux is spatially uniform.

How it hid: For homogeneous material, \(\psi_{\text{in}} = Q / \Sigt{}\) everywhere, so \(\Delta\psi = 0\). The scalar flux reduces to \(4\pi Q / \Sigt{}\) regardless of the weight. This is the fundamental degeneracy: weights are invisible to homogeneous tests.

What caught it: The first heterogeneous cross-verification against the CP solver revealed a 0.44 discrepancy. Re-deriving the weight formula from (13) identified the two missing factors.

Lesson: Always verify the transport sweep with a heterogeneous problem before trusting the homogeneous result. The weight formula should be derived from first principles, not guessed by analogy. See tests/l0_error_catalog.md ERR-019.

Eigenvalue Update

For reflective boundary conditions (zero leakage), the eigenvalue is:

(16)\[k_{\text{eff}} = \frac{\text{production}}{\text{absorption}} = \frac{\sum_i \sum_g (\nSigf{i,g} + 2 \Sigma_{2,i,g}^{\text{out}}) \phi_{i,g} \, A_i} {\sum_i \sum_g \Siga{i,g} \, \phi_{i,g} \, A_i}\]

where \(\Sigma_{2,i,g}^{\text{out}} = \sum_{g'} \Sigma_{2,g \to g'}\) is the total (n,2n) transfer out of group \(g\).

Derivation. Each (n,2n) reaction absorbs one neutron and produces two. In the eigenvalue balance:

  • Production (numerator): \(\nSigf{}\) (fission) + \(2\Sigma_2^{\text{out}}\) (two neutrons from each (n,2n))

  • Absorption (denominator): \(\Siga{} = \Sigma_c + \Sigma_f + \Sigma_L + \Sigma_2^{\text{out}}\) (capture + fission + leakage + (n,2n) removal)

The (n,2n) appears in BOTH numerator (2×, production) and denominator (1×, removal). The net contribution per (n,2n) reaction is +1 neutron.

When \(\Sigma_2 = 0\), (16) reduces to the standard \(k = \nSigf{}\phi A / \Siga{}\phi A\).

Implemented in MOCSolver.compute_keff(). Verified by test_moc_verification.py::TestL0N2nReaction::test_n2n_1g_analytical_keff.

Power Iteration

The MOC solver plugs into the generic power_iteration() loop via the EigenvalueSolver protocol:

  1. Fission source (MOCSolver.compute_fission_source()): \(F_{i,g} = \chi_{i,g} \cdot (\nSigf{i} \cdot \phi_i) / k\)

  2. Transport sweep (MOCSolver.solve_fixed_source()): Build \(Q\) from (8), sweep all tracks (forward + backward), accumulate \(\Delta\psi\) into delta_phi, update \(\phi\) from (15).

    The method performs n_inner_sweeps transport sweeps per outer iteration to converge the boundary angular fluxes. Within each sweep, the scattering source is updated from the latest flux.

    Boundary flux persistence: The angular fluxes at track entry/exit points (_fwd_bflux, _bwd_bflux) persist between outer iterations. This allows the reflective BCs to converge progressively without explicit cyclic tracking.

  3. Eigenvalue update (MOCSolver.compute_keff()): (16).

  4. Convergence (MOCSolver.converged()): \(|\Delta k| < \texttt{keff\_tol}\) and \(\|\Delta\phi\| / \|\phi\| < \texttt{flux\_tol}\).

The power iteration algorithm is shared with all ORPHEUS eigenvalue solvers; see The Power Iteration Algorithm in the homogeneous theory for the general formulation.

Cross-Section Data Layout

Cross sections flow through two paths:

  1. Mixture — per-material arrays for \(\Sigt{}\), \(\Siga{}\), \(\nSigf{}\), \(\chi\); sparse matrices for \(\Sigma_s\) (SigS[0]) and \(\Sigma_2\) (Sig2).

  2. Per-region assembly in MOCSolver.__init__() — loops over MOCMesh.region_mat_ids to build (n_regions, ng) arrays.

Indexing: sig_t[i, g] is the total cross section in FSR \(i\), group \(g\). Group 0 = fastest; group \(G-1\) = thermal. Region 0 = innermost annulus; region \(N-1\) = square border.

Scattering convention: SigS[g_from, g_to]. The in-scatter source uses the transpose: Q += SigS.T @ phi (same as CP, SN, MC).

Verification

102 Tests Across Four Levels

Level

File(s)

Tests

Coverage

L0

test_moc_quadrature.py

24

Weight sums, TY values, shapes, validation

L0

test_moc_ray_tracing.py

20

Ray-circle intersection, region ID, segments, volume, links

L0

test_moc_verification.py

27

Single-track attenuation, equilibrium flux, fission-only, (n,2n), scatter isolation, geometric invariants, protocol compliance, volume tracking, boundary conditions

L1

test_moc.py

6

Homogeneous eigenvalue (1G/2G/4G), heterogeneous (slow)

L1

test_moc_properties.py

4

Particle balance, positivity, flux consistency, thermal depression

L1

test_moc_verification.py

13

Eigenvalue + flux ratio, heterogeneous + monotonicity, particle balance (2G), flux positivity (all groups), material sensitivity

L2

test_moc_verification.py

4

Ray spacing convergence, azimuthal convergence, polar convergence

XV

test_moc_verification.py

2

MOC vs CP cross-verification (slow)

Run the full suite (excluding slow tests):

pytest tests/test_moc_quadrature.py tests/test_moc_ray_tracing.py \
       tests/test_moc.py tests/test_moc_properties.py \
       tests/test_moc_verification.py -v -k "not slow"

Homogeneous Infinite Medium

For homogeneous geometry with reflective BCs, the flux is spatially flat and \(k_{\text{eff}} = \lambda_{\max}(A^{-1}F)\).

Groups

\(k_\infty\)

Error (8 azi, TY-3)

Tolerance

1

1.5000

\(< 10^{-15}\) (exact)

\(< 10^{-4}\)

2

1.8750

\(\sim 10^{-6}\)

\(< 10^{-4}\)

4

1.4878

\(\sim 10^{-6}\)

\(< 10^{-4}\)

The 1-group case is exact because \(k = \nu\Sigma_f / \Sigma_a\) is independent of the angular flux distribution. The multi-group cases converge through the power iteration.

Why Homogeneous Verification Is Necessary But Insufficient

For a homogeneous medium, the flat-source approximation is exact: the source is genuinely spatially flat. Consequences:

  1. \(\Delta\psi = 0\) everywhere (all boundary fluxes equal \(Q / \Sigt{}\)).

  2. Angular integration weights cancel in the eigenvalue ratio.

  3. The transport sweep contributes nothing — \(\phi = 4\pi Q / \Sigt{}\).

This means homogeneous tests cannot detect:

  • Wrong angular integration weights (ERR-019)

  • Wrong scattering convention (SigS vs SigS.T)

  • Wrong boundary condition linking

  • Wrong ray spacing / volume conservation

Rule: Every transport solver must be verified on heterogeneous multi-group problems. Homogeneous success gives false confidence.

Heterogeneous Cross-Verification

For a 2-region cylindrical pin cell (fuel + coolant), the MOC eigenvalue is compared against the CP solver on the same geometry:

Groups

MOC \(k_{\text{eff}}\)

CP \(k_{\text{eff}}\)

Gap

2

0.6078

0.6072

0.001

4

0.5399

0.5394

0.001

The ~0.1% gap is consistent with the white-BC (CP) vs reflective-BC (MOC) approximation difference. This is verified by test_moc_verification.py::TestXVCrossVerification.

Convergence Properties

All convergence studies use a 2-region 2G pin cell (fuel + coolant, \(r_{\text{fuel}} = 0.5\) cm, pitch = 2.0 cm) with materials A (fissile) and B (scatterer) from the verification library.

Ray spacing convergence (\(N_\varphi = 32\), TY-3):

\(t_s\) (cm)

Tracks

\(k_{\text{eff}}\)

\(|\Delta k|\)

Ratio

0.16

528

0.608256

0.08

1036

0.608021

\(2.36 \times 10^{-4}\)

0.04

2048

0.608129

\(1.09 \times 10^{-4}\)

2.2

0.02

4092

0.608038

\(9.10 \times 10^{-5}\)

1.2

0.01

8168

0.608029

\(9.23 \times 10^{-6}\)

9.9

The eigenvalue converges with decreasing ray spacing. The convergence is not perfectly monotonic (non-monotone behaviour between 0.08 and 0.04 reflects the discrete nature of track placement on the annular boundary), but the refinement from 0.02 to 0.01 shows a clear order-of-magnitude improvement.

Azimuthal convergence (\(t_s = 0.02\) cm, TY-3):

\(N_\varphi\)

Tracks

\(k_{\text{eff}}\)

\(|\Delta k|\)

4

524

0.605797

8

1028

0.607311

\(1.51 \times 10^{-3}\)

16

2048

0.607959

\(6.47 \times 10^{-4}\)

32

4092

0.608038

\(7.96 \times 10^{-5}\)

64

8188

0.608057

\(1.89 \times 10^{-5}\)

The azimuthal convergence shows a clean reduction factor of ~2–8× per doubling of \(N_\varphi\). The error at \(N_\varphi = 64\) is dominated by the ray spacing, not the angular discretisation.

Polar convergence (\(N_\varphi = 16\), \(t_s = 0.03\) cm):

\(N_p\)

\(k_{\text{eff}}\)

\(|\Delta k|\) from TY-1

Quadrature points

TY-1

0.607538

1 per half-space

TY-2

0.607632

\(9.4 \times 10^{-5}\)

2 per half-space

TY-3

0.607790

\(2.5 \times 10^{-4}\)

3 per half-space

The polar convergence is weak: TY-1 to TY-3 gives only a \(\sim 2.5 \times 10^{-4}\) shift. This is because the TY quadrature is specifically optimised for the Bickley-function integrals that arise in MOC, so even TY-1 (a single polar angle) gives a surprisingly good approximation. Recommendation: TY-3 is the default; increasing to TY-2 or TY-1 is acceptable for quick exploratory runs.

Investigation History: ERR-019 (Weight Factor)

This section documents the development history of the scalar flux update formula to prevent future sessions from repeating the same error.

Symptoms

  1. All three homogeneous eigenvalue tests passed to machine precision (1G exact, 2G/4G \(\sim 10^{-6}\)).

  2. First heterogeneous test (2-region, 1G): MOC gave \(k = 1.344\), CP gave \(k = 0.902\).

  3. The discrepancy was 0.44 — far too large to be a discretisation effect.

The Root Cause

The delta_phi accumulation weight was omega_a * omega_p * ts, missing two factors:

  1. \(4\pi\) — the angular integration normalisation from \(\int_{4\pi} d\Omega = 4\pi\).

  2. \(\sin\theta_p\) — the 2-D to 3-D projection factor arising from the segment-averaged angular flux formula (14).

Why It Was Invisible to Homogeneous Tests

For homogeneous material with converged boundary fluxes:

\[\psi_{\text{in}} = \frac{Q}{\Sigt{}} \implies \Delta\psi = 0 \implies \delta\phi = 0\]

The scalar flux reduces to \(\phi = 4\pi Q / \Sigt{}\) regardless of the weight factor. This is the fundamental degeneracy: weights are invisible when \(\Delta\psi = 0\).

What Caught It

The first heterogeneous cross-verification against the CP solver showed a 0.44 discrepancy. Re-deriving the weight from (13) through (15) identified the two missing factors. After the fix, MOC agreed with CP to within \(1.2 \times 10^{-3}\) (consistent with the white-BC vs reflective-BC approximation difference).

Lesson

The angular integration weight in MOC contains factors that cancel for spatially uniform solutions:

  • \(4\pi\) cancels because the isotropic source already divides by \(4\pi\)

  • \(\sin\theta_p\) cancels because \(\Delta\psi = 0\)

Always verify the weight formula against a heterogeneous problem before trusting the homogeneous result. Derive the weight from first principles and verify each factor independently.

Design Decisions and Alternatives

Why Flat-Source (Not Linear-Source)

The flat-source approximation assumes \(Q(\mathbf{r}) = Q_i\) within each FSR. The alternative is the linear-source approximation, which expands the source as \(Q(\mathbf{r}) = Q_i + \nabla Q_i \cdot (\mathbf{r} - \mathbf{r}_i)\), requiring a spatial gradient per region.

Trade-off: Linear-source MOC converges at \(O(h^3)\) vs \(O(h^2)\) for flat-source (where \(h\) is the characteristic width of an FSR), allowing coarser meshes for the same accuracy. However, it requires gradient estimation per region per group per iteration, which doubles the storage and significantly complicates the implementation.

Decision: For a pin cell with ~3–20 FSRs (annular rings), the flat-source approximation with fine ray spacing is sufficient. The region boundaries are defined by material interfaces (fuel/clad/coolant), and the flat-source is exact within each material zone. Linear-source would only help for coarse radial subdivision within a material zone (e.g., subdividing the fuel pellet into fewer rings).

When to revisit: If the solver is extended to multi-assembly or whole-core geometry (where FSR counts explode), linear-source becomes essential for performance.

Why TY Quadrature (Not Gauss-Legendre)

The TY quadrature is optimised for the Bickley function \(\text{Ki}_3(\tau)\) that appears when the flat-source MOC angular flux is integrated over polar angle. In contrast, Gauss-Legendre (GL) quadrature optimises polynomial integration, which is not the relevant function form.

For the same number of polar points, TY gives ~2–5× smaller angular error than GL ([Yamamoto2007] Table 3). TY-3 (3 points) matches GL with 12–16 points on typical LWR pin-cell problems.

Connection to CP: The collision probability method uses the same Bickley function \(\text{Ki}_4(\tau)\) in its cylindrical kernel. The CP method integrates over all rays analytically (via the \(y\)-quadrature in CPMesh), while MOC traces individual rays and sums numerically. Both require accurate approximation of Bickley integrals; CP handles this with dense Gauss-Legendre quadrature over the impact parameter, while MOC uses TY quadrature over the polar angle. The TY approach is more efficient for MOC because each polar angle contributes to all tracks simultaneously, whereas CP evaluates the kernel per-quadrature-point per-region pair.

Why Inverse Wigner-Seitz (Not Direct Square)

An alternative to the inverse Wigner-Seitz approach would be to define the MOC geometry directly on the square cell (as the old pedagogical solver did with MoCGeometry, now removed). This would require a separate geometry class for the square cell + annuli combination.

Advantage of the current approach: The same Mesh1D created by pwr_pin_equivalent() is reused by CP, SN, MC, and now MOC. The pitch is encoded in the outer Wigner-Seitz radius, and each solver’s augmented geometry recovers whatever it needs. No duplication.

Limitation: The Wigner-Seitz approximation is exact for the cell area but not for the shape. Corner regions of the square cell (beyond the Wigner-Seitz circle but within the square) have slightly different neutronics than the cylindrical approximation suggests. However, for reflective BC pin-cell calculations, this is standard practice and the error is negligible for typical LWR geometries.

Gotchas and Subtleties

Void Regions (\(\Sigt{} = 0\))

If a region has zero total cross section (vacuum), the attenuation formula (6) degenerates: \(\psi^{\text{out}} = \psi^{\text{in}}\) (no interaction), \(\Delta\psi = 0\), \(Q/\Sigt{}\) is undefined. The code guards against this: segments in void regions are skipped (no contribution to delta_phi), and the scalar flux is set to zero.

In the ORPHEUS verification library, no test material has \(\Sigt{} = 0\), so this code path is exercised only by the protocol compliance test test_moc_verification.py::TestL0ProtocolCompliance.

Effective Spacing vs Requested Spacing

The requested ray spacing \(t_s\) is adjusted per azimuthal angle so that rays exactly tile the cell width: \(t_s^{\text{eff}} = (t_{\max} - t_{\min}) / n_{\text{rays}}\) (see (11)). The effective spacing is always \(\leq t_s\) but may differ slightly between angles.

This means the actual spatial resolution varies with azimuthal angle. Near-horizontal rays (\(\varphi \approx 0\)) see a narrower perpendicular extent and thus fewer rays than near-diagonal rays (\(\varphi \approx \pi/4\)). The azimuthal weight \(\omega_m^a = 1/N_\varphi\) treats all angles equally, so the spatial integration is slightly less accurate for some angles. This is standard practice and does not affect convergence order.

Flat-Source Breakdown Regime

The flat-source approximation breaks down when the source varies significantly within a single FSR. This occurs when:

  • An FSR is optically thick: \(\Sigt{} \cdot \ell \gg 1\) (strong spatial attenuation within the region)

  • The flux gradient is steep (near material interfaces)

For a typical LWR pin cell with fuel radius ~0.5 cm and thermal \(\Sigt{} \sim 0.5\) cm-1, the optical thickness of the fuel region is \(\sim 0.25\) (optically thin), so the flat-source approximation is excellent. For fast-spectrum systems with \(\Sigt{} \sim 0.1\) cm-1, the optical thickness is even smaller.

The flat-source error can be reduced by subdividing regions into thinner annuli (more cells in the Mesh1D). The CP solver’s pwr_pin_equivalent() default of 10 fuel + 3 clad + 7 coolant sub-cells is adequate for most applications.

Open Improvements

See 03.Method.Of.Characteristics/IMPROVEMENTS.md for the full tracker. Key open items:

  • MC-20260406-006: Sphinx theory chapter (this document) — move from IMPL to DONE

  • MC-20260406-007: Vectorise the Python transport sweep (5-deep loop) for performance

  • MC-20260406-008: Tighten heterogeneous test tolerances with regenerated Richardson references

  • MC-20260406-009: Cyclic track linking for single-sweep BC convergence

  • MC-20260406-010: CMFD coarse-mesh acceleration

References

[Boyd2014] (1,2)

W.R.D. Boyd, S. Shaner, L. Li, B. Forget, K. Smith, “The OpenMOC Method of Characteristics Neutral Particle Transport Code,” Annals of Nuclear Energy, 68, 43–52, 2014.

[Yamamoto2007] (1,2,3,4,5)

A. Yamamoto, M. Tabuchi, N. Sugimura, T. Ushio, M. Mori, “Derivation of Optimum Polar Angle Quadrature Set for the Method of Characteristics Based on Approximation Error for the Bickley Function,” Journal of Nuclear Science and Technology, 44(2), 129–136, 2007.

[Askew1972]

J.R. Askew, “A Characteristics Formulation of the Neutron Transport Equation in Complicated Geometries,” AEEW-M 1108, Winfrith, 1972.

[KnottYamamoto2010]

D. Knott and A. Yamamoto, “Lattice Physics Computations,” Chapter 9 in Handbook of Nuclear Engineering, Springer, 2010.