From ae2b61b1767c7ccd49d26ebdfae7e9e765753ec3 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 17 Dec 2018 20:28:15 -0500 Subject: [PATCH 1/9] differentiate burnable mateirals which have multiple instances --- openmc/deplete/operator.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index a0b2417c8..601f60051 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -185,8 +185,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,6 +218,10 @@ class Operator(TransportOperator): Nuclides in order of how they'll appear in the simulation. """ + + # Automatically distribute burnable materials + self._differentiate_burnable_mats() + burnable_mats = set() model_nuclides = set() volume = OrderedDict() From 86600e6c4a64b6afac4f2a0f55e3cd1f63017667 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 17 Dec 2018 23:21:31 -0500 Subject: [PATCH 2/9] power density can be speicifed in depletion --- openmc/deplete/integrator/cecm.py | 19 ++++++++++++++++--- openmc/deplete/integrator/predictor.py | 20 +++++++++++++++++--- openmc/deplete/operator.py | 5 +++++ openmc/material.py | 15 +++++++++++++++ 4 files changed, 53 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 0b17b49b5..43425e1b1 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -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 RuntimeError( + "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) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index a45017813..bf6e2553e 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -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 RuntimeError( + "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) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 601f60051..ac995bc85 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -96,6 +96,8 @@ 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 @@ -226,6 +228,8 @@ class Operator(TransportOperator): 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(): @@ -236,6 +240,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.get_heavy_metal_mass() # Make sure there are burnable materials if not burnable_mats: diff --git a/openmc/material.py b/openmc/material.py index 612a7dfb6..3657436d9 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -747,6 +747,21 @@ class Material(IDManagerMixin): raise ValueError("Volume must be set in order to determine mass.") return self.volume*self.get_mass_density(nuclide) + def get_heavy_metal_mass(self): + """Return mass of heavy metal nuclides. + + Returns + ------- + float + Mass in [g] + + """ + mass = 0.0 + for nuc, _, _ in self._nuclides: + if openmc.data.zam(nuc)[0] > 90: + mass += self.get_mass(nuc) + return mass + def clone(self, memo=None): """Create a copy of this material with a new unique ID. From 493341dbe1e44ef71044191e53b51aa2dac14e85 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 18 Dec 2018 10:46:25 -0500 Subject: [PATCH 3/9] heavy metal: atomic number >= 90 --- openmc/material.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 3657436d9..c7187d53e 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -758,7 +758,7 @@ class Material(IDManagerMixin): """ mass = 0.0 for nuc, _, _ in self._nuclides: - if openmc.data.zam(nuc)[0] > 90: + if openmc.data.zam(nuc)[0] >= 90: mass += self.get_mass(nuc) return mass From 8ae7c783b587e7eb73a65741fabd83fec92b1979 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 19 Dec 2018 11:40:15 -0500 Subject: [PATCH 4/9] print a warning if the maximum neutron energy is < 20 MeV, close #1138 --- src/input_xml.F90 | 9 ++++++--- src/nuclide_header.F90 | 4 ++-- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index dab02f52b..716a89361 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2258,6 +2258,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 @@ -2272,9 +2274,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() diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 296a4a6a2..fe35f7293 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 From c05bb28b7afd4c43a2d3d4438c76f86130d16a04 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 19 Dec 2018 12:04:51 -0500 Subject: [PATCH 5/9] fix zero timing of initialization and eclapsed time --- include/openmc/timer.h | 4 ++-- src/finalize.cpp | 10 ++++++---- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/include/openmc/timer.h b/include/openmc/timer.h index 8fbcee46b..ac8c81e0c 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -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 start_; //!< starting point for clock - double elapsed_ {0.0}; //!< elasped time in [s] + double elapsed_ {0.0}; //!< elapsed time in [s] }; //============================================================================== diff --git a/src/finalize.cpp b/src/finalize.cpp index 5885470c1..84b68e82e 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -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; From 6e1114dcabbc1f4e4c8d0185c38c66b249cb83e8 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 20 Dec 2018 09:10:02 -0500 Subject: [PATCH 6/9] address review comments --- openmc/deplete/integrator/cecm.py | 2 +- openmc/deplete/integrator/predictor.py | 2 +- openmc/deplete/operator.py | 13 ++++++++++--- openmc/material.py | 1 + src/input_xml.F90 | 2 +- 5 files changed, 14 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 43425e1b1..ac519cf87 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -47,7 +47,7 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True): """ if power is None: if power_density is None: - raise RuntimeError( + raise ValueError( "Neither power nor power density was specified.") if not isinstance(power_density, Iterable): power = power_density*operator.heavy_metal diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index bf6e2553e..969144f68 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -43,7 +43,7 @@ def predictor(operator, timesteps, power=None, power_density=None, """ if power is None: if power_density is None: - raise RuntimeError( + raise ValueError( "Neither power nor power density was specified.") if not isinstance(power_density, Iterable): power = power_density*operator.heavy_metal diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index ac995bc85..4dc164c2d 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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 ---------- @@ -102,13 +104,17 @@ class Operator(TransportOperator): 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=True): 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 @@ -221,8 +227,9 @@ class Operator(TransportOperator): """ - # Automatically distribute burnable materials - self._differentiate_burnable_mats() + if self.diff_burnable_mats: + # Automatically distribute burnable materials + self._differentiate_burnable_mats() burnable_mats = set() model_nuclides = set() diff --git a/openmc/material.py b/openmc/material.py index c7187d53e..08d314c95 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -747,6 +747,7 @@ class Material(IDManagerMixin): raise ValueError("Volume must be set in order to determine mass.") return self.volume*self.get_mass_density(nuclide) + @property def get_heavy_metal_mass(self): """Return mass of heavy metal nuclides. diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 716a89361..0ce836348 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2258,7 +2258,7 @@ 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 & + 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 From 2d613aab54da6463ce0b46f944ea7ce1a1244385 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 20 Dec 2018 09:41:37 -0500 Subject: [PATCH 7/9] remove heavy_metal_mass as we already have material.fissionable_mass --- openmc/deplete/operator.py | 2 +- openmc/material.py | 30 +++++------------------------- 2 files changed, 6 insertions(+), 26 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 4dc164c2d..37fbdd5e2 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -247,7 +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.get_heavy_metal_mass() + self.heavy_metal += mat.fissionable_mass # Make sure there are burnable materials if not burnable_mats: diff --git a/openmc/material.py b/openmc/material.py index 08d314c95..e75be28a0 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -249,15 +249,11 @@ class Material(IDManagerMixin): @property def fissionable_mass(self): - if self.volume is None: - raise ValueError("Volume must be set in order to determine mass.") - density = 0.0 - for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): - Z = openmc.data.zam(nuc)[0] - if Z >= 90: - density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ - / openmc.data.AVOGADRO - return density*self.volume + mass = 0.0 + for nuc, _, _ in self._nuclides: + if openmc.data.zam(nuc)[0] >= 90: + mass += self.get_mass(nuc) + return mass @classmethod def from_hdf5(cls, group): @@ -747,22 +743,6 @@ class Material(IDManagerMixin): raise ValueError("Volume must be set in order to determine mass.") return self.volume*self.get_mass_density(nuclide) - @property - def get_heavy_metal_mass(self): - """Return mass of heavy metal nuclides. - - Returns - ------- - float - Mass in [g] - - """ - mass = 0.0 - for nuc, _, _ in self._nuclides: - if openmc.data.zam(nuc)[0] >= 90: - mass += self.get_mass(nuc) - return mass - def clone(self, memo=None): """Create a copy of this material with a new unique ID. From 0a4a1e7e81b32b41bc837eb769a5da5c57fd6ce0 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 29 Dec 2018 09:01:09 -0500 Subject: [PATCH 8/9] do not differentiate material by default --- openmc/deplete/operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 37fbdd5e2..754e22f7d 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -109,7 +109,7 @@ class Operator(TransportOperator): """ def __init__(self, geometry, settings, chain_file=None, prev_results=None, - diff_burnable_mats=True): + diff_burnable_mats=False): super().__init__(chain_file) self.round_number = False self.settings = settings From 3846168e706b11b76fe9e518db613c64d24fccac Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 2 Jan 2019 10:39:23 -0500 Subject: [PATCH 9/9] revert the function calulcating fissionable mass --- openmc/material.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index e75be28a0..d7c93ebae 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -249,11 +249,15 @@ class Material(IDManagerMixin): @property def fissionable_mass(self): - mass = 0.0 - for nuc, _, _ in self._nuclides: - if openmc.data.zam(nuc)[0] >= 90: - mass += self.get_mass(nuc) - return mass + if self.volume is None: + raise ValueError("Volume must be set in order to determine mass.") + density = 0.0 + for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): + Z = openmc.data.zam(nuc)[0] + if Z >= 90: + density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ + / openmc.data.AVOGADRO + return density*self.volume @classmethod def from_hdf5(cls, group): @@ -715,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