Discrete Ordinates Solvers

Solver

Unified SN (Discrete Ordinates) eigenvalue solver.

Supports 1D (ny=1, GL quadrature) and 2D (Lebedev quadrature) with selectable inner solver strategy and scattering anisotropy order.

The transport sweep uses diamond-difference spatial discretization: - 1D: cumulative-product recurrence (~ms) - 2D: wavefront parallelism along anti-diagonals

Boundary conditions are reflective on all sides (infinite lattice).

See also

Discrete Ordinates Method (SN) — Key Facts, equations, gotchas.

class orpheus.sn.solver.SNResult(keff, keff_history, angular_flux, scalar_flux, geometry, quadrature, eg, elapsed_seconds)[source]

Bases: object

Results of an SN transport calculation.

The primary output is the angular flux (the direct solution of the SN equations). Scalar flux is derived by quadrature integration.

Parameters:
keff: float
keff_history: list[float]
angular_flux: ndarray
scalar_flux: ndarray
geometry: Mesh1D | Mesh2D
quadrature: AngularQuadrature
eg: ndarray
elapsed_seconds: float
class orpheus.sn.solver.SNSolver(materials, sn_mesh, inner_solver='source_iteration', scattering_order=0, keff_tol=1e-07, flux_tol=1e-06, max_inner=200, inner_tol=1e-08)[source]

Bases: object

Unified SN eigenvalue solver satisfying the EigenvalueSolver protocol.

Parameters:

  • materials (dict mapping material ID to Mixture.)

  • sn_mesh (SNMesh — augmented geometry (wraps Mesh1D or Mesh2D with) – precomputed streaming stencil).

  • inner_solver ("source_iteration" or "bicgstab".)

  • scattering_order (int — Legendre order for scattering (0 = P0).)

  • keff_tol (outer iteration convergence.)

  • flux_tol (outer iteration convergence.)

  • max_inner (inner iteration parameters.)

  • inner_tol (inner iteration parameters.)

sig_s: dict[int, list[ndarray]]
sig2: dict[int, ndarray]
sig_s0: dict[int, ndarray]
initial_flux_distribution()[source]

Initial scalar flux guess: ones(nx, ny, ng).

Return type:

ndarray

compute_fission_source(flux_distribution, keff)[source]

Fission source: χ · (νΣ_f · φ) / k.

Parameters:
Return type:

ndarray

solve_fixed_source(fission_source, flux_distribution)[source]

Solve the within-group transport equation for given fission source.

Returns updated scalar flux (nx, ny, ng).

Parameters:
Return type:

ndarray

compute_keff(flux_distribution)[source]

k = production / absorption (volume-weighted).

Parameters:

flux_distribution (ndarray)

Return type:

float

converged(keff, keff_old, flux_distribution, flux_old, iteration)[source]
Parameters:
Return type:

bool

orpheus.sn.solver.solve_sn(materials, mesh, quadrature, inner_solver='source_iteration', scattering_order=0, max_outer=500, keff_tol=1e-07, flux_tol=1e-06, max_inner=200, inner_tol=1e-08)[source]

Solve the multi-group SN eigenvalue problem.

Parameters:

  • materials (dict mapping material ID to Mixture.)

  • mesh (Mesh1D or Mesh2D (base geometry).)

  • quadrature (AngularQuadrature (GaussLegendre1D or LebedevSphere).)

  • inner_solver ("source_iteration" (default) or "bicgstab".)

  • scattering_order (Legendre order for scattering (0 = P0).)

  • max_outer (maximum outer (power) iterations.)

  • keff_tol (outer convergence.)

  • flux_tol (outer convergence.)

  • max_inner (inner solver parameters.)

  • inner_tol (inner solver parameters.)

Return type:

SNResult

Geometry

Augmented geometry for SN discrete ordinates transport.

SNMesh wraps a Mesh1D or Mesh2D and precomputes the coordinate-specific streaming stencil used by the transport sweep.

Three coordinate systems are supported: Cartesian (1D/2D), spherical (1D), and cylindrical (1D). Curvilinear geometries precompute angular redistribution coefficients (\(\alpha\)), the geometry factor \(\Delta A/w\), and Morel–Montry angular closure weights.

