diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 0b7faceb71..dd679a40a0 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -43,7 +43,7 @@ algorithms `_. SILEQIIntegrator Each of these classes expects a "transport operator" to be passed. OpenMC -provides The following classes implementing transpor operators: +provides the following transport operator classes: .. autosummary:: :toctree: generated @@ -55,7 +55,7 @@ provides The following classes implementing transpor operators: The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also have some knowledge of how nuclides transmute and decay. This is handled by the -:class:`Chain`. +:class:`Chain` class. Minimal Example --------------- diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 29c8114d93..229e90a4cf 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -21,7 +21,8 @@ material compositions over time. Each method appears as a different class. For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation using the CE/CM algorithm (deplete over a timestep using the middle-of-step reaction rates). An instance of :class:`~openmc.deplete.abc.TransportOperator` -is passed to one of these functions along with the timesteps and power level:: +is passed to one of these Integrator classes along with the timesteps and power +level:: power = 1200.0e6 # watts timesteps = [10.0, 10.0, 10.0] # days @@ -318,7 +319,7 @@ normalizing reaction rates: .. math:: - \phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \Sigma^f_i \cdot \rho_i)} + \phi = \frac{P}{V \cdot \sum_i (Q_i \cdot \sigma^f_i \cdot \n_i)} where :math:`\sum_i` is the sum over all nuclides :math:`i`. This equation makes the same assumptions and issues as discussed in diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index e3587f0527..3e2c37ece7 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -331,7 +331,7 @@ class NormalizationHelper(ABC): `fission_rates` for :meth:`update`. """ - def update(self, fission_rates, mat_index=None): + def update(self, fission_rates): """Update the normalization based on fission rates (only used for energy-based normalization) @@ -341,8 +341,6 @@ class NormalizationHelper(ABC): fission reaction rate for each isotope in the specified material. Should be ordered corresponding to initial ``rate_index`` used in :meth:`prepare` - mat_index : int - Material index """ @property diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 1e382576ee..3a7928a52d 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -426,7 +426,7 @@ class ChainFissionHelper(EnergyNormalizationHelper): self._fission_q_vector = fission_qs - def update(self, fission_rates, mat_index=None): + def update(self, fission_rates): """Update energy produced with fission rates in a material Parameters @@ -435,8 +435,6 @@ class ChainFissionHelper(EnergyNormalizationHelper): fission reaction rate for each isotope in the specified material. Should be ordered corresponding to initial ``rate_index`` used in :meth:`prepare` - mat_index : int - Unused """ self._energy += dot(fission_rates, self._fission_q_vector) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index b7956205a4..a2d89e441e 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -55,6 +55,10 @@ class IndependentOperator(OpenMCOperator): Dictionary of nuclides and their fission Q values [eV]. If not given, values will be pulled from the ``chain_file``. Only applicable if ``"normalization_mode" == "fission-q"``. + dilute_initial : float, optional + Initial atom density [atoms/cm^3] to add for nuclides that are zero + in initial condition to ensure they exist in the decay chain. + Only done for nuclides with reaction rates. reduce_chain : bool, optional If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the depletion chain up to ``reduce_chain_level``. @@ -108,8 +112,9 @@ class IndependentOperator(OpenMCOperator): micro_xs, chain_file, keff=None, - normalization_mode='source-rate', + normalization_mode='fission-q', fission_q=None, + dilute_initial=1.0e3, prev_results=None, reduce_chain=False, reduce_chain_level=None, @@ -135,6 +140,7 @@ class IndependentOperator(OpenMCOperator): chain_file, prev_results, fission_q=fission_q, + dilute_initial=dilute_initial, helper_kwargs=helper_kwargs, reduce_chain=reduce_chain, reduce_chain_level=reduce_chain_level) @@ -145,8 +151,9 @@ class IndependentOperator(OpenMCOperator): chain_file, nuc_units='atom/b-cm', keff=None, - normalization_mode='source-rate', + normalization_mode='fission-q', fission_q=None, + dilute_initial=1.0e3, prev_results=None, reduce_chain=False, reduce_chain_level=None, @@ -178,6 +185,10 @@ class IndependentOperator(OpenMCOperator): Dictionary of nuclides and their fission Q values [eV]. If not given, values will be pulled from the ``chain_file``. Only applicable if ``"normalization_mode" == "fission-q"``. + dilute_initial : float + Initial atom density [atoms/cm^3] to add for nuclides that + are zero in initial condition to ensure they exist in the decay + chain. Only done for nuclides with reaction rates. prev_results : Results, optional Results from a previous depletion calculation. reduce_chain : bool, optional @@ -202,6 +213,7 @@ class IndependentOperator(OpenMCOperator): keff=keff, normalization_mode=normalization_mode, fission_q=fission_q, + dilute_initial=dilute_initial, prev_results=prev_results, reduce_chain=reduce_chain, reduce_chain_level=reduce_chain_level, @@ -253,47 +265,6 @@ class IndependentOperator(OpenMCOperator): """Finds nuclides with cross section data""" return set(cross_sections.index) - class _IndependentNormalizationHelper(ChainFissionHelper): - """Class for calculating one-group flux based on a power. - - flux = Power / X, where X = volume * sum_i(Q_i * fission_micro_xs_i * density_i) - - Parameters - ---------- - op : openmc.deplete.IndependentOperator - Reference to the object encapsulating _IndependentNormalizationHelper. - We pass this so we don't have to duplicate :attr:`IndependentOperator.number`. - - """ - - def __init__(self, op): - rates = op.reaction_rates - self.nuc_ind_map = {ind: nuc for nuc, ind in rates.index_nuc.items()} - self._op = op - super().__init__() - - def update(self, fission_rates, mat_index=None): - """Update 'energy' produced with fission rates in a material. What - this actually calculates is the quantity X. - - Parameters - ---------- - fission_rates : numpy.ndarray - fission reaction rate for each isotope in the specified - material. Should be ordered corresponding to initial - ``rate_index`` used in :meth:`prepare` - mat_index : int - Material index - - """ - volume = self._op.number.get_mat_volume(mat_index) - densities = np.empty_like(fission_rates) - for i_nuc, nuc in self.nuc_ind_map.items(): - densities[i_nuc] = self._op.number.get_atom_density(mat_index, nuc) - fission_rates *= volume * densities - - super().update(fission_rates) - class _IndependentRateHelper(ReactionRateHelper): """Class for generating one-group reaction rates with flux and one-group cross sections. @@ -350,8 +321,10 @@ class IndependentOperator(OpenMCOperator): nuc = self.nuc_ind_map[i_nuc] rxn = self.rxn_ind_map[i_react] density = self._op.number.get_atom_density(mat_id, nuc) - self._results_cache[i_nuc, - i_react] = self._op.cross_sections[rxn][nuc] * density * volume + + # Sigma^j_i * V = sigma^j_i * rho * V + self._results_cache[i_nuc,i_react] = \ + self._op.cross_sections[rxn][nuc] * density * volume return self._results_cache @@ -370,7 +343,7 @@ class IndependentOperator(OpenMCOperator): self._rate_helper = self._IndependentRateHelper(self) if normalization_mode == "fission-q": - self._normalization_helper = self._IndependentNormalizationHelper(self) + self._normalization_helper = ChainFissionHelper() else: self._normalization_helper = SourceRateHelper() diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 982d0a55f8..ccbfab5389 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -1,7 +1,7 @@ """MicroXS module A pandas.DataFrame storing microscopic cross section data with -nuclides names as row indices and reaction names as column indices. +nuclide names as row indices and reaction names as column indices. """ import tempfile @@ -26,7 +26,7 @@ _valid_rxns.append('fission') class MicroXS(DataFrame): """Stores microscopic cross section data for use in - independent depletion. + transport-independent depletion. """ @classmethod @@ -35,7 +35,8 @@ class MicroXS(DataFrame): reaction_domain, chain_file, dilute_initial=1.0e3, - energy_bounds=(0, 20e6)): + energy_bounds=(0, 20e6), + run_kwargs=None): """Generate a one-group cross-section dataframe using OpenMC. Note that the ``openmc`` executable must be compiled. @@ -57,6 +58,8 @@ class MicroXS(DataFrame): Reaction names to tally energy_bound : 2-tuple of float, optional Bounds for the energy group. + run_kwargs : dict, optional + Keyword arguments for :meth:`openmc.model.Model.run()` Returns ------- @@ -89,7 +92,10 @@ class MicroXS(DataFrame): # create temporary run with tempfile.TemporaryDirectory() as temp_dir: - statepoint_path = model.run(cwd=temp_dir) + if run_kwargs is None: + run_kwargs = {} + run_kwargs.setdefault('cwd', temp_dir) + statepoint_path = model.run(**run_kwargs) with StatePoint(statepoint_path) as sp: for rx in xs: @@ -146,7 +152,7 @@ class MicroXS(DataFrame): for material in model.materials: if material.depletable: nuc_densities = material.get_nuclide_atom_densities() - dilute_density = 1.e-24 * dilute_initial + dilute_density = 1.0e-24 * dilute_initial material.set_density('sum') for nuc, density in nuc_densities.items(): material.remove_nuclide(nuc) diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 6e6f8dd8d4..ee569bf712 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -528,8 +528,7 @@ class OpenMCOperator(TransportOperator): # Accumulate energy from fission if fission_ind is not None: self._normalization_helper.update( - tally_rates[:, fission_ind], - mat_index=mat_index) + tally_rates[:, fission_ind]) # Divide by total number and store rates[i] = self._rate_helper.divide_by_adens(number) diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index b89bae088b..c2138b55c6 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -417,7 +417,16 @@ class Results(list): mat_id = str(mat.id) if mat_id in result.mat_to_ind: mat.volume = result.volume[mat_id] + + # Change density of all nuclides in material to atom/b-cm + atoms_per_barn_cm = mat.get_nuclide_atom_densities() + for nuc, value in atoms_per_barn_cm.items(): + mat.remove_nuclide(nuc) + mat.add_nuclide(nuc, value) mat.set_density('sum') + + # For nuclides in chain that have cross sections, replace + # density in original material with new density from results for nuc in result.nuc_to_ind: if nuc not in available_cross_sections: continue diff --git a/tests/micro_xs_simple.csv b/tests/micro_xs_simple.csv index 787ce74c39..71981c718c 100644 --- a/tests/micro_xs_simple.csv +++ b/tests/micro_xs_simple.csv @@ -1,13 +1,13 @@ nuclide,"(n,gamma)",fission -U234,22.231989822002465,0.49620744663749855 -U235,10.479008971197121,48.41787337164604 -U238,0.8673334105437324,0.10467880588762352 -U236,8.65171044607122,0.31948392400019293 -O16,7.497851000107524e-05,0.0 -O17,0.0004079227797153371,0.0 -I135,6.842395323713927,0.0 -Xe135,227463.86426990604,0.0 -Xe136,0.02317896034753588,0.0 +U234,22.231989822002454,0.4962074466374984 +U235,10.479008971197121,48.41787337164606 +U238,0.8673334105437321,0.1046788058876236 +U236,8.651710446071224,0.31948392400019293 +O16,7.497851000107522e-05,0.0 +O17,0.0004079227797153372,0.0 +I135,6.842395323713929,0.0 +Xe135,227463.8642699061,0.0 +Xe136,0.023178960347535887,0.0 Cs135,2.1721665580713623,0.0 -Gd157,12786.09939237018,0.0 +Gd157,12786.099392370175,0.0 Gd156,3.4006085445846983,0.0 diff --git a/tests/regression_tests/deplete_no_transport/test.py b/tests/regression_tests/deplete_no_transport/test.py index 95b842537c..f5ee4bf1db 100644 --- a/tests/regression_tests/deplete_no_transport/test.py +++ b/tests/regression_tests/deplete_no_transport/test.py @@ -20,6 +20,15 @@ def fuel(): return fuel +@pytest.fixture(scope="module") +def micro_xs(): + micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv' + return MicroXS.from_csv(micro_xs_file) + +@pytest.fixture(scope="module") +def chain_file(): + return Path(__file__).parents[2] / 'chain_simple.xml' + @pytest.mark.parametrize("multiproc, from_nuclides, normalization_mode, power, flux", [ (True, True,'source-rate', None, 1164719970082145.0), (False, True, 'source-rate', None, 1164719970082145.0), @@ -29,7 +38,15 @@ def fuel(): (False, False, 'source-rate', None, 1164719970082145.0), (True, False, 'fission-q', 174, None), (False, False, 'fission-q', 174, None)]) -def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclides, normalization_mode, power, flux): +def test_against_self(run_in_tmpdir, + fuel, + micro_xs, + chain_file, + multiproc, + from_nuclides, + normalization_mode, + power, + flux): """Transport free system test suite. Runs an OpenMC transport-free depletion calculation and verifies @@ -37,28 +54,22 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide """ # Create operator - micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv' - micro_xs = MicroXS.from_csv(micro_xs_file) - chain_file = Path(__file__).parents[2] / 'chain_simple.xml' - - if from_nuclides: - nuclides = {} - for nuc, dens in fuel.get_nuclide_atom_densities().items(): - nuclides[nuc] = dens - - op = IndependentOperator.from_nuclides( - fuel.volume, nuclides, micro_xs, chain_file, normalization_mode=normalization_mode) - - else: - op = IndependentOperator(openmc.Materials([fuel]), micro_xs, chain_file, normalization_mode=normalization_mode) + op = _create_operator(from_nuclides, + fuel, + micro_xs, + chain_file, + normalization_mode) # Power and timesteps - dt = [30] # single step + dt = [360] # single step # Perform simulation using the predictor algorithm openmc.deplete.pool.USE_MULTIPROCESSING = multiproc - openmc.deplete.PredictorIntegrator( - op, dt, power=power, source_rates=flux, timestep_units='d').integrate() + openmc.deplete.PredictorIntegrator(op, + dt, + power=power, + source_rates=flux, + timestep_units='s').integrate() # Get path to test and reference results path_test = op.output_dir / 'depletion_results.h5' @@ -73,9 +84,86 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide res_ref = openmc.deplete.Results(path_reference) # Assert same mats + _assert_same_mats(res_test, res_ref) + + tol = 1.0e-14 + _assert_atoms_equal(res_test, res_ref, tol) + _assert_reaction_rates_equal(res_test, res_ref, tol) + +@pytest.mark.parametrize("multiproc, dt, time_units, time_type, atom_tol, rx_tol ", [ + (True, 360, 's', 'minutes', 2.0e-3, 3.0e-2), + (False, 360, 's', 'minutes', 2.0e-3, 3.0e-2), + (True, 4, 'h', 'hours', 2.0e-3, 6.0e-2), + (False,4, 'h', 'hours', 2.0e-3, 6.0e-2), + (True, 5, 'd', 'days', 2.0e-3, 5.0e-2), + (False,5, 'd', 'days', 2.0e-3, 5.0e-2), + (True, 100, 'd', 'months', 4.0e-3, 9.0e-2), + (False, 100, 'd', 'months', 4.0e-3, 9.0e-2)]) +def test_against_coupled(run_in_tmpdir, + fuel, + micro_xs, + chain_file, + multiproc, + dt, + time_units, + time_type, + atom_tol, + rx_tol): + # Create operator + op = _create_operator(False, fuel, micro_xs, chain_file, 'fission-q') + + # Power and timesteps + dt = [dt] # single step + + # Perform simulation using the predictor algorithm + openmc.deplete.pool.USE_MULTIPROCESSING = multiproc + openmc.deplete.PredictorIntegrator( + op, dt, power=174, timestep_units=time_units).integrate() + + # Get path to test and reference results + path_test = op.output_dir / 'depletion_results.h5' + + ref_path = f'test_reference_coupled_{time_type}.h5' + path_reference = Path(__file__).with_name(ref_path) + + # Load the reference/test results + res_test = openmc.deplete.Results(path_test) + res_ref = openmc.deplete.Results(path_reference) + + # Assert same mats + _assert_same_mats(res_test, res_ref) + + _assert_atoms_equal(res_test, res_ref, atom_tol) + _assert_reaction_rates_equal(res_test, res_ref, rx_tol) + +def _create_operator(from_nuclides, + fuel, + micro_xs, + chain_file, + normalization_mode): + if from_nuclides: + nuclides = {} + for nuc, dens in fuel.get_nuclide_atom_densities().items(): + nuclides[nuc] = dens + + op = IndependentOperator.from_nuclides(fuel.volume, + nuclides, + micro_xs, + chain_file, + normalization_mode=normalization_mode) + + else: + op = IndependentOperator(openmc.Materials([fuel]), + micro_xs, + chain_file, + normalization_mode=normalization_mode) + + return op + +def _assert_same_mats(res_ref, res_test): for mat in res_ref[0].mat_to_ind: - assert mat in res_test[0].mat_to_ind, \ - "Material {} not in new results.".format(mat) + assert mat in res_test[0].mat_to_ind, \ + "Material {} not in new results.".format(mat) for nuc in res_ref[0].nuc_to_ind: assert nuc in res_test[0].nuc_to_ind, \ "Nuclide {} not in new results.".format(nuc) @@ -87,7 +175,7 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide assert nuc in res_ref[0].nuc_to_ind, \ "Nuclide {} not in old results.".format(nuc) - tol = 1.0e-6 +def _assert_atoms_equal(res_ref, res_test, tol): for mat in res_test[0].mat_to_ind: for nuc in res_test[0].nuc_to_ind: _, y_test = res_test.get_atoms(mat, nuc) @@ -104,3 +192,26 @@ def test_no_transport_from_nuclides(run_in_tmpdir, fuel, multiproc, from_nuclide assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format( mat, nuc, y_old, y_test) + +def _assert_reaction_rates_equal(res_ref, res_test, tol): + for reactions in res_test[0].rates: + for mat in reactions.index_mat: + for nuc in reactions.index_nuc: + for rx in reactions.index_rx: + y_test = res_test.get_reaction_rate(mat, nuc, rx)[1] / \ + res_test.get_atoms(mat, nuc)[1] + y_old = res_ref.get_reaction_rate(mat, nuc, rx)[1] / \ + res_ref.get_atoms(mat, nuc)[1] + + # Test each point + correct = True + for i, ref in enumerate(y_old): + if ref != y_test[i]: + if ref != 0.0: + correct = np.abs(y_test[i] - ref) / ref <= tol + else: + if y_test[i] != 0.0: + correct = False + + assert correct, "Discrepancy in mat {}, nuc {}, and rx {}\n{}\n{}".format( + mat, nuc, rx, y_old, y_test) diff --git a/tests/regression_tests/deplete_no_transport/test_reference_coupled_days.h5 b/tests/regression_tests/deplete_no_transport/test_reference_coupled_days.h5 new file mode 100644 index 0000000000..d97acffec6 Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_coupled_days.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_coupled_hours.h5 b/tests/regression_tests/deplete_no_transport/test_reference_coupled_hours.h5 new file mode 100644 index 0000000000..6d9bcc91ff Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_coupled_hours.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_coupled_minutes.h5 b/tests/regression_tests/deplete_no_transport/test_reference_coupled_minutes.h5 new file mode 100644 index 0000000000..1e22bf2b71 Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_coupled_minutes.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_coupled_months.h5 b/tests/regression_tests/deplete_no_transport/test_reference_coupled_months.h5 new file mode 100644 index 0000000000..43b45cb370 Binary files /dev/null and b/tests/regression_tests/deplete_no_transport/test_reference_coupled_months.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 b/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 index 82bf115673..826c5f3b01 100644 Binary files a/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 and b/tests/regression_tests/deplete_no_transport/test_reference_fission_q.h5 differ diff --git a/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 b/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 index 88ed73b176..f5161d3706 100644 Binary files a/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 and b/tests/regression_tests/deplete_no_transport/test_reference_source_rate.h5 differ