from .burnable_poison_rod import BurnablePoisonRod
from .fuel_pin import FuelPin
from .guide_tube import GuideTube
from .critical_leakage import CriticalLeakage
from .symmetry import Symmetry
from .. import PolarQuadrature
from ._ensleeve import (
_ensleeve_quarter,
_ensleeve_half_top,
_ensleeve_half_right,
_ensleeve_full,
)
from .equivalence_theory import compute_adf_cdf_from_cmfd, compute_adf_cdf_from_moc
from .._scarabee import (
borated_water,
Material,
CrossSection,
NDLibrary,
DepletionChain,
PinCellType,
Vector,
Direction,
Cartesian2D,
MOCDriver,
CMFD,
BoundaryCondition,
SimulationMode,
Legendre8,
YamamotoTabuchi6,
P1CriticalitySpectrum,
B1CriticalitySpectrum,
FundamentalModeCriticalitySpectrum,
ADF,
CDF,
DiffusionCrossSection,
DiffusionData,
FormFactors,
LeakageCorrections,
set_logging_level,
scarabee_log,
LogLevel,
)
from typing import Optional, List, Tuple, Union
import copy
from threading import Thread
import numpy as np
from scipy.optimize import curve_fit
[docs]class PWRAssembly:
"""
A PWRAssembly instance is responsible for performing all the lattice
calculations necessary to produce few-group cross sections for a
single PWR assembly.
Parameters
----------
shape : tuple of int
The number of pin cells in the full assembly along x and y.
pitch : float
The spacing between fuel pins.
ndl : NDLibrary
Nuclear data library used for the calculation.
cells : list of list of FuelPin or GuideTube
All of the cells which describe the assembly geometry. Should be
consistent with the symmetry argument.
moderator : dict
Parameters for defining the moderator. The boron concentration,
temperature, and pressure will be used to call
:py:func:`scarabee.borated_water`. Alternatively, the `material` key
can be used to provide a user specified moderator material. Default is
an empty dictionary (all default values). Acceptable keys are:
:boron-ppm: Moderator boron concentration in parts per million. Default
is 800 if not provided. (float)
:temperature: Moderator temperature in Kelvin. Default is 570 if not
provided. (float)
:pressure: Moderator pressure in MPa. Default is 15.5 if not provided.
(float)
:legendre-order: Maximum Legendre order to load for the moderator's
anisotropic scattering data. Default is 1 if not provided. (int)
:material: User defined material which will be used directly. The
scarabee.borated_water function will not be called.
(scarabee.Material)
symmetry : Symmetry
Symmetry of the fuel assembly. Default is Symmetry.Quarter.
independent_quadrant : bool
If symmetry is Symmetry.Quarter and this attribute is true, the quarter
assembly is treated as an independent assembly, with ADFs generated for
all sides and CDFs generated for all corners. Form factors are only
generated for the pins present in the quarter assembly. This improves
the modeling of asymmetric burnable absorber loadings in the nodal
calculation. Default value is False.
linear_power : float
Linear power density of the full assembly in kW/cm. This value should
not be reduced due to symmetry. Default is 42.
assembly_pitch : optional float
Spacing between fuel assemblies. If None, assembly pitch is calculated
from the shape and pitch. Default value is None.
spacer_grid_width : optional float
Width of the spacer grid material between pin cells.
Default value is None.
spacer_grid : optional Material
Material defining the spacer grid between pin cells.
grid_sleeve_width : optional float
Width of the grid sleeve around the assembly. Default value is None.
grid_sleeve : optional Material
Material defining the grid sleeve around the assembly.
Attributes
----------
shape : tuple of int
The number of pin cells in the full assembly along x and y.
pitch : float
The spacing between fuel pins.
symmetry : Symmetry
Symmetry of the fuel assembly.
independent_quadrant : bool
If symmetry is Symmetry.Quarter and this attribute is true, the quarter
assembly is treated as an independent assembly, with ADFs generated for
all sides and CDFs generated for all corners. Form factors are only
generated for the pins present in the quarter assembly. This improves
the modeling of asymmetric burnable absorber loadings in the nodal
calculation.
assembly_pitch : float
Spacing between fuel assemblies.
fuel_volume_fraction : float
Fraction of the assembly occupied by fuel.
moderator_volume_fraction : float
Fraction of the assembly occupied by moderator.
linear_power : float
Linear power density of the full assembly in kW/cm.
initial_heavy_metal_linear_mass : float
Initial linear density of heavy metal in the full assembly at the
beginning of life in units of kg/cm.
boron_ppm : float
Moderator boron concentration in parts per million.
moderator_temp : float
Moderator temperature in Kelvin.
moderator_pressure : float
Moderator pressure in MPa.
moderator_legendre_order : int
Maximum Legendre order to load for the moderator's anisotropic
scattering data.
moderator : Material
Material representing the assembly moderator.
spacer_grid_width : optional float
Width of the spacer grid material between pin cells.
spacer_grid : optional Material
Material defining the spacer grid between pin cells.
grid_sleeve_width : optional float
Width of the grid sleeve around the assembly.
grid_sleeve : optional Material
Material defining the grid sleeve around the assembly.
dancoff_moc_track_spacing : float
Spacing between tracks in the MOC calculations for determining Dancoff
corrections. Default value is 0.05 cm.
dancoff_moc_num_angles : int
Number of azimuthal angles in the MOC calculations for determining
Dancoff corrections. Default value is 32.
dancoff_flux_tolerance : float
Flux convergence tolerance for Dancoff correction calculations. Must be
in range (0., 1.E-2). Default value is 1.E-5.
moc_track_spacing : float
Spacing between tracks in the assembly MOC calculations. Default value
is 0.03 cm.
moc_num_angles : int
Number of azimuthal angles in the assembly MOC calculations. Default
value is 32.
moc_polar_quadrature : PolarQuadrature
The polar quadrature used for the MOC calculations. Default is
YamamotoTabuchi6 for isotropic calculations and Legendre8 for
anisotropic calculations.
flux_tolerance : float
Flux convergence tolerance for assembly calculations. Must be in range
(0., 1.E-2). Default value is 1.E-5.
keff_tolerance : float
Keff convergence tolerance for assembly calculations. Must be in range
(0., 1.E-2). Default value is 1.E-5.
anisotropic : bool
If False, isotropic scattering with the transport correction is used in
the assembly calculation. Otherwise, explicit anisotropic scattering is
modeled. Default value is False. **If you wish to turn anisotropic
scattering on, you must set this attribute before calling the solve
method for the first time.**
cmfd : bool
If True, Coarse Mesh Finite Difference diffusion will be used to
accelerate the assembly calculation. If True, a CMFD energy
condensation scheme must be provided (see cmfd_condensation_scheme).
Default value is True. **If you wish to turn CMFD off/on, you must set
this attribute before calling the solve method for the first time.**
cmfd_condensation_scheme : list of list of int
Energy condensation scheme for the CMFD acceleration of the assembly
calculation.
condensation_scheme : list of list of int
Energy condensation scheme to condense from the group structure of the
library to the few-groups used in the core solver.
prefer_moc_adf_cdf : bool
If True, the MOC results will be used to generate the ADFs and CDFs,
even if CMFD is being used. ADFs and CDFs that are generated in this
manner are typically less accurate than those generated with the CMFD
results. Default value is False.
leakage_corrections : bool
If True, the diffusion data contains few-group diffusion cross sections
which DO NOT have the critical leakage model applied. Instead, a
LeakageCorrections is created to update the few-group diffusion cross
sections based on the in-situ leakage from the nodal diffusion solver.
If False, the diffusion data does not have a LeakageCorrections instance
and the few-group cross sections DO have the critical leakage model
applied. Default is False.
leakage_model : CriticalLeakage
Model used to determine the critical leakage flux spectrum, also known
as the fundamental mode. Default method is homogeneous P1.
depletion_exposure_steps : optional ndarray
1D Numpy array of assembly burn-up exposure steps, in units of MWd/kg.
Default is None.
depletion_time_steps : optional ndarray
1D Numpy array of burn-up time steps, in units of days.
Default is None.
corrector_transport : bool
If True, the corrector depletion step performs a transport calculation
to obtain an updated solution for the flux based on the predictor
nuclide concentrations, performing two transport calculations per time
step. If False, the flux from the predictor transport calculation is
used to perform the corrector step, performing only one transport
calculation per time step. Default value is True.
exposures : ndarray
1D Numpy array of the total assembly burn-up exposures at which
material information is available, in units of MWd/kg. Default value
is an empty array before solve has been called.
times : ndarray
1D Numpy array of the total assembly burn-up times at which material
information is available, in units of days. Default value is an empty
array before solve has been called.
keff : float or ndarray
If depletion was not performed, this is a single float with keff for
the infinite assembly. If depletion was performed, this is a 1D Numpy
array for the values of keff at the tabulated burn-up exposures/times.
Default value is 1 before solve has been called.
depletion_data : optional DiffusionData or list of DiffusionData
If a single assembly calculation is being performed without depletion,
this attribute will the resulting DiffusionData instance based on the
few-group condensation_scheme attribute. If a depletion calculation is
performed, this attribute will be a list of DiffusionData instances,
with one for each burn-up point. Before solve has been called, this
attribute is None.
"""
[docs] def __init__(
self,
shape: Tuple[int, int],
pitch: float,
ndl: NDLibrary,
cells: List[List[Union[FuelPin, GuideTube]]],
moderator: dict = {},
symmetry: Symmetry = Symmetry.Quarter,
independent_quadrant: bool = False,
linear_power: float = 42.0,
assembly_pitch: Optional[float] = None,
spacer_grid_width: Optional[float] = None,
spacer_grid: Optional[Material] = None,
grid_sleeve_width: Optional[float] = None,
grid_sleeve: Optional[Material] = None,
):
self._ndl: NDLibrary = ndl
self._chain: DepletionChain = self._ndl.depletion_chain
self._symmetry: Symmetry = symmetry
self._independent_quadrant: bool = independent_quadrant
if self.symmetry != Symmetry.Quarter and self.independent_quadrant:
raise RuntimeError(
"Cannot have independent_quadrant set to True when not using quarter symmetry."
)
if len(shape) != 2:
raise ValueError("Shape must have 2 entries.")
if shape[0] <= 0 or shape[1] <= 0:
raise ValueError("Shape entries must be > 0.")
if shape[0] != shape[1]:
raise ValueError("Fuel assemblies must be square.")
self._shape = shape
# Compute the number of expected cells along each dimension [y][x]
# based on the symmetry of the problem being solved.
expx = self.shape[0]
expy = self.shape[1]
if self.symmetry == Symmetry.Half:
expy = (expy // 2) + (expy % 2)
elif self.symmetry == Symmetry.Quarter:
expy = (expy // 2) + (expy % 2)
expx = (expx // 2) + (expx % 2)
self._simulated_shape = (expx, expy)
self._initial_heavy_metal_linear_mass: float = 0.0
self._cells_set = False
self._cells: List[List[Union[FuelPin, GuideTube]]] = [[]]
self._set_cells(cells)
if pitch <= 0.5:
raise ValueError("Pitch must be > 0.5.")
self._pitch = pitch
# Compute assembly pitch
self._assembly_pitch: float = self.shape[0] * self.pitch
if assembly_pitch is not None:
try:
assembly_pitch = float(assembly_pitch)
except:
raise TypeError(
"The assembly_pitch argument must be convertible to a float."
)
if assembly_pitch <= 0.0:
raise RuntimeError("Assembly pitch must be > 0.")
elif assembly_pitch < self._assembly_pitch:
raise RuntimeError(
"Provided assembly pitch is smaller than that indicated by the pitch and shape."
)
self._assembly_pitch = assembly_pitch
# Can only compute the fuel volume fraction once the assembly pitch is known
self._fuel_volume_fraction: float = 0.0
self._compute_fuel_volume_fraction()
# Check the spacer grid parameters
self._spacer_grid_width: Optional[float] = None
self._spacer_grid: Optional[Material] = None
if spacer_grid is None and spacer_grid_width is not None:
raise RuntimeError("Spacer grid width is provided but material is not.")
elif spacer_grid is not None and spacer_grid_width is None:
raise RuntimeError("Spacer grid material is provided but width is not.")
else:
self._spacer_grid_width = spacer_grid_width
self._spacer_grid = spacer_grid
if self.spacer_grid_width is not None:
if self._spacer_grid_width <= 0.0:
raise ValueError("Spacer grid width must be > 0.")
elif self._spacer_grid_width >= self._pitch:
raise ValueError("Spacer grid width must be < pitch.")
# Check the grid sleeve parameters
self._grid_sleeve_width: Optional[float] = None
self._grid_sleeve: Optional[Material] = None
if grid_sleeve is None and grid_sleeve_width is not None:
raise RuntimeError("Grid sleeve width is provided but material is not.")
elif grid_sleeve is not None and grid_sleeve_width is None:
raise RuntimeError("Grid sleeve material is provided but width is not.")
else:
self._grid_sleeve_width = grid_sleeve_width
self._grid_sleeve = grid_sleeve
if self.grid_sleeve_width is not None:
if self._grid_sleeve_width <= 0.0:
raise ValueError("Grid sleeve width must be > 0.")
elif self._grid_sleeve_width >= 0.5 * (
self._assembly_pitch - self.shape[0] * self.pitch
):
raise ValueError(
"Grid sleeve width must be < the assembly gap width."
)
# Can only be computed after spacer grid and grid sleeve processing
self._moderator_volume_fraction: float = 0.0
self._compute_moderator_volume_fraction()
# Get moderator parameters from the dictionary.
self._boron_ppm = 800.0
self._moderator_temp = 570.0
self._moderator_pressure = 15.5
self._moderator_legendre_order = 1
for key in moderator:
if key == "boron-ppm":
self._boron_ppm = float(moderator[key])
if self._boron_ppm < 0.0:
raise ValueError("Boron concentration must be >= 0.")
elif key == "temperature":
self._moderator_temp = float(moderator[key])
if self._moderator_temp <= 0.0:
raise ValueError("Moderator temperature must be > 0.")
elif key == "pressure":
self._moderator_pressure = float(moderator[key])
if self._moderator_pressure <= 0.0:
raise ValueError("Moderator pressure must be > 0.")
elif key == "legendre-order":
self._moderator_legendre_order = int(moderator[key])
if self._moderator_legendre_order < 1:
raise ValueError("Moderator Legendre order must be >= 1.")
elif key == "material":
if isinstance(moderator[key], Material) == False:
raise TypeError(
"Provided moderator is not a scarabee.Material instance."
)
self._moderator = moderator[key]
else:
raise KeyError(f'Unknown key "{key}" provided in moderator dictionary.')
# If the moderator material wasn't directly provided, we make it from the info
try:
self._moderator
except AttributeError:
# Make material for borated water
self._moderator: Material = borated_water(
self.boron_ppm, self.moderator_temp, self.moderator_pressure, self._ndl
)
self._moderator.name = f"Moderator ({self.boron_ppm} ppm boron)"
self._moderator.max_legendre_order = self._moderator_legendre_order
if linear_power <= 0.0:
raise ValueError("Linear power must be > 0.")
self._linear_power = linear_power
# Set initial boundary conditions
self._x_min_bc = BoundaryCondition.Periodic
self._x_max_bc = BoundaryCondition.Periodic
self._y_min_bc = BoundaryCondition.Periodic
self._y_max_bc = BoundaryCondition.Periodic
if self.symmetry != Symmetry.Full:
self._y_min_bc = BoundaryCondition.Reflective
self._y_max_bc = BoundaryCondition.Reflective
if self.symmetry == Symmetry.Quarter:
self._x_min_bc = BoundaryCondition.Reflective
self._x_max_bc = BoundaryCondition.Reflective
# ======================================================================
# DANCOFF CORRECTION CALCULATION DATA
# ----------------------------------------------------------------------
# Make water xs for dancoff calculation
self._moderator_dancoff_xs: CrossSection = CrossSection(
np.array([self.moderator.potential_xs]),
np.array([self.moderator.potential_xs]),
np.array([[0.0]]),
"Moderator",
)
# Spacer grid and grid sleeve dancoff cross sections
self._spacer_grid_dancoff_xs: Optional[CrossSection] = None
self._grid_sleeve_dancoff_xs: Optional[CrossSection] = None
if self.spacer_grid is not None:
self._spacer_grid_dancoff_xs: CrossSection = CrossSection(
np.array([self.spacer_grid.potential_xs]),
np.array([self.spacer_grid.potential_xs]),
np.array([[0.0]]),
"Spacer Grid",
)
if self.grid_sleeve is not None:
self._grid_sleeve_dancoff_xs: CrossSection = CrossSection(
np.array([self.grid_sleeve.potential_xs]),
np.array([self.grid_sleeve.potential_xs]),
np.array([[0.0]]),
"Grid Sleeve",
)
# Isolated cell geometry for Dancoff correction calculations
self._isolated_dancoff_cells = []
self._isolated_dancoff_mocs = []
# Full geometry for Dancoff correction calculations
self._full_dancoff_cells = []
self._full_dancoff_geom = None
self._full_dancoff_moc = None
# Dancoff correction parameters
self._dancoff_isolation_scale = 10.0
self._dancoff_moc_track_spacing = 0.05
self._dancoff_moc_num_angles = 32
self._dancoff_flux_tolerance = 1.0e-5
self._fuel_dancoff_corrections = np.zeros(
(self._simulated_shape[1], self._simulated_shape[0])
)
self._clad_dancoff_corrections = np.zeros(
(self._simulated_shape[1], self._simulated_shape[0])
)
# ======================================================================
# TRANSPORT CALCULATION DATA
# ----------------------------------------------------------------------
# Make water xs for transport calculation
self._moderator_xs: CrossSection = self.moderator.dilution_xs(
self.moderator.size * [1.0e10], self._ndl
)
# Spacer grid and grid sleeve cross sections
self._spacer_grid_xs: Optional[CrossSection] = None
self._grid_sleeve_xs: Optional[CrossSection] = None
if self.spacer_grid is not None:
self._spacer_grid_xs: CrossSection = self.spacer_grid.dilution_xs(
self.spacer_grid.size * [1.0e10], self._ndl
)
if self.grid_sleeve is not None:
self._grid_sleeve_xs: CrossSection = self.grid_sleeve.dilution_xs(
self.grid_sleeve.size * [1.0e10], self._ndl
)
self._moc_track_spacing: float = 0.03
self._moc_num_angles: int = 64
# Default polar quadrature chosen later based on _anisotropic
self._moc_polar_quadrature: Optional[PolarQuadrature] = None
self._flux_tolerance: float = 1.0e-5
self._keff_tolerance: float = 1.0e-5
self._anisotropic: bool = False
self._cmfd: bool = True
self._prefer_moc_adf_cdf: bool = False
self._asmbly_cells = []
self._asmbly_geom: Optional[Cartesian2D] = None
self._asmbly_moc: Optional[MOCDriver] = None
self._leakage_corrections: bool = False
self._leakage_model: CriticalLeakage = CriticalLeakage.P1
self._infinite_flux_spectrum = (
None # To reset to infinite spectrum in MOC driver
)
# Depletion time steps in MWd/kg and depletion time steps in days.
# Initially starts as None (should be provided by user)
self._depletion_exposure_steps: Optional[np.ndarray] = None
self._depletion_time_steps: Optional[np.ndarray] = None
# Arrays for the depletion exposures (MWd/kg) and times (days)
self._exposures: np.ndarray = np.array([])
self._times: np.ndarray = np.array([])
# Indicates if 2 (True) or 1 (False) transport calculation should be
# performed per time step.
self._corrector_transport: bool = True
# Either a single value or list of values (for each depletion step)
self._keff: Union[float, List[float]] = 1.0
# Condensation scheme to make few-group cross sections
self._condensation_scheme: Optional[List[List[int]]] = (
self._ndl.condensation_scheme
)
# Condensation scheme to used for CMFD acceleration
self._cmfd_condensation_scheme: Optional[List[List[int]]] = (
self._ndl.cmfd_condensation_scheme
)
# If a single assembly calculation is performed, this attribute will
# contain a single DiffusionData instance. If a depletion calculation
# is performed, a DiffusionData instance will be generated for each
# burn-up step.
self._diffusion_data: Optional[Union[DiffusionData, List[DiffusionData]]] = None
self._form_factors: Optional[Union[FormFactors, List[FormFactors]]] = None
@property
def shape(self) -> Tuple[int, int]:
return self._shape
@property
def pitch(self) -> float:
return self._pitch
@property
def symmetry(self) -> Symmetry:
return self._symmetry
@property
def independent_quadrant(self) -> bool:
return self._independent_quadrant
@property
def assembly_pitch(self) -> float:
return self._assembly_pitch
@property
def fuel_volume_fraction(self) -> float:
return self._fuel_volume_fraction
@property
def moderator_volume_fraction(self) -> float:
return self._moderator_volume_fraction
@property
def spacer_grid_width(self) -> Optional[float]:
return self._spacer_grid_width
@property
def spacer_grid(self) -> Optional[Material]:
return self._spacer_grid
@property
def grid_sleeve_width(self) -> Optional[float]:
return self._grid_sleeve_width
@property
def grid_sleeve(self) -> Optional[Material]:
return self._grid_sleeve
@property
def linear_power(self) -> float:
return self._linear_power
@property
def initial_heavy_metal_linear_mass(self) -> float:
return self._initial_heavy_metal_linear_mass
@property
def boron_ppm(self) -> float:
return self._boron_ppm
@property
def moderator_temp(self) -> float:
return self._moderator_temp
@property
def moderator_pressure(self) -> float:
return self._moderator_pressure
@property
def moderator_legendre_order(self) -> int:
return self._moderator_legendre_order
@property
def moderator(self) -> Material:
return self._moderator
@property
def dancoff_moc_track_spacing(self) -> float:
return self._dancoff_moc_track_spacing
@dancoff_moc_track_spacing.setter
def dancoff_moc_track_spacing(self, dts: float) -> None:
if dts <= 0.0 or dts > 0.1:
raise ValueError("Dancoff track spacing must be in range (0, 0.1].")
self._dancoff_moc_track_spacing = dts
@property
def dancoff_moc_num_angles(self) -> int:
return self._dancoff_moc_num_angles
@dancoff_moc_num_angles.setter
def dancoff_moc_num_angles(self, dna: int) -> None:
if dna % 4 != 0:
raise ValueError(
"Number of angles for Dancoff correction calculation must be a multiple of 4."
)
if dna < 4:
raise ValueError(
"Number of angles for Dancoff correction calculation must be > 4."
)
self._dancoff_moc_num_angles = int(dna)
@property
def dancoff_flux_tolerance(self) -> float:
return self._dancoff_flux_tolerance
@dancoff_flux_tolerance.setter
def dancoff_flux_tolerance(self, tol: float) -> None:
try:
tol = float(tol)
except:
raise TypeError(
"Dancoff flux tolerance must be a floating point value in (0,1.E-2)."
)
if tol <= 0.0 or tol >= 1.0e-2:
raise ValueError("Dancoff flux tolerance must be in (0,1.E-2).")
self._dancoff_flux_tolerance = tol
@property
def moc_track_spacing(self) -> float:
return self._moc_track_spacing
@moc_track_spacing.setter
def moc_track_spacing(self, dts: float) -> None:
if dts <= 0.0 or dts > 0.1:
raise ValueError("Track spacing must be in range (0, 0.1].")
self._moc_track_spacing = dts
@property
def moc_num_angles(self) -> int:
return self._moc_num_angles
@moc_num_angles.setter
def moc_num_angles(self, dna: int) -> None:
if dna % 4 != 0:
raise ValueError(
"Number of angles for MOC calculation must be a multiple of 4."
)
if dna < 4:
raise ValueError("Number of angles for MOC calculation must be > 4.")
self._moc_num_angles = int(dna)
@property
def moc_polar_quadrature(self) -> Optional[PolarQuadrature]:
return self._moc_polar_quadrature
@moc_polar_quadrature.setter
def moc_polar_quadrature(self, pq: PolarQuadrature) -> None:
self._moc_polar_quadrature = pq
@property
def flux_tolerance(self) -> float:
return self._flux_tolerance
@flux_tolerance.setter
def flux_tolerance(self, tol: float) -> None:
try:
tol = float(tol)
except:
raise TypeError(
"Flux tolerance must be a floating point value in (0,1.E-2)."
)
if tol <= 0.0 or tol >= 1.0e-2:
raise ValueError("Flux tolerance must be in (0,1.E-2).")
self._flux_tolerance = tol
@property
def keff_tolerance(self) -> float:
return self._keff_tolerance
@keff_tolerance.setter
def keff_tolerance(self, tol: float) -> None:
try:
tol = float(tol)
except:
raise TypeError(
"Keff tolerance must be a floating point value in (0,1.E-2)."
)
if tol <= 0.0 or tol >= 1.0e-2:
raise ValueError("Keff tolerance must be in (0,1.E-2).")
self._keff_tolerance = tol
@property
def anisotropic(self) -> bool:
return self._anisotropic
@anisotropic.setter
def anisotropic(self, aniso: bool) -> None:
self._anisotropic = aniso
@property
def cmfd(self) -> bool:
return self._cmfd
@cmfd.setter
def cmfd(self, cmfd: bool) -> None:
self._cmfd = cmfd
@property
def cells(self) -> List[List[Union[FuelPin, GuideTube]]]:
return self._cells
@property
def leakage_corrections(self) -> bool:
return self._leakage_corrections
@leakage_corrections.setter
def leakage_corrections(self, lc: bool) -> None:
self._leakage_corrections = lc
@property
def leakage_model(self) -> CriticalLeakage:
return self._leakage_model
@leakage_model.setter
def leakage_model(self, clm: CriticalLeakage) -> None:
self._leakage_model = clm
@property
def depletion_exposure_steps(self) -> Optional[np.ndarray]:
return self._depletion_exposure_steps
@depletion_exposure_steps.setter
def depletion_exposure_steps(self, steps: Optional[np.ndarray]) -> None:
if steps is None:
self._depletion_exposure_steps = None
self._depletion_time_steps = None
return
if not isinstance(steps, np.ndarray):
raise TypeError("Depletion exposure steps must be a 1D Numpy array.")
if steps.ndim != 1:
raise ValueError("Depletion exposure steps must be a 1D Numpy array.")
if steps.size == 0:
raise ValueError(
"Depletion exposure steps must have at least one entry. Use None for no depletion."
)
for step in steps:
if step <= 0.0:
raise ValueError("Depletion exposure steps must be > 0.")
elif step > 5.0:
raise ValueError("Depletion exposure steps should be <= 5 MWd/kg.")
self._depletion_exposure_steps = steps
self._depletion_time_steps = self._depletion_exposure_steps.copy()
self._depletion_time_steps *= (
1.0e3 * (1.0 / self.linear_power) * self.initial_heavy_metal_linear_mass
)
@property
def depletion_time_steps(self) -> Optional[np.ndarray]:
return self._depletion_time_steps
@depletion_time_steps.setter
def depletion_time_steps(self, steps: Optional[np.ndarray]) -> None:
if steps is None:
self._depletion_time_steps = None
self._depletion_exposure_steps = None
return
if not isinstance(steps, np.ndarray):
raise TypeError("Depletion time steps must be a 1D Numpy array.")
if steps.ndim != 1:
raise ValueError("Depletion time steps must be a 1D Numpy array.")
if steps.size == 0:
raise ValueError(
"Depletion time steps must have at least one entry. Use None for no depletion."
)
for step in steps:
if step <= 0.0:
raise ValueError("Depletion time steps must be > 0.")
elif step > 100:
raise ValueError("Depletion time steps should be <= 100 days.")
self._depletion_time_steps = steps
self._depletion_exposure_steps = self._depletion_time_steps.copy()
self._depletion_exposure_steps /= (
1.0e3 * (1.0 / self.linear_power) * self.initial_heavy_metal_linear_mass
)
@property
def exposures(self) -> np.ndarray:
return self._exposures
@property
def times(self) -> np.ndarray:
return self._times
@property
def corrector_transport(self) -> bool:
return self._corrector_transport
@corrector_transport.setter
def corrector_transport(self, val: bool) -> None:
self._corrector_transport = val
@property
def keff(self) -> Union[float, List[float]]:
return self._keff
@property
def condensation_scheme(self) -> Optional[List[List[int]]]:
return self._condensation_scheme
@condensation_scheme.setter
def condensation_scheme(self, cs: Optional[List[List[int]]]) -> None:
if cs is None:
self._condensation_scheme = None
return
if not isinstance(cs, list):
raise TypeError("Condensation scheme must be a list of lists of ints.")
if len(cs) == 0:
raise TypeError("Condensation scheme must have at least one group.")
for G in range(len(cs)):
if not isinstance(cs[G], list):
raise TypeError("Condensation scheme must be a list of lists of ints.")
if len(cs[G]) != 2:
raise TypeError(
"Each entry in condensation scheme must have 2 entries."
)
try:
cs[G][0] = int(cs[G][0])
cs[G][1] = int(cs[G][1])
except:
raise TypeError(
f"Microgroup indices for macrogroup {G} are not convertible to ints."
)
if cs[G][1] < cs[G][0]:
raise ValueError(
f"The microgroup indices in macrogroup {G} are not ordered."
)
if G == 0 and cs[G][0] != 0:
raise ValueError(
"The first microgroup index of the 0 macrogroup must be 0."
)
elif G == len(cs) - 1 and cs[G][1] != self._ndl.ngroups - 1:
NG = self._ndl.ngroups - 1
raise ValueError(
f"The last microgroup index of the last macrogroup must be {NG}."
)
if G > 0:
if cs[G][0] != cs[G - 1][1] + 1:
raise ValueError("The condensation scheme is not continuous")
self._condensation_scheme = copy.deepcopy(cs)
@property
def prefer_moc_adf_cdf(self) -> bool:
return self._prefer_moc_adf_cdf
@prefer_moc_adf_cdf.setter
def prefer_moc_adf_cdf(self, value: bool) -> None:
self._prefer_moc_adf_cdf = value
@property
def cmfd_condensation_scheme(self) -> Optional[List[List[int]]]:
return self._cmfd_condensation_scheme
@cmfd_condensation_scheme.setter
def cmfd_condensation_scheme(self, cs: Optional[List[List[int]]]) -> None:
if cs is None:
self._cmfd_condensation_scheme = None
return
if not isinstance(cs, list):
raise TypeError("CMFD Condensation scheme must be a list of lists of ints.")
if len(cs) == 0:
raise TypeError("CMFD Condensation scheme must have at least one group.")
for G in range(len(cs)):
if not isinstance(cs[G], list):
raise TypeError(
"CMFD Condensation scheme must be a list of lists of ints."
)
if len(cs[G]) != 2:
raise TypeError(
"Each entry in CMFD condensation scheme must have 2 entries."
)
try:
cs[G][0] = int(cs[G][0])
cs[G][1] = int(cs[G][1])
except:
raise TypeError(
f"Microgroup indices for macrogroup {G} are not convertible to ints."
)
if cs[G][1] < cs[G][0]:
raise ValueError(
f"The microgroup indices in macrogroup {G} are not ordered."
)
if G == 0 and cs[G][0] != 0:
raise ValueError(
"The first microgroup index of the 0 macrogroup must be 0."
)
elif G == len(cs) - 1 and cs[G][1] != self._ndl.ngroups - 1:
NG = self._ndl.ngroups - 1
raise ValueError(
f"The last microgroup index of the last macrogroup must be {NG}."
)
if G > 0:
if cs[G][0] != cs[G - 1][1] + 1:
raise ValueError("The CMFD condensation scheme is not continuous")
self._cmfd_condensation_scheme = copy.deepcopy(cs)
@property
def diffusion_data(self) -> Optional[Union[DiffusionData, List[DiffusionData]]]:
return self._diffusion_data
@property
def form_factors(self) -> Optional[Union[FormFactors, List[FormFactors]]]:
return self._form_factors
def _set_cells(self, cells: List[List[Union[FuelPin, GuideTube]]]) -> None:
if len(cells) != self._simulated_shape[1]:
raise ValueError(
"Shape along y of cells list does not agree with assembly shape and symmetry."
)
for j in range(len(cells)):
if len(cells[j]) != self._simulated_shape[0]:
raise ValueError(
"Shape along x of cells list does not agree with assembly shape and symmetry."
)
# Make deep copies of all cells provided by the user. This makes sure
# that each duplicate FuelPin or GuideTube instance is unique.
self._cells = []
self._cells_set = False
self._initial_heavy_metal_linear_mass = 0.0
for j in range(len(cells)):
self._cells.append([])
for i in range(len(cells[j])):
self._cells[-1].append(copy.deepcopy(cells[j][i]))
if isinstance(cells[j][i], FuelPin):
lfm = cells[j][i].initial_fissionable_linear_mass
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
lfm *= 0.25
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
lfm *= 0.5
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
lfm *= 0.5
self._initial_heavy_metal_linear_mass += lfm
self._cells_set = True
if self.symmetry == Symmetry.Half:
self._initial_heavy_metal_linear_mass *= 2.0
elif self.symmetry == Symmetry.Quarter:
self._initial_heavy_metal_linear_mass *= 4.0
# Convert HM mass from g to kg
self._initial_heavy_metal_linear_mass *= 1.0e-3
def _compute_fuel_volume_fraction(self) -> None:
"""
Computes the fraction of the assembly which is fuel.
"""
sum_volume = 0.0
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if not isinstance(cell, FuelPin):
continue
pin_volume = np.pi * cell.fuel_radius * cell.fuel_radius
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
pin_volume *= 0.25
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
pin_volume *= 0.5
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
pin_volume *= 0.5
sum_volume += pin_volume
if self.symmetry == Symmetry.Quarter:
sum_volume *= 4.0
elif self.symmetry == Symmetry.Half:
sum_volume *= 2.0
self._fuel_volume_fraction = sum_volume / (
self.assembly_pitch * self.assembly_pitch
)
def _compute_moderator_volume_fraction(self) -> None:
"""
Computes the fraction of the assembly which is moderator.
"""
sum_volume = 0.0
pin_cell_volume = self.pitch**2.0
if self.spacer_grid_width is not None:
pin_cell_volume = (self.pitch - 2.0 * self.spacer_grid_width) ** 2.0
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
# Start with volume of moderator in a pin cell
cell_mod_volume = pin_cell_volume
if isinstance(cell, FuelPin):
# Remove volume of pin
cell_mod_volume -= np.pi * cell.clad_radius * cell.clad_radius
else:
# Guide tube
# First, remove volume of the tube
cell_mod_volume -= np.pi * (
cell.outer_radius * cell.outer_radius
- cell.inner_radius * cell.inner_radius
)
# Remove volume of the burnable poison rod (if present)
if isinstance(cell.fill, BurnablePoisonRod):
cell_mod_volume -= (
np.pi
* cell.fill.outer_clad_radius
* cell.fill.outer_clad_radius
)
# Add center volume back if the center is moderator !
if cell.fill.center is not None:
cell_mod_volume += (
np.pi
* cell.fill.center_radius
* cell.fill.center_radius
)
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
cell_mod_volume *= 0.25
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
cell_mod_volume *= 0.5
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
cell_mod_volume *= 0.5
sum_volume += cell_mod_volume
if self.symmetry == Symmetry.Quarter:
sum_volume *= 4.0
elif self.symmetry == Symmetry.Half:
sum_volume *= 2.0
# Add volume for the gap around the assembly
if self.grid_sleeve_width is None:
sum_volume += self.assembly_pitch**2.0 - (self.shape[0] * self.pitch) ** 2.0
else:
sum_volume += (
self.assembly_pitch**2.0
- (self.shape[0] * self.pitch + 2.0 * self.grid_sleeve_width) ** 2.0
)
self._moderator_volume_fraction = sum_volume / (
self.assembly_pitch * self.assembly_pitch
)
# ==========================================================================
# Interrogation Methods
[docs] def get_average_fuel_nuclide_density(self, t: int, nuclide: str) -> float:
"""
Computes the average density of a nuclide within all the fuel pellets
in the assembly at a single depletion time step.
Parameters
----------
t : int
Depletion time step index.
nuclide : str
Name of the nuclide.
Returns
-------
float
Average density of the nuclide at depletion time step t across all
the fuel pellets in units of atoms per barn-cm.
"""
sum_density = 0.0
sum_volume = 0.0
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if not isinstance(cell, FuelPin):
continue
pin_avg_density = cell.get_average_fuel_nuclide_density(t, nuclide)
pin_volume = np.pi * cell.fuel_radius * cell.fuel_radius
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
pin_volume *= 0.25
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
pin_volume *= 0.5
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
pin_volume *= 0.5
sum_density += pin_volume * pin_avg_density
sum_volume += pin_volume
return sum_density / sum_volume
# ==========================================================================
# Dancoff Correction Related Methods
def _init_dancoff_components(self) -> None:
scarabee_log(
LogLevel.Info, "Initializing Dancoff correction calculation components."
)
set_logging_level(LogLevel.Warning)
self._init_isolated_dancoff_components()
self._init_full_dancoff_components()
self._save_dancoff_fsr_indexes()
set_logging_level(LogLevel.Info)
def _init_isolated_dancoff_components(self) -> None:
# Isolated pitch
iso_pitch = self._dancoff_isolation_scale * self.pitch
for j in range(len(self.cells)):
self._isolated_dancoff_cells.append([])
self._isolated_dancoff_mocs.append([])
for i in range(len(self.cells[j])):
cell = None
geom = None
moc = None
x_min_bc = BoundaryCondition.Vacuum
x_max_bc = BoundaryCondition.Vacuum
y_min_bc = BoundaryCondition.Vacuum
y_max_bc = BoundaryCondition.Vacuum
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
0.5 * iso_pitch,
0.5 * iso_pitch,
PinCellType.I,
True,
)
geom = Cartesian2D([0.5 * iso_pitch], [0.5 * iso_pitch])
geom.set_tiles([cell])
x_min_bc = BoundaryCondition.Reflective
y_min_bc = BoundaryCondition.Reflective
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
0.5 * iso_pitch,
iso_pitch,
PinCellType.XP,
True,
)
geom = Cartesian2D([0.5 * iso_pitch], [iso_pitch])
geom.set_tiles([cell])
x_min_bc = BoundaryCondition.Reflective
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
iso_pitch,
0.5 * iso_pitch,
PinCellType.YP,
True,
)
geom = Cartesian2D([iso_pitch], [0.5 * iso_pitch])
geom.set_tiles([cell])
y_min_bc = BoundaryCondition.Reflective
# Otherwise, we just make the full cell
else:
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
iso_pitch,
iso_pitch,
PinCellType.Full,
True,
)
geom = Cartesian2D([iso_pitch], [iso_pitch])
geom.set_tiles([cell])
# Save to cells
self._isolated_dancoff_cells[-1].append(cell)
# Make the MOCDriver
moc = MOCDriver(geom)
moc.sim_mode = SimulationMode.FixedSource
moc.x_min_bc = x_min_bc
moc.x_max_bc = x_max_bc
moc.y_min_bc = y_min_bc
moc.y_max_bc = y_max_bc
# Generate tracks in serial as each call will run with threads
moc.generate_tracks(
self._dancoff_moc_num_angles,
self._dancoff_moc_track_spacing,
YamamotoTabuchi6(),
)
# Save the MOC
self._isolated_dancoff_mocs[-1].append(moc)
def _init_full_dancoff_components(self) -> None:
pitch = self.pitch
if self.spacer_grid_width is not None:
pitch -= 2.0 * self.spacer_grid_width
# Get all the cells
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
# Get the geometry bit from the cell
cell = None
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
0.5 * pitch,
0.5 * pitch,
PinCellType.I,
False,
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_quarter(
cell,
pitch,
self.spacer_grid_width,
self._spacer_grid_dancoff_xs,
)
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
0.5 * pitch,
pitch,
PinCellType.XP,
False,
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_half_right(
cell,
pitch,
self.spacer_grid_width,
self._spacer_grid_dancoff_xs,
)
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
pitch,
0.5 * pitch,
PinCellType.YP,
False,
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_half_top(
cell,
pitch,
self.spacer_grid_width,
self._spacer_grid_dancoff_xs,
)
# Otherwise, we just make the full cell
else:
cell = self.cells[j][i].make_dancoff_moc_cell(
self._moderator_dancoff_xs,
pitch,
pitch,
PinCellType.Full,
False,
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_full(
cell,
pitch,
self.spacer_grid_width,
self._spacer_grid_dancoff_xs,
)
# Save to cells
self._full_dancoff_cells.append(cell)
# Construct the Cartesian2D geometry
dx = self._simulated_shape[0] * [self.pitch]
dy = self._simulated_shape[1] * [self.pitch]
if self.symmetry != Symmetry.Full and self.shape[1] % 2 == 1:
dy[0] *= 0.5
if self.symmetry == Symmetry.Quarter and self.shape[0] % 2 == 1:
dx[0] *= 0.5
self._full_dancoff_geom = Cartesian2D(dx, dy)
self._full_dancoff_geom.set_tiles(self._full_dancoff_cells)
# Add gap if needed
gap_width = 0.5 * (self.assembly_pitch - self.shape[0] * self.pitch)
if self.grid_sleeve_width is not None:
gap_width -= self.grid_sleeve_width
if gap_width > 0.0:
pins_geom = self._full_dancoff_geom
if self.symmetry == Symmetry.Quarter:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_quarter(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_dancoff_xs,
)
self._full_dancoff_geom, _ = _ensleeve_quarter(
pins_geom, self.pitch, gap_width, self._moderator_dancoff_xs
)
elif self.symmetry == Symmetry.Half:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_half_top(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_dancoff_xs,
)
self._full_dancoff_geom, _ = _ensleeve_half_top(
pins_geom, self.pitch, gap_width, self._moderator_dancoff_xs
)
else:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_full(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_dancoff_xs,
)
self._full_dancoff_geom, _ = _ensleeve_full(
pins_geom, self.pitch, gap_width, self._moderator_dancoff_xs
)
# Construct the MOC
self._full_dancoff_moc = MOCDriver(self._full_dancoff_geom)
self._full_dancoff_moc.sim_mode = SimulationMode.FixedSource
self._full_dancoff_moc.x_min_bc = self._x_min_bc
self._full_dancoff_moc.x_max_bc = self._x_max_bc
self._full_dancoff_moc.y_min_bc = self._y_min_bc
self._full_dancoff_moc.y_max_bc = self._y_max_bc
self._full_dancoff_moc.generate_tracks(
self._dancoff_moc_num_angles,
self._dancoff_moc_track_spacing,
YamamotoTabuchi6(),
)
def _save_dancoff_fsr_indexes(self) -> None:
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
cell.populate_dancoff_fsr_indexes(isomoc, self._full_dancoff_moc)
def _dancoff_components_initialized(self) -> bool:
if self._full_dancoff_moc is None:
return False
return True
[docs] def set_dancoff_moderator_xs(self) -> None:
"""
Updates the moderator cross section for all Dancoff correction calculations.
"""
self._moderator_dancoff_xs.set(
CrossSection(
np.array([self.moderator.potential_xs]),
np.array([self.moderator.potential_xs]),
np.array([[0.0]]),
"Moderator",
)
)
[docs] def set_dancoff_spacer_grid_sleeve_xs(self) -> None:
"""
Updates the spacer grid and grid sleeve cross sections for all Dancoff
correction calculations.
"""
if self.spacer_grid is not None:
self._spacer_grid_dancoff_xs.set(
CrossSection(
np.array([self.spacer_grid.potential_xs]),
np.array([self.spacer_grid.potential_xs]),
np.array([[0.0]]),
"Spacer Grid",
)
)
if self.grid_sleeve is not None:
self._grid_sleeve_dancoff_xs.set(
CrossSection(
np.array([self.grid_sleeve.potential_xs]),
np.array([self.grid_sleeve.potential_xs]),
np.array([[0.0]]),
"Grid Sleeve",
)
)
[docs] def compute_fuel_dancoff_corrections(self) -> None:
"""
Recomputes all Dancoff corrections for the fuel regions in the problem,
using the most recent material definitions. All fuel is shelf-shielded
together, regardless of wether or not is is UO2 or MOX.
"""
scarabee_log(LogLevel.Info, "Computing Dancoff corrections for the fuel.")
set_logging_level(LogLevel.Warning)
if not self._dancoff_components_initialized():
raise RuntimeError(
"Dancoff calculation components have not been initialized."
)
# Set the xs and sources for all cells
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
isomoc.flux_tolerance = self.dancoff_flux_tolerance
cell.set_xs_for_fuel_dancoff_calculation()
cell.set_isolated_dancoff_fuel_sources(isomoc, self.moderator)
cell.set_full_dancoff_fuel_sources(
self._full_dancoff_moc, self.moderator
)
self._full_dancoff_moc.flux_tolerance = self.dancoff_flux_tolerance
# Solve all the MOCs in parallel
threads = []
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if isinstance(cell, FuelPin):
isomoc = self._isolated_dancoff_mocs[j][i]
threads.append(Thread(target=isomoc.solve))
threads[-1].start()
threads.append(Thread(target=self._full_dancoff_moc.solve))
threads[-1].start()
for t in threads:
t.join()
# Go through and let each cell compute the Dancoff correction if it holds
# a fuel pin.
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
if isinstance(cell, FuelPin):
C = cell.compute_fuel_dancoff_correction(
isomoc, self._full_dancoff_moc
)
self._fuel_dancoff_corrections[j, i] = C
set_logging_level(LogLevel.Info)
[docs] def compute_clad_dancoff_corrections(self) -> None:
"""
Recomputes all Dancoff corrections for the fuel pin cladding regions
in the problem, using the most recent material definitions. All
cladding is shelf-shielded together.
"""
scarabee_log(LogLevel.Info, "Computing Dancoff corrections for the cladding.")
set_logging_level(LogLevel.Warning)
if self._full_dancoff_moc is None:
raise RuntimeError(
"Dancoff calculation components have not been initialized."
)
# Set the xs and sources for all cells
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
isomoc.flux_tolerance = self.dancoff_flux_tolerance
cell.set_xs_for_clad_dancoff_calculation(self._ndl)
cell.set_isolated_dancoff_clad_sources(
isomoc, self.moderator, self._ndl
)
cell.set_full_dancoff_clad_sources(
self._full_dancoff_moc, self.moderator, self._ndl
)
self._full_dancoff_moc.flux_tolerance = self.dancoff_flux_tolerance
# Solve all the MOCs in parallel
threads = []
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
threads.append(Thread(target=isomoc.solve))
threads[-1].start()
threads.append(Thread(target=self._full_dancoff_moc.solve))
threads[-1].start()
for t in threads:
t.join()
# Go through and let each cell compute the Dancoff correction if it holds
# a fuel pin.
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
isomoc = self._isolated_dancoff_mocs[j][i]
C = cell.compute_clad_dancoff_correction(isomoc, self._full_dancoff_moc)
self._clad_dancoff_corrections[j, i] = C
set_logging_level(LogLevel.Info)
[docs] def apply_dancoff_corrections(self) -> None:
"""
Appends all fuel and cladding Dancoff corrections to the appropriate cell.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if isinstance(cell, FuelPin):
cell.append_fuel_dancoff_correction(
self._fuel_dancoff_corrections[j, i]
)
cell.append_clad_dancoff_correction(
self._clad_dancoff_corrections[j, i]
)
[docs] def self_shield_and_xs_update(self) -> None:
"""
Computes a new set of Dancoff corrections for the fuel and the
cladding. After, these are applied to all the cells in the problem.
"""
if not self._dancoff_components_initialized():
self._init_dancoff_components()
# Update the Dancoff cross sections held by the assembly
self.set_dancoff_moderator_xs()
self.set_dancoff_spacer_grid_sleeve_xs()
# Compute Dancoff corrections
self.compute_fuel_dancoff_corrections()
self.compute_clad_dancoff_corrections()
self.apply_dancoff_corrections()
# ==========================================================================
# Transport Calculation Related Methods
def _init_moc(self) -> None:
pitch = self.pitch
if self.spacer_grid_width is not None:
pitch -= 2.0 * self.spacer_grid_width
# Get all the cells
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
# Get the geometry bit from the cell
cell = None
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
cell = self.cells[j][i].make_moc_cell(
self._moderator_xs,
0.5 * pitch,
0.5 * pitch,
PinCellType.I,
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_quarter(
cell, pitch, self.spacer_grid_width, self._spacer_grid_xs
)
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
cell = self.cells[j][i].make_moc_cell(
self._moderator_xs, 0.5 * pitch, pitch, PinCellType.XP
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_half_right(
cell, pitch, self.spacer_grid_width, self._spacer_grid_xs
)
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
cell = self.cells[j][i].make_moc_cell(
self._moderator_xs, pitch, 0.5 * pitch, PinCellType.YP
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_half_top(
cell, pitch, self.spacer_grid_width, self._spacer_grid_xs
)
# Otherwise, we just make the full cell
else:
cell = self.cells[j][i].make_moc_cell(
self._moderator_xs, pitch, pitch, PinCellType.Full
)
if self.spacer_grid_width is not None:
cell, _ = _ensleeve_full(
cell, pitch, self.spacer_grid_width, self._spacer_grid_xs
)
# Save to cells
self._asmbly_cells.append(cell)
# Construct the Cartesian2D geometry
dx = self._simulated_shape[0] * [self.pitch]
dy = self._simulated_shape[1] * [self.pitch]
if self.symmetry != Symmetry.Full and self.shape[1] % 2 == 1:
dy[0] *= 0.5
if self.symmetry == Symmetry.Quarter and self.shape[0] % 2 == 1:
dx[0] *= 0.5
self._asmbly_geom = Cartesian2D(dx, dy)
self._asmbly_geom.set_tiles(self._asmbly_cells)
# Add gap if needed
gap_width = 0.5 * (self.assembly_pitch - self.shape[0] * self.pitch)
if self.grid_sleeve_width is not None:
gap_width -= self.grid_sleeve_width
if gap_width > 0.0:
pins_geom = self._asmbly_geom
if self.symmetry == Symmetry.Quarter:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_quarter(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_xs,
)
self._asmbly_geom, _ = _ensleeve_quarter(
pins_geom, self.pitch, gap_width, self._moderator_xs
)
elif self.symmetry == Symmetry.Half:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_half_top(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_xs,
)
self._asmbly_geom, _ = _ensleeve_half_top(
pins_geom, self.pitch, gap_width, self._moderator_xs
)
else:
if self.grid_sleeve_width is not None:
pins_geom, _ = _ensleeve_full(
pins_geom,
self.pitch,
self._grid_sleeve_width,
self._grid_sleeve_xs,
)
self._asmbly_geom, _ = _ensleeve_full(
pins_geom, self.pitch, gap_width, self._moderator_xs
)
# Construct the MOC
self._asmbly_moc = MOCDriver(self._asmbly_geom, anisotropic=self.anisotropic)
self._asmbly_moc.x_min_bc = self._x_min_bc
self._asmbly_moc.x_max_bc = self._x_max_bc
self._asmbly_moc.y_min_bc = self._y_min_bc
self._asmbly_moc.y_max_bc = self._y_max_bc
# Apply CMFD if turned on
if self.cmfd and self.cmfd_condensation_scheme is None:
raise RuntimeError(
"CMFD is enabled but no CMFD condensation scheme has been provided."
)
elif self.cmfd:
extra_width = 0.5 * (self.assembly_pitch - self.shape[0] * self.pitch)
dx_cmfd = self._simulated_shape[0] * [self.pitch]
dy_cmfd = self._simulated_shape[1] * [self.pitch]
if self.symmetry != Symmetry.Full and self.shape[1] % 2 == 1:
dy_cmfd[0] *= 0.5
if self.symmetry == Symmetry.Quarter and self.shape[0] % 2 == 1:
dx_cmfd[0] *= 0.5
if self.symmetry == Symmetry.Full:
dx_cmfd[0] += extra_width
dx_cmfd[-1] += extra_width
dy_cmfd[0] += extra_width
dy_cmfd[-1] += extra_width
elif self.symmetry == Symmetry.Half:
dx_cmfd[0] += extra_width
dx_cmfd[-1] += extra_width
dy_cmfd[-1] += extra_width
else:
dx_cmfd[-1] += extra_width
dy_cmfd[-1] += extra_width
self._asmbly_moc.cmfd = CMFD(
dx_cmfd, dy_cmfd, self.cmfd_condensation_scheme
)
if self.anisotropic and self.moc_polar_quadrature is None:
self.moc_polar_quadrature = Legendre8()
elif self.moc_polar_quadrature is None:
self.moc_polar_quadrature = YamamotoTabuchi6()
# Trace tracks
self._asmbly_moc.generate_tracks(
self.moc_num_angles,
self.moc_track_spacing,
self.moc_polar_quadrature,
)
self._save_fsr_indexes()
def _save_fsr_indexes(self) -> None:
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
cell.populate_fsr_indexes(self._asmbly_moc)
[docs] def plot(self) -> None:
"""
Launches the graphical geometry plotter for the assembly calculation.
"""
if self._asmbly_moc is None:
raise RuntimeError(
"Cannot launch plotter. MOCDriver has not been initialized."
)
self._asmbly_moc.plot()
[docs] def set_moderator_xs(self) -> None:
"""
Updates the moderator cross section for transport calculations.
"""
self._moderator_xs.set(
self.moderator.dilution_xs(self.moderator.size * [1.0e10], self._ndl)
)
[docs] def set_spacer_grid_sleeve_xs(self) -> None:
"""
Updates the spacer grid and grid sleeve cross sections for transport
calculations.
"""
if self.spacer_grid is not None:
self._spacer_grid_xs.set(
self.spacer_grid.dilution_xs(
self.spacer_grid.size * [1.0e10], self._ndl
)
)
if self.grid_sleeve is not None:
self._grid_sleeve_xs.set(
self.grid_sleeve.dilution_xs(
self.grid_sleeve.size * [1.0e10], self._ndl
)
)
[docs] def recompute_all_xs(self) -> None:
"""
Computes and applies all cross sections using the most recent material
information and Dancoff corrections.
"""
self.recompute_all_fuel_xs()
self.recompute_all_gap_xs()
self.recompute_all_clad_xs()
self.recompute_all_guide_tube_fill_xs()
self.set_moderator_xs()
self.set_spacer_grid_sleeve_xs()
[docs] def recompute_all_self_shielded_xs(self) -> None:
"""
Computes and applies all fuel and caldding cross sections using the
most recent material information and Dancoff corrections.
"""
self.recompute_all_fuel_xs()
self.recompute_all_clad_xs()
[docs] def recompute_all_fuel_xs(self) -> None:
"""
Computes and applies all fuel cross sections using the most recent
material information and Dancoff corrections.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if isinstance(cell, FuelPin):
cell.set_fuel_xs_for_depletion_step(-1, self._ndl)
[docs] def recompute_all_clad_xs(self) -> None:
"""
Computes and applies all cladding cross sections using the most recent
material information and Dancoff corrections.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
cell.set_clad_xs_for_depletion_step(-1, self._ndl)
[docs] def recompute_all_gap_xs(self) -> None:
"""
Computes and applies all gap cross sections using the most recent
material information.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if isinstance(cell, FuelPin):
cell.set_gap_xs(self._ndl)
[docs] def recompute_all_guide_tube_fill_xs(self) -> None:
"""
Computes and applies all cross sections for the fill objects of guide
tubes. These could be for burnable poison rods or control rods.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if isinstance(cell, GuideTube):
cell.set_fill_xs_for_depletion_step(-1, self._ndl)
[docs] def apply_leakage_model(self, scilent: bool = False) -> None:
"""
Applied the critical leakage model to the assembly, modifying the flux
in the MOC simulation directly.
Parameters
----------
scilent : bool
If True, nothing is written to the log. If False, the type of
spectrum calculation, as well as kinf and the buckling are logged.
Default value is False.
"""
# If no leakage, just return
if self.leakage_model == CriticalLeakage.NoLeakage:
return
if not scilent:
scarabee_log(LogLevel.Info, "")
self._infinite_flux_spectrum = self._asmbly_moc.homogenize_flux_spectrum()
homogenized_moc = self._asmbly_moc.homogenize()
if self.leakage_model == CriticalLeakage.P1:
if not scilent:
scarabee_log(
LogLevel.Info, "Performing P1 criticality spectrum calculation"
)
critical_spectrum = P1CriticalitySpectrum(homogenized_moc)
elif self.leakage_model == CriticalLeakage.B1:
if not scilent:
scarabee_log(
LogLevel.Info, "Performing B1 criticality spectrum calculation"
)
critical_spectrum = B1CriticalitySpectrum(homogenized_moc)
else:
if not scilent:
scarabee_log(
LogLevel.Info,
"Performing fundamental mode criticality spectrum calculation",
)
critical_spectrum = FundamentalModeCriticalitySpectrum(homogenized_moc)
self._asmbly_moc.apply_criticality_spectrum(critical_spectrum.flux)
if not scilent:
scarabee_log(
LogLevel.Info, "Kinf : {:.5f}".format(critical_spectrum.k_inf)
)
scarabee_log(
LogLevel.Info, "Buckling: {:.5f}".format(critical_spectrum.buckling)
)
[docs] def apply_infinite_spectrum(self) -> None:
"""
Undoes the critical flux spectrum adjustment to the MOCDriver.
This permits subsequent transport calcualtions to converge much faster.
"""
if self.leakage_model == CriticalLeakage.NoLeakage:
return
if self._infinite_flux_spectrum is not None and self._asmbly_moc is not None:
self._asmbly_moc.apply_criticality_spectrum(self._infinite_flux_spectrum)
[docs] def obtain_flux_spectra(self) -> None:
"""
Computes the average flux spectrum for material regions whcih are
depleted. This includes each fuel ring in fuel pins and the poison in
burnable poison rods.
"""
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
self.cells[j][i].obtain_flux_spectra(self._asmbly_moc)
[docs] def normalize_flux_to_power(self) -> None:
"""
Normalizes the flux spectra based on the specified linear power density
for the assembly. It assumes that all power comes from fission.
"""
assembly_linear_power = 0.0
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if not isinstance(cell, FuelPin):
continue
pin_linear_power = cell.compute_pin_linear_power(self._ndl)
# First check for quarter symmetry and being corner pin
if (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
and i == 0
):
pin_linear_power *= 0.25
# Next, check for being on the side with a half pin in quarter symmetry
elif (
self.symmetry == Symmetry.Quarter
and self.shape[0] % 2 == 1
and i == 0
):
pin_linear_power *= 0.5
# Next, check for being on the bottom row with a half pin
elif (
self.symmetry != Symmetry.Full
and self.shape[1] % 2 == 1
and j == self._simulated_shape[1] - 1
):
pin_linear_power *= 0.5
assembly_linear_power += pin_linear_power
if self.symmetry == Symmetry.Half:
assembly_linear_power *= 2.0
elif self.symmetry == Symmetry.Quarter:
assembly_linear_power *= 4.0
# Compute normalization factor.
# Multiply by 10^3 to go from kW to W on linear power
f = 1.0e3 * self.linear_power / assembly_linear_power
# Normalize flux spectra
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i].normalize_flux_spectrum(f)
def _compute_form_factors(self) -> FormFactors:
"""
Computes the one group pin power form factors for the full assembly.
Returns
-------
FormFactors
A FormFactors object for the pin power reconstruction. First index
is y (from high to low) and second index is x (from low to high).
"""
if self.symmetry == Symmetry.Quarter and self.independent_quadrant:
ff = np.zeros((self._simulated_shape[1], self._simulated_shape[0]))
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if not isinstance(cell, FuelPin):
continue
# Pin Power
ff[j, i] = cell.compute_pin_linear_power(self._ndl)
mean_ff = np.mean(ff)
ff /= mean_ff
# Create the widths arrays
x_widths = np.zeros(len(self.cells[0]))
y_widths = np.zeros(len(self.cells))
x_widths[:] = self.pitch
y_widths[:] = self.pitch
if self.shape[0] % 2 == 1 and self.shape[1] % 2 == 1:
x_widths[0] *= 0.5
y_widths[0] *= 0.5
gap_width = 0.5 * (self.assembly_pitch - self.shape[0] * self.pitch)
if gap_width > 0.0:
x_widths[-1] += gap_width
y_widths[-1] += gap_width
return FormFactors(ff, x_widths, y_widths)
else:
ff = np.zeros((self.shape[1], self.shape[0]))
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
if not isinstance(cell, FuelPin):
continue
# Pin Power
pp = cell.compute_pin_linear_power(self._ndl)
if self.symmetry == Symmetry.Full:
ff[j, i] = pp
elif self.symmetry == Symmetry.Half:
ff[j, i] = pp
ff[-(j + 1), i] = pp
elif self.symmetry == Symmetry.Quarter:
ox = self.shape[1] // 2
ff[j, i + ox] = pp
ff[j, -(i + ox + 1)] = pp
ff[-(j + 1), i + ox] = pp
ff[-(j + 1), -(i + ox + 1)] = pp
mean_ff = np.mean(ff)
ff /= mean_ff
# Create the widths arrays
x_widths = np.zeros(self.shape[1])
y_widths = np.zeros(self.shape[0])
x_widths[:] = self.pitch
y_widths[:] = self.pitch
gap_width = 0.5 * (self.assembly_pitch - self.shape[0] * self.pitch)
if gap_width > 0.0:
x_widths[0] += gap_width
y_widths[0] += gap_width
x_widths[-1] += gap_width
y_widths[-1] += gap_width
return FormFactors(ff, x_widths, y_widths)
def _compute_few_group_xs(self) -> DiffusionCrossSection:
"""
Computes the few-group diffusion cross sections for the problem.
Returns
-------
DiffusionCrossSection
Few-group diffusion cross sections for the assembly.
Raises
------
RuntimeError
If the condensation_scheme attribute is not set.
"""
if self.condensation_scheme is None:
raise RuntimeError("Energy condensation scheme not set.")
# According to Smith, one should do energy condensation on the
# diffusion coefficients, and not on the transport cross sections which
# one could then use to make diffusion coefficients [1]. This is in
# contradiction to Lattice Physics Computations which states that
# either method is acceptable [2]. In light of these comments, I have
# chosen to go with Smith's recommendation of performing energy
# condensation on the diffusion coefficients.
# Get homogenized cross sections for assembly in MOC group structure
homog_xs = self._asmbly_moc.homogenize()
# Obtain the flux spectrum for condensation
if self.leakage_corrections or self.leakage_model == CriticalLeakage.NoLeakage:
flux_spectrum = self._asmbly_moc.homogenize_flux_spectrum()
elif self.leakage_model == CriticalLeakage.P1:
flux_spectrum = P1CriticalitySpectrum(homog_xs).flux
elif self.leakage_model == CriticalLeakage.B1:
flux_spectrum = B1CriticalitySpectrum(homog_xs).flux
else:
flux_spectrum = FundamentalModeCriticalitySpectrum(homog_xs).flux
# Convert xs to diffusion xs, then condense
diff_xs = homog_xs.diffusion_xs()
return diff_xs.condense(self.condensation_scheme, flux_spectrum)
def _compute_leakage_corrections(self) -> LeakageCorrections:
"""
Computes the coefficients used by the nodal diffusion solver to update
the few-group cross sections based on the in-situ leakage in the full
core simulation.
"""
if self.condensation_scheme is None:
raise RuntimeError("Energy condensation scheme not set.")
if self.leakage_model == CriticalLeakage.NoLeakage:
raise RuntimeError(
"Cannot generate leakage corrections without a leakage model."
)
# Fitting function
def fnc(xdata, C):
return C * xdata
# Number of few-groups
NG = len(self.condensation_scheme)
# initialize the leakage correction object
lc = LeakageCorrections(NG)
# Get list of bucklings, and generate all few-group cross section sets
NB = 11 # Must be odd !
B2s = np.linspace(start=-0.01, stop=0.01, num=NB)
xss = []
homog_xs = self._asmbly_moc.homogenize()
diff_xs = homog_xs.diffusion_xs()
for B2 in B2s:
if self.leakage_model == CriticalLeakage.P1:
flux_spectrum = P1CriticalitySpectrum(homog_xs, B2).flux
elif self.leakage_model == CriticalLeakage.B1:
flux_spectrum = B1CriticalitySpectrum(homog_xs, B2).flux
else:
flux_spectrum = FundamentalModeCriticalitySpectrum(homog_xs, B2).flux
xss.append(diff_xs.condense(self.condensation_scheme, flux_spectrum))
xs_ref = xss[int(NB / 2)]
D = np.zeros(NB)
Ea = np.zeros(NB)
Ef = np.zeros(NB)
vEf = np.zeros(NB)
Es = np.zeros(NB)
for G_in in range(NG):
D.fill(0.0)
Ea.fill(0.0)
Ef.fill(0.0)
vEf.fill(0.0)
# For this group, fill all the buckling dependent data
for b in range(NB):
D[b] = xss[b].D(G_in)
Ea[b] = xss[b].Ea(G_in)
Ef[b] = xss[b].Ef(G_in)
vEf[b] = xss[b].vEf(G_in)
# Compute the leakage to loss ratio
LRr = (D * B2s) / (xs_ref.Er(G_in))
# Compute fractional changes
f_D = (D - xs_ref.D(G_in)) / xs_ref.D(G_in)
f_Ea = (Ea - xs_ref.Ea(G_in)) / xs_ref.Ea(G_in)
f_Ef = (Ef - xs_ref.Ef(G_in)) / xs_ref.Ef(G_in)
f_vEf = (vEf - xs_ref.vEf(G_in)) / xs_ref.vEf(G_in)
# Fit leakage correction coefficients
C_D, _ = curve_fit(fnc, LRr, f_D)
C_Ea, _ = curve_fit(fnc, LRr, f_Ea)
C_Ef, _ = curve_fit(fnc, LRr, f_Ef)
C_vEf, _ = curve_fit(fnc, LRr, f_vEf)
# Store computed values
lc.set_D(G_in, C_D)
lc.set_Ea(G_in, C_Ea)
lc.set_Ef(G_in, C_Ef)
lc.set_vEf(G_in, C_vEf)
# Do same for each down-scattering transition
for G_out in range(G_in + 1, NG):
Es.fill(0.0)
for b in range(NB):
Es[b] = xss[b].Es(G_in, G_out)
f_Es = (Es - xs_ref.Es(G_in, G_out)) / xs_ref.Es(G_in, G_out)
C_Es, _ = curve_fit(fnc, LRr, f_Es)
lc.set_Es(G_in, G_out, C_Es)
return lc
def _compute_diffusion_data_and_form_factors(
self,
) -> Tuple[DiffusionData, FormFactors]:
"""
Computes the nodal diffusion data for the assembly along with the form
factors for pin power reconstruction.
If self.leakage_corrections is True, the coefficients used to update
the few-group cross sections in the nodal diffusion solver based on the
in-situ leakge will be computed. In this case, the few-group cross
sections DO NOT have the leakage model applied. Otherwise, the leakage
correction coefficients are not generated and the few-group cross
sections are leakage corrected.
Returns
-------
DiffusionData
Few-group diffusion cross sections and discontinuity factors.
FormFactors
Form factors to be used for pin power reconstruction.
"""
if self.condensation_scheme is None:
raise RuntimeError("Energy condensation scheme not set.")
diff_xs = self._compute_few_group_xs()
ff = self._compute_form_factors()
if (
self.cmfd
and self.cmfd_condensation_scheme is not None
and not self.prefer_moc_adf_cdf
):
# If we are using CMFD, we should use the currents to generate the
# ADFs and CDFs, as it is far more accurate, as explained by Smith [1].
adf, cdf = compute_adf_cdf_from_cmfd(
self._asmbly_moc,
self.condensation_scheme,
self.symmetry,
self.independent_quadrant,
)
else:
if not self.prefer_moc_adf_cdf:
mssg = "CMFD is not on. Accurate discontinuity factors cannot be computed without CMFD."
scarabee_log(LogLevel.Warning, mssg)
adf, cdf = compute_adf_cdf_from_moc(
self._asmbly_moc,
self.condensation_scheme,
self.symmetry,
self.independent_quadrant,
)
diff_data = DiffusionData(diff_xs, adf, cdf)
if self.leakage_corrections and self.leakage_model == CriticalLeakage.NoLeakage:
scarabee_log(
LogLevel.Warning,
"Leakage corrections were requested, but no critical leakage model is provided.",
)
scarabee_log(
LogLevel.Warning,
"Leakage correction coefficients WILL NOT be generated.",
)
elif (
self.leakage_corrections and self.leakage_model != CriticalLeakage.NoLeakage
):
diff_data.leakage_corrections = self._compute_leakage_corrections()
return diff_data, ff
def _run_assembly_calculation(
self,
self_shield: bool,
apply_dancoff_corrections: bool = False,
transport: bool = True,
) -> None:
"""
Runs a single MOC calculation, applies critical leakage model, obtains
the flux spectra for each fuel region, and then normalize the flux
based on the linear assembly power. Self-shielding can be performed
if desired. If not, the last computed Dancoff corrections can be
applied to all fuel pins (to run a time step without explicit
self-shielding step, using previously known Dancoff corrections).
Paramters
---------
self_shield : bool
If True, self-shielding is performed for the fuel and cladding.
apply_dancoff_corrections : bool, default False
If self_shield is False and this option is True, the previously
obtained Dancoff corrections are applied to all cells.
transport : bool, default True
If True, the MOC calculation is performed. Otherwise, the MOC
calculation is not performed, but the leakage correction and flux
normalization are.
"""
if self_shield:
# If we want self-shielding, do that stuff
self.self_shield_and_xs_update()
elif apply_dancoff_corrections:
# Sets dancoff corrections, even if self-shielding wasn't performed
self.apply_dancoff_corrections()
self.recompute_all_xs()
if self._asmbly_moc is None:
self._init_moc()
self._asmbly_moc.flux_tolerance = self.flux_tolerance
self._asmbly_moc.keff_tolerance = self.keff_tolerance
if transport:
set_logging_level(LogLevel.Warning)
self._asmbly_moc.solve()
set_logging_level(LogLevel.Info)
scarabee_log(LogLevel.Info, "")
scarabee_log(LogLevel.Info, "Kinf: {:.5f}".format(self._asmbly_moc.keff))
self.apply_leakage_model(scilent=self.leakage_corrections)
self.obtain_flux_spectra()
self.normalize_flux_to_power()
self.apply_infinite_spectrum()
def _predict_depletion(self, dt: float, dtm1: Optional[float]) -> None:
# Do all depletions in parallel
threads = []
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
threads.append(
Thread(
target=cell.predict_depletion,
args=(self._chain, self._ndl, dt, dtm1),
)
)
threads[-1].start()
for t in threads:
t.join()
def _correct_depletion(self, dt: float, dtm1: Optional[float]) -> None:
# Do all depletions in parallel
threads = []
for j in range(len(self.cells)):
for i in range(len(self.cells[j])):
cell = self.cells[j][i]
threads.append(
Thread(
target=cell.correct_depletion,
args=(self._chain, self._ndl, dt, dtm1),
)
)
threads[-1].start()
for t in threads:
t.join()
def _run_depletion_steps(self) -> None:
if self._chain is None:
raise RuntimeError(
"No depletion chain is present. Cannot run depletion calculation."
)
self._keff = np.zeros(self.depletion_exposure_steps.size + 1)
self._exposures = np.zeros(self.depletion_exposure_steps.size + 1)
self._times = np.zeros(self.depletion_exposure_steps.size + 1)
self._diffusion_data = []
self._form_factors = []
for t, dt in enumerate(self.depletion_time_steps):
if t > 0:
self._exposures[t] = (
self._exposures[t - 1] + self.depletion_exposure_steps[t - 1]
)
self._times[t] = self._times[t - 1] + self.depletion_time_steps[t - 1]
# Don't write separating line at first iteration
scarabee_log(LogLevel.Info, 60 * "-")
scarabee_log(LogLevel.Info, "Running Time Step {:}".format(t))
scarabee_log(
LogLevel.Info, "Exposure: {:.3E} MWd/kg".format(self._exposures[t])
)
scarabee_log(LogLevel.Info, "Time : {:.3E} days".format(self._times[t]))
scarabee_log(LogLevel.Info, "")
# Convert days to seconds
dt_sec = dt * 60.0 * 60.0 * 24.0
dtm1_sec = None
if t > 0:
dtm1_sec = self.depletion_time_steps[t - 1] * 60.0 * 60.0 * 24.0
scarabee_log(LogLevel.Info, "Predictor:")
# Run initial calcualtion for this time step
self._run_assembly_calculation(True)
scarabee_log(LogLevel.Info, "")
self._keff[t] = self._asmbly_moc.keff
dd, ff = self._compute_diffusion_data_and_form_factors()
self._diffusion_data.append(dd)
self._form_factors.append(ff)
# Predic isotopes at midpoint of step
self._predict_depletion(dt_sec, dtm1_sec)
scarabee_log(LogLevel.Info, "Corrector:")
# Run the a new transport calcualtion to get rates
self._run_assembly_calculation(False, transport=self.corrector_transport)
# Do correction step for isotopes
self._correct_depletion(dt_sec, dtm1_sec)
scarabee_log(LogLevel.Info, "")
# Run last step at the end to get keff for our final material compositions
scarabee_log(LogLevel.Info, "")
scarabee_log(LogLevel.Info, 60 * "-")
self._exposures[-1] = self._exposures[-2] + self.depletion_exposure_steps[-1]
self._times[-1] = self._times[-2] + dt
scarabee_log(
LogLevel.Info, "Running Time Step {:}".format(self._times.size - 1)
)
scarabee_log(
LogLevel.Info, "Exposure: {:.3E} MWd/kg".format(self._exposures[-1])
)
scarabee_log(LogLevel.Info, "Time : {:.3E} days".format(self._times[-1]))
scarabee_log(LogLevel.Info, "")
self._run_assembly_calculation(True)
self._keff[-1] = self._asmbly_moc.keff
dd, ff = self._compute_diffusion_data_and_form_factors()
self._diffusion_data.append(dd)
self._form_factors.append(ff)
scarabee_log(LogLevel.Info, "")
[docs] def solve(self) -> None:
"""
Solves the assembly problem. If depletion_exposure_steps or
depletion_time_steps are None, then a single k-eigenvalue problem will
be run. Otherwise, the depletion calculation will be performed.
"""
if self.depletion_exposure_steps is None:
# Single one-off calulcation
self._run_assembly_calculation(True)
self._keff = self._asmbly_moc.keff
self._diffusion_data, self._form_factors = (
self._compute_diffusion_data_and_form_factors()
)
else:
# Run depletion steps
self._run_depletion_steps()
# REFERENCES
#
# [1] K. S. Smith, “Nodal diffusion methods and lattice physics data in LWR
# analyses: Understanding numerous subtle details,” Prog Nucl Energ,
# vol. 101, pp. 360–369, 2017, doi: 10.1016/j.pnucene.2017.06.013.
#
# [2] D. Knott and A. Yamamoto, "Lattice Physics Computations" in
# Handbook of Nuclear Engineering, 2010, p 1226.
