crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1141 from liangjg/assm_depletion

Browse source 
Differentiate material internally and allow for power density in depletion
This commit is contained in:
Paul Romano 2019-01-02 13:52:46 -06:00 • committed by GitHub
commit b7cbaaa808
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 89 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

4
include/openmc/timer.h
View file

@ -37,7 +37,7 @@ public:
Timer() {};
//! Start running the timer
void start ();
void start();
//! Get total elapsed time in seconds
//! \return Elapsed time in [s]
@ -52,7 +52,7 @@ public:
private:
bool running_ {false}; //!< is timer running?
std::chrono::time_point<clock> start_; //!< starting point for clock
double elapsed_ {0.0}; //!< elasped time in [s]
double elapsed_ {0.0}; //!< elapsed time in [s]
};
//==============================================================================

19
openmc/deplete/integrator/cecm.py
View file

@ -7,7 +7,7 @@ from .cram import deplete
from ..results import Results
def cecm(operator, timesteps, power, print_out=True):
def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
r"""Deplete using the CE/CM algorithm.
Implements the second order `CE/CM predictor-corrector algorithm
@ -31,16 +31,29 @@ def cecm(operator, timesteps, power, print_out=True):
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2].
actually an area in [cm^2]. Either `power` or `power_density` must be
specified.
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by initial
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)

20
openmc/deplete/integrator/predictor.py
View file

@ -7,7 +7,8 @@ from .cram import deplete
from ..results import Results
def predictor(operator, timesteps, power, print_out=True):
def predictor(operator, timesteps, power=None, power_density=None,
print_out=True):
r"""Deplete using a first-order predictor algorithm.
Implements the first-order predictor algorithm. This algorithm is
@ -26,16 +27,29 @@ def predictor(operator, timesteps, power, print_out=True):
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]. Note that values are not cumulative.
power : float or iterable of float
power : float or iterable of float, optional
Power of the reactor in [W]. A single value indicates that the power is
constant over all timesteps. An iterable indicates potentially different
power levels for each timestep. For a 2D problem, the power can be given
in [W/cm] as long as the "volume" assigned to a depletion material is
actually an area in [cm^2].
actually an area in [cm^2]. Either `power` or `power_density` must be
specified.
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by initial
heavy metal inventory to get total power if `power` is not speficied.
print_out : bool, optional
Whether or not to print out time.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)

40
openmc/deplete/operator.py
View file

@ -70,6 +70,8 @@ class Operator(TransportOperator):
Results from a previous depletion calculation. If this argument is
specified, the depletion calculation will start from the latest state
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances
Attributes
----------
@ -96,17 +98,23 @@ class Operator(TransportOperator):
Reaction rates from the last operator step.
burnable_mats : list of str
All burnable material IDs
heavy_metal : float
Initial heavy metal inventory
local_mats : list of str
All burnable material IDs being managed by a single process
prev_res : ResultsList
Results from a previous depletion calculation
diff_burnable_mats : bool
Whether to differentiate burnable materials with multiple instances
"""
def __init__(self, geometry, settings, chain_file=None, prev_results=None):
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
diff_burnable_mats=False):
super().__init__(chain_file)
self.round_number = False
self.settings = settings
self.geometry = geometry
self.diff_burnable_mats = diff_burnable_mats
if prev_results != None:
# Reload volumes into geometry
@ -185,8 +193,28 @@ class Operator(TransportOperator):
return copy.deepcopy(op_result)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material
"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.geometry.get_all_materials().values()
if mat.depletable and mat.num_instances > 1])
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats and cell.num_instances > 1:
cell.fill = [cell.fill.clone()
for i in range(cell.num_instances)]
def _get_burnable_mats(self):
"""Determine depletable materials, volumes, and nuclids
"""Determine depletable materials, volumes, and nuclides
Returns
-------
@ -198,10 +226,17 @@ class Operator(TransportOperator):
Nuclides in order of how they'll appear in the simulation.
"""
if self.diff_burnable_mats:
# Automatically distribute burnable materials
self._differentiate_burnable_mats()
burnable_mats = set()
model_nuclides = set()
volume = OrderedDict()
self.heavy_metal = 0.0
# Iterate once through the geometry to get dictionaries
for mat in self.geometry.get_all_materials().values():
for nuclide in mat.get_nuclides():
@ -212,6 +247,7 @@ class Operator(TransportOperator):
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
volume[str(mat.id)] = mat.volume
self.heavy_metal += mat.fissionable_mass
# Make sure there are burnable materials
if not burnable_mats:

4
openmc/material.py
View file

@ -719,9 +719,9 @@ class Material(IDManagerMixin):
"""
mass_density = 0.0
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
if nuclide is None or nuclide == nuc:
density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
/ openmc.data.AVOGADRO
mass_density += density_i
return mass_density

10
src/finalize.cpp
View file

@ -23,6 +23,10 @@ int openmc_finalize()
// Clear results
openmc_reset();
// Reset timers
reset_timers();
reset_timers_f();
// Reset global variables
settings::assume_separate = false;
settings::check_overlaps = false;
@ -115,10 +119,6 @@ int openmc_reset()
simulation::k_abs_tra = 0.0;
simulation::k_sum = {0.0, 0.0};
// Reset timers
reset_timers();
reset_timers_f();
return 0;
}
@ -126,6 +126,8 @@ int openmc_hard_reset()
{
// Reset all tallies and timers
openmc_reset();
reset_timers();
reset_timers_f();
// Reset total generations and keff guess
simulation::keff = 1.0;

9
src/input_xml.F90
View file

@ -2208,6 +2208,8 @@ contains
call write_message("Maximum neutron transport energy: " // &
trim(to_str(energy_max(NEUTRON))) // " eV for " // &
trim(adjustl(nuclides(i) % name)), 7)
if (master .and. energy_max(NEUTRON) < 20.e6) call warning("Maximum &
&neutron energy is below 20 MeV. This may bias the results.")
exit
end if
end if
@ -2222,9 +2224,10 @@ contains
exit
end if
end do
if (.not. mp_found) call warning("Windowed multipole functionality is &
&turned on, but no multipole libraries were found. Make sure that &
&windowed multipole data is present in your cross_sections.xml file.")
if (master .and. .not. mp_found) call warning("Windowed multipole &
&functionality is turned on, but no multipole libraries were found. &
&Make sure that windowed multipole data is present in your &
&cross_sections.xml file.")
end if
call already_read % clear()

4
src/nuclide_header.F90
View file

@ -685,8 +685,8 @@ contains
subroutine nuclide_init_grid(this, E_min, E_max, M)
class(Nuclide), intent(inout) :: this
real(8), intent(in) :: E_min ! Minimum energy in MeV
real(8), intent(in) :: E_max ! Maximum energy in MeV
real(8), intent(in) :: E_min ! Minimum energy in eV
real(8), intent(in) :: E_max ! Maximum energy in eV
integer, intent(in) :: M ! Number of equally log-spaced bins
integer :: i, j, k ! Loop indices