class orpheus.sn.geometry.SNMesh(mesh, quadrature)[source]

Bases: object

Augmented geometry for the discrete ordinates method.

Wraps a Mesh1D or Mesh2D and precomputes the streaming stencil (diamond-difference coefficients that depend only on geometry and angular quadrature, not on cross sections).

For Cartesian geometry the stencil stores:

  • streaming_x[n, i] = \(2|\mu_{x,n}| / \Delta x_i\)

  • streaming_y[n, j] = \(2|\mu_{y,n}| / \Delta y_j\)

For future curvilinear geometries, additional curvature terms (\(\alpha_n / r_i\)) will be stored in self.curvature.

Parameters:

  • mesh (Mesh1D or Mesh2D) – Base geometry.

  • quadrature (AngularQuadrature) – Angular quadrature (Gauss–Legendre, Lebedev, etc.).

property volumes: ndarray

Cell volumes, shape (nx, ny).

property is_1d: bool

True if this is a 1-D mesh (ny == 1).

Angular Quadrature

Angular quadrature for SN transport.

Provides a protocol and two implementations: - GaussLegendre1D: for 1D slab problems (mu on [-1,1], weights sum to 2) - LebedevSphere: for 2D/3D problems (directions on unit sphere, weights sum to 4π)

The solver uses 1/sum(weights) as the isotropic normalization factor, making it quadrature-agnostic.

class orpheus.sn.quadrature.AngularQuadrature(*args, **kwargs)[source]

Bases: Protocol

Contract for any angular quadrature usable by the SN solver.

mu_x: ndarray
mu_y: ndarray
weights: ndarray
N: int
reflection_index(axis)[source]

Index array: ref[n] = partner of ordinate n reflected in axis.

Parameters:

axis (str)

Return type:

ndarray

spherical_harmonics(L)[source]

(N, L+1, 2L+1) real spherical harmonics Y[n, l, l+m].

Parameters:

L (int)

Return type:

ndarray

class orpheus.sn.quadrature.GaussLegendre1D(mu_x, mu_y, weights, N)[source]

Bases: object

Gauss-Legendre quadrature on [-1, 1] for 1D slab transport.

mu_x = GL points, mu_y = 0. Weights sum to 2.

Parameters:
mu_x: ndarray
mu_y: ndarray
weights: ndarray
N: int
property mu: ndarray

Alias for mu_x (the 1D direction cosines).

classmethod create(n_ordinates=16)[source]

Build N-point GL quadrature (must be even for SN).

Parameters:

n_ordinates (int)

Return type:

GaussLegendre1D

reflection_index(axis)[source]

GL is symmetric: partner of i is N-1-i for x-reflection.

Parameters:

axis (str)

Return type:

ndarray

spherical_harmonics(L)[source]

1D harmonics: Y_0^0=1, Y_1^0=μ_x (only x-component in 1D).

Parameters:

L (int)

Return type:

ndarray

class orpheus.sn.quadrature.LebedevSphere(mu_x, mu_y, mu_z, weights, N, _ref_x, _ref_y, _ref_z)[source]

Bases: object

Lebedev quadrature on the unit sphere for 2D/3D transport.

Weights sum to 4π. Directions cover the full sphere.

Parameters:
mu_x: ndarray
mu_y: ndarray
mu_z: ndarray
weights: ndarray
N: int
classmethod create(order=17)[source]

Build Lebedev quadrature from scipy.

Parameters:

order (int)

Return type:

LebedevSphere

reflection_index(axis)[source]
Parameters:

axis (str)

Return type:

ndarray

spherical_harmonics(L)[source]

(N, L+1, 2L+1) real spherical harmonics for Lebedev ordinates.

Parameters:

L (int)

Return type:

ndarray

class orpheus.sn.quadrature.LevelSymmetricSN(mu_x, mu_y, mu_z, weights, N, _ref_x, _ref_y, _ref_z, n_levels, level_indices, level_mu)[source]

Bases: object

Level-symmetric S_N quadrature on the unit sphere.

