Method of Characteristics Solvers

2D Method of Characteristics (MoC) neutron transport solver for a pin cell.

Solves the multi-group transport equation on a 2-D pin-cell geometry (concentric annuli inside a square lattice cell) with reflective boundary conditions. Uses the flat-source MOC with a product angular quadrature (azimuthal x Tabuchi-Yamamoto polar) and exact ray tracing through circular boundaries.

Reference: Boyd et al. (2014) “The OpenMOC Method of Characteristics”, Ann. Nucl. Energy 68, 43-52.

See also

Method of Characteristics (MOC) — Key Facts, equations, gotchas.

class orpheus.moc.solver.MoCResult(keff, keff_history, flux_per_material, scalar_flux, moc_mesh, eg, elapsed_seconds)[source]

Bases: object

Results of a Method of Characteristics calculation.

Parameters:
keff: float
keff_history: list[float]
flux_per_material: dict[int, ndarray]
scalar_flux: ndarray
moc_mesh: MOCMesh
eg: ndarray
elapsed_seconds: float
property flux_fuel: ndarray
property flux_clad: ndarray
property flux_cool: ndarray
orpheus.moc.solver.solve_moc(materials, mesh=None, n_azi=16, n_polar=3, ray_spacing=0.05, max_outer=500, keff_tol=1e-06, flux_tol=1e-05, n_inner_sweeps=15)[source]

Run the 2D Method of Characteristics transport calculation.

Parameters:

  • materials (dict[int, Mixture]) – Macroscopic cross sections keyed by material ID.

  • mesh (Mesh1D, optional) – Cylindrical 1-D Wigner-Seitz mesh. Defaults to the standard PWR pin cell via pwr_pin_equivalent().

  • n_azi (int) – Number of azimuthal angles in [0, pi).

  • n_polar (int) – Number of TY polar angles per half-space (1, 2, or 3).

  • ray_spacing (float) – Perpendicular spacing between parallel rays (cm).

  • max_outer (int) – Maximum number of outer (power) iterations.

  • keff_tol (float) – Convergence tolerances.

  • flux_tol (float) – Convergence tolerances.

  • n_inner_sweeps (int) – Number of transport sweeps per outer iteration.

Return type:

MoCResult

Angular quadrature for the Method of Characteristics.

Product quadrature: uniform azimuthal angles x Tabuchi-Yamamoto (TY) polar angles. The TY quadrature is optimised for the Bickley-function integrals that arise when the 2-D MOC flat-source solution is integrated over polar angle (Yamamoto et al., J. Nucl. Sci. Technol. 44(2), 2007).

class orpheus.moc.quadrature.MOCQuadrature(n_azi, n_polar, phi, omega_azi, sin_polar, omega_polar)[source]

Bases: object

Product quadrature for 2-D Method of Characteristics.

Azimuthal: n_azi uniform angles in [0, pi). Each azimuthal angle defines a family of parallel tracks; the supplementary angle (phi + pi) is the same physical track traversed in the opposite direction.

Polar: Tabuchi-Yamamoto quadrature with n_polar angles per half-space. The polar angle enters the MOC only through the effective 3-D optical thickness tau / sin(theta_p).

Parameters:
n_azi

Number of azimuthal angles in [0, pi).

Type:

int

n_polar

Number of polar angles per half-space (1, 2, or 3).

Type:

int

phi

Azimuthal angles in [0, pi) (radians).

Type:

ndarray, shape (n_azi,)

omega_azi

Azimuthal weights, each = 1 / n_azi, summing to 1.

Type:

ndarray, shape (n_azi,)

sin_polar

sin(theta_p) for each TY polar angle.

Type:

ndarray, shape (n_polar,)

omega_polar

TY polar weights (sum = 0.5 for one hemisphere).

Type:

ndarray, shape (n_polar,)

n_azi: int
n_polar: int
phi: ndarray
omega_azi: ndarray
sin_polar: ndarray
omega_polar: ndarray
classmethod create(n_azi=16, n_polar=3)[source]

Build an MOC product quadrature.

Parameters:

  • n_azi (int) – Number of azimuthal angles in [0, pi). Must be >= 2.

  • n_polar (int) – Number of TY polar angles per half-space (1, 2, or 3).

Return type:

MOCQuadrature

Augmented geometry for the Method of Characteristics.

