From 3e2ddceaac8484503c57fef4248af50399067b96 Mon Sep 17 00:00:00 2001 From: josh Date: Tue, 11 Oct 2022 03:03:50 +0000 Subject: [PATCH 1/7] Add tolerance for interp temps outside of bounds --- docs/source/io_formats/settings.rst | 11 +++++++-- openmc/settings.py | 16 +++++++------ src/cell.cpp | 4 ++-- src/nuclide.cpp | 35 +++++++++++++++++++++++++++++ src/thermal.cpp | 35 +++++++++++++++++++++-------- 5 files changed, 81 insertions(+), 20 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1a8f63716f..7a794bf003 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -832,7 +832,9 @@ cell, the nearest temperature at which cross sections are given is to be applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of "interpolation" indicates that cross sections are to be linear-linear interpolated between temperatures at which nuclear data are present (see -:ref:`temperature_treatment`). +:ref:`temperature_treatment`). With the "interpolation" method, temperatures +outside of the bounds of the nuclear data may be accepted, provided they still +fall within the tolerance (see :ref:`temperature_tolerance`). *Default*: "nearest" @@ -871,7 +873,12 @@ The ```` element specifies a tolerance in Kelvin that is to be applied when the "nearest" temperature method is used. For example, if a cell temperature is 340 K and the tolerance is 15 K, then the closest temperature in the range of 325 K to 355 K will be used to evaluate cross -sections. +sections. If the ```` is "interpolation", the tolerance +specified applies to cell temperatures outside of the data bounds. For example, +If a cell is specified at 695K, a tolerance of 15K and data only available at +700K and 1000K, the cell's cross sections will be evaluated at 700K, since +desired temperature of 695K is within the tolerance of the actual data despite +not being bounded on both sides. *Default*: 10 K diff --git a/openmc/settings.py b/openmc/settings.py index ad5a9b28ae..9d9b421ca3 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -186,13 +186,15 @@ class Settings: 'default', 'method', 'range', 'tolerance', and 'multipole'. The value for 'default' should be a float representing the default temperature in Kelvin. The value for 'method' should be 'nearest' or 'interpolation'. - If the method is 'nearest', 'tolerance' indicates a range of temperature - within which cross sections may be used. The value for 'range' should be - a pair a minimum and maximum temperatures which are used to indicate - that cross sections be loaded at all temperatures within the - range. 'multipole' is a boolean indicating whether or not the windowed - multipole method should be used to evaluate resolved resonance cross - sections. + If the method is 'nearest', 'tolerance' indicates a range of + temperature within which cross sections may be used. If the method is + 'interpolation', 'tolerance' indicates the range of temperatures outside + of the available cross section temperatures where cross sections will + evaluate to the nearer bound. The value for 'range' should be a pair a + minimum and maximum temperatures which are used to indicate that cross + sections be loaded at all temperatures within the range. 'multipole' is + a boolean indicating whether or not the windowed multipole method should + be used to evaluate resolved resonance cross sections. trace : tuple or list Show detailed information about a single particle, indicated by three integers: the batch number, generation number, and particle number diff --git a/src/cell.cpp b/src/cell.cpp index 4ff35498c6..d69f79a3e7 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -334,10 +334,10 @@ double Cell::temperature(int32_t instance) const void Cell::set_temperature(double T, int32_t instance, bool set_contained) { if (settings::temperature_method == TemperatureMethod::INTERPOLATION) { - if (T < data::temperature_min) { + if (T < data::temperature_min - settings::temperature_tolerance) { throw std::runtime_error {"Temperature is below minimum temperature at " "which data is available."}; - } else if (T > data::temperature_max) { + } else if (T > data::temperature_max + settings::temperature_tolerance) { throw std::runtime_error {"Temperature is above maximum temperature at " "which data is available."}; } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 564e27d446..4c1706509b 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -169,6 +169,22 @@ Nuclide::Nuclide(hid_t group, const vector& temperature) } if (!found_pair) { + // If no pairs found, check if the desired temperature falls just + // outside of data + if (T_desired - temps_available.front() <= + -settings::temperature_tolerance) { + if (!contains(temps_to_read, temps_available.front())) { + temps_to_read.push_back(std::round(temps_available.front())); + } + break; + } + if (T_desired - temps_available.back() <= + settings::temperature_tolerance) { + if (!contains(temps_to_read, temps_available.back())) { + temps_to_read.push_back(std::round(temps_available.back())); + } + break; + } fatal_error( "Nuclear data library does not contain cross sections for " + name_ + " at temperatures that bound " + std::to_string(T_desired) + " K."); @@ -646,6 +662,16 @@ void Nuclide::calculate_xs( } break; case TemperatureMethod::INTERPOLATION: + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i_temp = 0; + break; + } + if (kT > kTs_.back()) { + i_temp = kTs_.size() - 1; + break; + } + // Find temperatures that bound the actual temperature for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) { if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) @@ -969,6 +995,15 @@ std::pair Nuclide::find_temperature(double T) const } break; case TemperatureMethod::INTERPOLATION: + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i_temp = 0; + break; + } + if (kT > kTs_.back()) { + i_temp = kTs_.size() - 1; + break; + } // Find temperatures that bound the actual temperature while (kTs_[i_temp + 1] < kT && i_temp + 1 < n - 1) ++i_temp; diff --git a/src/thermal.cpp b/src/thermal.cpp index 1a9592d0af..66e1fb0095 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -117,6 +117,16 @@ ThermalScattering::ThermalScattering( } } if (!found) { + // If no pairs found, check if the desired temperature falls within + // bounds' tolerance + if (T - temps_available[0] <= -settings::temperature_tolerance) { + temps_to_read.push_back(std::round(temps_available[0])); + break; + } + if (T - temps_available[n - 1] <= settings::temperature_tolerance) { + temps_to_read.push_back(std::round(temps_available[n - 1])); + break; + } fatal_error( fmt::format("Nuclear data library does not contain cross " "sections for {} at temperatures that bound {} K.", @@ -159,20 +169,27 @@ void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, auto n = kTs_.size(); if (n > 1) { - // Find temperatures that bound the actual temperature - while (kTs_[i + 1] < kT && i + 1 < n - 1) - ++i; - if (settings::temperature_method == TemperatureMethod::NEAREST) { + while (kTs_[i + 1] < kT && i + 1 < n - 1) + ++i; // Pick closer of two bounding temperatures if (kT - kTs_[i] > kTs_[i + 1] - kT) ++i; - } else { - // Randomly sample between temperature i and i+1 - double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]); - if (f > prn(seed)) - ++i; + // If current kT outside of the bounds of available, snap to the bound + if (kT < kTs_.front()) { + i = 0; + } else if (kT > kTs_.back()) { + i = kTs_.size() - 1; + } else { + // Find temperatures that bound the actual temperature + while (kTs_[i + 1] < kT && i + 1 < n - 1) + ++i; + // Randomly sample between temperature i and i+1 + double f = (kT - kTs_[i]) / (kTs_[i + 1] - kTs_[i]); + if (f > prn(seed)) + ++i; + } } } From b0a300aedcd453b46f530547afbb6280c4b564f0 Mon Sep 17 00:00:00 2001 From: josh Date: Thu, 13 Oct 2022 21:24:25 +0000 Subject: [PATCH 2/7] add to existing temp interp test, update doccs --- docs/source/io_formats/settings.rst | 8 +- tests/unit_tests/test_temp_interp.py | 121 ++++++++++++++++++++++++--- 2 files changed, 113 insertions(+), 16 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 7a794bf003..1a29e00ee7 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -875,10 +875,10 @@ cell temperature is 340 K and the tolerance is 15 K, then the closest temperature in the range of 325 K to 355 K will be used to evaluate cross sections. If the ```` is "interpolation", the tolerance specified applies to cell temperatures outside of the data bounds. For example, -If a cell is specified at 695K, a tolerance of 15K and data only available at -700K and 1000K, the cell's cross sections will be evaluated at 700K, since -desired temperature of 695K is within the tolerance of the actual data despite -not being bounded on both sides. +if a cell is specified at 695K, a tolerance of 15K and data is only available +at 700K and 1000K, the cell's cross sections will be evaluated at 700K, since +the desired temperature of 695K is within the tolerance of the actual data +despite not being bounded on both sides. *Default*: 10 K diff --git a/tests/unit_tests/test_temp_interp.py b/tests/unit_tests/test_temp_interp.py index f87d529619..c64ebe46ad 100644 --- a/tests/unit_tests/test_temp_interp.py +++ b/tests/unit_tests/test_temp_interp.py @@ -86,6 +86,70 @@ def make_fake_cross_section(): lib.export_to_xml('cross_sections_fake.xml') +def fake_thermal_scattering(): + """Create a fake thermal scattering library for U-235 at 294K and 600K + """ + fake_tsl = openmc.data.ThermalScattering("c_U_fake", 1.9968, 4.9, [0.0253]) + fake_tsl.nuclides = ['U235'] + + # Create elastic reaction + bragg_edges = [0.00370672, 0.00494229] + factors = [0.00375735, 0.01386287] + coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors) + incoherent_xs_294 = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672]) + elastic_xs_base = openmc.data.Sum((coherent_xs, incoherent_xs_294)) + elastic_xs = {'294K': elastic_xs_base, '600K': elastic_xs_base} + coherent_dist = openmc.data.CoherentElasticAE(coherent_xs) + incoherent_dist_294 = openmc.data.IncoherentElasticAEDiscrete([ + [-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6] + ]) + incoherent_dist_600 = openmc.data.IncoherentElasticAEDiscrete([ + [-0.1, -0.2, 0.2, 0.1], [-0.1, -0.2, 0.2, 0.1] + ]) + elastic_dist = { + '294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_294), + '600K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_600) + } + fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) + + # Create inelastic reaction + inelastic_xs = { + '294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35]), + '600K': openmc.data.Tabulated1D([1.0e-2, 10], [1.4, 5]) + } + breakpoints = [3] + interpolation = [2] + energy = [1.0e-5, 4.3e-2, 4.9] + energy_out = [ + openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]), + ] + for eout in energy_out: + eout.normalize() + eout.c = eout.cdf() + discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8) + discrete.c = discrete.cdf()[1:] + mu = [[discrete]*4]*3 + dist = openmc.data.IncoherentInelasticAE( + breakpoints, interpolation, energy, energy_out, mu) + inelastic_dist = {'294K': dist, '600K': dist} + inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist) + fake_tsl.inelastic = inelastic + + return fake_tsl + + +def edit_fake_cross_sections(): + """Edit the test cross sections xml to include fake thermal scattering data + """ + lib = openmc.data.DataLibrary.from_xml("cross_sections_fake.xml") + c_U_fake = fake_thermal_scattering() + c_U_fake.export_to_hdf5("c_U_fake.h5") + lib.register_file("c_U_fake.h5") + lib.export_to_xml("cross_sections_fake.xml") + + @pytest.fixture(scope='module') def model(tmp_path_factory): tmp_path = tmp_path_factory.mktemp("temp_interp") @@ -119,21 +183,23 @@ def model(tmp_path_factory): @pytest.mark.parametrize( - ["method", "temperature", "fission_expected"], + ["method", "temperature", "fission_expected", "tolerance"], [ - ("nearest", 300.0, 0.5), - ("nearest", 600.0, 1.0), - ("nearest", 900.0, 0.5), - ("interpolation", 360.0, 0.6), - ("interpolation", 450.0, 0.75), - ("interpolation", 540.0, 0.9), - ("interpolation", 660.0, 0.9), - ("interpolation", 750.0, 0.75), - ("interpolation", 840.0, 0.6), + ("nearest", 300.0, 0.5, 10), + ("nearest", 600.0, 1.0, 10), + ("nearest", 900.0, 0.5, 10), + ("interpolation", 360.0, 0.6, 10), + ("interpolation", 450.0, 0.75, 10), + ("interpolation", 540.0, 0.9, 10), + ("interpolation", 660.0, 0.9, 10), + ("interpolation", 750.0, 0.75, 10), + ("interpolation", 840.0, 0.6, 10), + ("interpolation", 295.0, 0.5, 10), + ("interpolation", 990.0, 0.5, 100), ] ) -def test_interpolation(model, method, temperature, fission_expected): - model.settings.temperature = {'method': method, 'default': temperature} +def test_interpolation(model, method, temperature, fission_expected, tolerance): + model.settings.temperature = {'method': method, 'default': temperature, "tolerance": tolerance} sp_filename = model.run() with openmc.StatePoint(sp_filename) as sp: t = sp.tallies[model.tallies[0].id] @@ -152,3 +218,34 @@ def test_interpolation(model, method, temperature, fission_expected): assert k.n == pytest.approx(nu*fission_expected) else: assert abs(k.n - nu*fission_expected) <= 3*k.s + + +def test_temperature_interpolation_tolerance(model): + """Test applying global and cell temperatures with thermal scattering libraries + """ + edit_fake_cross_sections() + model.materials[0].add_s_alpha_beta("c_U_fake") + + # Default k-effective, using the thermal scattering data's minimum available temperature + model.settings.temperature = {'method': "nearest", 'default': 294, "tolerance": 50} + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + default_k = sp.keff.n + + # Get k-effective with temperature below the minimum but in interpolation mode + model.settings.temperature = {'method': "interpolation", 'default': 255, "tolerance": 50} + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + interpolated_k = sp.keff.n + + # Get the k-effective with the temperature applied to the cell, instead of globally + model.settings.temperature = {'method': "interpolation", 'default': 500, "tolerance": 50} + for cell in model.geometry.get_all_cells().values(): + cell.temperature = 275 + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + cell_k = sp.keff.n + + # All calculated k-effectives should be equal + assert default_k == pytest.approx(interpolated_k) + assert interpolated_k == pytest.approx(cell_k) From 7cb6883647b16227eb4bcdd26ab082734cca6459 Mon Sep 17 00:00:00 2001 From: josh Date: Sun, 16 Oct 2022 03:57:11 +0000 Subject: [PATCH 3/7] output bounding temp for cell temp out of bounds --- src/cell.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index d69f79a3e7..30663e04d1 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -334,12 +334,12 @@ double Cell::temperature(int32_t instance) const void Cell::set_temperature(double T, int32_t instance, bool set_contained) { if (settings::temperature_method == TemperatureMethod::INTERPOLATION) { - if (T < data::temperature_min - settings::temperature_tolerance) { - throw std::runtime_error {"Temperature is below minimum temperature at " - "which data is available."}; - } else if (T > data::temperature_max + settings::temperature_tolerance) { - throw std::runtime_error {"Temperature is above maximum temperature at " - "which data is available."}; + if (T < (data::temperature_min - settings::temperature_tolerance)) { + throw std::runtime_error {fmt::format("Temperature of {} K is below minimum temperature at " + "which data is available of {} K.", T, data::temperature_min)}; + } else if (T > (data::temperature_max + settings::temperature_tolerance)) { + throw std::runtime_error {fmt::format("Temperature of {} K is above maximum temperature at " + "which data is available of {} K.", T, data::temperature_max)}; } } From fc0e2c9bd06a42f436a81cca5d79e4460b8455f1 Mon Sep 17 00:00:00 2001 From: josh Date: Sun, 30 Oct 2022 03:28:14 +0000 Subject: [PATCH 4/7] undo formatting change of settings file --- openmc/settings.py | 715 +++++++++++++++++++++------------------------ 1 file changed, 340 insertions(+), 375 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index e6561ffb93..f13b233de9 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -17,14 +17,14 @@ from openmc.checkvalue import PathLike class RunMode(Enum): - EIGENVALUE = "eigenvalue" - FIXED_SOURCE = "fixed source" - PLOT = "plot" - VOLUME = "volume" - PARTICLE_RESTART = "particle restart" + EIGENVALUE = 'eigenvalue' + FIXED_SOURCE = 'fixed source' + PLOT = 'plot' + VOLUME = 'volume' + PARTICLE_RESTART = 'particle restart' -_RES_SCAT_METHODS = ["dbrc", "rvs"] +_RES_SCAT_METHODS = ['dbrc', 'rvs'] class Settings: @@ -190,7 +190,7 @@ class Settings: temperature within which cross sections may be used. If the method is 'interpolation', 'tolerance' indicates the range of temperatures outside of the available cross section temperatures where cross sections will - evaluate to the nearer bound. The value for 'range' should be a pair a + evaluate to the nearer bound. The value for 'range' should be a pair of minimum and maximum temperatures which are used to indicate that cross sections be loaded at all temperatures within the range. 'multipole' is a boolean indicating whether or not the windowed multipole method should @@ -245,7 +245,7 @@ class Settings: self._max_order = None # Source subelement - self._source = cv.CheckedList(Source, "source distributions") + self._source = cv.CheckedList(Source, 'source distributions') self._confidence_intervals = None self._electron_treatment = None @@ -290,8 +290,7 @@ class Settings: self._resonance_scattering = {} self._volume_calculations = cv.CheckedList( - VolumeCalculation, "volume calculations" - ) + VolumeCalculation, 'volume calculations') self._create_fission_neutrons = None self._delayed_photon_scaling = None @@ -301,7 +300,7 @@ class Settings: self._event_based = None self._max_particles_in_flight = None self._write_initial_source = None - self._weight_windows = cv.CheckedList(WeightWindows, "weight windows") + self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows') self._weight_windows_on = None self._max_splits = None self._max_tracks = None @@ -499,109 +498,103 @@ class Settings: @run_mode.setter def run_mode(self, run_mode: str): - cv.check_value("run mode", run_mode, {x.value for x in RunMode}) + cv.check_value('run mode', run_mode, {x.value for x in RunMode}) for mode in RunMode: if mode.value == run_mode: self._run_mode = mode @batches.setter def batches(self, batches: int): - cv.check_type("batches", batches, Integral) - cv.check_greater_than("batches", batches, 0) + cv.check_type('batches', batches, Integral) + cv.check_greater_than('batches', batches, 0) self._batches = batches @generations_per_batch.setter def generations_per_batch(self, generations_per_batch: int): - cv.check_type("generations per patch", generations_per_batch, Integral) - cv.check_greater_than("generations per batch", generations_per_batch, 0) + cv.check_type('generations per patch', generations_per_batch, Integral) + cv.check_greater_than('generations per batch', generations_per_batch, 0) self._generations_per_batch = generations_per_batch @inactive.setter def inactive(self, inactive: int): - cv.check_type("inactive batches", inactive, Integral) - cv.check_greater_than("inactive batches", inactive, 0, True) + cv.check_type('inactive batches', inactive, Integral) + cv.check_greater_than('inactive batches', inactive, 0, True) self._inactive = inactive @max_lost_particles.setter def max_lost_particles(self, max_lost_particles: int): - cv.check_type("max_lost_particles", max_lost_particles, Integral) - cv.check_greater_than("max_lost_particles", max_lost_particles, 0) + cv.check_type('max_lost_particles', max_lost_particles, Integral) + cv.check_greater_than('max_lost_particles', max_lost_particles, 0) self._max_lost_particles = max_lost_particles @rel_max_lost_particles.setter def rel_max_lost_particles(self, rel_max_lost_particles: float): - cv.check_type("rel_max_lost_particles", rel_max_lost_particles, Real) - cv.check_greater_than("rel_max_lost_particles", rel_max_lost_particles, 0) - cv.check_less_than("rel_max_lost_particles", rel_max_lost_particles, 1) + cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real) + cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0) + cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1) self._rel_max_lost_particles = rel_max_lost_particles @particles.setter def particles(self, particles: int): - cv.check_type("particles", particles, Integral) - cv.check_greater_than("particles", particles, 0) + cv.check_type('particles', particles, Integral) + cv.check_greater_than('particles', particles, 0) self._particles = particles @keff_trigger.setter def keff_trigger(self, keff_trigger: dict): if not isinstance(keff_trigger, dict): - msg = ( - f'Unable to set a trigger on keff from "{keff_trigger}" ' - "which is not a Python dictionary" - ) + msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ + 'which is not a Python dictionary' raise ValueError(msg) - elif "type" not in keff_trigger: - msg = ( - f'Unable to set a trigger on keff from "{keff_trigger}" ' - 'which does not have a "type" key' - ) + elif 'type' not in keff_trigger: + msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ + 'which does not have a "type" key' raise ValueError(msg) - elif keff_trigger["type"] not in ["variance", "std_dev", "rel_err"]: - msg = "Unable to set a trigger on keff with " 'type "{0}"'.format( - keff_trigger["type"] - ) + elif keff_trigger['type'] not in ['variance', 'std_dev', 'rel_err']: + msg = 'Unable to set a trigger on keff with ' \ + 'type "{0}"'.format(keff_trigger['type']) raise ValueError(msg) - elif "threshold" not in keff_trigger: - msg = ( - f'Unable to set a trigger on keff from "{keff_trigger}" ' - 'which does not have a "threshold" key' - ) + elif 'threshold' not in keff_trigger: + msg = f'Unable to set a trigger on keff from "{keff_trigger}" ' \ + 'which does not have a "threshold" key' raise ValueError(msg) - elif not isinstance(keff_trigger["threshold"], Real): - msg = "Unable to set a trigger on keff with " 'threshold "{0}"'.format( - keff_trigger["threshold"] - ) + elif not isinstance(keff_trigger['threshold'], Real): + msg = 'Unable to set a trigger on keff with ' \ + 'threshold "{0}"'.format(keff_trigger['threshold']) raise ValueError(msg) self._keff_trigger = keff_trigger @energy_mode.setter def energy_mode(self, energy_mode: str): - cv.check_value("energy mode", energy_mode, ["continuous-energy", "multi-group"]) + cv.check_value('energy mode', energy_mode, + ['continuous-energy', 'multi-group']) self._energy_mode = energy_mode @max_order.setter def max_order(self, max_order: Optional[int]): if max_order is not None: - cv.check_type("maximum scattering order", max_order, Integral) - cv.check_greater_than("maximum scattering order", max_order, 0, True) + cv.check_type('maximum scattering order', max_order, Integral) + cv.check_greater_than('maximum scattering order', max_order, 0, + True) self._max_order = max_order @source.setter def source(self, source: Union[Source, typing.Iterable[Source]]): if not isinstance(source, MutableSequence): source = [source] - self._source = cv.CheckedList(Source, "source distributions", source) + self._source = cv.CheckedList(Source, 'source distributions', source) @output.setter def output(self, output: dict): - cv.check_type("output", output, Mapping) + cv.check_type('output', output, Mapping) for key, value in output.items(): - cv.check_value("output key", key, ("summary", "tallies", "path")) - if key in ("summary", "tallies"): + cv.check_value('output key', key, ('summary', 'tallies', 'path')) + if key in ('summary', 'tallies'): cv.check_type(f"output['{key}']", value, bool) else: cv.check_type("output['path']", value, str) @@ -609,257 +602,245 @@ class Settings: @verbosity.setter def verbosity(self, verbosity: int): - cv.check_type("verbosity", verbosity, Integral) - cv.check_greater_than("verbosity", verbosity, 1, True) - cv.check_less_than("verbosity", verbosity, 10, True) + cv.check_type('verbosity', verbosity, Integral) + cv.check_greater_than('verbosity', verbosity, 1, True) + cv.check_less_than('verbosity', verbosity, 10, True) self._verbosity = verbosity @sourcepoint.setter def sourcepoint(self, sourcepoint: dict): - cv.check_type("sourcepoint options", sourcepoint, Mapping) + cv.check_type('sourcepoint options', sourcepoint, Mapping) for key, value in sourcepoint.items(): - if key == "batches": - cv.check_type("sourcepoint batches", value, Iterable, Integral) + if key == 'batches': + cv.check_type('sourcepoint batches', value, Iterable, Integral) for batch in value: - cv.check_greater_than("sourcepoint batch", batch, 0) - elif key == "separate": - cv.check_type("sourcepoint separate", value, bool) - elif key == "write": - cv.check_type("sourcepoint write", value, bool) - elif key == "overwrite": - cv.check_type("sourcepoint overwrite", value, bool) + cv.check_greater_than('sourcepoint batch', batch, 0) + elif key == 'separate': + cv.check_type('sourcepoint separate', value, bool) + elif key == 'write': + cv.check_type('sourcepoint write', value, bool) + elif key == 'overwrite': + cv.check_type('sourcepoint overwrite', value, bool) else: - raise ValueError( - f"Unknown key '{key}' encountered when " - "setting sourcepoint options." - ) + raise ValueError(f"Unknown key '{key}' encountered when " + "setting sourcepoint options.") self._sourcepoint = sourcepoint @statepoint.setter def statepoint(self, statepoint: dict): - cv.check_type("statepoint options", statepoint, Mapping) + cv.check_type('statepoint options', statepoint, Mapping) for key, value in statepoint.items(): - if key == "batches": - cv.check_type("statepoint batches", value, Iterable, Integral) + if key == 'batches': + cv.check_type('statepoint batches', value, Iterable, Integral) for batch in value: - cv.check_greater_than("statepoint batch", batch, 0) + cv.check_greater_than('statepoint batch', batch, 0) else: - raise ValueError( - f"Unknown key '{key}' encountered when " - "setting statepoint options." - ) + raise ValueError(f"Unknown key '{key}' encountered when " + "setting statepoint options.") self._statepoint = statepoint @surf_source_read.setter def surf_source_read(self, surf_source_read: dict): - cv.check_type("surface source reading options", surf_source_read, Mapping) + cv.check_type('surface source reading options', surf_source_read, Mapping) for key, value in surf_source_read.items(): - cv.check_value("surface source reading key", key, ("path")) - if key == "path": - cv.check_type("path to surface source file", value, str) + cv.check_value('surface source reading key', key, + ('path')) + if key == 'path': + cv.check_type('path to surface source file', value, str) self._surf_source_read = surf_source_read @surf_source_write.setter def surf_source_write(self, surf_source_write: dict): - cv.check_type("surface source writing options", surf_source_write, Mapping) + cv.check_type('surface source writing options', surf_source_write, Mapping) for key, value in surf_source_write.items(): - cv.check_value( - "surface source writing key", key, ("surface_ids", "max_particles") - ) - if key == "surface_ids": - cv.check_type( - "surface ids for source banking", value, Iterable, Integral - ) + cv.check_value('surface source writing key', key, + ('surface_ids', 'max_particles')) + if key == 'surface_ids': + cv.check_type('surface ids for source banking', value, + Iterable, Integral) for surf_id in value: - cv.check_greater_than("surface id for source banking", surf_id, 0) - elif key == "max_particles": - cv.check_type( - "maximum particle banks on surfaces per process", value, Integral - ) - cv.check_greater_than( - "maximum particle banks on surfaces per process", value, 0 - ) + cv.check_greater_than('surface id for source banking', + surf_id, 0) + elif key == 'max_particles': + cv.check_type('maximum particle banks on surfaces per process', + value, Integral) + cv.check_greater_than('maximum particle banks on surfaces per process', + value, 0) self._surf_source_write = surf_source_write @confidence_intervals.setter def confidence_intervals(self, confidence_intervals: bool): - cv.check_type("confidence interval", confidence_intervals, bool) + cv.check_type('confidence interval', confidence_intervals, bool) self._confidence_intervals = confidence_intervals @electron_treatment.setter def electron_treatment(self, electron_treatment: str): - cv.check_value("electron treatment", electron_treatment, ["led", "ttb"]) + cv.check_value('electron treatment', electron_treatment, ['led', 'ttb']) self._electron_treatment = electron_treatment @photon_transport.setter def photon_transport(self, photon_transport: bool): - cv.check_type("photon transport", photon_transport, bool) + cv.check_type('photon transport', photon_transport, bool) self._photon_transport = photon_transport @ptables.setter def ptables(self, ptables: bool): - cv.check_type("probability tables", ptables, bool) + cv.check_type('probability tables', ptables, bool) self._ptables = ptables @seed.setter def seed(self, seed: int): - cv.check_type("random number generator seed", seed, Integral) - cv.check_greater_than("random number generator seed", seed, 0) + cv.check_type('random number generator seed', seed, Integral) + cv.check_greater_than('random number generator seed', seed, 0) self._seed = seed @survival_biasing.setter def survival_biasing(self, survival_biasing: bool): - cv.check_type("survival biasing", survival_biasing, bool) + cv.check_type('survival biasing', survival_biasing, bool) self._survival_biasing = survival_biasing @cutoff.setter def cutoff(self, cutoff: dict): if not isinstance(cutoff, Mapping): - msg = ( - f'Unable to set cutoff from "{cutoff}" which is not a ' - "Python dictionary" - ) + msg = f'Unable to set cutoff from "{cutoff}" which is not a '\ + 'Python dictionary' raise ValueError(msg) for key in cutoff: - if key == "weight": - cv.check_type("weight cutoff", cutoff[key], Real) - cv.check_greater_than("weight cutoff", cutoff[key], 0.0) - elif key == "weight_avg": - cv.check_type("average survival weight", cutoff[key], Real) - cv.check_greater_than("average survival weight", cutoff[key], 0.0) - elif key in [ - "energy_neutron", - "energy_photon", - "energy_electron", - "energy_positron", - ]: - cv.check_type("energy cutoff", cutoff[key], Real) - cv.check_greater_than("energy cutoff", cutoff[key], 0.0) + if key == 'weight': + cv.check_type('weight cutoff', cutoff[key], Real) + cv.check_greater_than('weight cutoff', cutoff[key], 0.0) + elif key == 'weight_avg': + cv.check_type('average survival weight', cutoff[key], Real) + cv.check_greater_than('average survival weight', + cutoff[key], 0.0) + elif key in ['energy_neutron', 'energy_photon', 'energy_electron', + 'energy_positron']: + cv.check_type('energy cutoff', cutoff[key], Real) + cv.check_greater_than('energy cutoff', cutoff[key], 0.0) else: - msg = ( - f'Unable to set cutoff to "{key}" which is unsupported ' "by OpenMC" - ) + msg = f'Unable to set cutoff to "{key}" which is unsupported ' \ + 'by OpenMC' self._cutoff = cutoff @entropy_mesh.setter def entropy_mesh(self, entropy: RegularMesh): - cv.check_type("entropy mesh", entropy, RegularMesh) + cv.check_type('entropy mesh', entropy, RegularMesh) self._entropy_mesh = entropy @trigger_active.setter def trigger_active(self, trigger_active: bool): - cv.check_type("trigger active", trigger_active, bool) + cv.check_type('trigger active', trigger_active, bool) self._trigger_active = trigger_active @trigger_max_batches.setter def trigger_max_batches(self, trigger_max_batches: int): - cv.check_type("trigger maximum batches", trigger_max_batches, Integral) - cv.check_greater_than("trigger maximum batches", trigger_max_batches, 0) + cv.check_type('trigger maximum batches', trigger_max_batches, Integral) + cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0) self._trigger_max_batches = trigger_max_batches @trigger_batch_interval.setter def trigger_batch_interval(self, trigger_batch_interval: int): - cv.check_type("trigger batch interval", trigger_batch_interval, Integral) - cv.check_greater_than("trigger batch interval", trigger_batch_interval, 0) + cv.check_type('trigger batch interval', trigger_batch_interval, Integral) + cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0) self._trigger_batch_interval = trigger_batch_interval @no_reduce.setter def no_reduce(self, no_reduce: bool): - cv.check_type("no reduction option", no_reduce, bool) + cv.check_type('no reduction option', no_reduce, bool) self._no_reduce = no_reduce @tabular_legendre.setter def tabular_legendre(self, tabular_legendre: dict): - cv.check_type("tabular_legendre settings", tabular_legendre, Mapping) + cv.check_type('tabular_legendre settings', tabular_legendre, Mapping) for key, value in tabular_legendre.items(): - cv.check_value("tabular_legendre key", key, ["enable", "num_points"]) - if key == "enable": - cv.check_type("enable tabular_legendre", value, bool) - elif key == "num_points": - cv.check_type("num_points tabular_legendre", value, Integral) - cv.check_greater_than("num_points tabular_legendre", value, 0) + cv.check_value('tabular_legendre key', key, + ['enable', 'num_points']) + if key == 'enable': + cv.check_type('enable tabular_legendre', value, bool) + elif key == 'num_points': + cv.check_type('num_points tabular_legendre', value, Integral) + cv.check_greater_than('num_points tabular_legendre', value, 0) self._tabular_legendre = tabular_legendre @temperature.setter def temperature(self, temperature: dict): - cv.check_type("temperature settings", temperature, Mapping) + cv.check_type('temperature settings', temperature, Mapping) for key, value in temperature.items(): - cv.check_value( - "temperature key", - key, - ["default", "method", "tolerance", "multipole", "range"], - ) - if key == "default": - cv.check_type("default temperature", value, Real) - elif key == "method": - cv.check_value( - "temperature method", value, ["nearest", "interpolation"] - ) - elif key == "tolerance": - cv.check_type("temperature tolerance", value, Real) - elif key == "multipole": - cv.check_type("temperature multipole", value, bool) - elif key == "range": - cv.check_length("temperature range", value, 2) + cv.check_value('temperature key', key, + ['default', 'method', 'tolerance', 'multipole', + 'range']) + if key == 'default': + cv.check_type('default temperature', value, Real) + elif key == 'method': + cv.check_value('temperature method', value, + ['nearest', 'interpolation']) + elif key == 'tolerance': + cv.check_type('temperature tolerance', value, Real) + elif key == 'multipole': + cv.check_type('temperature multipole', value, bool) + elif key == 'range': + cv.check_length('temperature range', value, 2) for T in value: - cv.check_type("temperature", T, Real) + cv.check_type('temperature', T, Real) self._temperature = temperature @trace.setter def trace(self, trace: Iterable): - cv.check_type("trace", trace, Iterable, Integral) - cv.check_length("trace", trace, 3) - cv.check_greater_than("trace batch", trace[0], 0) - cv.check_greater_than("trace generation", trace[1], 0) - cv.check_greater_than("trace particle", trace[2], 0) + cv.check_type('trace', trace, Iterable, Integral) + cv.check_length('trace', trace, 3) + cv.check_greater_than('trace batch', trace[0], 0) + cv.check_greater_than('trace generation', trace[1], 0) + cv.check_greater_than('trace particle', trace[2], 0) self._trace = trace @track.setter def track(self, track: typing.Iterable[typing.Iterable[int]]): - cv.check_type("track", track, Iterable) + cv.check_type('track', track, Iterable) for t in track: if len(t) != 3: msg = f'Unable to set the track to "{t}" since its length is not 3' raise ValueError(msg) - cv.check_greater_than("track batch", t[0], 0) - cv.check_greater_than("track generation", t[1], 0) - cv.check_greater_than("track particle", t[2], 0) - cv.check_type("track batch", t[0], Integral) - cv.check_type("track generation", t[1], Integral) - cv.check_type("track particle", t[2], Integral) + cv.check_greater_than('track batch', t[0], 0) + cv.check_greater_than('track generation', t[1], 0) + cv.check_greater_than('track particle', t[2], 0) + cv.check_type('track batch', t[0], Integral) + cv.check_type('track generation', t[1], Integral) + cv.check_type('track particle', t[2], Integral) self._track = track @ufs_mesh.setter def ufs_mesh(self, ufs_mesh: RegularMesh): - cv.check_type("UFS mesh", ufs_mesh, RegularMesh) - cv.check_length("UFS mesh dimension", ufs_mesh.dimension, 3) - cv.check_length("UFS mesh lower-left corner", ufs_mesh.lower_left, 3) - cv.check_length("UFS mesh upper-right corner", ufs_mesh.upper_right, 3) + cv.check_type('UFS mesh', ufs_mesh, RegularMesh) + cv.check_length('UFS mesh dimension', ufs_mesh.dimension, 3) + cv.check_length('UFS mesh lower-left corner', ufs_mesh.lower_left, 3) + cv.check_length('UFS mesh upper-right corner', ufs_mesh.upper_right, 3) self._ufs_mesh = ufs_mesh @resonance_scattering.setter def resonance_scattering(self, res: dict): - cv.check_type("resonance scattering settings", res, Mapping) - keys = ("enable", "method", "energy_min", "energy_max", "nuclides") + cv.check_type('resonance scattering settings', res, Mapping) + keys = ('enable', 'method', 'energy_min', 'energy_max', 'nuclides') for key, value in res.items(): - cv.check_value("resonance scattering dictionary key", key, keys) - if key == "enable": - cv.check_type("resonance scattering enable", value, bool) - elif key == "method": - cv.check_value("resonance scattering method", value, _RES_SCAT_METHODS) - elif key == "energy_min": - name = "resonance scattering minimum energy" + cv.check_value('resonance scattering dictionary key', key, keys) + if key == 'enable': + cv.check_type('resonance scattering enable', value, bool) + elif key == 'method': + cv.check_value('resonance scattering method', value, + _RES_SCAT_METHODS) + elif key == 'energy_min': + name = 'resonance scattering minimum energy' cv.check_type(name, value, Real) cv.check_greater_than(name, value, 0) - elif key == "energy_max": - name = "resonance scattering minimum energy" + elif key == 'energy_max': + name = 'resonance scattering minimum energy' cv.check_type(name, value, Real) cv.check_greater_than(name, value, 0) - elif key == "nuclides": - cv.check_type("resonance scattering nuclides", value, Iterable, str) + elif key == 'nuclides': + cv.check_type('resonance scattering nuclides', value, + Iterable, str) self._resonance_scattering = res @volume_calculations.setter @@ -869,69 +850,67 @@ class Settings: if not isinstance(vol_calcs, MutableSequence): vol_calcs = [vol_calcs] self._volume_calculations = cv.CheckedList( - VolumeCalculation, "stochastic volume calculations", vol_calcs - ) + VolumeCalculation, 'stochastic volume calculations', vol_calcs) @create_fission_neutrons.setter def create_fission_neutrons(self, create_fission_neutrons: bool): - cv.check_type("Whether create fission neutrons", create_fission_neutrons, bool) + cv.check_type('Whether create fission neutrons', + create_fission_neutrons, bool) self._create_fission_neutrons = create_fission_neutrons @delayed_photon_scaling.setter def delayed_photon_scaling(self, value: bool): - cv.check_type("delayed photon scaling", value, bool) + cv.check_type('delayed photon scaling', value, bool) self._delayed_photon_scaling = value @event_based.setter def event_based(self, value: bool): - cv.check_type("event based", value, bool) + cv.check_type('event based', value, bool) self._event_based = value @max_particles_in_flight.setter def max_particles_in_flight(self, value: int): - cv.check_type("max particles in flight", value, Integral) - cv.check_greater_than("max particles in flight", value, 0) + cv.check_type('max particles in flight', value, Integral) + cv.check_greater_than('max particles in flight', value, 0) self._max_particles_in_flight = value @material_cell_offsets.setter def material_cell_offsets(self, value: bool): - cv.check_type("material cell offsets", value, bool) + cv.check_type('material cell offsets', value, bool) self._material_cell_offsets = value @log_grid_bins.setter def log_grid_bins(self, log_grid_bins: int): - cv.check_type("log grid bins", log_grid_bins, Real) - cv.check_greater_than("log grid bins", log_grid_bins, 0) + cv.check_type('log grid bins', log_grid_bins, Real) + cv.check_greater_than('log grid bins', log_grid_bins, 0) self._log_grid_bins = log_grid_bins @write_initial_source.setter def write_initial_source(self, value: bool): - cv.check_type("write initial source", value, bool) + cv.check_type('write initial source', value, bool) self._write_initial_source = value @weight_windows.setter - def weight_windows( - self, value: Union[WeightWindows, typing.Iterable[WeightWindows]] - ): + def weight_windows(self, value: Union[WeightWindows, typing.Iterable[WeightWindows]]): if not isinstance(value, MutableSequence): value = [value] - self._weight_windows = cv.CheckedList(WeightWindows, "weight windows", value) + self._weight_windows = cv.CheckedList(WeightWindows, 'weight windows', value) @weight_windows_on.setter def weight_windows_on(self, value): - cv.check_type("weight windows on", value, bool) + cv.check_type('weight windows on', value, bool) self._weight_windows_on = value @max_splits.setter def max_splits(self, value: int): - cv.check_type("maximum particle splits", value, Integral) - cv.check_greater_than("max particle splits", value, 0) + cv.check_type('maximum particle splits', value, Integral) + cv.check_greater_than('max particle splits', value, 0) self._max_splits = value @max_tracks.setter def max_tracks(self, value: int): - cv.check_type("maximum particle tracks", value, Integral) - cv.check_greater_than("maximum particle tracks", value, 0, True) + cv.check_type('maximum particle tracks', value, Integral) + cv.check_greater_than('maximum particle tracks', value, 0, True) self._max_tracks = value def _create_run_mode_subelement(self, root): @@ -998,7 +977,7 @@ class Settings: element = ET.SubElement(root, "output") for key, value in sorted(self._output.items()): subelement = ET.SubElement(element, key) - if key in ("summary", "tallies"): + if key in ('summary', 'tallies'): subelement.text = str(value).lower() else: subelement.text = value @@ -1011,49 +990,50 @@ class Settings: def _create_statepoint_subelement(self, root): if self._statepoint: element = ET.SubElement(root, "state_point") - if "batches" in self._statepoint: + if 'batches' in self._statepoint: subelement = ET.SubElement(element, "batches") - subelement.text = " ".join(str(x) for x in self._statepoint["batches"]) + subelement.text = ' '.join( + str(x) for x in self._statepoint['batches']) def _create_sourcepoint_subelement(self, root): if self._sourcepoint: element = ET.SubElement(root, "source_point") - if "batches" in self._sourcepoint: + if 'batches' in self._sourcepoint: subelement = ET.SubElement(element, "batches") - subelement.text = " ".join(str(x) for x in self._sourcepoint["batches"]) + subelement.text = ' '.join( + str(x) for x in self._sourcepoint['batches']) - if "separate" in self._sourcepoint: + if 'separate' in self._sourcepoint: subelement = ET.SubElement(element, "separate") - subelement.text = str(self._sourcepoint["separate"]).lower() + subelement.text = str(self._sourcepoint['separate']).lower() - if "write" in self._sourcepoint: + if 'write' in self._sourcepoint: subelement = ET.SubElement(element, "write") - subelement.text = str(self._sourcepoint["write"]).lower() + subelement.text = str(self._sourcepoint['write']).lower() # Overwrite latest subelement - if "overwrite" in self._sourcepoint: + if 'overwrite' in self._sourcepoint: subelement = ET.SubElement(element, "overwrite_latest") - subelement.text = str(self._sourcepoint["overwrite"]).lower() + subelement.text = str(self._sourcepoint['overwrite']).lower() def _create_surf_source_read_subelement(self, root): if self._surf_source_read: element = ET.SubElement(root, "surf_source_read") - if "path" in self._surf_source_read: + if 'path' in self._surf_source_read: subelement = ET.SubElement(element, "path") - subelement.text = self._surf_source_read["path"] + subelement.text = self._surf_source_read['path'] def _create_surf_source_write_subelement(self, root): if self._surf_source_write: element = ET.SubElement(root, "surf_source_write") - if "surface_ids" in self._surf_source_write: + if 'surface_ids' in self._surf_source_write: subelement = ET.SubElement(element, "surface_ids") - subelement.text = " ".join( - str(x) for x in self._surf_source_write["surface_ids"] - ) - if "max_particles" in self._surf_source_write: + subelement.text = ' '.join( + str(x) for x in self._surf_source_write['surface_ids']) + if 'max_particles' in self._surf_source_write: subelement = ET.SubElement(element, "max_particles") - subelement.text = str(self._surf_source_write["max_particles"]) + subelement.text = str(self._surf_source_write['max_particles']) def _create_confidence_intervals(self, root): if self._confidence_intervals is not None: @@ -1097,14 +1077,12 @@ class Settings: # use default heuristic for entropy mesh if not set by user if self.entropy_mesh.dimension is None: if self.particles is None: - raise RuntimeError( - "Number of particles must be set in order to " - "use entropy mesh dimension heuristic" - ) + raise RuntimeError("Number of particles must be set in order to " \ + "use entropy mesh dimension heuristic") else: - n = ceil((self.particles / 20.0) ** (1.0 / 3.0)) + n = ceil((self.particles / 20.0)**(1.0 / 3.0)) d = len(self.entropy_mesh.lower_left) - self.entropy_mesh.dimension = (n,) * d + self.entropy_mesh.dimension = (n,)*d # See if a element already exists -- if not, add it path = f"./mesh[@id='{self.entropy_mesh.id}']" @@ -1137,10 +1115,10 @@ class Settings: if self.tabular_legendre: element = ET.SubElement(root, "tabular_legendre") subelement = ET.SubElement(element, "enable") - subelement.text = str(self._tabular_legendre["enable"]).lower() - if "num_points" in self._tabular_legendre: + subelement.text = str(self._tabular_legendre['enable']).lower() + if 'num_points' in self._tabular_legendre: subelement = ET.SubElement(element, "num_points") - subelement.text = str(self._tabular_legendre["num_points"]) + subelement.text = str(self._tabular_legendre['num_points']) def _create_temperature_subelements(self, root): if self.temperature: @@ -1148,20 +1126,20 @@ class Settings: element = ET.SubElement(root, f"temperature_{key}") if isinstance(value, bool): element.text = str(value).lower() - elif key == "range": - element.text = " ".join(str(T) for T in value) + elif key == 'range': + element.text = ' '.join(str(T) for T in value) else: element.text = str(value) def _create_trace_subelement(self, root): if self._trace is not None: element = ET.SubElement(root, "trace") - element.text = " ".join(map(str, self._trace)) + element.text = ' '.join(map(str, self._trace)) def _create_track_subelement(self, root): if self._track is not None: element = ET.SubElement(root, "track") - element.text = " ".join(map(str, itertools.chain(*self._track))) + element.text = ' '.join(map(str, itertools.chain(*self._track))) def _create_ufs_mesh_subelement(self, root): if self.ufs_mesh is not None: @@ -1176,22 +1154,22 @@ class Settings: def _create_resonance_scattering_subelement(self, root): res = self.resonance_scattering if res: - elem = ET.SubElement(root, "resonance_scattering") - if "enable" in res: - subelem = ET.SubElement(elem, "enable") - subelem.text = str(res["enable"]).lower() - if "method" in res: - subelem = ET.SubElement(elem, "method") - subelem.text = res["method"] - if "energy_min" in res: - subelem = ET.SubElement(elem, "energy_min") - subelem.text = str(res["energy_min"]) - if "energy_max" in res: - subelem = ET.SubElement(elem, "energy_max") - subelem.text = str(res["energy_max"]) - if "nuclides" in res: - subelem = ET.SubElement(elem, "nuclides") - subelem.text = " ".join(res["nuclides"]) + elem = ET.SubElement(root, 'resonance_scattering') + if 'enable' in res: + subelem = ET.SubElement(elem, 'enable') + subelem.text = str(res['enable']).lower() + if 'method' in res: + subelem = ET.SubElement(elem, 'method') + subelem.text = res['method'] + if 'energy_min' in res: + subelem = ET.SubElement(elem, 'energy_min') + subelem.text = str(res['energy_min']) + if 'energy_max' in res: + subelem = ET.SubElement(elem, 'energy_max') + subelem.text = str(res['energy_max']) + if 'nuclides' in res: + subelem = ET.SubElement(elem, 'nuclides') + subelem.text = ' '.join(res['nuclides']) def _create_create_fission_neutrons_subelement(self, root): if self._create_fission_neutrons is not None: @@ -1253,7 +1231,7 @@ class Settings: elem.text = str(self._max_tracks) def _eigenvalue_from_xml_element(self, root): - elem = root.find("eigenvalue") + elem = root.find('eigenvalue') if elem is not None: self._run_mode_from_xml_element(elem) self._particles_from_xml_element(elem) @@ -1264,168 +1242,161 @@ class Settings: self._generations_per_batch_from_xml_element(elem) def _run_mode_from_xml_element(self, root): - text = get_text(root, "run_mode") + text = get_text(root, 'run_mode') if text is not None: self.run_mode = text def _particles_from_xml_element(self, root): - text = get_text(root, "particles") + text = get_text(root, 'particles') if text is not None: self.particles = int(text) def _batches_from_xml_element(self, root): - text = get_text(root, "batches") + text = get_text(root, 'batches') if text is not None: self.batches = int(text) def _inactive_from_xml_element(self, root): - text = get_text(root, "inactive") + text = get_text(root, 'inactive') if text is not None: self.inactive = int(text) def _max_lost_particles_from_xml_element(self, root): - text = get_text(root, "max_lost_particles") + text = get_text(root, 'max_lost_particles') if text is not None: self.max_lost_particles = int(text) def _rel_max_lost_particles_from_xml_element(self, root): - text = get_text(root, "rel_max_lost_particles") + text = get_text(root, 'rel_max_lost_particles') if text is not None: self.rel_max_lost_particles = float(text) def _generations_per_batch_from_xml_element(self, root): - text = get_text(root, "generations_per_batch") + text = get_text(root, 'generations_per_batch') if text is not None: self.generations_per_batch = int(text) def _keff_trigger_from_xml_element(self, root): - elem = root.find("keff_trigger") + elem = root.find('keff_trigger') if elem is not None: - trigger = get_text(elem, "type") - threshold = float(get_text(elem, "threshold")) - self.keff_trigger = {"type": trigger, "threshold": threshold} + trigger = get_text(elem, 'type') + threshold = float(get_text(elem, 'threshold')) + self.keff_trigger = {'type': trigger, 'threshold': threshold} def _source_from_xml_element(self, root): - for elem in root.findall("source"): + for elem in root.findall('source'): self.source.append(Source.from_xml_element(elem)) def _volume_calcs_from_xml_element(self, root): volume_elems = root.findall("volume_calc") if volume_elems: - self.volume_calculations = [ - VolumeCalculation.from_xml_element(elem) for elem in volume_elems - ] + self.volume_calculations = [VolumeCalculation.from_xml_element(elem) + for elem in volume_elems] def _output_from_xml_element(self, root): - elem = root.find("output") + elem = root.find('output') if elem is not None: self.output = {} - for key in ("summary", "tallies", "path"): + for key in ('summary', 'tallies', 'path'): value = get_text(elem, key) if value is not None: - if key in ("summary", "tallies"): - value = value in ("true", "1") + if key in ('summary', 'tallies'): + value = value in ('true', '1') self.output[key] = value def _statepoint_from_xml_element(self, root): - elem = root.find("state_point") + elem = root.find('state_point') if elem is not None: - text = get_text(elem, "batches") + text = get_text(elem, 'batches') if text is not None: - self.statepoint["batches"] = [int(x) for x in text.split()] + self.statepoint['batches'] = [int(x) for x in text.split()] def _sourcepoint_from_xml_element(self, root): - elem = root.find("source_point") + elem = root.find('source_point') if elem is not None: - for key in ("separate", "write", "overwrite_latest", "batches"): + for key in ('separate', 'write', 'overwrite_latest', 'batches'): value = get_text(elem, key) if value is not None: - if key in ("separate", "write"): - value = value in ("true", "1") - elif key == "overwrite_latest": - value = value in ("true", "1") - key = "overwrite" + if key in ('separate', 'write'): + value = value in ('true', '1') + elif key == 'overwrite_latest': + value = value in ('true', '1') + key = 'overwrite' else: value = [int(x) for x in value.split()] self.sourcepoint[key] = value def _surf_source_read_from_xml_element(self, root): - elem = root.find("surf_source_read") + elem = root.find('surf_source_read') if elem is not None: - value = get_text(elem, "path") + value = get_text(elem, 'path') if value is not None: - self.surf_source_read["path"] = value + self.surf_source_read['path'] = value def _surf_source_write_from_xml_element(self, root): - elem = root.find("surf_source_write") + elem = root.find('surf_source_write') if elem is not None: - for key in ("surface_ids", "max_particles"): + for key in ('surface_ids', 'max_particles'): value = get_text(elem, key) if value is not None: - if key == "surface_ids": + if key == 'surface_ids': value = [int(x) for x in value.split()] - elif key in ("max_particles"): + elif key in ('max_particles'): value = int(value) self.surf_source_write[key] = value def _confidence_intervals_from_xml_element(self, root): - text = get_text(root, "confidence_intervals") + text = get_text(root, 'confidence_intervals') if text is not None: - self.confidence_intervals = text in ("true", "1") + self.confidence_intervals = text in ('true', '1') def _electron_treatment_from_xml_element(self, root): - text = get_text(root, "electron_treatment") + text = get_text(root, 'electron_treatment') if text is not None: self.electron_treatment = text def _energy_mode_from_xml_element(self, root): - text = get_text(root, "energy_mode") + text = get_text(root, 'energy_mode') if text is not None: self.energy_mode = text def _max_order_from_xml_element(self, root): - text = get_text(root, "max_order") + text = get_text(root, 'max_order') if text is not None: self.max_order = int(text) def _photon_transport_from_xml_element(self, root): - text = get_text(root, "photon_transport") + text = get_text(root, 'photon_transport') if text is not None: - self.photon_transport = text in ("true", "1") + self.photon_transport = text in ('true', '1') def _ptables_from_xml_element(self, root): - text = get_text(root, "ptables") + text = get_text(root, 'ptables') if text is not None: - self.ptables = text in ("true", "1") + self.ptables = text in ('true', '1') def _seed_from_xml_element(self, root): - text = get_text(root, "seed") + text = get_text(root, 'seed') if text is not None: self.seed = int(text) def _survival_biasing_from_xml_element(self, root): - text = get_text(root, "survival_biasing") + text = get_text(root, 'survival_biasing') if text is not None: - self.survival_biasing = text in ("true", "1") + self.survival_biasing = text in ('true', '1') def _cutoff_from_xml_element(self, root): - elem = root.find("cutoff") + elem = root.find('cutoff') if elem is not None: self.cutoff = {} - for key in ( - "energy_neutron", - "energy_photon", - "energy_electron", - "energy_positron", - "weight", - "weight_avg", - ): + for key in ('energy_neutron', 'energy_photon', 'energy_electron', + 'energy_positron', 'weight', 'weight_avg'): value = get_text(elem, key) if value is not None: self.cutoff[key] = float(value) def _entropy_mesh_from_xml_element(self, root): - text = get_text(root, "entropy_mesh") + text = get_text(root, 'entropy_mesh') if text is not None: path = f"./mesh[@id='{int(text)}']" elem = root.find(path) @@ -1433,65 +1404,65 @@ class Settings: self.entropy_mesh = RegularMesh.from_xml_element(elem) def _trigger_from_xml_element(self, root): - elem = root.find("trigger") + elem = root.find('trigger') if elem is not None: - self.trigger_active = get_text(elem, "active") in ("true", "1") - text = get_text(elem, "max_batches") + self.trigger_active = get_text(elem, 'active') in ('true', '1') + text = get_text(elem, 'max_batches') if text is not None: self.trigger_max_batches = int(text) - text = get_text(elem, "batch_interval") + text = get_text(elem, 'batch_interval') if text is not None: self.trigger_batch_interval = int(text) def _no_reduce_from_xml_element(self, root): - text = get_text(root, "no_reduce") + text = get_text(root, 'no_reduce') if text is not None: - self.no_reduce = text in ("true", "1") + self.no_reduce = text in ('true', '1') def _verbosity_from_xml_element(self, root): - text = get_text(root, "verbosity") + text = get_text(root, 'verbosity') if text is not None: self.verbosity = int(text) def _tabular_legendre_from_xml_element(self, root): - elem = root.find("tabular_legendre") + elem = root.find('tabular_legendre') if elem is not None: - text = get_text(elem, "enable") - self.tabular_legendre["enable"] = text in ("true", "1") - text = get_text(elem, "num_points") + text = get_text(elem, 'enable') + self.tabular_legendre['enable'] = text in ('true', '1') + text = get_text(elem, 'num_points') if text is not None: - self.tabular_legendre["num_points"] = int(text) + self.tabular_legendre['num_points'] = int(text) def _temperature_from_xml_element(self, root): - text = get_text(root, "temperature_default") + text = get_text(root, 'temperature_default') if text is not None: - self.temperature["default"] = float(text) - text = get_text(root, "temperature_tolerance") + self.temperature['default'] = float(text) + text = get_text(root, 'temperature_tolerance') if text is not None: - self.temperature["tolerance"] = float(text) - text = get_text(root, "temperature_method") + self.temperature['tolerance'] = float(text) + text = get_text(root, 'temperature_method') if text is not None: - self.temperature["method"] = text - text = get_text(root, "temperature_range") + self.temperature['method'] = text + text = get_text(root, 'temperature_range') if text is not None: - self.temperature["range"] = [float(x) for x in text.split()] - text = get_text(root, "temperature_multipole") + self.temperature['range'] = [float(x) for x in text.split()] + text = get_text(root, 'temperature_multipole') if text is not None: - self.temperature["multipole"] = text in ("true", "1") + self.temperature['multipole'] = text in ('true', '1') def _trace_from_xml_element(self, root): - text = get_text(root, "trace") + text = get_text(root, 'trace') if text is not None: self.trace = [int(x) for x in text.split()] def _track_from_xml_element(self, root): - text = get_text(root, "track") + text = get_text(root, 'track') if text is not None: values = [int(x) for x in text.split()] self.track = list(zip(values[::3], values[1::3], values[2::3])) def _ufs_mesh_from_xml_element(self, root): - text = get_text(root, "ufs_mesh") + text = get_text(root, 'ufs_mesh') if text is not None: path = f"./mesh[@id='{int(text)}']" elem = root.find(path) @@ -1499,82 +1470,80 @@ class Settings: self.ufs_mesh = RegularMesh.from_xml_element(elem) def _resonance_scattering_from_xml_element(self, root): - elem = root.find("resonance_scattering") + elem = root.find('resonance_scattering') if elem is not None: - keys = ("enable", "method", "energy_min", "energy_max", "nuclides") + keys = ('enable', 'method', 'energy_min', 'energy_max', 'nuclides') for key in keys: value = get_text(elem, key) if value is not None: - if key == "enable": - value = value in ("true", "1") - elif key in ("energy_min", "energy_max"): + if key == 'enable': + value = value in ('true', '1') + elif key in ('energy_min', 'energy_max'): value = float(value) - elif key == "nuclides": + elif key == 'nuclides': value = value.split() self.resonance_scattering[key] = value def _create_fission_neutrons_from_xml_element(self, root): - text = get_text(root, "create_fission_neutrons") + text = get_text(root, 'create_fission_neutrons') if text is not None: - self.create_fission_neutrons = text in ("true", "1") + self.create_fission_neutrons = text in ('true', '1') def _delayed_photon_scaling_from_xml_element(self, root): - text = get_text(root, "delayed_photon_scaling") + text = get_text(root, 'delayed_photon_scaling') if text is not None: - self.delayed_photon_scaling = text in ("true", "1") + self.delayed_photon_scaling = text in ('true', '1') def _event_based_from_xml_element(self, root): - text = get_text(root, "event_based") + text = get_text(root, 'event_based') if text is not None: - self.event_based = text in ("true", "1") + self.event_based = text in ('true', '1') def _max_particles_in_flight_from_xml_element(self, root): - text = get_text(root, "max_particles_in_flight") + text = get_text(root, 'max_particles_in_flight') if text is not None: self.max_particles_in_flight = int(text) def _material_cell_offsets_from_xml_element(self, root): - text = get_text(root, "material_cell_offsets") + text = get_text(root, 'material_cell_offsets') if text is not None: - self.material_cell_offsets = text in ("true", "1") + self.material_cell_offsets = text in ('true', '1') def _log_grid_bins_from_xml_element(self, root): - text = get_text(root, "log_grid_bins") + text = get_text(root, 'log_grid_bins') if text is not None: self.log_grid_bins = int(text) def _write_initial_source_from_xml_element(self, root): - text = get_text(root, "write_initial_source") + text = get_text(root, 'write_initial_source') if text is not None: - self.write_initial_source = text in ("true", "1") + self.write_initial_source = text in ('true', '1') def _weight_windows_from_xml_element(self, root): - for elem in root.findall("weight_windows"): + for elem in root.findall('weight_windows'): ww = WeightWindows.from_xml_element(elem, root) self.weight_windows.append(ww) - text = get_text(root, "weight_windows_on") + text = get_text(root, 'weight_windows_on') if text is not None: - self.weight_windows_on = text in ("true", "1") + self.weight_windows_on = text in ('true', '1') def _max_splits_from_xml_element(self, root): - text = get_text(root, "max_splits") + text = get_text(root, 'max_splits') if text is not None: self.max_splits = int(text) def _max_tracks_from_xml_element(self, root): - text = get_text(root, "max_tracks") + text = get_text(root, 'max_tracks') if text is not None: self.max_tracks = int(text) - def export_to_xml(self, path: PathLike = "settings.xml"): + def export_to_xml(self, path: PathLike = 'settings.xml'): """Export simulation settings to an XML file. - Parameters ---------- path : str Path to file to write. Defaults to 'settings.xml'. - """ # Reset xml element tree @@ -1631,29 +1600,25 @@ class Settings: # Check if path is a directory p = Path(path) if p.is_dir(): - p /= "settings.xml" + p /= 'settings.xml' # Write the XML Tree to the settings.xml file reorder_attributes(root_element) # TODO: Remove when support is Python 3.8+ tree = ET.ElementTree(root_element) - tree.write(str(p), xml_declaration=True, encoding="utf-8") + tree.write(str(p), xml_declaration=True, encoding='utf-8') @classmethod - def from_xml(cls, path: PathLike = "settings.xml"): + def from_xml(cls, path: PathLike = 'settings.xml'): """Generate settings from XML file - .. versionadded:: 0.13.0 - Parameters ---------- path : str, optional Path to settings XML file - Returns ------- openmc.Settings Settings object - """ tree = ET.parse(path) root = tree.getroot() @@ -1707,4 +1672,4 @@ class Settings: # TODO: Get volume calculations - return settings + return settings \ No newline at end of file From f6a7b61c4a64f18165608c4114570e6ca8b061ec Mon Sep 17 00:00:00 2001 From: josh Date: Sun, 30 Oct 2022 03:31:26 +0000 Subject: [PATCH 5/7] address a few more formatting issues --- openmc/settings.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/settings.py b/openmc/settings.py index f13b233de9..e48c0d3ad2 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -1540,10 +1540,12 @@ class Settings: def export_to_xml(self, path: PathLike = 'settings.xml'): """Export simulation settings to an XML file. + Parameters ---------- path : str Path to file to write. Defaults to 'settings.xml'. + """ # Reset xml element tree @@ -1610,15 +1612,19 @@ class Settings: @classmethod def from_xml(cls, path: PathLike = 'settings.xml'): """Generate settings from XML file + .. versionadded:: 0.13.0 + Parameters ---------- path : str, optional Path to settings XML file + Returns ------- openmc.Settings Settings object + """ tree = ET.parse(path) root = tree.getroot() @@ -1672,4 +1678,4 @@ class Settings: # TODO: Get volume calculations - return settings \ No newline at end of file + return settings From 81b88cc018b17156bf08213a0e1348517ec8d428 Mon Sep 17 00:00:00 2001 From: josh Date: Wed, 2 Nov 2022 02:23:37 +0000 Subject: [PATCH 6/7] Clean up comparisons and lump test cross sections --- src/nuclide.cpp | 6 ++--- src/thermal.cpp | 6 +++-- tests/unit_tests/test_temp_interp.py | 35 ++++++++++------------------ 3 files changed, 19 insertions(+), 28 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 4c1706509b..134e2d6b9a 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -171,14 +171,14 @@ Nuclide::Nuclide(hid_t group, const vector& temperature) if (!found_pair) { // If no pairs found, check if the desired temperature falls just // outside of data - if (T_desired - temps_available.front() <= - -settings::temperature_tolerance) { + if (std::abs(T_desired - temps_available.front()) <= + settings::temperature_tolerance) { if (!contains(temps_to_read, temps_available.front())) { temps_to_read.push_back(std::round(temps_available.front())); } break; } - if (T_desired - temps_available.back() <= + if (std::abs(T_desired - temps_available.back()) <= settings::temperature_tolerance) { if (!contains(temps_to_read, temps_available.back())) { temps_to_read.push_back(std::round(temps_available.back())); diff --git a/src/thermal.cpp b/src/thermal.cpp index 66e1fb0095..b54f07b362 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -119,11 +119,13 @@ ThermalScattering::ThermalScattering( if (!found) { // If no pairs found, check if the desired temperature falls within // bounds' tolerance - if (T - temps_available[0] <= -settings::temperature_tolerance) { + if (std::abs(T - temps_available[0]) <= + settings::temperature_tolerance) { temps_to_read.push_back(std::round(temps_available[0])); break; } - if (T - temps_available[n - 1] <= settings::temperature_tolerance) { + if (std::abs(T - temps_available[n - 1]) <= + settings::temperature_tolerance) { temps_to_read.push_back(std::round(temps_available[n - 1])); break; } diff --git a/tests/unit_tests/test_temp_interp.py b/tests/unit_tests/test_temp_interp.py index c64ebe46ad..a28aafb6f8 100644 --- a/tests/unit_tests/test_temp_interp.py +++ b/tests/unit_tests/test_temp_interp.py @@ -10,7 +10,7 @@ import pytest def make_fake_cross_section(): - """Create fake U235 nuclide + """Create fake U235 nuclide with a fake thermal scattering library attached This nuclide is designed to have k_inf=1 at 300 K, k_inf=2 at 600 K, and k_inf=1 at 900 K. The absorption cross section is also constant with @@ -81,16 +81,9 @@ def make_fake_cross_section(): # Export HDF5 file u235_fake.export_to_hdf5('U235_fake.h5', 'w') - lib = openmc.data.DataLibrary() - lib.register_file('U235_fake.h5') - lib.export_to_xml('cross_sections_fake.xml') - - -def fake_thermal_scattering(): - """Create a fake thermal scattering library for U-235 at 294K and 600K - """ - fake_tsl = openmc.data.ThermalScattering("c_U_fake", 1.9968, 4.9, [0.0253]) - fake_tsl.nuclides = ['U235'] + # Create a fake thermal scattering library attached to the fake U235 data + c_U_fake = openmc.data.ThermalScattering("c_U_fake", 1.9968, 4.9, [0.0253]) + c_U_fake.nuclides = ['U235'] # Create elastic reaction bragg_edges = [0.00370672, 0.00494229] @@ -110,7 +103,7 @@ def fake_thermal_scattering(): '294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_294), '600K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist_600) } - fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) + c_U_fake.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) # Create inelastic reaction inelastic_xs = { @@ -135,19 +128,16 @@ def fake_thermal_scattering(): breakpoints, interpolation, energy, energy_out, mu) inelastic_dist = {'294K': dist, '600K': dist} inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist) - fake_tsl.inelastic = inelastic + c_U_fake.inelastic = inelastic - return fake_tsl - - -def edit_fake_cross_sections(): - """Edit the test cross sections xml to include fake thermal scattering data - """ - lib = openmc.data.DataLibrary.from_xml("cross_sections_fake.xml") - c_U_fake = fake_thermal_scattering() + # Export HDF5 file c_U_fake.export_to_hdf5("c_U_fake.h5") + + # Create a data library of the fake nuclide and its thermal scattering data + lib = openmc.data.DataLibrary() + lib.register_file('U235_fake.h5') lib.register_file("c_U_fake.h5") - lib.export_to_xml("cross_sections_fake.xml") + lib.export_to_xml('cross_sections_fake.xml') @pytest.fixture(scope='module') @@ -223,7 +213,6 @@ def test_interpolation(model, method, temperature, fission_expected, tolerance): def test_temperature_interpolation_tolerance(model): """Test applying global and cell temperatures with thermal scattering libraries """ - edit_fake_cross_sections() model.materials[0].add_s_alpha_beta("c_U_fake") # Default k-effective, using the thermal scattering data's minimum available temperature From dc5f2746d3a2c0355402e3141f55dd0cee2ce4d6 Mon Sep 17 00:00:00 2001 From: josh Date: Mon, 14 Nov 2022 20:36:11 +0000 Subject: [PATCH 7/7] avoid double reading thermal scattering libraries --- src/thermal.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/thermal.cpp b/src/thermal.cpp index b54f07b362..982bf226a9 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -119,20 +119,22 @@ ThermalScattering::ThermalScattering( if (!found) { // If no pairs found, check if the desired temperature falls within // bounds' tolerance - if (std::abs(T - temps_available[0]) <= - settings::temperature_tolerance) { + if (std::abs(T - temps_available[0]) <= settings::temperature_tolerance){ + if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[0])) == + temps_to_read.end()) { temps_to_read.push_back(std::round(temps_available[0])); - break; - } - if (std::abs(T - temps_available[n - 1]) <= - settings::temperature_tolerance) { + }} + else if (std::abs(T - temps_available[n - 1]) <= settings::temperature_tolerance){ + if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temps_available[n - 1])) == + temps_to_read.end()){ temps_to_read.push_back(std::round(temps_available[n - 1])); - break; - } - fatal_error( + }} + else { + fatal_error( fmt::format("Nuclear data library does not contain cross " "sections for {} at temperatures that bound {} K.", name_, std::round(T))); + } } } }