Standard triangular quadrature with N/2 μ-levels per hemisphere. Provides the level_indices structure required by the cylindrical SN sweep for azimuthal redistribution.

Weights sum to 4π.

Parameters:
mu_x: ndarray
mu_y: ndarray
mu_z: ndarray
weights: ndarray
N: int
n_levels: int
level_indices: list[ndarray]
level_mu: ndarray
classmethod create(sn_order=4)[source]

Build S_N level-symmetric quadrature of given order.

Parameters:

sn_order (int)

Return type:

LevelSymmetricSN

reflection_index(axis)[source]
Parameters:

axis (str)

Return type:

ndarray

spherical_harmonics(L)[source]
Parameters:

L (int)

Return type:

ndarray

class orpheus.sn.quadrature.ProductQuadrature(mu_x, mu_y, mu_z, weights, N, _ref_x, _ref_y, _ref_z, n_levels, level_indices, level_mu)[source]

Bases: object

Product quadrature: Gauss-Legendre(μ) × equispaced(φ).

The polar angle θ is discretised via Gauss-Legendre on μ = cos θ ∈ [-1, 1]. The azimuthal angle φ is discretised uniformly on [0, 2π).

Direction cosines: - μ_z = μ (axial, = cos θ) - μ_x = η = sin(θ) cos(φ) (radial for cylindrical) - μ_y = ξ = sin(θ) sin(φ) (azimuthal for cylindrical)

Weights: w = w_GL(μ) · (2π / n_phi) — sum to 4π.

Provides level_indices for the cylindrical sweep.

Parameters:
mu_x: ndarray
mu_y: ndarray
mu_z: ndarray
weights: ndarray
N: int
n_levels: int
level_indices: list[ndarray]
level_mu: ndarray
classmethod create(n_mu=8, n_phi=8)[source]

Build product quadrature with n_mu GL points and n_phi azimuthal points.

Parameters:

  • n_mu (int) – Number of Gauss-Legendre points in μ (polar).

  • n_phi (int) – Number of equispaced points in φ (azimuthal).

Return type:

ProductQuadrature

reflection_index(axis)[source]
Parameters:

axis (str)

Return type:

ndarray

spherical_harmonics(L)[source]
Parameters:

L (int)

Return type:

ndarray

Transport Sweep

Diamond-difference transport sweep for 1D and 2D Cartesian meshes.

Two paths, automatically dispatched: - 1D cumprod: for ny=1, mu_y=0 (Gauss-Legendre). Solves the

recurrence via cumulative products — O(nc) numpy ops, ~ms.

  • 2D wavefront: for general 2D. Sweeps cells along anti-diagonals (i+j=const), vectorized within each diagonal.

Both paths use the precomputed streaming stencil from SNMesh to avoid redundant per-ordinate per-cell divisions.

orpheus.sn.sweep.transport_sweep(Q, sig_t, sn_mesh, psi_bc, Q_aniso=None)[source]

Perform one full diamond-difference transport sweep.

Parameters:

  • Q ((nx, ny, ng) isotropic source density.)

  • sig_t ((nx, ny, ng) total macroscopic cross section.)

  • sn_mesh (SNMesh — augmented geometry with precomputed stencil.)

  • psi_bc (mutable dict storing persistent boundary fluxes) – for reflective BCs between outer iterations.

  • Q_aniso ((N, nx, ny, ng) per-ordinate anisotropic source (P1+) – scattering). None for isotropic-only (P0).

Returns:

  • angular_flux ((N, nx, ny, ng) angular flux per ordinate.)

  • scalar_flux ((nx, ny, ng) = Σ_n w_n ψ_n.)

Return type:

tuple[np.ndarray, np.ndarray]

Direct Transport Operator

Direct transport operator for Krylov inner solves.

Provides the explicit operator T: ψ → T·ψ via finite differences, used by the bicgstab inner solver path in SNSolver.

Two geometries are supported:

  • Cartesian 2D — T = μ_x ∂/∂x + μ_y ∂/∂y + Σ_t

  • Spherical 1D — T = μ (A ∂/∂r)/V + (α ∂/∂μ)/V + Σ_t