MOCMesh wraps a cylindrical Mesh1D and a MOCQuadrature, precomputing all ray-tracing data (tracks, segments, reflective boundary links) for a pin-cell with concentric annuli inside a square lattice cell.

Inverse Wigner-Seitz: the Mesh1D is built with geometry.factories.pwr_pin_equivalent(), whose outer edge is the Wigner-Seitz radius r_cell = pitch / sqrt(pi). MOCMesh recovers the pitch and reinterprets the outermost annular region as the square border bounded by the cell walls.

class orpheus.moc.geometry.Segment(region_id, length)[source]

Bases: object

One segment of a ray through a flat-source region.

Parameters:
region_id: int
length: float
class orpheus.moc.geometry.Track(segments, azi_index, entry_point, exit_point, entry_surface, exit_surface, fwd_link, fwd_link_fwd, bwd_link, bwd_link_fwd)[source]

Bases: object

A single characteristic ray across the unit cell.

Each track supports two traversal directions:

  • Forward — direction (cos phi, sin phi) with sin >= 0. Enters at entry_point (entry_surface), exits at exit_point (exit_surface).

  • Backward — direction (-cos phi, -sin phi). Enters at exit_point, exits at entry_point.

Reflective BC links describe where outgoing angular flux goes:

  • fwd_link / fwd_link_fwd: when the forward sweep exits, its outgoing flux feeds into track fwd_link. If fwd_link_fwd is True the flux enters that track’s forward entry; otherwise its backward entry.

  • bwd_link / bwd_link_fwd: same for the backward sweep.

Parameters:
segments: tuple[Segment, ...]
azi_index: int
entry_point: tuple[float, float]
exit_point: tuple[float, float]
entry_surface: int
exit_surface: int
fwd_link: int
fwd_link_fwd: bool
bwd_link: int
bwd_link_fwd: bool
class orpheus.moc.geometry.MOCMesh(mesh, quadrature, ray_spacing=0.05)[source]

Bases: object

Augmented geometry for the Method of Characteristics.

Wraps a cylindrical Mesh1D (Wigner-Seitz pin cell) and an MOCQuadrature, precomputing all ray-tracing data for a 2-D pin-cell transport calculation.

The outermost Mesh1D region is reinterpreted as the square border via the inverse Wigner-Seitz transformation: pitch = edges[-1] * sqrt(pi).

Parameters:

  • mesh (Mesh1D)

  • quadrature (MOCQuadrature)

  • ray_spacing (float)

pitch: float
n_regions: int
radii: ndarray
region_mat_ids: ndarray
region_areas: ndarray
tracks: list[Track]
tracks_per_azi: list[list[int]]
effective_spacing(azi_index)[source]

Effective perpendicular ray spacing for azimuthal angle index.

Parameters:

azi_index (int)

Return type:

float

MOC transport solver satisfying the EigenvalueSolver protocol.

Solves the multi-group neutron transport equation via the Method of Characteristics on a 2-D pin-cell geometry with reflective boundary conditions. The flat-source approximation is used within each FSR, and the angular discretisation uses a product quadrature (azimuthal x Tabuchi-Yamamoto polar).

Reference: Boyd et al. (2014) “The OpenMOC Method of Characteristics Neutral Particle Transport Code”, Ann. Nucl. Energy 68, 43-52.

class orpheus.moc.core.MOCSolver(moc_mesh, materials, keff_tol=1e-06, flux_tol=1e-05, n_inner_sweeps=15)[source]

Bases: object

2-D Method of Characteristics eigenvalue solver.

Satisfies the EigenvalueSolver protocol.

Flux representation: (n_regions, ng) — scalar flux per flat-source region and energy group.

Parameters:
initial_flux_distribution()[source]
Return type:

ndarray

compute_fission_source(flux_distribution, keff)[source]

Fission source: chi * (SigP . phi) / keff.

Parameters:
Return type:

ndarray

solve_fixed_source(fission_source, flux_distribution)[source]

MOC transport sweep with inner scattering iterations.

Each inner sweep: 1. Build isotropic source Q from fission + scattering + (n,2n) 2. Sweep all tracks (forward + backward), accumulate delta_phi 3. Update scalar flux from delta_phi (Boyd Eq. 45)

Parameters:
Return type:

ndarray

compute_keff(flux_distribution)[source]

k_eff = production / absorption (zero leakage, reflective BCs).

Parameters:

flux_distribution (ndarray)

Return type:

float

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

bool