Monte Carlo Neutron Transport

Monte Carlo neutron transport solver for a PWR pin cell.

Simulates neutron histories in a 2D unit cell using Woodcock delta tracking, analog absorption with fission weight adjustment, Russian roulette, and splitting. Returns keff with statistical uncertainty.

Architecture (MT-20260406-008):

• Particle / Neutron — per-particle state (extensible to gamma transport).

• NeutronBank — population of neutrons with array-backed storage for performance, clean API for population control.

• Extracted functions: _precompute_xs(), _random_walk(), _russian_roulette(), _split_heavy().

• solve_monte_carlo() — orchestrator (~60 lines).

Port of MATLAB monteCarloPWR.m.

See also

Monte Carlo Neutron Transport — Key Facts, equations, gotchas.

class orpheus.mc.solver.MCGeometry(*args, **kwargs)

Bases: Protocol

Geometry interface for Monte Carlo delta-tracking.

Implementations must provide: - material_id_at(x, y) — return the material ID at position (x, y) - pitch — unit cell side length for periodic boundary conditions

The delta-tracking algorithm only needs the material at the collision point, not distance-to-surface. This makes the interface simple and extensible to future CSG implementations.

pitch: float

material_id_at(x, y)

Return the material ID at position (x, y).

Parameters:

• x (float)

• y (float)

Return type:

int

class orpheus.mc.solver.ConcentricPinCell(radii, mat_ids, pitch)

Bases: object

Concentric cylindrical pin cell in a square lattice.

Parameters:

• radii (outer radius of each annular region (innermost first).)

• mat_ids (material ID for each annulus (same length as radii).) – Regions outside the outermost radius get the last material ID.

• pitch (side length of the square unit cell (cm).)

radii: list[float]

mat_ids: list[int]

pitch: float

material_id_at(x, y)

Return material ID at (x, y) based on distance from pin center.

Parameters:

• x (float)

• y (float)

Return type:

int

classmethod default_pwr(pitch=3.6)

Default PWR geometry matching the original MATLAB layout.

Fuel: r < 0.9 cm (material 2) Clad: 0.9 ≤ r < 1.1 cm (material 1) Cool: r ≥ 1.1 cm (material 0)

Parameters:

pitch (float)

Return type:

ConcentricPinCell

class orpheus.mc.solver.SlabPinCell(boundaries, mat_ids, pitch)

Bases: object

1D slab geometry embedded in a 2D square cell (for MATLAB compatibility).

Determines material based on x-coordinate only. Matches the original MATLAB monteCarloPWR.m hard-coded regions.

Parameters:

• boundaries (list of x-coordinates separating regions.)

• mat_ids (material ID for each region (len = len(boundaries) + 1).) – mat_ids[0] is for x < boundaries[0], mat_ids[-1] for x > boundaries[-1].

• pitch (unit cell side length (cm).)

material_id_at(x, y)

Parameters:

• x (float)

• y (float)

Return type:

int

classmethod default_pwr(pitch=3.6)

Original MATLAB slab geometry: cool|clad|fuel|clad|cool.

Parameters:

pitch (float)

Return type:

SlabPinCell

class orpheus.mc.solver.MCMesh(mesh, pitch)

Bases: object

Augmented geometry for Monte Carlo delta-tracking.

Wraps a Mesh1D and provides point-wise material lookup, following the same pattern as CPMesh and SNMesh.

Supported coordinate systems:

• Cartesian — 1-D slab embedded in a 2-D square cell. Material determined by x-coordinate only.

• Cylindrical — concentric annuli (Wigner-Seitz pin cell). Material determined by radial distance from cell centre.

Parameters:

• mesh (Mesh1D) – Base geometry (from geometry.factories).

• pitch (float) – Unit cell side length (cm) for periodic boundary conditions.

material_id_at(x, y)

Return the material ID at position (x, y).

Parameters:

• x (float)

• y (float)

Return type:

int

class orpheus.mc.solver.Particle(x, y, weight, alive=True)

Bases: object

Base class for transported particles.

Designed for future extension to gamma transport or other particle types that share position and weight but differ in energy treatment.

Parameters:

• x (float)

• y (float)

• weight (float)

• alive (bool)

x: float

y: float

weight: float

alive: bool = True

class orpheus.mc.solver.Neutron(x, y, weight, alive=True, group=0)

Bases: Particle

A neutron with discrete energy group.

Parameters:

• x (float)

• y (float)

• weight (float)

• alive (bool)

• group (int)

group: int = 0

class orpheus.mc.solver.NeutronBank(x, y, weight, group, n)

Bases: object

Array-backed population of neutrons.

Stores parallel arrays for performance-critical inner loops while providing a clean population-control API (normalize, compact, grow).

Parameters:

• x (ndarray — positions.)

• y (ndarray — positions.)

• weight (ndarray — statistical weights.)

• group (ndarray — energy group indices.)

• n (int — number of live neutrons (rest is buffer).)

x

y

weight

group

n

classmethod initialize(n_neutrons, pitch, chi_cum, ng, rng, buffer_factor=4)

Create an initial neutron population.

Positions uniform in the cell, weights = 1, groups sampled from the fission spectrum.

Parameters:

• n_neutrons (int)

• pitch (float)

• chi_cum (ndarray)

• ng (int)

• rng (Generator)

• buffer_factor (int)

Return type:

NeutronBank

normalize_weights(target)

Scale live weights so they sum to target.

Parameters:

target (float)

Return type:

None

save_start_weights()

Return a copy of live weights (for roulette and keff).

Return type:

ndarray

compact()

Remove dead neutrons (weight = 0) by compacting arrays.

Return type:

None

grow(n_extra)

Extend arrays to accommodate n_extra more particles.

Parameters:

n_extra (int)

Return type:

None

class orpheus.mc.solver.MCParams(n_neutrons=100, n_inactive=100, n_active=2000, pitch=3.6, seed=None, geometry=None)

Bases: object

Monte Carlo simulation parameters.

Parameters:

• n_neutrons (int)

• n_inactive (int)

• n_active (int)

• pitch (float)

• seed (int | None)

• geometry (MCGeometry | None)

n_neutrons: int = 100

n_inactive: int = 100

n_active: int = 2000

pitch: float = 3.6

seed: int | None = None

geometry: MCGeometry | None = None

class orpheus.mc.solver.MCResult(keff, sigma, keff_history, sigma_history, flux_per_lethargy, eg_mid, elapsed_seconds)

Bases: object

Results of a Monte Carlo calculation.

Parameters:

• keff (float)

• sigma (float)

• keff_history (ndarray)

• sigma_history (ndarray)

• flux_per_lethargy (ndarray)

• eg_mid (ndarray)

• elapsed_seconds (float)

keff: float

sigma: float

keff_history: ndarray

sigma_history: ndarray

flux_per_lethargy: ndarray

eg_mid: ndarray

elapsed_seconds: float

orpheus.mc.solver.solve_monte_carlo(materials, params=None)

Run Monte Carlo neutron transport simulation.

Parameters:

• materials (dict mapping material ID (0=cool, 1=clad, 2=fuel) to Mixture.)

• params (MCParams (default: 100 neutrons, 100 inactive + 2000 active).)

Return type:

MCResult