The sweep-based solver (source iteration) inverts T implicitly via diamond-difference sweeps. This module forms T explicitly so that scipy’s Krylov solvers (BiCGSTAB, GMRES) can solve T·ψ = b directly.

class orpheus.sn.operator.EquationMap(n_eq, n_unknowns, ordinate, ix, iy)[source]

Bases: object

Mapping between 1D solution vector and 4D angular flux.

Parameters:
n_eq: int
n_unknowns: int
ordinate: ndarray
ix: ndarray
iy: ndarray
orpheus.sn.operator.build_equation_map(nx, ny, quad, ng)[source]

Identify which (ordinate, cell) combos are unknowns.

Filter: mu_z >= 0 (upper hemisphere), and NOT incoming at reflective boundaries (those are determined by reflection).

Parameters:

  • nx (int)

  • ny (int)

  • quad (AngularQuadrature)

  • ng (int)

Return type:

EquationMap

orpheus.sn.operator.solution_to_angular_flux(solution, eq_map, quad, nx, ny, ng)[source]

Convert 1D solution vector to 4D angular flux (ng, N, nx, ny).

Applies z-reflection and reflective BCs to fill the full array from the reduced set of unknowns.

Parameters:

  • solution (ndarray)

  • eq_map (EquationMap)

  • quad (AngularQuadrature)

  • nx (int)

  • ny (int)

  • ng (int)

Return type:

ndarray

orpheus.sn.operator.angular_flux_to_scalar(fi, quad, nx, ny, ng)[source]

Integrate angular flux to scalar flux: φ = Σ w_n ψ_n.

Parameters:
Return type:

ndarray

orpheus.sn.operator.transport_operator_matvec(solution, eq_map, quad, sig_t, nx, ny, ng, dx, dy)[source]

Apply the streaming + collision operator T·ψ.

Parameters:

  • solution ((n_unknowns,) flattened angular flux vector.)

  • sig_t ((nx, ny, ng) total cross section.)

  • eq_map (EquationMap)

  • quad (AngularQuadrature)

  • nx (int)

  • ny (int)

  • ng (int)

  • dx (ndarray)

  • dy (ndarray)

Return type:

(n_unknowns,) result of T applied to the angular flux.

orpheus.sn.operator.build_transport_linear_operator(eq_map, quad, sig_t, nx, ny, ng, dx, dy)[source]

Build scipy LinearOperator for T = μ·∇ + Σ_t.

Parameters:
Return type:

LinearOperator

orpheus.sn.operator.build_rhs(fission_source, scalar_flux, eq_map, quad, sig_s, sig2, mat_map, nx, ny, ng, scattering_order=0, angular_flux=None)[source]

Build the RHS source vector for T·ψ = b.

All isotropic sources are divided by sum(weights) — the angular normalization for the discrete angular flux equation.

For Pn scattering (scattering_order > 0), the scattering source is per-ordinate using Legendre moments of the angular flux:

qS(n) = Σ_l (2l+1) · Σ_s^l^T @ [Σ_m fiL_lm · Y_lm(n)] / sum_w

Parameters:

  • fission_source ((nx, ny, ng) — already divided by sum(w) by the caller.)

  • scalar_flux ((nx, ny, ng) — current scalar flux for scattering.)

  • sig_s (dict[mat_id → list of (ng, ng)] Legendre scattering matrices.)

  • sig2 (dict[mat_id → (ng, ng)] (n,2n) matrices.)

  • scattering_order (Legendre order L (0 = P0 isotropic).)

  • angular_flux ((ng, N, nx, ny) angular flux for computing Legendre) – moments. Required if scattering_order > 0.

  • eq_map (EquationMap)

  • quad (AngularQuadrature)

  • mat_map (ndarray)

  • nx (int)

  • ny (int)

  • ng (int)

Return type:

(n_unknowns,) RHS vector.

orpheus.sn.operator.build_equation_map_spherical(nx, quad, ng)[source]

Equation map for spherical 1D: all (ordinate, cell) pairs except incoming directions at the outer reflective boundary.

