1-D Two-Group Diffusion Solver

1D two-group neutron diffusion solver for a PWR subassembly.

Solves the 1D neutron diffusion equation with vacuum boundary conditions using a finite-difference discretization and BiCGSTAB inner iterations.

Port of MATLAB CORE1D.m + funCORE1D.m.

See also

Verification Suite — diffusion verification (buckling eigenvalue, Richardson).

class orpheus.diffusion.solver.CoreGeometry(bot_refl_height=50.0, fuel_height=300.0, top_refl_height=50.0, f2f=21.61, dz=5.0)[source]

Bases: object

1D geometry for a PWR subassembly.

Parameters:
bot_refl_height: float = 50.0
fuel_height: float = 300.0
top_refl_height: float = 50.0
f2f: float = 21.61
dz: float = 5.0
property n_refl_bot: int
property n_fuel: int
property n_refl_top: int
property n_cells: int
property n_faces: int
property az: float

Cross-sectional area of one subassembly (cm^2).

property dv: float

Volume of one cell (cm^3).

class orpheus.diffusion.solver.TwoGroupXS(transport, absorption, fission, production, chi, scattering)[source]

Bases: object

Two-group macroscopic cross sections for a material.

Parameters:
transport: ndarray
absorption: ndarray
fission: ndarray
production: ndarray
chi: ndarray
scattering: ndarray
class orpheus.diffusion.solver.DiffusionResult(keff, flux, current, z_cells, z_faces, geometry)[source]

Bases: object

Results of a 1D diffusion calculation.

Parameters:
keff: float
flux: ndarray
current: ndarray
z_cells: ndarray
z_faces: ndarray
geometry: CoreGeometry
class orpheus.diffusion.solver.DiffusionSolver(geom, reflector_xs, fuel_xs, *, errtol=1e-06, maxiter_inner=2000)[source]

Bases: object

1D two-group diffusion eigenvalue solver (EigenvalueSolver protocol).

Wraps the finite-difference diffusion operator with vacuum boundary conditions and BiCGSTAB inner solves into the generic power iteration framework.

Parameters:

  • geom (CoreGeometry)

  • reflector_xs (TwoGroupXS)

  • fuel_xs (TwoGroupXS)

  • errtol (float)

  • maxiter_inner (int)

initial_flux_distribution()[source]

Return flat initial guess, shape (ng, nc).

Return type:

ndarray

compute_fission_source(flux_distribution, keff)[source]

Fission RHS: Q_f = chi * sum_g(nu_sig_f * phi * dV) / keff.

Standard eigenvalue formulation: A*phi = (1/k)*F*phi.

Parameters:
Return type:

ndarray

solve_fixed_source(fission_source, flux_distribution)[source]

Solve A*phi = fission_source via BiCGSTAB.

The solution is conditioned by dividing out max(|phi|) to prevent overflow in the absence of per-iteration power normalization.

Parameters:
Return type:

ndarray

compute_keff(flux_distribution)[source]

k_eff = production / (absorption + leakage).

Parameters:

flux_distribution (ndarray)

Return type:

float

converged(keff, keff_old, flux_distribution, flux_old, iteration)[source]

Check keff change and relative solution change.

Parameters:
Return type:

bool

orpheus.diffusion.solver.solve_diffusion_1d(geom=None, reflector_xs=None, fuel_xs=None, power_W=17675200.0, errtol=1e-06, maxiter=2000)[source]

Solve the 1D two-group diffusion eigenvalue problem.

Parameters:

  • geom (CoreGeometry (default: standard PWR SA).)

  • reflector_xs (TwoGroupXS (default: from MATLAB CORE1D.m).)

  • fuel_xs (TwoGroupXS (default: from MATLAB CORE1D.m).)

  • power_W (float — target power for flux normalization (W).)

  • errtol (float)

  • maxiter (int)

Return type:

DiffusionResult