Parameters:

  • nx (int)

  • quad (AngularQuadrature)

  • ng (int)

Return type:

EquationMap

orpheus.sn.operator.solution_to_angular_flux_spherical(solution, eq_map, quad, nx, ng)[source]

Convert 1D solution vector to angular flux array (ng, N, nx, 1).

Applies reflective BC at the outer boundary.

Parameters:

  • solution (ndarray)

  • eq_map (EquationMap)

  • quad (AngularQuadrature)

  • nx (int)

  • ng (int)

Return type:

ndarray

orpheus.sn.operator.transport_operator_matvec_spherical(solution, eq_map, quad, sig_t, nx, ng, face_areas, volumes, alpha_half, redist_dAw, tau_mm)[source]

Apply the spherical transport operator T·ψ.

\[(T\psi)_{n,i} = \frac{\mu_n}{V_i} \bigl[A_{i+\frac12}\psi_{i+\frac12} - A_{i-\frac12}\psi_{i-\frac12}\bigr] + \frac{\Delta A_i}{w_n V_i} \bigl[\alpha_{n+\frac12}\psi_{n+\frac12} - \alpha_{n-\frac12}\psi_{n-\frac12}\bigr] + \Sigma_t \psi_{n,i}\]

The \(\Delta A / w\) geometry factor (redist_dAw, precomputed in SNMesh) ensures per-ordinate flat-flux consistency (Bailey et al. 2009).

Face fluxes are approximated by arithmetic averages of cell-centre values.

Parameters:
Return type:

ndarray

orpheus.sn.operator.build_transport_linear_operator_spherical(eq_map, quad, sig_t, nx, ng, face_areas, volumes, alpha_half, redist_dAw, tau_mm)[source]

Build scipy LinearOperator for spherical T.

Parameters:
Return type:

LinearOperator

orpheus.sn.operator.build_rhs_spherical(fission_source, scalar_flux, eq_map, quad, sig_s, sig2, mat_map, nx, ng, scattering_order=0, angular_flux=None)[source]

Build the RHS source vector for spherical T·ψ = b.

Same structure as Cartesian build_rhs but with spherical equation map (no y-direction, no z-reflection filtering).

Parameters:
Return type:

ndarray

orpheus.sn.operator.build_equation_map_cylindrical(nx, quad, ng)

Equation map for spherical 1D: all (ordinate, cell) pairs except incoming directions at the outer reflective boundary.

Parameters:

  • nx (int)

  • quad (AngularQuadrature)

  • ng (int)

Return type:

EquationMap

orpheus.sn.operator.solution_to_angular_flux_cylindrical(solution, eq_map, quad, nx, ng)

Convert 1D solution vector to angular flux array (ng, N, nx, 1).

Applies reflective BC at the outer boundary.

Parameters:

  • solution (ndarray)

  • eq_map (EquationMap)

  • quad (AngularQuadrature)

  • nx (int)

  • ng (int)

Return type:

ndarray

orpheus.sn.operator.transport_operator_matvec_cylindrical(solution, eq_map, quad, sig_t, nx, ng, face_areas, volumes, alpha_per_level, redist_dAw_per_level, tau_mm_per_level)[source]

Apply the cylindrical transport operator T·ψ.

Per-level azimuthal redistribution with geometry-weighted \(\Delta A / w\) factor and Morel–Montry angular closure.

Parameters:
Return type:

ndarray

orpheus.sn.operator.build_transport_linear_operator_cylindrical(eq_map, quad, sig_t, nx, ng, face_areas, volumes, alpha_per_level, redist_dAw_per_level, tau_mm_per_level)[source]

Build scipy LinearOperator for cylindrical T.

Parameters:
Return type:

LinearOperator

orpheus.sn.operator.build_rhs_cylindrical(fission_source, scalar_flux, eq_map, quad, sig_s, sig2, mat_map, nx, ng, scattering_order=0, angular_flux=None)

Build the RHS source vector for spherical T·ψ = b.

Same structure as Cartesian build_rhs but with spherical equation map (no y-direction, no z-reflection filtering).

Parameters:
Return type:

ndarray