From ce54673f1d3afc4e725ed34e2040f3e47b9e2956 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 16 Mar 2019 22:01:41 -0400 Subject: [PATCH 01/61] include photon heating number in the library (from ACE) --- openmc/data/photon.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 47bdbf1d9..444136b31 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -42,6 +42,7 @@ _REACTION_NAME = { 516: 'Total pair production', 517: 'Pair production, nuclear field', 522: 'Photoelectric absorption', + 525: 'Heating', 526: 'Electro-atomic scattering', 527: 'Electro-atomic bremsstrahlung', 528: 'Electro-atomic excitation', @@ -446,9 +447,13 @@ class IncidentPhoton(EqualityMixin): # Read each reaction data = cls(Z) - for mt in (502, 504, 515, 522): + for mt in (502, 504, 515, 522, 525): data.reactions[mt] = PhotonReaction.from_ace(ace, mt) + # Get heating cross sections [eV*barn] from factors [MeV per collision] + data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in + (502, 504, 515, 522)]) * EV_PER_MEV + # Compton profiles n_shell = ace.nxs[5] if n_shell != 0: @@ -644,6 +649,11 @@ class IncidentPhoton(EqualityMixin): photoelec_group = group.create_group('photoelectric') photoelec_group.create_dataset('xs', data=self[522].xs(union_grid)) + # Write heating cross section + if 525 in self: + heat_group = group.create_group('heating') + heat_group.create_dataset('xs', data=self[525].xs(union_grid)) + # Write photoionization cross sections shell_group = group.create_group('subshells') designators = [] @@ -885,14 +895,18 @@ class PhotonReaction(EqualityMixin): elif mt == 522: # Photoelectric idx = ace.jxs[1] + 3*n + elif mt == 525: + # Heating + idx = ace.jxs[5] else: raise ValueError('ACE photoatomic cross sections do not have ' 'data for MT={}.'.format(mt)) # Store cross section xs = ace.xss[idx : idx+n].copy() - nonzero = (xs != 0.0) - xs[nonzero] = np.exp(xs[nonzero]) + if mt in (502, 504, 515, 522): + nonzero = (xs != 0.0) + xs[nonzero] = np.exp(xs[nonzero]) rx.xs = Tabulated1D(energy, xs, [n], [5]) # Get form factors for incoherent/coherent scattering From 0a18d85a1310e3eaed02db89f9140c4969e224d2 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 18 Mar 2019 13:56:12 -0400 Subject: [PATCH 02/61] read in photon heating number --- include/openmc/photon.h | 1 + src/photon.cpp | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 902d24232..0fbdbc0c9 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -67,6 +67,7 @@ public: xt::xtensor pair_production_total_; xt::xtensor pair_production_electron_; xt::xtensor pair_production_nuclear_; + xt::xtensor heating_; // Form factors Tabulated1D incoherent_form_factor_; diff --git a/src/photon.cpp b/src/photon.cpp index 79345b90e..86bf4390a 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -96,6 +96,11 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) read_dataset(rgroup, "xs", photoelectric_total_); close_group(rgroup); + // Read heating + rgroup = open_group(group, "heating"); + read_dataset(rgroup, "xs", heating_); + close_group(rgroup); + // Read subshell photoionization cross section and atomic relaxation data rgroup = open_group(group, "subshells"); std::vector designators; @@ -280,6 +285,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) xt::log(photoelectric_total_), -500.0); pair_production_total_ = xt::where(pair_production_total_ > 0.0, xt::log(pair_production_total_), -500.0); + heating_ = xt::where(heating_ > 0.0, xt::log(heating_), -500.0); } void PhotonInteraction::compton_scatter(double alpha, bool doppler, From 7ae1b024ed044c3e098cdba17d4e5f536e1dd779 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 18 Mar 2019 19:20:51 -0400 Subject: [PATCH 03/61] do heating tallies for both neutron and photon --- include/openmc/constants.h | 3 ++- openmc/capi/tally.py | 2 +- src/output.cpp | 1 + src/photon.cpp | 10 ++++++--- src/reaction.cpp | 2 ++ src/tallies/tally.cpp | 5 +++-- src/tallies/tally_scoring.cpp | 41 +++++++++++++++++++++++++++++++++++ 7 files changed, 57 insertions(+), 7 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index d9307abe2..dd9ecf211 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -230,7 +230,7 @@ constexpr int N_XD {204}; constexpr int N_XT {205}; constexpr int N_X3HE {206}; constexpr int N_XA {207}; -constexpr int HEATING {301}; +constexpr int NEUTRON_HEATING {301}; constexpr int DAMAGE_ENERGY {444}; constexpr int COHERENT {502}; constexpr int INCOHERENT {504}; @@ -366,6 +366,7 @@ constexpr int SCORE_INVERSE_VELOCITY {-13}; // flux-weighted inverse velocity constexpr int SCORE_FISS_Q_PROMPT {-14}; // prompt fission Q-value constexpr int SCORE_FISS_Q_RECOV {-15}; // recoverable fission Q-value constexpr int SCORE_DECAY_RATE {-16}; // delayed neutron precursor decay rate +constexpr int SCORE_HEATING {-17}; // nuclear heating (neutron or photon) // Tally map bin finding constexpr int NO_BIN_FOUND {-1}; diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index eec3a2b66..9529a31f2 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -89,7 +89,7 @@ _SCORES = { -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission', -9: 'current', -10: 'events', -11: 'delayed-nu-fission', -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt', - -15: 'fission-q-recoverable', -16: 'decay-rate' + -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'heating' } _ESTIMATORS = { 1: 'analog', 2: 'tracklength', 3: 'collision' diff --git a/src/output.cpp b/src/output.cpp index df0976d94..0c7ebcbda 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -618,6 +618,7 @@ const std::unordered_map score_names = { {SCORE_FISS_Q_PROMPT, "Prompt fission power"}, {SCORE_FISS_Q_RECOV, "Recoverable fission power"}, {SCORE_CURRENT, "Current"}, + {SCORE_HEATING, "Heating"}, }; //! Create an ASCII output file showing all tally results. diff --git a/src/photon.cpp b/src/photon.cpp index 86bf4390a..20000c4a9 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -97,9 +97,13 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) close_group(rgroup); // Read heating - rgroup = open_group(group, "heating"); - read_dataset(rgroup, "xs", heating_); - close_group(rgroup); + if (object_exists(group, "heating")) { + rgroup = open_group(group, "heating"); + read_dataset(rgroup, "xs", heating_); + close_group(rgroup); + } else { + heating_ = xt::zeros_like(energy_); + } // Read subshell photoionization cross section and atomic relaxation data rgroup = open_group(group, "subshells"); diff --git a/src/reaction.cpp b/src/reaction.cpp index d2405d899..222822230 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -121,6 +121,8 @@ std::string reaction_name(int mt) return "fission-q-prompt"; } else if (mt == SCORE_FISS_Q_RECOV) { return "fission-q-recoverable"; + } else if (mt == SCORE_HEATING) { + return "heating"; // Normal ENDF-based reactions } else if (mt == TOTAL_XS) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 682db3540..d26725db2 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -109,6 +109,9 @@ score_str_to_int(std::string score_str) if (score_str == "fission-q-recoverable") return SCORE_FISS_Q_RECOV; + if (score_str == "heating") + return SCORE_HEATING; + if (score_str == "current") return SCORE_CURRENT; @@ -207,8 +210,6 @@ score_str_to_int(std::string score_str) return N_X3HE; if (score_str == "(n,Xa)" || score_str == "He4-production") return N_XA; - if (score_str == "heating") - return HEATING; if (score_str == "damage-energy") return DAMAGE_ENERGY; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 7e235fa45..d574ed109 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -7,6 +7,7 @@ #include "openmc/material.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/photon.h" #include "openmc/reaction_product.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -1148,6 +1149,46 @@ score_general_ce(Particle* p, int i_tally, int start_index, break; + case SCORE_HEATING: + if (p->type_ == Particle::Type::neutron) { + score_bin = NEUTRON_HEATING; + // No break here, continue to run the default neutron MT scoring + } else if (p->type_ == Particle::Type::photon) { + if (tally.estimator_ != ESTIMATOR_ANALOG) { + // Calculate photon heating cross section on-the-fly + score = 0.; + if (i_nuclide >= 0) { + // Find the element corresponding to the nuclide + auto name = data::nuclides[i_nuclide]->name_; + int pos = name.find_first_of("0123456789"); + std::string element = name.substr(0, pos); + int i_element = data::element_map[element]; + auto& heating {data::elements[i_element].heating_}; + auto i_grid = simulation::micro_photon_xs[i_element].index_grid; + auto f = simulation::micro_photon_xs[i_element].interp_factor; + score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; + } else { + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_element = material.element_[i]; + auto atom_density = material.atom_density_(i); + auto& heating {data::elements[i_element].heating_}; + auto i_grid = simulation::micro_photon_xs[i_element].index_grid; + auto f = simulation::micro_photon_xs[i_element].interp_factor; + score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; + } + } + } + } + break; + } else { + // Nuclear heating of any other particles not implmented + break; + } + default: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need From 2af3af48dde89fc862c5a7e225397474e171c850 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 18 Mar 2019 21:57:46 -0400 Subject: [PATCH 04/61] heating tally specification --- docs/source/usersguide/tallies.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 2ced8acf0..6ab4af84c 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -256,9 +256,9 @@ The following tables show all valid scores: |inverse-velocity |The flux-weighted inverse velocity where the | | |velocity is in units of centimeters per second. | +----------------------+---------------------------------------------------+ - |heating |Total neutron heating in units of eV per source | - | |particle. This corresponds to MT=301 produced by | - | |NJOY's HEATR module. | + |heating |Total nuclear heating in units of eV per source | + | |particle. For neutrons, this corresponds to MT=301 | + | |produced by NJOY's HEATR module. | +----------------------+---------------------------------------------------+ |kappa-fission |The recoverable energy production rate due to | | |fission. The recoverable energy is defined as the | From 25d5a67fa17bee94cfe1adc59b5e666b7aa65a74 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 19 Mar 2019 13:41:59 -0400 Subject: [PATCH 05/61] read photon data from hdf5 --- openmc/data/photon.py | 144 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 143 insertions(+), 1 deletion(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 444136b31..9d453929c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -33,6 +33,7 @@ def _subshell(i): else: return _SUBSHELLS[i - 1] +_SUBSHELL_MT = {s: i + 534 for i, s in enumerate(_SUBSHELLS)} _REACTION_NAME = { 501: 'Total photon interaction', @@ -324,10 +325,31 @@ class AtomicRelaxation(EqualityMixin): # Return instance of class return cls(binding_energy, num_electrons, transitions) + @classmethod + def from_hdf5(cls, group_or_filename): + """Generate atomic relaxation data from HDF5 group + + Parameters + ---------- + group_or_filename : h5py.Group or str + HDF5 group containing interaction data. If given as a string, it is + assumed to be the filename for the HDF5 file, and the first group is + used to read from. + + Returns + ------- + openmc.data.AtomicRelaxation + Atomic relaxation data + + """ + + # Return instance of class + return None + #return cls(binding_energy, num_electrons, transitions) + def to_hdf5(self, group): raise NotImplementedError - class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -717,6 +739,126 @@ class IncidentPhoton(EqualityMixin): else: brem_group.create_dataset(key, data=value) + @classmethod + def from_hdf5(cls, group_or_filename): + """Generate photon reaction from an HDF5 group + + Parameters + ---------- + group_or_filename : h5py.Group or str + HDF5 group containing interaction data. If given as a string, it is + assumed to be the filename for the HDF5 file, and the first group is + used to read from. + + Returns + ------- + openmc.data.IncidentPhoton + Photon interaction data + + """ + if isinstance(group_or_filename, h5py.Group): + group = group_or_filename + else: + h5file = h5py.File(str(group_or_filename), 'r') + + # Make sure version matches + if 'version' in h5file.attrs: + major, minor = h5file.attrs['version'] + # For now all versions of HDF5 data can be read + else: + raise IOError( + 'HDF5 data does not indicate a version. Your installation of ' + 'the OpenMC Python API expects version {}.x data.' + .format(HDF5_VERSION_MAJOR)) + + group = list(h5file.values())[0] + + atomic_number = Z = group.attrs['Z'] + data = cls(Z) + + # Read energy grid + data.energy = group['energy'].value + n = data.energy.size + + # Read coherent scattering cross section + rgroup = group['coherent'] + rx = PhotonReaction(502) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + if 'anomalous_real' in rgroup: + rx.anomalous_real = Tabulated1D.from_hdf5(rgroup['anomalous_real']) + if 'anomalous_imag' in rgroup: + rx.anomalous_imag = Tabulated1D.from_hdf5(rgroup['anomalous_imag']) + data.reactions[502] = rx + + # Read incoherent scattering cross section + rgroup = group['incoherent'] + rx = PhotonReaction(504) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + if 'scattering_factor' in rgroup: + dset = rgroup['scattering_factor'] + rx.scattering_factor = Tabulated1D.from_hdf5(dset) + data.reactions[504] = rx + + # Read electron-field pair production cross section + if 'pair_production_electron' in group: + rgoup = group['pair_production_eletron'] + data.reactions[515] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + + # Read nuclear-field pair production cross section + if 'pair_production_nuclear' in group: + rgoup = group['pair_production_nuclear'] + data.reactions[517] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + + # Read photoelectric cross section + if 'photoelectric' in group: + rgoup = group['photoelectric'] + data.reactions[522] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + + # Read heating cross section + if 'heating' in group: + rgoup = group['heating'] + data.reactions[525] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + + # Read photoionization cross sections and atomic relaxation + rgroup = group['subshells'] + data.atomic_relaxation = AtomicRelaxation.from_hdf5(rgroup) + designators = rgroup.attrs['designators'] + n_shell = designators.size + for d in designators: + mt = _SUBSHELL_MT[d] + sub_group = rgroup[d] + xs = sub_group['xs'].value + threshold_idx = sub_group['xs'].attrs['threshold_idx'] + data.reactions[mt] = Tabulated1D(energy[threshold_idx:], xs, + [len(xs)], [5]) + + # Read Compton profiles + if 'compton_profiles' in group: + rgroup = group['compton_profiles'] + data.compton_profile['num_electrons'] = rgroup['num_electrons'].value + data.compton_profile['binding_energy'] = rgroup['binding_energy'].value + + # Get electron momentum values + pz = rgroup['pz'].value + J = rgroup['J'].value + if n_shell != J.shape[0]: + raise ValueError("'J' array shape is not consistent with the " + "number of shells") + if pz.size != J.shape[1]: + raise ValueError("'J' array shape is not consistent with the " + "'pz' array shape") + data.compton_profile['J'] = [Tabulated1D(pz, J[k]) for k in n_shell] + + # Read bremsstrahlung + if 'bremsstrahlung' in group: + rgroup = group['bremsstrahlung'] + data.bremsstrahlung['I'] = rgoup.attrs['I'] + for key in ('dcs', 'electron_energy', 'ionization_energy', + 'num_electrons', 'photon_energy'): + data.bremsstrahlung[key] = rgroup[key].value + + return data + def _add_bremsstrahlung(self): """Add the data used in the thick-target bremsstrahlung approximation From 4d80367ae429286f20afc6df93838212a2ffc8a9 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 19 Mar 2019 17:41:53 -0400 Subject: [PATCH 06/61] from_hdf5 completed --- openmc/data/photon.py | 100 ++++++++++++++++++++++-------------------- 1 file changed, 53 insertions(+), 47 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 9d453929c..ea0e0cb92 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -325,28 +325,6 @@ class AtomicRelaxation(EqualityMixin): # Return instance of class return cls(binding_energy, num_electrons, transitions) - @classmethod - def from_hdf5(cls, group_or_filename): - """Generate atomic relaxation data from HDF5 group - - Parameters - ---------- - group_or_filename : h5py.Group or str - HDF5 group containing interaction data. If given as a string, it is - assumed to be the filename for the HDF5 file, and the first group is - used to read from. - - Returns - ------- - openmc.data.AtomicRelaxation - Atomic relaxation data - - """ - - # Return instance of class - return None - #return cls(binding_energy, num_electrons, transitions) - def to_hdf5(self, group): raise NotImplementedError @@ -773,11 +751,11 @@ class IncidentPhoton(EqualityMixin): group = list(h5file.values())[0] - atomic_number = Z = group.attrs['Z'] + Z = group.attrs['Z'] data = cls(Z) # Read energy grid - data.energy = group['energy'].value + data.energy = energy= group['energy'].value n = data.energy.size # Read coherent scattering cross section @@ -801,58 +779,86 @@ class IncidentPhoton(EqualityMixin): # Read electron-field pair production cross section if 'pair_production_electron' in group: - rgoup = group['pair_production_eletron'] - data.reactions[515] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + rgroup = group['pair_production_electron'] + rx = PhotonReaction(515) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + data.reactions[515] = rx # Read nuclear-field pair production cross section if 'pair_production_nuclear' in group: - rgoup = group['pair_production_nuclear'] - data.reactions[517] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + rgroup = group['pair_production_nuclear'] + rx = PhotonReaction(517) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + data.reactions[517] = rx # Read photoelectric cross section if 'photoelectric' in group: - rgoup = group['photoelectric'] - data.reactions[522] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + rgroup = group['photoelectric'] + rx = PhotonReaction(522) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + data.reactions[522] = rx # Read heating cross section if 'heating' in group: - rgoup = group['heating'] - data.reactions[525] = Tabulated1D(energy, rgoup['xs'].value, [n], [5]) + rgroup = group['heating'] + rx = PhotonReaction(525) + rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) + data.reactions[525] = rx # Read photoionization cross sections and atomic relaxation rgroup = group['subshells'] - data.atomic_relaxation = AtomicRelaxation.from_hdf5(rgroup) - designators = rgroup.attrs['designators'] - n_shell = designators.size - for d in designators: - mt = _SUBSHELL_MT[d] - sub_group = rgroup[d] + designators = [i.decode() for i in rgroup.attrs['designators']] + binding_energy = {} + num_electrons = {} + transitions = {} + columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] + for shell in designators: + mt = _SUBSHELL_MT[shell] + sub_group = rgroup[shell] xs = sub_group['xs'].value threshold_idx = sub_group['xs'].attrs['threshold_idx'] - data.reactions[mt] = Tabulated1D(energy[threshold_idx:], xs, - [len(xs)], [5]) + rx = PhotonReaction(mt) + rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5]) + data.reactions[mt] = rx + + # Read subshell binding energy and number of electrons + binding_energy[shell] = sub_group.attrs['binding_energy'] + num_electrons[shell] = sub_group.attrs['num_electrons'] + + # Read transition data + if 'transitions' in sub_group: + t_value = sub_group['transitions'].value + records = [] + secondaries = [_subshell(int(i)) for i in t_value[:, 0]] + for i, s in enumerate(secondaries): + records.append((s, t_value[i, 1], t_value[i, 2], + t_value[i, 3])) + transitions[shell] = pd.DataFrame.from_records(records, + columns=columns) + + if binding_energy: + data.atomic_relaxation = AtomicRelaxation(binding_energy, + num_electrons, transitions) # Read Compton profiles if 'compton_profiles' in group: rgroup = group['compton_profiles'] - data.compton_profile['num_electrons'] = rgroup['num_electrons'].value - data.compton_profile['binding_energy'] = rgroup['binding_energy'].value + profile = data.compton_profiles + profile['num_electrons'] = rgroup['num_electrons'].value + profile['binding_energy'] = rgroup['binding_energy'].value # Get electron momentum values pz = rgroup['pz'].value J = rgroup['J'].value - if n_shell != J.shape[0]: - raise ValueError("'J' array shape is not consistent with the " - "number of shells") if pz.size != J.shape[1]: raise ValueError("'J' array shape is not consistent with the " "'pz' array shape") - data.compton_profile['J'] = [Tabulated1D(pz, J[k]) for k in n_shell] + profile['J'] = [Tabulated1D(pz, Jk) for Jk in J] # Read bremsstrahlung if 'bremsstrahlung' in group: rgroup = group['bremsstrahlung'] - data.bremsstrahlung['I'] = rgoup.attrs['I'] + data.bremsstrahlung['I'] = rgroup.attrs['I'] for key in ('dcs', 'electron_energy', 'ionization_energy', 'num_electrons', 'photon_energy'): data.bremsstrahlung[key] = rgroup[key].value From bb30f78b0d0646e50c252194a5b4cf1058cc5753 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 09:44:40 -0400 Subject: [PATCH 07/61] using pd.DataFrame to replace the strings --- openmc/data/photon.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index ea0e0cb92..9c046c525 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -595,6 +595,9 @@ class IncidentPhoton(EqualityMixin): mode : {'r', r+', 'w', 'x', 'a'} Mode that is used to open the HDF5 file. This is the second argument to the :class:`h5py.File` constructor. + libver : {'earliest', 'latest'} + Compatibility mode for the HDF5 file. 'latest' will produce files + that are less backwards compatible but have performance benefits. """ # Open file and write version @@ -755,8 +758,8 @@ class IncidentPhoton(EqualityMixin): data = cls(Z) # Read energy grid - data.energy = energy= group['energy'].value - n = data.energy.size + energy= group['energy'].value + n = energy.size # Read coherent scattering cross section rgroup = group['coherent'] @@ -811,6 +814,8 @@ class IncidentPhoton(EqualityMixin): binding_energy = {} num_electrons = {} transitions = {} + shell_values = _SUBSHELLS.copy() + shell_values.insert(0, None) columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] for shell in designators: mt = _SUBSHELL_MT[shell] @@ -827,18 +832,16 @@ class IncidentPhoton(EqualityMixin): # Read transition data if 'transitions' in sub_group: - t_value = sub_group['transitions'].value - records = [] - secondaries = [_subshell(int(i)) for i in t_value[:, 0]] - for i, s in enumerate(secondaries): - records.append((s, t_value[i, 1], t_value[i, 2], - t_value[i, 3])) - transitions[shell] = pd.DataFrame.from_records(records, - columns=columns) + df = pd.DataFrame(sub_group['transitions'].value, + columns=columns) + # Replace float indexes back to subshell strings + df[columns[:2]] = df[columns[:2]].replace( + np.arange(float(len(shell_values))), shell_values) + transitions[shell] = df if binding_energy: data.atomic_relaxation = AtomicRelaxation(binding_energy, - num_electrons, transitions) + num_electrons, transitions) # Read Compton profiles if 'compton_profiles' in group: From 69e1a9d30d1e7adc4a4b5ea33fe20f8ffcb45701 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 10:42:06 -0400 Subject: [PATCH 08/61] Add from_hdf5 for atomic relaxation data --- openmc/data/photon.py | 359 +++++++++++++++++++++++------------------- 1 file changed, 200 insertions(+), 159 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 9c046c525..6c66efc01 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -325,8 +325,69 @@ class AtomicRelaxation(EqualityMixin): # Return instance of class return cls(binding_energy, num_electrons, transitions) + @classmethod + def from_hdf5(cls, group): + """Generate atomic relaxation data from an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.AtomicRelaxation + Atomic relaxation data + + """ + # Create data dictionaries + binding_energy = {} + num_electrons = {} + transitions = {} + + designators = [s.decode() for s in group.attrs['designators']] + shell_values = [None] + _SUBSHELLS + columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] + for shell in designators: + # Shell group + sub_group = group[shell] + + # Read subshell binding energy and number of electrons + if 'binding_energy' in sub_group.attrs: + binding_energy[shell] = sub_group.attrs['binding_energy'] + if 'num_electrons' in sub_group.attrs: + num_electrons[shell] = sub_group.attrs['num_electrons'] + + # Read transition data + if 'transitions' in sub_group: + df = pd.DataFrame(sub_group['transitions'].value, + columns=columns) + # Replace float indexes back to subshell strings + df[columns[:2]] = df[columns[:2]].replace( + np.arange(float(len(shell_values))), shell_values) + transitions[shell] = df + + return cls(binding_energy, num_electrons, transitions) + def to_hdf5(self, group): - raise NotImplementedError + """Write atomic relaxation data to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + + group.attrs['mt'] = self.mt + if self.mt in REACTION_NAME: + group.attrs['label'] = np.string_(REACTION_NAME[self.mt]) + else: + group.attrs['label'] = np.string_(self.mt) + group.attrs['Q_value'] = self.q_value + group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 + group.attrs['redundant'] = 1 if self.redundant else 0 + class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -585,141 +646,6 @@ class IncidentPhoton(EqualityMixin): return data - def export_to_hdf5(self, path, mode='a', libver='earliest'): - """Export incident photon data to an HDF5 file. - - Parameters - ---------- - path : str - Path to write HDF5 file to - mode : {'r', r+', 'w', 'x', 'a'} - Mode that is used to open the HDF5 file. This is the second argument - to the :class:`h5py.File` constructor. - libver : {'earliest', 'latest'} - Compatibility mode for the HDF5 file. 'latest' will produce files - that are less backwards compatible but have performance benefits. - - """ - # Open file and write version - f = h5py.File(str(path), mode, libver=libver) - f.attrs['filetype'] = np.string_('data_photon') - if 'version' not in f.attrs: - f.attrs['version'] = np.array(HDF5_VERSION) - - group = f.create_group(self.name) - group.attrs['Z'] = Z = self.atomic_number - - # Determine union energy grid - union_grid = np.array([]) - for rx in self: - union_grid = np.union1d(union_grid, rx.xs.x) - group.create_dataset('energy', data=union_grid) - - # Write coherent scattering cross section - rx = self.reactions[502] - coh_group = group.create_group('coherent') - coh_group.create_dataset('xs', data=rx.xs(union_grid)) - if rx.scattering_factor is not None: - # Create integrated form factor - ff = deepcopy(rx.scattering_factor) - ff.x *= ff.x - ff.y *= ff.y/Z**2 - int_ff = Tabulated1D(ff.x, ff.integral()) - int_ff.to_hdf5(coh_group, 'integrated_scattering_factor') - if rx.anomalous_real is not None: - rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real') - if rx.anomalous_imag is not None: - rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag') - - # Write incoherent scattering cross section - rx = self[504] - incoh_group = group.create_group('incoherent') - incoh_group.create_dataset('xs', data=rx.xs(union_grid)) - if rx.scattering_factor is not None: - rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor') - - # Write electron-field pair production cross section - if 515 in self: - pair_group = group.create_group('pair_production_electron') - pair_group.create_dataset('xs', data=self[515].xs(union_grid)) - - # Write nuclear-field pair production cross section - if 517 in self: - pair_group = group.create_group('pair_production_nuclear') - pair_group.create_dataset('xs', data=self[517].xs(union_grid)) - - # Write photoelectric cross section - photoelec_group = group.create_group('photoelectric') - photoelec_group.create_dataset('xs', data=self[522].xs(union_grid)) - - # Write heating cross section - if 525 in self: - heat_group = group.create_group('heating') - heat_group.create_dataset('xs', data=self[525].xs(union_grid)) - - # Write photoionization cross sections - shell_group = group.create_group('subshells') - designators = [] - for mt, rx in self.reactions.items(): - if mt >= 534 and mt <= 572: - # Get name of subshell - shell = _SUBSHELLS[mt - 534] - designators.append(shell) - sub_group = shell_group.create_group(shell) - - if self.atomic_relaxation is not None: - relax = self.atomic_relaxation - # Write subshell binding energy and number of electrons - sub_group.attrs['binding_energy'] = relax.binding_energy[shell] - sub_group.attrs['num_electrons'] = relax.num_electrons[shell] - - # Write transition data with replacements - if shell in relax.transitions: - shell_values = _SUBSHELLS.copy() - shell_values.insert(0, None) - df = relax.transitions[shell].replace( - shell_values, range(len(shell_values))) - sub_group.create_dataset( - 'transitions', data=df.values.astype(float)) - - # Determine threshold - threshold = rx.xs.x[0] - idx = np.searchsorted(union_grid, threshold, side='right') - 1 - - # Interpolate cross section onto union grid and write - photoionization = rx.xs(union_grid[idx:]) - sub_group.create_dataset('xs', data=photoionization) - assert len(union_grid) == len(photoionization) + idx - sub_group['xs'].attrs['threshold_idx'] = idx - - shell_group.attrs['designators'] = np.array(designators, dtype='S') - - # Write Compton profiles - if self.compton_profiles: - compton_group = group.create_group('compton_profiles') - - profile = self.compton_profiles - compton_group.create_dataset('num_electrons', - data=profile['num_electrons']) - compton_group.create_dataset('binding_energy', - data=profile['binding_energy']) - - # Get electron momentum values - compton_group.create_dataset('pz', data=profile['J'][0].x) - - # Create/write 2D array of profiles - J = np.array([Jk.y for Jk in profile['J']]) - compton_group.create_dataset('J', data=J) - - # Write bremsstrahlung - if self.bremsstrahlung: - brem_group = group.create_group('bremsstrahlung') - for key, value in self.bremsstrahlung.items(): - if key == 'I': - brem_group.attrs[key] = value - else: - brem_group.create_dataset(key, data=value) - @classmethod def from_hdf5(cls, group_or_filename): """Generate photon reaction from an HDF5 group @@ -808,15 +734,9 @@ class IncidentPhoton(EqualityMixin): rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) data.reactions[525] = rx - # Read photoionization cross sections and atomic relaxation rgroup = group['subshells'] + # Read photoionization cross sections designators = [i.decode() for i in rgroup.attrs['designators']] - binding_energy = {} - num_electrons = {} - transitions = {} - shell_values = _SUBSHELLS.copy() - shell_values.insert(0, None) - columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] for shell in designators: mt = _SUBSHELL_MT[shell] sub_group = rgroup[shell] @@ -826,22 +746,8 @@ class IncidentPhoton(EqualityMixin): rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5]) data.reactions[mt] = rx - # Read subshell binding energy and number of electrons - binding_energy[shell] = sub_group.attrs['binding_energy'] - num_electrons[shell] = sub_group.attrs['num_electrons'] - - # Read transition data - if 'transitions' in sub_group: - df = pd.DataFrame(sub_group['transitions'].value, - columns=columns) - # Replace float indexes back to subshell strings - df[columns[:2]] = df[columns[:2]].replace( - np.arange(float(len(shell_values))), shell_values) - transitions[shell] = df - - if binding_energy: - data.atomic_relaxation = AtomicRelaxation(binding_energy, - num_electrons, transitions) + # Read atomic relaxation + data.atomic_relaxation = AtomicRelaxation.from_hdf5(rgroup) # Read Compton profiles if 'compton_profiles' in group: @@ -868,6 +774,141 @@ class IncidentPhoton(EqualityMixin): return data + def export_to_hdf5(self, path, mode='a', libver='earliest'): + """Export incident photon data to an HDF5 file. + + Parameters + ---------- + path : str + Path to write HDF5 file to + mode : {'r', r+', 'w', 'x', 'a'} + Mode that is used to open the HDF5 file. This is the second argument + to the :class:`h5py.File` constructor. + libver : {'earliest', 'latest'} + Compatibility mode for the HDF5 file. 'latest' will produce files + that are less backwards compatible but have performance benefits. + + """ + # Open file and write version + f = h5py.File(str(path), mode, libver=libver) + f.attrs['filetype'] = np.string_('data_photon') + if 'version' not in f.attrs: + f.attrs['version'] = np.array(HDF5_VERSION) + + group = f.create_group(self.name) + group.attrs['Z'] = Z = self.atomic_number + + # Determine union energy grid + union_grid = np.array([]) + for rx in self: + union_grid = np.union1d(union_grid, rx.xs.x) + group.create_dataset('energy', data=union_grid) + + # Write coherent scattering cross section + rx = self.reactions[502] + coh_group = group.create_group('coherent') + coh_group.create_dataset('xs', data=rx.xs(union_grid)) + if rx.scattering_factor is not None: + # Create integrated form factor + ff = deepcopy(rx.scattering_factor) + ff.x *= ff.x + ff.y *= ff.y/Z**2 + int_ff = Tabulated1D(ff.x, ff.integral()) + int_ff.to_hdf5(coh_group, 'integrated_scattering_factor') + if rx.anomalous_real is not None: + rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real') + if rx.anomalous_imag is not None: + rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag') + + # Write incoherent scattering cross section + rx = self[504] + incoh_group = group.create_group('incoherent') + incoh_group.create_dataset('xs', data=rx.xs(union_grid)) + if rx.scattering_factor is not None: + rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor') + + # Write electron-field pair production cross section + if 515 in self: + pair_group = group.create_group('pair_production_electron') + pair_group.create_dataset('xs', data=self[515].xs(union_grid)) + + # Write nuclear-field pair production cross section + if 517 in self: + pair_group = group.create_group('pair_production_nuclear') + pair_group.create_dataset('xs', data=self[517].xs(union_grid)) + + # Write photoelectric cross section + photoelec_group = group.create_group('photoelectric') + photoelec_group.create_dataset('xs', data=self[522].xs(union_grid)) + + # Write heating cross section + if 525 in self: + heat_group = group.create_group('heating') + heat_group.create_dataset('xs', data=self[525].xs(union_grid)) + + # Write photoionization cross sections + shell_group = group.create_group('subshells') + designators = [] + shell_values = [None] + _SUBSHELLS + for mt, rx in self.reactions.items(): + if mt >= 534 and mt <= 572: + # Get name of subshell + shell = _SUBSHELLS[mt - 534] + designators.append(shell) + sub_group = shell_group.create_group(shell) + + if self.atomic_relaxation is not None: + relax = self.atomic_relaxation + # Write subshell binding energy and number of electrons + sub_group.attrs['binding_energy'] = relax.binding_energy[shell] + sub_group.attrs['num_electrons'] = relax.num_electrons[shell] + + # Write transition data with replacements + if shell in relax.transitions: + df = relax.transitions[shell].replace( + shell_values, range(len(shell_values))) + sub_group.create_dataset( + 'transitions', data=df.values.astype(float)) + + # Determine threshold + threshold = rx.xs.x[0] + idx = np.searchsorted(union_grid, threshold, side='right') - 1 + + # Interpolate cross section onto union grid and write + photoionization = rx.xs(union_grid[idx:]) + sub_group.create_dataset('xs', data=photoionization) + assert len(union_grid) == len(photoionization) + idx + sub_group['xs'].attrs['threshold_idx'] = idx + + shell_group.attrs['designators'] = np.array(designators, dtype='S') + + # Write Compton profiles + if self.compton_profiles: + compton_group = group.create_group('compton_profiles') + + profile = self.compton_profiles + compton_group.create_dataset('num_electrons', + data=profile['num_electrons']) + compton_group.create_dataset('binding_energy', + data=profile['binding_energy']) + + # Get electron momentum values + compton_group.create_dataset('pz', data=profile['J'][0].x) + + # Create/write 2D array of profiles + J = np.array([Jk.y for Jk in profile['J']]) + compton_group.create_dataset('J', data=J) + + # Write bremsstrahlung + if self.bremsstrahlung: + brem_group = group.create_group('bremsstrahlung') + for key, value in self.bremsstrahlung.items(): + if key == 'I': + brem_group.attrs[key] = value + else: + brem_group.create_dataset(key, data=value) + + def _add_bremsstrahlung(self): """Add the data used in the thick-target bremsstrahlung approximation From 06610823376f295197f56fbb7d0d42faafcbb040 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 11:01:06 -0400 Subject: [PATCH 09/61] added to_hdf5 for atomic relaxation --- openmc/data/photon.py | 38 ++++++++++++++++---------------------- 1 file changed, 16 insertions(+), 22 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 6c66efc01..2722d0324 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -369,24 +369,28 @@ class AtomicRelaxation(EqualityMixin): return cls(binding_energy, num_electrons, transitions) - def to_hdf5(self, group): + def to_hdf5(self, group, shell): """Write atomic relaxation data to an HDF5 group Parameters ---------- group : h5py.Group HDF5 group to write to + shell : str + The subshell to write data for """ - group.attrs['mt'] = self.mt - if self.mt in REACTION_NAME: - group.attrs['label'] = np.string_(REACTION_NAME[self.mt]) - else: - group.attrs['label'] = np.string_(self.mt) - group.attrs['Q_value'] = self.q_value - group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 - group.attrs['redundant'] = 1 if self.redundant else 0 + # Write subshell binding energy and number of electrons + group.attrs['binding_energy'] = self.binding_energy[shell] + group.attrs['num_electrons'] = self.num_electrons[shell] + + # Write transition data with replacements + if shell in self.transitions: + shell_values = [None] + _SUBSHELLS + df = self.transitions[shell].replace( + shell_values, range(len(shell_values))) + group.create_dataset('transitions', data=df.values.astype(float)) class IncidentPhoton(EqualityMixin): @@ -849,7 +853,6 @@ class IncidentPhoton(EqualityMixin): # Write photoionization cross sections shell_group = group.create_group('subshells') designators = [] - shell_values = [None] + _SUBSHELLS for mt, rx in self.reactions.items(): if mt >= 534 and mt <= 572: # Get name of subshell @@ -857,18 +860,9 @@ class IncidentPhoton(EqualityMixin): designators.append(shell) sub_group = shell_group.create_group(shell) - if self.atomic_relaxation is not None: - relax = self.atomic_relaxation - # Write subshell binding energy and number of electrons - sub_group.attrs['binding_energy'] = relax.binding_energy[shell] - sub_group.attrs['num_electrons'] = relax.num_electrons[shell] - - # Write transition data with replacements - if shell in relax.transitions: - df = relax.transitions[shell].replace( - shell_values, range(len(shell_values))) - sub_group.create_dataset( - 'transitions', data=df.values.astype(float)) + # Write atomic relaxation + if shell in self.atomic_relaxation.subshells: + self.atomic_relaxation.to_hdf5(sub_group, shell) # Determine threshold threshold = rx.xs.x[0] From 36641b10291aac2fdea83e0d385809cd2983176f Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 13:23:44 -0400 Subject: [PATCH 10/61] create from_hdf5 in photon reactions to simplify photon export --- openmc/data/photon.py | 222 +++++++++++++++++++++--------------------- 1 file changed, 112 insertions(+), 110 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 2722d0324..663c6a2b7 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -36,56 +36,56 @@ def _subshell(i): _SUBSHELL_MT = {s: i + 534 for i, s in enumerate(_SUBSHELLS)} _REACTION_NAME = { - 501: 'Total photon interaction', - 502: 'Photon coherent scattering', - 504: 'Photon incoherent scattering', - 515: 'Pair production, electron field', - 516: 'Total pair production', - 517: 'Pair production, nuclear field', - 522: 'Photoelectric absorption', - 525: 'Heating', - 526: 'Electro-atomic scattering', - 527: 'Electro-atomic bremsstrahlung', - 528: 'Electro-atomic excitation', - 534: 'K (1s1/2) subshell photoelectric', - 535: 'L1 (2s1/2) subshell photoelectric', - 536: 'L2 (2p1/2) subshell photoelectric', - 537: 'L3 (2p3/2) subshell photoelectric', - 538: 'M1 (3s1/2) subshell photoelectric', - 539: 'M2 (3p1/2) subshell photoelectric', - 540: 'M3 (3p3/2) subshell photoelectric', - 541: 'M4 (3d3/2) subshell photoelectric', - 542: 'M5 (3d5/2) subshell photoelectric', - 543: 'N1 (4s1/2) subshell photoelectric', - 544: 'N2 (4p1/2) subshell photoelectric', - 545: 'N3 (4p3/2) subshell photoelectric', - 546: 'N4 (4d3/2) subshell photoelectric', - 547: 'N5 (4d5/2) subshell photoelectric', - 548: 'N6 (4f5/2) subshell photoelectric', - 549: 'N7 (4f7/2) subshell photoelectric', - 550: 'O1 (5s1/2) subshell photoelectric', - 551: 'O2 (5p1/2) subshell photoelectric', - 552: 'O3 (5p3/2) subshell photoelectric', - 553: 'O4 (5d3/2) subshell photoelectric', - 554: 'O5 (5d5/2) subshell photoelectric', - 555: 'O6 (5f5/2) subshell photoelectric', - 556: 'O7 (5f7/2) subshell photoelectric', - 557: 'O8 (5g7/2) subshell photoelectric', - 558: 'O9 (5g9/2) subshell photoelectric', - 559: 'P1 (6s1/2) subshell photoelectric', - 560: 'P2 (6p1/2) subshell photoelectric', - 561: 'P3 (6p3/2) subshell photoelectric', - 562: 'P4 (6d3/2) subshell photoelectric', - 563: 'P5 (6d5/2) subshell photoelectric', - 564: 'P6 (6f5/2) subshell photoelectric', - 565: 'P7 (6f7/2) subshell photoelectric', - 566: 'P8 (6g7/2) subshell photoelectric', - 567: 'P9 (6g9/2) subshell photoelectric', - 568: 'P10 (6h9/2) subshell photoelectric', - 569: 'P11 (6h11/2) subshell photoelectric', - 570: 'Q1 (7s1/2) subshell photoelectric', - 571: 'Q2 (7p1/2) subshell photoelectric', - 572: 'Q3 (7p3/2) subshell photoelectric' + 501: ('Total photon interaction', 'total'), + 502: ('Photon coherent scattering', 'coherent'), + 504: ('Photon incoherent scattering', 'incoherent'), + 515: ('Pair production, electron field', 'pair_production_electron'), + 516: ('Total pair production', 'pair_production_total'), + 517: ('Pair production, nuclear field', 'pair_production_nuclear'), + 522: ('Photoelectric absorption', 'photoelectric'), + 525: ('Heating', 'heating'), + 526: ('Electro-atomic scattering', 'electro_atomic_scat'), + 527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'), + 528: ('Electro-atomic excitation', 'electro_atomic_excit'), + 534: ('K (1s1/2) subshell photoelectric', 'K'), + 535: ('L1 (2s1/2) subshell photoelectric', 'L1'), + 536: ('L2 (2p1/2) subshell photoelectric', 'L2'), + 537: ('L3 (2p3/2) subshell photoelectric', 'L3'), + 538: ('M1 (3s1/2) subshell photoelectric', 'M1'), + 539: ('M2 (3p1/2) subshell photoelectric', 'M2'), + 540: ('M3 (3p3/2) subshell photoelectric', 'M3'), + 541: ('M4 (3d3/2) subshell photoelectric', 'M4'), + 542: ('M5 (3d5/2) subshell photoelectric', 'M5'), + 543: ('N1 (4s1/2) subshell photoelectric', 'N1'), + 544: ('N2 (4p1/2) subshell photoelectric', 'N2'), + 545: ('N3 (4p3/2) subshell photoelectric', 'N3'), + 546: ('N4 (4d3/2) subshell photoelectric', 'N4'), + 547: ('N5 (4d5/2) subshell photoelectric', 'N5'), + 548: ('N6 (4f5/2) subshell photoelectric', 'N6'), + 549: ('N7 (4f7/2) subshell photoelectric', 'N7'), + 550: ('O1 (5s1/2) subshell photoelectric', 'O1'), + 551: ('O2 (5p1/2) subshell photoelectric', 'O2'), + 552: ('O3 (5p3/2) subshell photoelectric', 'O3'), + 553: ('O4 (5d3/2) subshell photoelectric', 'O4'), + 554: ('O5 (5d5/2) subshell photoelectric', 'O5'), + 555: ('O6 (5f5/2) subshell photoelectric', 'O6'), + 556: ('O7 (5f7/2) subshell photoelectric', 'O7'), + 557: ('O8 (5g7/2) subshell photoelectric', 'O8'), + 558: ('O9 (5g9/2) subshell photoelectric', 'O9'), + 559: ('P1 (6s1/2) subshell photoelectric', 'P1'), + 560: ('P2 (6p1/2) subshell photoelectric', 'P2'), + 561: ('P3 (6p3/2) subshell photoelectric', 'P3'), + 562: ('P4 (6d3/2) subshell photoelectric', 'P4'), + 563: ('P5 (6d5/2) subshell photoelectric', 'P5'), + 564: ('P6 (6f5/2) subshell photoelectric', 'P6'), + 565: ('P7 (6f7/2) subshell photoelectric', 'P7'), + 566: ('P8 (6g7/2) subshell photoelectric', 'P8'), + 567: ('P9 (6g9/2) subshell photoelectric', 'P9'), + 568: ('P10 (6h9/2) subshell photoelectric', 'P10'), + 569: ('P11 (6h11/2) subshell photoelectric', 'P11'), + 570: ('Q1 (7s1/2) subshell photoelectric', 'Q1'), + 571: ('Q2 (7p1/2) subshell photoelectric', 'Q2'), + 572: ('Q3 (7p3/2) subshell photoelectric', 'Q3') } # Compton profiles are read from a pre-generated HDF5 file when they are first @@ -678,8 +678,8 @@ class IncidentPhoton(EqualityMixin): # For now all versions of HDF5 data can be read else: raise IOError( - 'HDF5 data does not indicate a version. Your installation of ' - 'the OpenMC Python API expects version {}.x data.' + 'HDF5 data does not indicate a version. Your installation ' + 'of the OpenMC Python API expects version {}.x data.' .format(HDF5_VERSION_MAJOR)) group = list(h5file.values())[0] @@ -689,69 +689,24 @@ class IncidentPhoton(EqualityMixin): # Read energy grid energy= group['energy'].value - n = energy.size - # Read coherent scattering cross section - rgroup = group['coherent'] - rx = PhotonReaction(502) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - if 'anomalous_real' in rgroup: - rx.anomalous_real = Tabulated1D.from_hdf5(rgroup['anomalous_real']) - if 'anomalous_imag' in rgroup: - rx.anomalous_imag = Tabulated1D.from_hdf5(rgroup['anomalous_imag']) - data.reactions[502] = rx + # Read cross section data + for mt, (name, key) in _REACTION_NAME.items(): + if key in group: + rgroup = group[key] + elif key in group['subshells']: + rgroup = group['subshells'][key] + else: + continue - # Read incoherent scattering cross section - rgroup = group['incoherent'] - rx = PhotonReaction(504) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - if 'scattering_factor' in rgroup: - dset = rgroup['scattering_factor'] - rx.scattering_factor = Tabulated1D.from_hdf5(dset) - data.reactions[504] = rx + data.reactions[mt] = PhotonReaction.from_hdf5(rgroup, mt, energy) - # Read electron-field pair production cross section - if 'pair_production_electron' in group: - rgroup = group['pair_production_electron'] - rx = PhotonReaction(515) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - data.reactions[515] = rx - - # Read nuclear-field pair production cross section - if 'pair_production_nuclear' in group: - rgroup = group['pair_production_nuclear'] - rx = PhotonReaction(517) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - data.reactions[517] = rx - - # Read photoelectric cross section - if 'photoelectric' in group: - rgroup = group['photoelectric'] - rx = PhotonReaction(522) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - data.reactions[522] = rx - - # Read heating cross section - if 'heating' in group: - rgroup = group['heating'] - rx = PhotonReaction(525) - rx.xs = Tabulated1D(energy, rgroup['xs'].value, [n], [5]) - data.reactions[525] = rx - - rgroup = group['subshells'] - # Read photoionization cross sections - designators = [i.decode() for i in rgroup.attrs['designators']] - for shell in designators: - mt = _SUBSHELL_MT[shell] - sub_group = rgroup[shell] - xs = sub_group['xs'].value - threshold_idx = sub_group['xs'].attrs['threshold_idx'] - rx = PhotonReaction(mt) - rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5]) - data.reactions[mt] = rx + # Check for necessary reactions + for mt in [502, 504, 522]: + assert mt in data, "reaction {} not found".format(mt) # Read atomic relaxation - data.atomic_relaxation = AtomicRelaxation.from_hdf5(rgroup) + data.atomic_relaxation = AtomicRelaxation.from_hdf5(group['subshells']) # Read Compton profiles if 'compton_profiles' in group: @@ -1182,3 +1137,50 @@ class PhotonReaction(EqualityMixin): params, rx.anomalous_imag = get_tab1_record(file_obj) return rx + + @classmethod + def from_hdf5(cls, group, mt, energy): + """Generate photon reaction from an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + mt : int + The MT value of the reaction to get data for + energy : Iterable of float + arrays of energies at which cross sections are tabulated at. + + Returns + ------- + openmc.data.PhotonReaction + Photon reaction data + + """ + # Create instance + rx = cls(mt) + + # Cross sections + xs = group['xs'].value + # Replace zero elements to small non-zero to enable log-log + xs[xs == 0.0] = np.exp(-500.0) + + # Threshold + threshold_idx = 0 + if 'threshold_idx' in group['xs'].attrs: + threshold_idx = group['xs'].attrs['threshold_idx'] + + # Store + rx.xs = Tabulated1D(energy[threshold_idx:], xs, [len(xs)], [5]) + + # Check for anomalous scattering factor + if 'anomalous_real' in group: + rx.anomalous_real = Tabulated1D.from_hdf5(group['anomalous_real']) + if 'anomalous_imag' in group: + rx.anomalous_imag = Tabulated1D.from_hdf5(group['anomalous_imag']) + + # Check for factors / scattering functions + if 'scattering_factor' in group: + rx.scattering_factor = Tabulated1D.from_hdf5(group['scattering_factor']) + + return rx From 4eee743ec50c109cac8f8323e2d4191cb4cfc985 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 14:18:57 -0400 Subject: [PATCH 11/61] added to_hdf5 in photon reaction to simplify export --- openmc/data/photon.py | 120 ++++++++++++++++++++---------------------- 1 file changed, 58 insertions(+), 62 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 663c6a2b7..73c48527d 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -33,8 +33,6 @@ def _subshell(i): else: return _SUBSHELLS[i - 1] -_SUBSHELL_MT = {s: i + 534 for i, s in enumerate(_SUBSHELLS)} - _REACTION_NAME = { 501: ('Total photon interaction', 'total'), 502: ('Photon coherent scattering', 'coherent'), @@ -763,71 +761,25 @@ class IncidentPhoton(EqualityMixin): union_grid = np.union1d(union_grid, rx.xs.x) group.create_dataset('energy', data=union_grid) - # Write coherent scattering cross section - rx = self.reactions[502] - coh_group = group.create_group('coherent') - coh_group.create_dataset('xs', data=rx.xs(union_grid)) - if rx.scattering_factor is not None: - # Create integrated form factor - ff = deepcopy(rx.scattering_factor) - ff.x *= ff.x - ff.y *= ff.y/Z**2 - int_ff = Tabulated1D(ff.x, ff.integral()) - int_ff.to_hdf5(coh_group, 'integrated_scattering_factor') - if rx.anomalous_real is not None: - rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real') - if rx.anomalous_imag is not None: - rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag') - - # Write incoherent scattering cross section - rx = self[504] - incoh_group = group.create_group('incoherent') - incoh_group.create_dataset('xs', data=rx.xs(union_grid)) - if rx.scattering_factor is not None: - rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor') - - # Write electron-field pair production cross section - if 515 in self: - pair_group = group.create_group('pair_production_electron') - pair_group.create_dataset('xs', data=self[515].xs(union_grid)) - - # Write nuclear-field pair production cross section - if 517 in self: - pair_group = group.create_group('pair_production_nuclear') - pair_group.create_dataset('xs', data=self[517].xs(union_grid)) - - # Write photoelectric cross section - photoelec_group = group.create_group('photoelectric') - photoelec_group.create_dataset('xs', data=self[522].xs(union_grid)) - - # Write heating cross section - if 525 in self: - heat_group = group.create_group('heating') - heat_group.create_dataset('xs', data=self[525].xs(union_grid)) - - # Write photoionization cross sections + # Write cross sections shell_group = group.create_group('subshells') designators = [] for mt, rx in self.reactions.items(): - if mt >= 534 and mt <= 572: - # Get name of subshell - shell = _SUBSHELLS[mt - 534] - designators.append(shell) - sub_group = shell_group.create_group(shell) + name, key = _REACTION_NAME[mt] + if mt in [502, 504, 515, 517, 522, 525]: + sub_group = group.create_group(key) + elif mt >= 534 and mt <= 572: + # Subshell + designators.append(key) + sub_group = shell_group.create_group(key) # Write atomic relaxation - if shell in self.atomic_relaxation.subshells: - self.atomic_relaxation.to_hdf5(sub_group, shell) + if key in self.atomic_relaxation.subshells: + self.atomic_relaxation.to_hdf5(sub_group, key) + else: + continue - # Determine threshold - threshold = rx.xs.x[0] - idx = np.searchsorted(union_grid, threshold, side='right') - 1 - - # Interpolate cross section onto union grid and write - photoionization = rx.xs(union_grid[idx:]) - sub_group.create_dataset('xs', data=photoionization) - assert len(union_grid) == len(photoionization) + idx - sub_group['xs'].attrs['threshold_idx'] = idx + rx.to_hdf5(sub_group, union_grid, Z) shell_group.attrs['designators'] = np.array(designators, dtype='S') @@ -1149,7 +1101,7 @@ class PhotonReaction(EqualityMixin): mt : int The MT value of the reaction to get data for energy : Iterable of float - arrays of energies at which cross sections are tabulated at. + arrays of energies at which cross sections are tabulated at Returns ------- @@ -1184,3 +1136,47 @@ class PhotonReaction(EqualityMixin): rx.scattering_factor = Tabulated1D.from_hdf5(group['scattering_factor']) return rx + + def to_hdf5(self, group, energy, Z): + """Write photon reaction to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + energy : Iterable of float + arrays of energies at which cross sections are tabulated at + Z : int + atomic number + + """ + + # Write cross sections + if self.mt >= 534 and self.mt <= 572: + # Determine threshold + threshold = self.xs.x[0] + idx = np.searchsorted(energy, threshold, side='right') - 1 + + # Interpolate cross section onto union grid and write + photoionization = self.xs(energy[idx:]) + group.create_dataset('xs', data=photoionization) + assert len(energy) == len(photoionization) + idx + group['xs'].attrs['threshold_idx'] = idx + else: + group.create_dataset('xs', data=self.xs(energy)) + + # Write scattering factor + if self.scattering_factor is not None: + if self.mt == 502: + # Create integrated form factor + ff = deepcopy(self.scattering_factor) + ff.x *= ff.x + ff.y *= ff.y/Z**2 + int_ff = Tabulated1D(ff.x, ff.integral()) + int_ff.to_hdf5(group, 'integrated_scattering_factor') + else: + self.scattering_factor.to_hdf5(group, 'scattering_factor') + if self.anomalous_real is not None: + self.anomalous_real.to_hdf5(group, 'anomalous_real') + if self.anomalous_imag is not None: + self.anomalous_imag.to_hdf5(group, 'anomalous_imag') From a2375910ac7d8fd474f7d9a45400e459661a9ba9 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 15:18:13 -0400 Subject: [PATCH 12/61] update heating tally to allow for analog estimator --- src/tallies/tally_scoring.cpp | 99 ++++++++++++++++++++++++++++++++++- 1 file changed, 97 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index d574ed109..9fd31925f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1152,9 +1152,104 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_HEATING: if (p->type_ == Particle::Type::neutron) { score_bin = NEUTRON_HEATING; - // No break here, continue to run the default neutron MT scoring + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Calculate material heating cross section + double macro_heating = 0.; + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = simulation::micro_xs[j_nuclide].index_grid; + auto f = simulation::micro_xs[j_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + if (i_grid >= xs.threshold) { + macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density; + } + } + } + // All events score to an heating rate bin + if (settings::survival_biasing) { + // We need to account for the fact that some weight was already + // absorbed + score = p->wgt_last_ + p->wgt_absorb_; + } else { + score = p->wgt_last_; + } + score *= macro_heating * flux / simulation::material_xs.total; + } else { + // Calculate neutron heating cross section on-the-fly + score = 0.; + if (i_nuclide >= 0) { + const auto& nuc {*data::nuclides[i_nuclide]}; + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = simulation::micro_xs[i_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = simulation::micro_xs[i_nuclide].index_grid; + auto f = simulation::micro_xs[i_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + if (i_grid >= xs.threshold) { + score = ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux; + } + } + } else { + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = simulation::micro_xs[j_nuclide].index_grid; + auto f = simulation::micro_xs[j_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + if (i_grid >= xs.threshold) { + score += ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density + * flux; + } + } + } + } + } + } } else if (p->type_ == Particle::Type::photon) { - if (tally.estimator_ != ESTIMATOR_ANALOG) { + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // Calculate material heating cross section + double macro_heating = 0.; + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_element = material.element_[i]; + auto atom_density = material.atom_density_(i); + auto& heating {data::elements[i_element].heating_}; + auto i_grid = simulation::micro_photon_xs[i_element].index_grid; + auto f = simulation::micro_photon_xs[i_element].interp_factor; + macro_heating += std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density; + } + // All events score to an heating rate bin + if (settings::survival_biasing) { + // We need to account for the fact that some weight was already + // absorbed + score = p->wgt_last_ + p->wgt_absorb_; + } else { + score = p->wgt_last_; + } + score *= macro_heating * flux / simulation::material_xs.total; + } else { // Calculate photon heating cross section on-the-fly score = 0.; if (i_nuclide >= 0) { From 90c70f910ad55512f2bbcf9b24e7c0bb59cf60cd Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 20 Mar 2019 21:47:14 -0400 Subject: [PATCH 13/61] first result looks good --- openmc/data/photon.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 73c48527d..f8400b736 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -1174,8 +1174,7 @@ class PhotonReaction(EqualityMixin): ff.y *= ff.y/Z**2 int_ff = Tabulated1D(ff.x, ff.integral()) int_ff.to_hdf5(group, 'integrated_scattering_factor') - else: - self.scattering_factor.to_hdf5(group, 'scattering_factor') + self.scattering_factor.to_hdf5(group, 'scattering_factor') if self.anomalous_real is not None: self.anomalous_real.to_hdf5(group, 'anomalous_real') if self.anomalous_imag is not None: From 21888e6cb56939283ffa1a12435ee10eda69aa8b Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 21 Mar 2019 11:38:49 -0400 Subject: [PATCH 14/61] adjust heating number when reading from ACE --- openmc/data/photon.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index f8400b736..4e266bc9d 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -513,10 +513,6 @@ class IncidentPhoton(EqualityMixin): for mt in (502, 504, 515, 522, 525): data.reactions[mt] = PhotonReaction.from_ace(ace, mt) - # Get heating cross sections [eV*barn] from factors [MeV per collision] - data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in - (502, 504, 515, 522)]) * EV_PER_MEV - # Compton profiles n_shell = ace.nxs[5] if n_shell != 0: @@ -997,7 +993,17 @@ class PhotonReaction(EqualityMixin): # Store cross section xs = ace.xss[idx : idx+n].copy() - if mt in (502, 504, 515, 522): + if mt == 525: + # Get heating cross sections [eV barn] from factors [MeV/collision] + # by multiplying the total xs (sum of (502, 504, 515, 522)) + idx = ace.jxs[1] + n + total_xs = ace.xss[idx : idx + 4*n].copy() + total_xs = total_xs.reshape((4, n)) + nonzero = (total_xs != 0.0) + total_xs[nonzero] = np.exp(total_xs[nonzero]) + total_xs = total_xs.sum(axis=0) + xs *= total_xs * EV_PER_MEV + else: nonzero = (xs != 0.0) xs[nonzero] = np.exp(xs[nonzero]) rx.xs = Tabulated1D(energy, xs, [n], [5]) From acf4b3616c3c87e75fc50bcc5c87377e0f6ccbc6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 21 Mar 2019 11:45:57 -0400 Subject: [PATCH 15/61] fix photon reaction __repr__ --- openmc/data/photon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 4e266bc9d..b2edebacd 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -904,7 +904,7 @@ class PhotonReaction(EqualityMixin): def __repr__(self): if self.mt in _REACTION_NAME: return "".format( - self.mt, _REACTION_NAME[self.mt]) + self.mt, _REACTION_NAME[self.mt][0]) else: return "".format(self.mt) From e1a575b4d02ddd383f6cae4278e66a7496016065 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 21 Mar 2019 16:53:20 -0400 Subject: [PATCH 16/61] fix zero heating for analog estimator tally --- src/tallies/tally_scoring.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 9fd31925f..820b51401 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1278,11 +1278,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } } - break; - } else { - // Nuclear heating of any other particles not implmented - break; } + // Nuclear heating of any other particles not implmented + break; default: if (tally.estimator_ == ESTIMATOR_ANALOG) { From dfdb07471815c797d0d1d7189529c61faeb3852e Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 23 Mar 2019 15:24:49 -0400 Subject: [PATCH 17/61] currently only the heating factors (not heating cross sections) can be extracted from ace as there is no 517 cross section in photoatomic library --- openmc/data/photon.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index b2edebacd..57f68510c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -513,6 +513,12 @@ class IncidentPhoton(EqualityMixin): for mt in (502, 504, 515, 522, 525): data.reactions[mt] = PhotonReaction.from_ace(ace, mt) + # Get heating cross sections [eV*barn] from factors [eV per collision] + # by multiplying with total xs (sum of (502, 504, 515, 517, 522)) + # but currently we cannot do this right as we miss 517) + data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in + (502, 504, 515, 522)]) + # Compton profiles n_shell = ace.nxs[5] if n_shell != 0: @@ -994,15 +1000,8 @@ class PhotonReaction(EqualityMixin): # Store cross section xs = ace.xss[idx : idx+n].copy() if mt == 525: - # Get heating cross sections [eV barn] from factors [MeV/collision] - # by multiplying the total xs (sum of (502, 504, 515, 522)) - idx = ace.jxs[1] + n - total_xs = ace.xss[idx : idx + 4*n].copy() - total_xs = total_xs.reshape((4, n)) - nonzero = (total_xs != 0.0) - total_xs[nonzero] = np.exp(total_xs[nonzero]) - total_xs = total_xs.sum(axis=0) - xs *= total_xs * EV_PER_MEV + # Get heating factors in [eV per collision] + xs *= EV_PER_MEV else: nonzero = (xs != 0.0) xs[nonzero] = np.exp(xs[nonzero]) From cdde6b55b33e22b346819444644d824780330ece Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 26 Mar 2019 12:48:12 -0400 Subject: [PATCH 18/61] generate photon heating number from ENDF --- openmc/data/photon.py | 55 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 57f68510c..78b48593e 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -9,6 +9,8 @@ import h5py import numpy as np import pandas as pd from scipy.interpolate import CubicSpline +from scipy.integrate import quad +from warnings import warn from openmc.mixin import EqualityMixin import openmc.checkvalue as cv @@ -19,13 +21,20 @@ from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record from .function import Tabulated1D +# Constants +MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy +PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms +FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant +CM_PER_ANGSTROM = 1.0e-8 +RE = CM_PER_ANGSTROM * PLANCK_C / (2.0 * np.pi * FINE_STRUCTURE * MASS_ELECTRON_EV) + +# Electron subshell labels _SUBSHELLS = ['K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', 'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3', 'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11', 'Q1', 'Q2', 'Q3'] - # Helper function to map designator to subshell string or None def _subshell(i): if i == 0: @@ -648,6 +657,9 @@ class IncidentPhoton(EqualityMixin): # Add bremsstrahlung DCS data data._add_bremsstrahlung() + # Add heating cross sections + #data._compute_heating() + return data @classmethod @@ -876,6 +888,47 @@ class IncidentPhoton(EqualityMixin): self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy'] self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number]) + def _compute_heating(self): + """Compute heating cross sections + + """ + # Determine union energy grid + energy = np.array([]) + for mt in (504, 515, 517, 522): + if mt not in self: + warn('Cross sections for MT={} do not exist. The total heating ' + 'cross sections may be wrong'.format(mt)) + else: + energy = np.union1d(energy, self[mt].xs.x) + + heating_xs = np.zeros_like(energy) + + if 504 in self: + rx = self[504] + def xs_mu(mu, E): + k = E/MASS_ELECTRON_EV + k_p = k/(1.0 + k*(1.0 - mu)) + x = E * np.sqrt(0.5*(1.0 - mu)) / PLANCK_C + return np.pi * RE*RE * k_p/k * k_p/k * (k_p + 1.0/k_p + + mu*mu - 1.0) * rx.scattering_factor(x) + def xs_mu_E(mu, E): + E_p = E/(1.0 + E/MASS_ELECTRON_EV*(1-mu)) + return xs_mu(mu, E) * E_p + e_new = np.array([quad(xs_mu_E, -1, 1, args=(e,))[0]/quad(xs_mu, -1, 1, args=(e,))[0] for e in energy]) + heating_xs += (energy - e_new)*rx.xs(energy) + + if 515 in self: + heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy) + + if 517 in self: + heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[517].xs(energy) + + if 522 in self: + heating_xs += (energy - 0.)*self[522].xs(energy) + + h_rx = PhotonReaction(525) + h_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5]) + self.reactions[525] = h_rx class PhotonReaction(EqualityMixin): """Photon-induced reaction From cfb11882a74bfa8a1df2adde82e382a933c556d2 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 26 Mar 2019 13:06:22 -0400 Subject: [PATCH 19/61] update tally_scoring --- src/tallies/tally_scoring.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 820b51401..f249bc121 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1163,10 +1163,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + auto i_temp = p->neutron_xs_[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[j_nuclide].index_grid; - auto f = simulation::micro_xs[j_nuclide].interp_factor; + auto i_grid = p->neutron_xs_[j_nuclide].index_grid; + auto f = p->neutron_xs_[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1182,7 +1182,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= macro_heating * flux / simulation::material_xs.total; + score *= macro_heating * flux / p->macro_xs_.total; } else { // Calculate neutron heating cross section on-the-fly score = 0.; @@ -1191,10 +1191,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = simulation::micro_xs[i_nuclide].index_temp; + auto i_temp = p->neutron_xs_[i_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[i_nuclide].index_grid; - auto f = simulation::micro_xs[i_nuclide].interp_factor; + auto i_grid = p->neutron_xs_[i_nuclide].index_grid; + auto f = p->neutron_xs_[i_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score = ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1211,10 +1211,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = simulation::micro_xs[j_nuclide].index_temp; + auto i_temp = p->neutron_xs_[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[j_nuclide].index_grid; - auto f = simulation::micro_xs[j_nuclide].interp_factor; + auto i_grid = p->neutron_xs_[j_nuclide].index_grid; + auto f = p->neutron_xs_[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score += ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1248,7 +1248,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= macro_heating * flux / simulation::material_xs.total; + score *= macro_heating * flux / p->macro_xs_.total; } else { // Calculate photon heating cross section on-the-fly score = 0.; From deda8ce0767cc89e53f9b3fac7ee6e1882758f39 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 27 Mar 2019 15:26:42 -0400 Subject: [PATCH 20/61] tally photon heating with energy lost for analog estimator --- src/tallies/tally_scoring.cpp | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index f249bc121..b0bf8f448 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1228,19 +1228,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } else if (p->type_ == Particle::Type::photon) { if (tally.estimator_ == ESTIMATOR_ANALOG) { - // Calculate material heating cross section - double macro_heating = 0.; - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto i_element = material.element_[i]; - auto atom_density = material.atom_density_(i); - auto& heating {data::elements[i_element].heating_}; - auto i_grid = simulation::micro_photon_xs[i_element].index_grid; - auto f = simulation::micro_photon_xs[i_element].interp_factor; - macro_heating += std::exp(heating(i_grid) + f * (heating(i_grid+1) - - heating(i_grid))) * atom_density; - } - // All events score to an heating rate bin + // Use photon energy lost for heating deposition if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed @@ -1248,7 +1236,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= macro_heating * flux / p->macro_xs_.total; + score *= (E - p->E_) * flux; } else { // Calculate photon heating cross section on-the-fly score = 0.; From f5cff7baefeecd259e8d0714382c036443d88cbf Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 28 Mar 2019 16:58:46 -0400 Subject: [PATCH 21/61] generate heating number --- openmc/data/photon.py | 11 ++++++----- src/tallies/tally_scoring.cpp | 8 ++++---- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 78b48593e..fd6a8fd17 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -523,8 +523,7 @@ class IncidentPhoton(EqualityMixin): data.reactions[mt] = PhotonReaction.from_ace(ace, mt) # Get heating cross sections [eV*barn] from factors [eV per collision] - # by multiplying with total xs (sum of (502, 504, 515, 517, 522)) - # but currently we cannot do this right as we miss 517) + # by multiplying with total xs (sum of (502, 504, 515, 522)) data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in (502, 504, 515, 522)]) @@ -658,7 +657,7 @@ class IncidentPhoton(EqualityMixin): data._add_bremsstrahlung() # Add heating cross sections - #data._compute_heating() + data._compute_heating() return data @@ -909,12 +908,14 @@ class IncidentPhoton(EqualityMixin): k = E/MASS_ELECTRON_EV k_p = k/(1.0 + k*(1.0 - mu)) x = E * np.sqrt(0.5*(1.0 - mu)) / PLANCK_C - return np.pi * RE*RE * k_p/k * k_p/k * (k_p + 1.0/k_p + + return np.pi*RE*RE*k_p/k*k_p/k*(k_p + 1.0/k_p + mu*mu - 1.0) * rx.scattering_factor(x) def xs_mu_E(mu, E): E_p = E/(1.0 + E/MASS_ELECTRON_EV*(1-mu)) return xs_mu(mu, E) * E_p - e_new = np.array([quad(xs_mu_E, -1, 1, args=(e,))[0]/quad(xs_mu, -1, 1, args=(e,))[0] for e in energy]) + e_new = np.array([quad(xs_mu_E, -1, 1, args=(e,), epsabs=0.0, + epsrel=1.0e-3)[0] / quad(xs_mu, -1, 1, args=(e,), + epsabs=0.0, epsrel=1.0e-3)[0] for e in energy]) heating_xs += (energy - e_new)*rx.xs(energy) if 515 in self: diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index b0bf8f448..8fb85db5d 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1247,8 +1247,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, std::string element = name.substr(0, pos); int i_element = data::element_map[element]; auto& heating {data::elements[i_element].heating_}; - auto i_grid = simulation::micro_photon_xs[i_element].index_grid; - auto f = simulation::micro_photon_xs[i_element].interp_factor; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - heating(i_grid))) * atom_density * flux; } else { @@ -1258,8 +1258,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto i_element = material.element_[i]; auto atom_density = material.atom_density_(i); auto& heating {data::elements[i_element].heating_}; - auto i_grid = simulation::micro_photon_xs[i_element].index_grid; - auto f = simulation::micro_photon_xs[i_element].interp_factor; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - heating(i_grid))) * atom_density * flux; } From ebf4482d5b481a5c5d538e9d71c1864d218aaf6c Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 28 Mar 2019 16:59:45 -0400 Subject: [PATCH 22/61] fix photon analog --- src/physics.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/physics.cpp b/src/physics.cpp index 1cfd31658..258cea924 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -56,7 +56,6 @@ void collision(Particle* p) if (p->E_ < settings::energy_cutoff[type]) { p->alive_ = false; p->wgt_ = 0.0; - p->wgt_last_ = 0.0; } // Display information about collision From a034fb2e767cc6e76aad321ca4763a27b67fcc27 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 4 Apr 2019 11:22:03 -0400 Subject: [PATCH 23/61] added a data member in particle to track # of secondary particles --- include/openmc/constants.h | 5 ++++- include/openmc/particle.h | 3 ++- src/particle.cpp | 10 +++++++++- src/physics.cpp | 10 +++++++++- 4 files changed, 24 insertions(+), 4 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index dd9ecf211..a6397b2a7 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -105,9 +105,10 @@ constexpr std::array SUBSHELLS = { "Q1", "Q2", "Q3" }; -// Void material +// Void material and nuclide // TODO: refactor and remove constexpr int MATERIAL_VOID {-1}; +constexpr int NUCLIDE_VOID {-1}; // ============================================================================ // CROSS SECTION RELATED CONSTANTS @@ -127,6 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2}; // Reaction types // TODO: Convert to enum +constexpr int NONE_REACTION {0}; constexpr int TOTAL_XS {1}; constexpr int ELASTIC {2}; constexpr int N_NONELASTIC {3}; @@ -344,6 +346,7 @@ constexpr int ESTIMATOR_COLLISION {3}; // TODO: Convert to enum constexpr int EVENT_SURFACE {-2}; constexpr int EVENT_LATTICE {-1}; +constexpr int EVENT_KILL {0}; constexpr int EVENT_SCATTER {1}; constexpr int EVENT_ABSORB {2}; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 8ec53de44..fcaf209a0 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -184,7 +184,7 @@ public: //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] //! \param type Particle type - void create_secondary(Direction u, double E, Type type) const; + void create_secondary(Direction u, double E, Type type); //! initialize from a source site // @@ -261,6 +261,7 @@ public: // Post-collision physical data int n_bank_ {0}; //!< number of fission sites banked + int n_bank_second_ {0}; //!< number of secondary particles double wgt_bank_ {0.0}; //!< weight of fission sites banked int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission //!< sites banked diff --git a/src/particle.cpp b/src/particle.cpp index 4e0066b20..42889a981 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -73,7 +73,7 @@ Particle::clear() } void -Particle::create_secondary(Direction u, double E, Type type) const +Particle::create_secondary(Direction u, double E, Type type) { simulation::secondary_bank.emplace_back(); @@ -83,6 +83,8 @@ Particle::create_secondary(Direction u, double E, Type type) const bank.r = this->r(); bank.u = u; bank.E = settings::run_CE ? E : g_; + + n_bank_second_ += 1; } void @@ -157,6 +159,11 @@ Particle::transport() u_last_ = this->u(); r_last_ = this->r(); + // Reset event variables + event_ = EVENT_KILL; + event_nuclide_ = NUCLIDE_VOID; + event_mt_ = NONE_REACTION; + // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles @@ -309,6 +316,7 @@ Particle::transport() // Reset banked weight during collision n_bank_ = 0; + n_bank_second_ = 0; wgt_bank_ = 0.0; for (int& v : n_delayed_bank_) v = 0; diff --git a/src/physics.cpp b/src/physics.cpp index 67af6dc42..018677a54 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -61,7 +61,9 @@ void collision(Particle* p) // Display information about collision if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; - if (p->type_ == Particle::Type::neutron) { + if (p->event_ == EVENT_KILL) { + msg << " " << " Killed. Energy = " << p->E_ << " eV."; + } else if (p->type_ == Particle::Type::neutron) { msg << " " << reaction_name(p->event_mt_) << " with " << data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else { @@ -226,6 +228,7 @@ void sample_photon_reaction(Particle* p) if (prob > cutoff) { double mu = element.rayleigh_scatter(alpha); p->u() = rotate_angle(p->u(), mu, nullptr); + p->event_ = EVENT_SCATTER; p->event_mt_ = COHERENT; return; } @@ -268,6 +271,7 @@ void sample_photon_reaction(Particle* p) phi += PI; p->E_ = alpha_out*MASS_ELECTRON_EV; p->u() = rotate_angle(p->u(), mu, &phi); + p->event_ = EVENT_SCATTER; p->event_mt_ = INCOHERENT; return; } @@ -319,6 +323,7 @@ void sample_photon_reaction(Particle* p) // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons element.atomic_relaxation(shell, *p); + p->event_ = EVENT_ABSORB; p->event_mt_ = 533 + shell.index_subshell; p->alive_ = false; p->E_ = 0.0; @@ -344,6 +349,7 @@ void sample_photon_reaction(Particle* p) u = rotate_angle(p->u(), mu_positron, nullptr); p->create_secondary(u, E_positron, Particle::Type::positron); + p->event_ = EVENT_ABSORB; p->event_mt_ = PAIR_PROD; p->alive_ = false; p->E_ = 0.0; @@ -361,6 +367,7 @@ void sample_electron_reaction(Particle* p) p->E_ = 0.0; p->alive_ = false; + p->event_ = EVENT_ABSORB; } void sample_positron_reaction(Particle* p) @@ -386,6 +393,7 @@ void sample_positron_reaction(Particle* p) p->E_ = 0.0; p->alive_ = false; + p->event_ = EVENT_ABSORB; } int sample_nuclide(const Particle* p) From c91aa58e8e4f50fa6fe92e8cb3b0b9d52447640d Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 4 Apr 2019 15:26:05 -0400 Subject: [PATCH 24/61] added analog heating tally for non-neutron particles --- src/tallies/tally_scoring.cpp | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 8fb85db5d..3a759496b 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1151,7 +1151,6 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_HEATING: if (p->type_ == Particle::Type::neutron) { - score_bin = NEUTRON_HEATING; if (tally.estimator_ == ESTIMATOR_ANALOG) { // Calculate material heating cross section double macro_heating = 0.; @@ -1160,7 +1159,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[score_bin]; + auto m = nuc.reaction_index_[NEUTRON_HEATING]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; auto i_temp = p->neutron_xs_[j_nuclide].index_temp; @@ -1188,7 +1187,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = 0.; if (i_nuclide >= 0) { const auto& nuc {*data::nuclides[i_nuclide]}; - auto m = nuc.reaction_index_[score_bin]; + auto m = nuc.reaction_index_[NEUTRON_HEATING]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; auto i_temp = p->neutron_xs_[i_nuclide].index_temp; @@ -1208,7 +1207,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[score_bin]; + auto m = nuc.reaction_index_[NEUTRON_HEATING]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; auto i_temp = p->neutron_xs_[j_nuclide].index_temp; @@ -1226,18 +1225,19 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } } - } else if (p->type_ == Particle::Type::photon) { + } else { if (tally.estimator_ == ESTIMATOR_ANALOG) { - // Use photon energy lost for heating deposition - if (settings::survival_biasing) { - // We need to account for the fact that some weight was already - // absorbed - score = p->wgt_last_ + p->wgt_absorb_; - } else { - score = p->wgt_last_; + // Score direct energy deposition in the collision + score = E - p->E_; + // We need to substract the energy of the secondary particles since + // they will be transported individually later + for (auto i = 0; i < p->n_bank_second_; ++i) { + auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i; + const auto& bank = simulation::secondary_bank[i_bank]; + score -= bank.E; } - score *= (E - p->E_) * flux; - } else { + score *= p->wgt_last_ * flux; + } else if (p->type_ == Particle::Type::photon) { // Calculate photon heating cross section on-the-fly score = 0.; if (i_nuclide >= 0) { @@ -1267,7 +1267,6 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } } - // Nuclear heating of any other particles not implmented break; default: From ec49c32a619dad1ee348931082982fa44f28e5e6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 4 Apr 2019 16:02:36 -0400 Subject: [PATCH 25/61] fix positron and electron collision message --- src/physics.cpp | 6 ++-- src/tallies/tally.cpp | 3 +- src/tallies/tally_scoring.cpp | 52 ++++++++++++++++++----------------- 3 files changed, 33 insertions(+), 28 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 018677a54..e8611dbd7 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -62,13 +62,15 @@ void collision(Particle* p) if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; if (p->event_ == EVENT_KILL) { - msg << " " << " Killed. Energy = " << p->E_ << " eV."; + msg << " Killed. Energy = " << p->E_ << " eV."; } else if (p->type_ == Particle::Type::neutron) { msg << " " << reaction_name(p->event_mt_) << " with " << data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; - } else { + } else if (p->type_ == Particle::Type::photon) { msg << " " << reaction_name(p->event_mt_) << " with " << data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV."; + } else { + msg << " Disappeared. Energy = " << p->E_ << " eV."; } write_message(msg, 1); } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index d6ccff0f1..2cb0945ab 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -9,6 +9,7 @@ #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/particle.h" +#include "openmc/reaction.h" #include "openmc/reaction_product.h" #include "openmc/settings.h" #include "openmc/simulation.h" @@ -772,7 +773,7 @@ void read_tallies_xml() case SCORE_PROMPT_NU_FISSION: case SCORE_DECAY_RATE: warning("Particle filter is not used with photon transport" - " on and " + std::to_string(score) + " score."); + " on and " + reaction_name(score) + " score."); break; } } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 3a759496b..ea046ed04 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1237,31 +1237,33 @@ score_general_ce(Particle* p, int i_tally, int start_index, score -= bank.E; } score *= p->wgt_last_ * flux; - } else if (p->type_ == Particle::Type::photon) { - // Calculate photon heating cross section on-the-fly - score = 0.; - if (i_nuclide >= 0) { - // Find the element corresponding to the nuclide - auto name = data::nuclides[i_nuclide]->name_; - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(0, pos); - int i_element = data::element_map[element]; - auto& heating {data::elements[i_element].heating_}; - auto i_grid = p->photon_xs_[i_element].index_grid; - auto f = p->photon_xs_[i_element].interp_factor; - score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - - heating(i_grid))) * atom_density * flux; - } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto i_element = material.element_[i]; - auto atom_density = material.atom_density_(i); - auto& heating {data::elements[i_element].heating_}; - auto i_grid = p->photon_xs_[i_element].index_grid; - auto f = p->photon_xs_[i_element].interp_factor; - score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - - heating(i_grid))) * atom_density * flux; + } else { + if (p->type_ == Particle::Type::photon) { + // Calculate photon heating cross section on-the-fly + score = 0.; + if (i_nuclide >= 0) { + // Find the element corresponding to the nuclide + auto name = data::nuclides[i_nuclide]->name_; + int pos = name.find_first_of("0123456789"); + std::string element = name.substr(0, pos); + int i_element = data::element_map[element]; + auto& heating {data::elements[i_element].heating_}; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; + score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; + } else { + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_element = material.element_[i]; + auto atom_density = material.atom_density_(i); + auto& heating {data::elements[i_element].heating_}; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; + score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; + } } } } From f4541ab1c8a58f032803589f88fc08127eb2d595 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 5 Apr 2019 22:18:22 -0400 Subject: [PATCH 26/61] photon heating tally with analog estimator works now --- src/tallies/tally_scoring.cpp | 61 ++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index ea046ed04..b2e0761dd 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1150,6 +1150,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_HEATING: + score = 0.; if (p->type_ == Particle::Type::neutron) { if (tally.estimator_ == ESTIMATOR_ANALOG) { // Calculate material heating cross section @@ -1184,7 +1185,6 @@ score_general_ce(Particle* p, int i_tally, int start_index, score *= macro_heating * flux / p->macro_xs_.total; } else { // Calculate neutron heating cross section on-the-fly - score = 0.; if (i_nuclide >= 0) { const auto& nuc {*data::nuclides[i_nuclide]}; auto m = nuc.reaction_index_[NEUTRON_HEATING]; @@ -1225,7 +1225,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, } } } - } else { + } else if (p->type_ == Particle::Type::photon) { if (tally.estimator_ == ESTIMATOR_ANALOG) { // Score direct energy deposition in the collision score = E - p->E_; @@ -1234,36 +1234,39 @@ score_general_ce(Particle* p, int i_tally, int start_index, for (auto i = 0; i < p->n_bank_second_; ++i) { auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i; const auto& bank = simulation::secondary_bank[i_bank]; - score -= bank.E; + if (bank.particle == Particle::Type::photon || + bank.particle == Particle::Type::neutron) { + score -= bank.E; + } else if (bank.particle == Particle::Type::positron) { + // Annihilation of the positron will produce two new photons + score -= 2*MASS_ELECTRON_EV; + } } score *= p->wgt_last_ * flux; } else { - if (p->type_ == Particle::Type::photon) { - // Calculate photon heating cross section on-the-fly - score = 0.; - if (i_nuclide >= 0) { - // Find the element corresponding to the nuclide - auto name = data::nuclides[i_nuclide]->name_; - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(0, pos); - int i_element = data::element_map[element]; - auto& heating {data::elements[i_element].heating_}; - auto i_grid = p->photon_xs_[i_element].index_grid; - auto f = p->photon_xs_[i_element].interp_factor; - score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - - heating(i_grid))) * atom_density * flux; - } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto i_element = material.element_[i]; - auto atom_density = material.atom_density_(i); - auto& heating {data::elements[i_element].heating_}; - auto i_grid = p->photon_xs_[i_element].index_grid; - auto f = p->photon_xs_[i_element].interp_factor; - score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - - heating(i_grid))) * atom_density * flux; - } + // Calculate photon heating cross section on-the-fly + if (i_nuclide >= 0) { + // Find the element corresponding to the nuclide + auto name = data::nuclides[i_nuclide]->name_; + int pos = name.find_first_of("0123456789"); + std::string element = name.substr(0, pos); + int i_element = data::element_map[element]; + auto& heating {data::elements[i_element].heating_}; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; + score = std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; + } else { + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_element = material.element_[i]; + auto atom_density = material.atom_density_(i); + auto& heating {data::elements[i_element].heating_}; + auto i_grid = p->photon_xs_[i_element].index_grid; + auto f = p->photon_xs_[i_element].interp_factor; + score += std::exp(heating(i_grid) + f * (heating(i_grid+1) - + heating(i_grid))) * atom_density * flux; } } } From 5db6db53f11ebd2af2e85159d623a8d59f4744cc Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 5 Apr 2019 22:18:51 -0400 Subject: [PATCH 27/61] fix photon heating generation from endf --- openmc/data/photon.py | 59 ++++++++++++++++++++++++++++++------------- 1 file changed, 41 insertions(+), 18 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index fb40b0a6b..974b9323c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -26,7 +26,8 @@ MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant CM_PER_ANGSTROM = 1.0e-8 -RE = CM_PER_ANGSTROM * PLANCK_C / (2.0 * np.pi * FINE_STRUCTURE * MASS_ELECTRON_EV) +# classical electron radius in cm +R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * np.pi * FINE_STRUCTURE * MASS_ELECTRON_EV) # Electron subshell labels _SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', @@ -961,10 +962,27 @@ class IncidentPhoton(EqualityMixin): self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number]) def _compute_heating(self): - """Compute heating cross sections + """Compute heating cross sections (KERMA) + + Photon energy is deposited as energy loss in three reactions: + incoherent scattering, photoelectric effect and pair production. + The point-wise heating cross section is calculated as: + + .. math:: + \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PE, PP \right\} + + \overline{E}_I (E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} + + \overline{E}_{PE} &= 0 + + \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV + + The differential cross section for incoherent scattering can be + found in the theory manual. """ - # Determine union energy grid + + # Determine a union energy grid energy = np.array([]) for mt in (504, 515, 517, 522): if mt not in self: @@ -977,19 +995,24 @@ class IncidentPhoton(EqualityMixin): if 504 in self: rx = self[504] - def xs_mu(mu, E): - k = E/MASS_ELECTRON_EV - k_p = k/(1.0 + k*(1.0 - mu)) - x = E * np.sqrt(0.5*(1.0 - mu)) / PLANCK_C - return np.pi*RE*RE*k_p/k*k_p/k*(k_p + 1.0/k_p + + def dsigma_dmu(mu, E): + k = E / MASS_ELECTRON_EV + kout = k / (1.0 + k * (1.0 - mu)) + x = E * np.sqrt(0.5 * (1.0 - mu)) / PLANCK_C + return np.pi * R0**2 * (kout/k)**2 * (kout/k + k/kout + mu*mu - 1.0) * rx.scattering_factor(x) - def xs_mu_E(mu, E): - E_p = E/(1.0 + E/MASS_ELECTRON_EV*(1-mu)) - return xs_mu(mu, E) * E_p - e_new = np.array([quad(xs_mu_E, -1, 1, args=(e,), epsabs=0.0, - epsrel=1.0e-3)[0] / quad(xs_mu, -1, 1, args=(e,), - epsabs=0.0, epsrel=1.0e-3)[0] for e in energy]) - heating_xs += (energy - e_new)*rx.xs(energy) + def dsigma_dmu_times_eout(mu, E): + Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1 - mu)) + return dsigma_dmu(mu, E) * Eout + def average_eout(E): + integral_sigma = quad(dsigma_dmu, -1.0, 1.0, + args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] + integral_sigma_e = quad(dsigma_dmu_times_eout, -1.0, 1.0, + args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] + return integral_sigma_e / integral_sigma + + e_out = np.vectorize(lambda x: average_eout(x))(energy) + heating_xs += (energy - e_out) * rx.xs(energy) if 515 in self: heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy) @@ -1000,9 +1023,9 @@ class IncidentPhoton(EqualityMixin): if 522 in self: heating_xs += (energy - 0.)*self[522].xs(energy) - h_rx = PhotonReaction(525) - h_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5]) - self.reactions[525] = h_rx + heat_rx = PhotonReaction(525) + heat_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5]) + self.reactions[525] = heat_rx class PhotonReaction(EqualityMixin): """Photon-induced reaction From 63527fe9087982b3b34bbd2e76ec19e666552392 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 6 Apr 2019 08:02:07 -0400 Subject: [PATCH 28/61] heating analog tally for specific nuclide --- src/tallies/tally_scoring.cpp | 72 ++++++++++++++++++++++++----------- 1 file changed, 49 insertions(+), 23 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index b2e0761dd..e7cde6bd7 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1153,28 +1153,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = 0.; if (p->type_ == Particle::Type::neutron) { if (tally.estimator_ == ESTIMATOR_ANALOG) { - // Calculate material heating cross section - double macro_heating = 0.; - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); - const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[NEUTRON_HEATING]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[j_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[j_nuclide].index_grid; - auto f = p->neutron_xs_[j_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density; - } - } - } - // All events score to an heating rate bin + // All events score to a heating tally bin. We actually use a + // collision estimator in place of an analog one since there is no + // reaction-wise heating cross section if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed @@ -1182,7 +1163,52 @@ score_general_ce(Particle* p, int i_tally, int start_index, } else { score = p->wgt_last_; } - score *= macro_heating * flux / p->macro_xs_.total; + if (i_nuclide >= 0) { + // Calculate nuclide heating cross section + double macro_heating = 0.; + const auto& nuc {*data::nuclides[i_nuclide]}; + auto m = nuc.reaction_index_[NEUTRON_HEATING]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = p->neutron_xs_[i_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = p->neutron_xs_[i_nuclide].index_grid; + auto f = p->neutron_xs_[i_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + if (i_grid >= xs.threshold) { + macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density; + } + } + score *= macro_heating * flux / p->neutron_xs_[i_nuclide].total; + } else { + if (p->material_ != MATERIAL_VOID) { + // Calculate material heating cross section + double macro_heating = 0.; + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const auto& nuc {*data::nuclides[j_nuclide]}; + auto m = nuc.reaction_index_[NEUTRON_HEATING]; + if (m == C_NONE) continue; + const auto& rxn {*nuc.reactions_[m]}; + auto i_temp = p->neutron_xs_[j_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + auto i_grid = p->neutron_xs_[j_nuclide].index_grid; + auto f = p->neutron_xs_[j_nuclide].interp_factor; + const auto& xs {rxn.xs_[i_temp]}; + if (i_grid >= xs.threshold) { + macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density; + } + } + } + score *= macro_heating * flux / p->macro_xs_.total; + } else { + score = 0.; + } + } } else { // Calculate neutron heating cross section on-the-fly if (i_nuclide >= 0) { From 7536de564480e592678de7173ea1496353acbcb6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 10 Apr 2019 13:46:44 -0400 Subject: [PATCH 29/61] Account for fluorescent photons release in heating cross section computing --- openmc/data/photon.py | 76 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 15 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 974b9323c..aa01eac9b 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -161,6 +161,7 @@ class AtomicRelaxation(EqualityMixin): self.binding_energy = binding_energy self.num_electrons = num_electrons self.transitions = transitions + self._e_fluorescence = {} @property def binding_energy(self): @@ -390,6 +391,41 @@ class AtomicRelaxation(EqualityMixin): _SUBSHELLS, range(len(_SUBSHELLS))) group.create_dataset('transitions', data=df.values.astype(float)) + def energy_fluorescence(self, shell): + """Compute expected energy of fluorescent photons for the shell + """ + + if shell not in self.binding_energy: + raise KeyError('Invalid shell {}.'.format(shell)) + + if shell in self._e_fluorescence: + # Already computed + return self._e_fluorescence[shell] + else: + e = 0.0 + if shell not in self.transitions or self.transitions[shell].empty: + e = self.binding_energy[shell] + else: + df = self.transitions[shell] + for index, row in df.iterrows(): + e_row = 0.0 + primary = row['secondary'] + secondary = row['tertiary'] + if secondary is None: + # Fluorescent photon release in radiative transition + e_row += row['energy (eV)'] + else: + # Fill pole left by auger eletron + e_row += self.energy_fluorescence(secondary) + + # Fill photonelectron pole + e_row += self.energy_fluorescence(primary) + + # Expected fluorescent photon energy + e += e_row * row['probability'] + + self._e_fluorescence[shell] = e + return e class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -965,20 +1001,20 @@ class IncidentPhoton(EqualityMixin): """Compute heating cross sections (KERMA) Photon energy is deposited as energy loss in three reactions: - incoherent scattering, photoelectric effect and pair production. + incoherent scattering, pair production and photoelectric effect. The point-wise heating cross section is calculated as: .. math:: - \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PE, PP \right\} + \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PP, PE \right\} \overline{E}_I (E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} - \overline{E}_{PE} &= 0 - \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV - The differential cross section for incoherent scattering can be - found in the theory manual. + \overline{E}_{PE} &= E_{fluorescent photons} + + The differential cross section representation for incoherent + scattering can be found in the theory manual. """ @@ -987,41 +1023,51 @@ class IncidentPhoton(EqualityMixin): for mt in (504, 515, 517, 522): if mt not in self: warn('Cross sections for MT={} do not exist. The total heating ' - 'cross sections may be wrong'.format(mt)) + 'to be calculated is probably incorrect'.format(mt)) else: energy = np.union1d(energy, self[mt].xs.x) heating_xs = np.zeros_like(energy) + # Incoherent scattering if 504 in self: rx = self[504] - def dsigma_dmu(mu, E): + def _dsigma_dmu(mu, E): k = E / MASS_ELECTRON_EV kout = k / (1.0 + k * (1.0 - mu)) x = E * np.sqrt(0.5 * (1.0 - mu)) / PLANCK_C return np.pi * R0**2 * (kout/k)**2 * (kout/k + k/kout + mu*mu - 1.0) * rx.scattering_factor(x) - def dsigma_dmu_times_eout(mu, E): + def _eout_dsigma_dmu(mu, E): Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1 - mu)) - return dsigma_dmu(mu, E) * Eout - def average_eout(E): - integral_sigma = quad(dsigma_dmu, -1.0, 1.0, + return Eout * _dsigma_dmu(mu, E) + def _eout_average(E): + #TODO: the integrals are currently time consuming + integral_sigma = quad(_dsigma_dmu, -1.0, 1.0, args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] - integral_sigma_e = quad(dsigma_dmu_times_eout, -1.0, 1.0, + integral_sigma_e = quad(_eout_dsigma_dmu, -1.0, 1.0, args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] return integral_sigma_e / integral_sigma - e_out = np.vectorize(lambda x: average_eout(x))(energy) + e_out = np.vectorize(lambda x: _eout_average(x))(energy) heating_xs += (energy - e_out) * rx.xs(energy) + # Pair production, electron field if 515 in self: heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy) + # Pair production, nuclear field if 517 in self: heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[517].xs(energy) + # Photoelectric effect if 522 in self: - heating_xs += (energy - 0.)*self[522].xs(energy) + # Account for fluorescent photons + for mt, rx in self.reactions.items(): + if mt >= 534 and mt <= 572: + shell = _REACTION_NAME[mt][1] + e_f = self.atomic_relaxation.energy_fluorescence(shell) + heating_xs += (energy - e_f) * rx.xs(energy) heat_rx = PhotonReaction(525) heat_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5]) From e7ab0b704c859e1caaaba7d6052bb2890eb3a530 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 10 Apr 2019 16:18:26 -0400 Subject: [PATCH 30/61] fix tallying of nculide-specific heating for neutrons --- src/tallies/tally_scoring.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e7cde6bd7..578aa004c 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1177,7 +1177,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density; + + f * xs.value[i_grid-xs.threshold+1]); } } score *= macro_heating * flux / p->neutron_xs_[i_nuclide].total; @@ -2130,12 +2130,11 @@ void score_analog_tally_ce(Particle* p) auto i_nuclide = tally.nuclides_[i]; // Tally this event in the present nuclide bin if that bin represents - // the event nuclide or the total material. Note that the i_nuclide - // and flux arguments for score_general are not used for analog - // tallies. + // the event nuclide or the total material. Note that the atomic + // density argument for score_general is not used for analog tallies. if (i_nuclide == p->event_nuclide_ || i_nuclide == -1) score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, - -1, -1., filter_weight); + i_nuclide, -1., filter_weight); } } else { From a795de2db59f58e14d676a40163306de810bd39c Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 10 Apr 2019 16:43:25 -0400 Subject: [PATCH 31/61] save event nuclide rather than event element so that it is compatible with nuclide tally filters --- include/openmc/physics.h | 2 +- src/physics.cpp | 21 +++++++++++++-------- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 107380eb8..f0e6dde01 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -44,7 +44,7 @@ void sample_positron_reaction(Particle* p); //! //! \param[in] p Particle //! \return Index in the data::nuclides vector -int sample_nuclide(const Particle* p); +int sample_nuclide(Particle* p); //! Determine the average total, prompt, and delayed neutrons produced from //! fission and creates appropriate bank sites. diff --git a/src/physics.cpp b/src/physics.cpp index e8611dbd7..9b4a1a9c5 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -68,7 +68,7 @@ void collision(Particle* p) data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else if (p->type_ == Particle::Type::photon) { msg << " " << reaction_name(p->event_mt_) << " with " << - data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV."; + data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else { msg << " Disappeared. Energy = " << p->E_ << " eV."; } @@ -81,9 +81,6 @@ void sample_neutron_reaction(Particle* p) // Sample a nuclide within the material int i_nuclide = sample_nuclide(p); - // Save which nuclide particle had collision with - p->event_nuclide_ = i_nuclide; - // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not // change when sampling fission sites. The following block handles all @@ -213,7 +210,6 @@ void sample_photon_reaction(Particle* p) // Sample element within material int i_element = sample_element(p); - p->event_nuclide_ = i_element; const auto& micro {p->photon_xs_[i_element]}; const auto& element {data::elements[i_element]}; @@ -398,7 +394,7 @@ void sample_positron_reaction(Particle* p) p->event_ = EVENT_ABSORB; } -int sample_nuclide(const Particle* p) +int sample_nuclide(Particle* p) { // Sample cumulative distribution function double cutoff = prn() * p->macro_xs_.total; @@ -415,7 +411,11 @@ int sample_nuclide(const Particle* p) // Increment probability to compare to cutoff prob += atom_density * p->neutron_xs_[i_nuclide].total; - if (prob >= cutoff) return i_nuclide; + if (prob >= cutoff) { + // Save which nuclide particle had collision with + p->event_nuclide_ = i_nuclide; + return i_nuclide; + } } // If we reach here, no nuclide was sampled @@ -442,7 +442,12 @@ int sample_element(Particle* p) // Increment probability to compare to cutoff prob += sigma; - if (prob > cutoff) return i_element; + if (prob > cutoff) { + // Save which nuclide particle had collision with + p->event_nuclide_ = mat->nuclide_[i]; + + return i_element; + } } // If we made it here, no element was sampled From 57fcb078d325c53480c244c4274b8a1cee44196f Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 11 Apr 2019 12:14:18 -0400 Subject: [PATCH 32/61] refactoring --- include/openmc/constants.h | 4 ++-- include/openmc/particle.h | 2 +- include/openmc/physics.h | 2 +- src/particle.cpp | 4 ++-- src/physics.cpp | 13 ++++++------- 5 files changed, 12 insertions(+), 13 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a6397b2a7..49f980cc2 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -108,7 +108,7 @@ constexpr std::array SUBSHELLS = { // Void material and nuclide // TODO: refactor and remove constexpr int MATERIAL_VOID {-1}; -constexpr int NUCLIDE_VOID {-1}; +constexpr int NUCLIDE_NONE {-1}; // ============================================================================ // CROSS SECTION RELATED CONSTANTS @@ -128,7 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2}; // Reaction types // TODO: Convert to enum -constexpr int NONE_REACTION {0}; +constexpr int REACTION_NONE {0}; constexpr int TOTAL_XS {1}; constexpr int ELASTIC {2}; constexpr int N_NONELASTIC {3}; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index fcaf209a0..62923c11f 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -261,7 +261,7 @@ public: // Post-collision physical data int n_bank_ {0}; //!< number of fission sites banked - int n_bank_second_ {0}; //!< number of secondary particles + int n_bank_second_ {0}; //!< number of secondary particles banked double wgt_bank_ {0.0}; //!< weight of fission sites banked int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission //!< sites banked diff --git a/include/openmc/physics.h b/include/openmc/physics.h index f0e6dde01..107380eb8 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -44,7 +44,7 @@ void sample_positron_reaction(Particle* p); //! //! \param[in] p Particle //! \return Index in the data::nuclides vector -int sample_nuclide(Particle* p); +int sample_nuclide(const Particle* p); //! Determine the average total, prompt, and delayed neutrons produced from //! fission and creates appropriate bank sites. diff --git a/src/particle.cpp b/src/particle.cpp index 42889a981..cbc6217ca 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -161,8 +161,8 @@ Particle::transport() // Reset event variables event_ = EVENT_KILL; - event_nuclide_ = NUCLIDE_VOID; - event_mt_ = NONE_REACTION; + event_nuclide_ = NUCLIDE_NONE; + event_mt_ = REACTION_NONE; // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the diff --git a/src/physics.cpp b/src/physics.cpp index 9b4a1a9c5..cd29e97c2 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -81,6 +81,9 @@ void sample_neutron_reaction(Particle* p) // Sample a nuclide within the material int i_nuclide = sample_nuclide(p); + // Save which nuclide particle had collision with + p->event_nuclide_ = i_nuclide; + // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not // change when sampling fission sites. The following block handles all @@ -394,7 +397,7 @@ void sample_positron_reaction(Particle* p) p->event_ = EVENT_ABSORB; } -int sample_nuclide(Particle* p) +int sample_nuclide(const Particle* p) { // Sample cumulative distribution function double cutoff = prn() * p->macro_xs_.total; @@ -411,11 +414,7 @@ int sample_nuclide(Particle* p) // Increment probability to compare to cutoff prob += atom_density * p->neutron_xs_[i_nuclide].total; - if (prob >= cutoff) { - // Save which nuclide particle had collision with - p->event_nuclide_ = i_nuclide; - return i_nuclide; - } + if (prob >= cutoff) return i_nuclide; } // If we reach here, no nuclide was sampled @@ -443,7 +442,7 @@ int sample_element(Particle* p) // Increment probability to compare to cutoff prob += sigma; if (prob > cutoff) { - // Save which nuclide particle had collision with + // Save which nuclide particle had collision with for tally purpose p->event_nuclide_ = mat->nuclide_[i]; return i_element; From 85b1af458ab15240249dd18a92a48c91abd66345 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 11 Apr 2019 13:59:11 -0400 Subject: [PATCH 33/61] add to_element(), a function to convert nuclide name to element name --- include/openmc/string_utils.h | 2 ++ src/cross_sections.cpp | 3 +-- src/material.cpp | 3 +-- src/physics.cpp | 3 ++- src/string_utils.cpp | 6 ++++++ src/tallies/tally_scoring.cpp | 4 ++-- 6 files changed, 14 insertions(+), 7 deletions(-) diff --git a/include/openmc/string_utils.h b/include/openmc/string_utils.h index 35c25e548..db7635d0a 100644 --- a/include/openmc/string_utils.h +++ b/include/openmc/string_utils.h @@ -10,6 +10,8 @@ std::string& strtrim(std::string& s); char* strtrim(char* c_str); +std::string to_element(std::string const& name); + void to_lower(std::string& str); int word_count(std::string const& str); diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 2216b37da..5d8d18c21 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -239,8 +239,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, already_read.insert(name); // Check if elemental data has been read, if needed - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(0, pos); + std::string element = to_element(name); if (settings::photon_transport) { if (already_read.find(element) == already_read.end()) { // Read photon interaction data from HDF5 photon library diff --git a/src/material.cpp b/src/material.cpp index de0d7e5de..74d777fbd 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -227,8 +227,7 @@ Material::Material(pugi::xml_node node) // If the corresponding element hasn't been encountered yet and photon // transport will be used, we need to add its symbol to the element_dict if (settings::photon_transport) { - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(0, pos); + std::string element = to_element(name); // Make sure photon cross section data is available LibraryKey key {Library::Type::photon, element}; diff --git a/src/physics.cpp b/src/physics.cpp index cd29e97c2..0eb8c494e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -17,6 +17,7 @@ #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/string_utils.h" #include "openmc/thermal.h" #include "openmc/tallies/tally.h" @@ -68,7 +69,7 @@ void collision(Particle* p) data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else if (p->type_ == Particle::Type::photon) { msg << " " << reaction_name(p->event_mt_) << " with " << - data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; + to_element(data::nuclides[p->event_nuclide_]->name_) << ". Energy = " << p->E_ << " eV."; } else { msg << " Disappeared. Energy = " << p->E_ << " eV."; } diff --git a/src/string_utils.cpp b/src/string_utils.cpp index c1a203e48..40a20f232 100644 --- a/src/string_utils.cpp +++ b/src/string_utils.cpp @@ -26,6 +26,12 @@ char* strtrim(char* c_str) } +std::string to_element(std::string const& name) { + int pos = name.find_first_of("0123456789"); + return name.substr(0, pos); +} + + void to_lower(std::string& str) { for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 578aa004c..74cf9dcd4 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -12,6 +12,7 @@ #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/string_utils.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_delayedgroup.h" @@ -1274,8 +1275,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, if (i_nuclide >= 0) { // Find the element corresponding to the nuclide auto name = data::nuclides[i_nuclide]->name_; - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(0, pos); + std::string element = to_element(name); int i_element = data::element_map[element]; auto& heating {data::elements[i_element].heating_}; auto i_grid = p->photon_xs_[i_element].index_grid; From 39a38f9dca25358bf2568e8e4973a657c55562fc Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 11 Apr 2019 14:14:05 -0400 Subject: [PATCH 34/61] update document for photon heating --- docs/source/io_formats/nuclear_data.rst | 4 ++++ docs/source/usersguide/tallies.rst | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 81e4f4bcd..957bfa3a3 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -165,6 +165,10 @@ Incident Photon Data - **J** (*double[][]*) -- Compton profile for each subshell in units of :math:`\hbar / (me^2)` +**//heating/** + +:Datasets: - **xs** (*double[]*) -- Total heating cross section in [b-eV] + **//incoherent/** :Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b] diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 6ab4af84c..f005b5060 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -258,7 +258,10 @@ The following tables show all valid scores: +----------------------+---------------------------------------------------+ |heating |Total nuclear heating in units of eV per source | | |particle. For neutrons, this corresponds to MT=301 | - | |produced by NJOY's HEATR module. | + | |produced by NJOY's HEATR module while for photons, | + | |this is tallied from either direct photon energy | + | |deposition (analog estimator) or pre-generated | + | |photon heating number. | +----------------------+---------------------------------------------------+ |kappa-fission |The recoverable energy production rate due to | | |fission. The recoverable energy is defined as the | From 1969767c91ec98fad6a7c2e24e61aa9ba1ba7a99 Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 11 Apr 2019 15:57:39 -0400 Subject: [PATCH 35/61] fix formula --- openmc/data/photon.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 141838e49..cd1d8c05c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -549,8 +549,8 @@ class IncidentPhoton(EqualityMixin): for mt in (502, 504, 515, 522, 525): data.reactions[mt] = PhotonReaction.from_ace(ace, mt) - # Get heating cross sections [eV*barn] from factors [eV per collision] - # by multiplying with total xs (sum of (502, 504, 515, 522)) + # Get heating cross sections [eV-barn] from factors [eV per collision] + # by multiplying with total xs data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in (502, 504, 515, 522)]) @@ -849,7 +849,6 @@ class IncidentPhoton(EqualityMixin): else: brem_group.create_dataset(key, data=value) - def _add_bremsstrahlung(self): """Add the data used in the thick-target bremsstrahlung approximation @@ -915,20 +914,20 @@ class IncidentPhoton(EqualityMixin): self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number]) def _compute_heating(self): - """Compute heating cross sections (KERMA) + r"""Compute heating cross sections (KERMA) Photon energy is deposited as energy loss in three reactions: incoherent scattering, pair production and photoelectric effect. The point-wise heating cross section is calculated as: .. math:: - \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PP, PE \right\} + \sigma_{Hx}(E) &= (E - \overline{E}_x(E)) \cdot \sigma_x(E), x \in \left\{I, PP, PE \right\} - \overline{E}_I (E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} + \overline{E}_I(E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV - \overline{E}_{PE} &= E_{fluorescent photons} + \overline{E}_{PE} &= E(\text{fluorescent photons}) The differential cross section representation for incoherent scattering can be found in the theory manual. From acc69ea6cb6891076fcc9424b4f1f7a0e087d1a3 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 12 Apr 2019 20:38:36 -0400 Subject: [PATCH 36/61] updated tests to include photon heating --- .../photon_production/inputs_true.dat | 9 ++++-- .../photon_production/results_true.dat | 18 +++++++++++ .../photon_production/test.py | 31 ++++++++++--------- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 1d43c8913..f82b4a6a7 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -49,17 +49,20 @@ 2 - total + Al27 total + total heating tracklength 2 - total + Al27 total + total heating collision 2 - total + Al27 total + total heating analog diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 7df7e5dbc..a7725576c 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -2,11 +2,29 @@ tally 1: 9.403000E-01 8.841641E-01 tally 2: +1.819886E-01 +3.311985E-02 +1.960159E+05 +3.842224E+10 8.281718E-01 6.858685E-01 +1.960159E+05 +3.842224E+10 tally 3: +1.799069E-01 +3.236650E-02 +1.945225E+05 +3.783901E+10 8.242000E-01 6.793056E-01 +1.945225E+05 +3.783901E+10 tally 4: +2.308000E-01 +5.326864E-02 +1.988563E+05 +3.954384E+10 8.242000E-01 6.793056E-01 +1.988612E+05 +3.954578E+10 diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 47215fff3..c6a1fc390 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -51,20 +51,23 @@ class SourceTestHarness(PyAPITestHarness): current_tally = openmc.Tally() current_tally.filters = [surface_filter, particle_filter] current_tally.scores = ['current'] - total_tally_tracklength = openmc.Tally() - total_tally_tracklength.filters = [particle_filter] - total_tally_tracklength.scores = ['total'] - total_tally_tracklength.estimator = 'tracklength' - total_tally_collision = openmc.Tally() - total_tally_collision.filters = [particle_filter] - total_tally_collision.scores = ['total'] - total_tally_collision.estimator = 'collision' - total_tally_analog = openmc.Tally() - total_tally_analog.filters = [particle_filter] - total_tally_analog.scores = ['total'] - total_tally_analog.estimator = 'analog' - tallies = openmc.Tallies([current_tally, total_tally_tracklength, - total_tally_collision, total_tally_analog]) + tally_tracklength = openmc.Tally() + tally_tracklength.filters = [particle_filter] + tally_tracklength.scores = ['total', 'heating'] + tally_tracklength.nuclides = ['Al27', 'total'] + tally_tracklength.estimator = 'tracklength' + tally_collision = openmc.Tally() + tally_collision.filters = [particle_filter] + tally_collision.scores = ['total', 'heating'] + tally_collision.nuclides = ['Al27', 'total'] + tally_collision.estimator = 'collision' + tally_analog = openmc.Tally() + tally_analog.filters = [particle_filter] + tally_analog.scores = ['total', 'heating'] + tally_analog.nuclides = ['Al27', 'total'] + tally_analog.estimator = 'analog' + tallies = openmc.Tallies([current_tally, tally_tracklength, + tally_collision, tally_analog]) tallies.export_to_xml() def _get_results(self): From 435b21eecb92d759f212954942f158404f0b1cea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 17 Apr 2019 10:40:11 -0500 Subject: [PATCH 37/61] Use fixed-tolerance Gaussian quadrature for integrals --- openmc/data/photon.py | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index cd1d8c05c..1de84721e 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -3,13 +3,14 @@ from collections.abc import Mapping, Callable from copy import deepcopy from io import StringIO from numbers import Integral, Real +from math import pi, sqrt import os import h5py import numpy as np import pandas as pd from scipy.interpolate import CubicSpline -from scipy.integrate import quad +from scipy.integrate import quadrature from warnings import warn from openmc.mixin import EqualityMixin @@ -27,7 +28,7 @@ PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant CM_PER_ANGSTROM = 1.0e-8 # classical electron radius in cm -R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * np.pi * FINE_STRUCTURE * MASS_ELECTRON_EV) +R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * pi * FINE_STRUCTURE * MASS_ELECTRON_EV) # Electron subshell labels _SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', @@ -948,24 +949,26 @@ class IncidentPhoton(EqualityMixin): # Incoherent scattering if 504 in self: rx = self[504] - def _dsigma_dmu(mu, E): + + def dsigma_dmu(mu, E): k = E / MASS_ELECTRON_EV - kout = k / (1.0 + k * (1.0 - mu)) - x = E * np.sqrt(0.5 * (1.0 - mu)) / PLANCK_C - return np.pi * R0**2 * (kout/k)**2 * (kout/k + k/kout + + krat = 1.0 / (1.0 + k * (1.0 - mu)) + x = E * sqrt(0.5 * (1.0 - mu)) / PLANCK_C + return pi * R0*R0 * krat*krat * (krat + 1/krat + mu*mu - 1.0) * rx.scattering_factor(x) - def _eout_dsigma_dmu(mu, E): - Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1 - mu)) + + def eout_dsigma_dmu(mu, E): + Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1.0 - mu)) return Eout * _dsigma_dmu(mu, E) - def _eout_average(E): - #TODO: the integrals are currently time consuming - integral_sigma = quad(_dsigma_dmu, -1.0, 1.0, - args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] - integral_sigma_e = quad(_eout_dsigma_dmu, -1.0, 1.0, - args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] + + def eout_average(E): + integral_sigma = quadrature(_dsigma_dmu, -1.0, 1.0, + (E,), rtol=1e-3, vec_func=False)[0] + integral_sigma_e = quadrature(_eout_dsigma_dmu, -1.0, 1.0, + (E,), rtol=1e-3, vec_func=False)[0] return integral_sigma_e / integral_sigma - e_out = np.vectorize(lambda x: _eout_average(x))(energy) + e_out = np.vectorize(_eout_average)(energy) heating_xs += (energy - e_out) * rx.xs(energy) # Pair production, electron field From f90e8dd3dccceb20fea9f6d299afd0378be1dccb Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 17 Apr 2019 14:11:55 -0400 Subject: [PATCH 38/61] address @paulromano's comments --- include/openmc/string_utils.h | 2 +- openmc/data/photon.py | 90 ++++++++++++++++++----------------- src/string_utils.cpp | 2 +- 3 files changed, 49 insertions(+), 45 deletions(-) diff --git a/include/openmc/string_utils.h b/include/openmc/string_utils.h index db7635d0a..a60fd06a8 100644 --- a/include/openmc/string_utils.h +++ b/include/openmc/string_utils.h @@ -10,7 +10,7 @@ std::string& strtrim(std::string& s); char* strtrim(char* c_str); -std::string to_element(std::string const& name); +std::string to_element(const std::string& name); void to_lower(std::string& str); diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 1de84721e..da576b5c8 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -10,8 +10,7 @@ import h5py import numpy as np import pandas as pd from scipy.interpolate import CubicSpline -from scipy.integrate import quadrature -from warnings import warn +from scipy.integrate import quad from openmc.mixin import EqualityMixin import openmc.checkvalue as cv @@ -394,6 +393,17 @@ class AtomicRelaxation(EqualityMixin): def energy_fluorescence(self, shell): """Compute expected energy of fluorescent photons for the shell + + Parameters + ---------- + shell : str + The subshell to compute + + Returns + ------- + float + Energy of fluorescent photons + """ if shell not in self.binding_energy: @@ -402,31 +412,28 @@ class AtomicRelaxation(EqualityMixin): if shell in self._e_fluorescence: # Already computed return self._e_fluorescence[shell] + e = 0.0 + if shell not in self.transitions or self.transitions[shell].empty: + e = self.binding_energy[shell] else: - e = 0.0 - if shell not in self.transitions or self.transitions[shell].empty: - e = self.binding_energy[shell] - else: - df = self.transitions[shell] - for index, row in df.iterrows(): - e_row = 0.0 - primary = row['secondary'] - secondary = row['tertiary'] - if secondary is None: - # Fluorescent photon release in radiative transition - e_row += row['energy (eV)'] - else: - # Fill pole left by auger eletron - e_row += self.energy_fluorescence(secondary) + df = self.transitions[shell] + for primary, secondary, energy, prob in df.itertuples(index=False): + e_row = 0.0 + if secondary is None: + # Fluorescent photon release in radiative transition + e_row += energy + else: + # Fill the hole left by auger electron + e_row += self.energy_fluorescence(secondary) - # Fill photonelectron pole - e_row += self.energy_fluorescence(primary) + # Fill the photoelectron hole + e_row += self.energy_fluorescence(primary) - # Expected fluorescent photon energy - e += e_row * row['probability'] + # Expected fluorescent photon energy + e += e_row * prob - self._e_fluorescence[shell] = e - return e + self._e_fluorescence[shell] = e + return e class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -917,31 +924,28 @@ class IncidentPhoton(EqualityMixin): def _compute_heating(self): r"""Compute heating cross sections (KERMA) - Photon energy is deposited as energy loss in three reactions: - incoherent scattering, pair production and photoelectric effect. - The point-wise heating cross section is calculated as: + Photon energy is deposited as energy loss in three reactions: + incoherent scattering, pair production and photoelectric effect. + The point-wise heating cross section is calculated as: - .. math:: - \sigma_{Hx}(E) &= (E - \overline{E}_x(E)) \cdot \sigma_x(E), x \in \left\{I, PP, PE \right\} + .. math:: + \sigma_{Hx}(E) &= (E - \overline{E}_x(E)) \cdot \sigma_x(E), x \in \left\{I, PP, PE \right\} - \overline{E}_I(E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} + \overline{E}_I(E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} - \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV + \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV - \overline{E}_{PE} &= E(\text{fluorescent photons}) + \overline{E}_{PE} &= E(\text{fluorescent photons}) - The differential cross section representation for incoherent - scattering can be found in the theory manual. + The differential cross section representation for incoherent + scattering can be found in the theory manual. """ # Determine a union energy grid energy = np.array([]) for mt in (504, 515, 517, 522): - if mt not in self: - warn('Cross sections for MT={} do not exist. The total heating ' - 'to be calculated is probably incorrect'.format(mt)) - else: + if mt in self: energy = np.union1d(energy, self[mt].xs.x) heating_xs = np.zeros_like(energy) @@ -959,16 +963,16 @@ class IncidentPhoton(EqualityMixin): def eout_dsigma_dmu(mu, E): Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1.0 - mu)) - return Eout * _dsigma_dmu(mu, E) + return Eout * dsigma_dmu(mu, E) def eout_average(E): - integral_sigma = quadrature(_dsigma_dmu, -1.0, 1.0, - (E,), rtol=1e-3, vec_func=False)[0] - integral_sigma_e = quadrature(_eout_dsigma_dmu, -1.0, 1.0, - (E,), rtol=1e-3, vec_func=False)[0] + integral_sigma = quad(dsigma_dmu, -1.0, 1.0, + args=(E,), epsabs=0.0, epsrel=1e-3)[0] + integral_sigma_e = quad(eout_dsigma_dmu, -1.0, 1.0, + args=(E,), epsabs=0.0, epsrel=1e-3)[0] return integral_sigma_e / integral_sigma - e_out = np.vectorize(_eout_average)(energy) + e_out = np.vectorize(eout_average)(energy) heating_xs += (energy - e_out) * rx.xs(energy) # Pair production, electron field diff --git a/src/string_utils.cpp b/src/string_utils.cpp index 40a20f232..614e38767 100644 --- a/src/string_utils.cpp +++ b/src/string_utils.cpp @@ -26,7 +26,7 @@ char* strtrim(char* c_str) } -std::string to_element(std::string const& name) { +std::string to_element(const std::string& name) { int pos = name.find_first_of("0123456789"); return name.substr(0, pos); } From 8f48ecf7dec4ee6a36b431b965ce23c91429f096 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 18 Apr 2019 08:22:40 -0500 Subject: [PATCH 39/61] Add method for generating settings from an XML file --- openmc/mesh.py | 37 +++++ openmc/settings.py | 255 +++++++++++++++++++++++++++++++++-- openmc/source.py | 76 ++++++++++- openmc/stats/multivariate.py | 130 ++++++++++++++++++ openmc/stats/univariate.py | 155 +++++++++++++++++++++ 5 files changed, 641 insertions(+), 12 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index ba07ec178..a30b333fb 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -246,6 +246,43 @@ class Mesh(IDManagerMixin): return element + @classmethod + def from_xml_element(cls, elem): + """Generate mesh from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.Mesh + Mesh generated from XML element + + """ + mesh_id = int(elem.get('id')) + mesh = cls(mesh_id) + mesh.type = elem.get('type') + + dimension = elem.findtext('dimension') + if dimension is not None: + mesh.dimension = [int(x) for x in dimension.split()] + + lower_left = elem.findtext('lower_left') + if lower_left is not None: + mesh.lower_left = [float(x) for x in lower_left.split()] + + upper_right = elem.findtext('upper_right') + if upper_right is not None: + mesh.upper_right = [float(x) for x in upper_right.split()] + + width = elem.findtext('width') + if width is not None: + mesh.width = [float(x) for x in width.split()] + + return mesh + def build_cells(self, bc=['reflective'] * 6): """Generates a lattice of universes with the same dimensionality as the mesh object. The individual cells/universes produced diff --git a/openmc/settings.py b/openmc/settings.py index f348bc6cb..a7ec08679 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -693,29 +693,29 @@ class Settings(object): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode - def _create_batches_subelement(self, run_mode_element): + def _create_batches_subelement(self, root): if self._batches is not None: - element = ET.SubElement(run_mode_element, "batches") + element = ET.SubElement(root, "batches") element.text = str(self._batches) - def _create_generations_per_batch_subelement(self, run_mode_element): + def _create_generations_per_batch_subelement(self, root): if self._generations_per_batch is not None: - element = ET.SubElement(run_mode_element, "generations_per_batch") + element = ET.SubElement(root, "generations_per_batch") element.text = str(self._generations_per_batch) - def _create_inactive_subelement(self, run_mode_element): + def _create_inactive_subelement(self, root): if self._inactive is not None: - element = ET.SubElement(run_mode_element, "inactive") + element = ET.SubElement(root, "inactive") element.text = str(self._inactive) - def _create_particles_subelement(self, run_mode_element): + def _create_particles_subelement(self, root): if self._particles is not None: - element = ET.SubElement(run_mode_element, "particles") + element = ET.SubElement(root, "particles") element.text = str(self._particles) - def _create_keff_trigger_subelement(self, run_mode_element): + def _create_keff_trigger_subelement(self, root): if self._keff_trigger is not None: - element = ET.SubElement(run_mode_element, "keff_trigger") + element = ET.SubElement(root, "keff_trigger") for key in self._keff_trigger: subelement = ET.SubElement(element, key) @@ -985,3 +985,238 @@ class Settings(object): # Write the XML Tree to the settings.xml file tree = ET.ElementTree(root_element) tree.write(str(p), xml_declaration=True, encoding='utf-8') + + @classmethod + def from_xml(cls, path='settings.xml'): + """Generate settings from XML file + + Parameters + ---------- + path : str, optional + Path to settings XML file + + Returns + ------- + openmc.Settings + Settings object + + """ + tree = ET.parse(path) + root = tree.getroot() + + settings = cls() + + # Get the run mode + elem = root.find('run_mode') + if elem is not None: + settings.run_mode = elem.text + + # Get number of particles + elem = root.find('particles') + if elem is not None: + settings.particles = int(elem.text) + + # Get number of batches + elem = root.find('batches') + if elem is not None: + settings.batches = int(elem.text) + + # Get number of inactive batches + elem = root.find('inactive') + if elem is not None: + settings.inactive = int(elem.text) + + # Get number of generations per batch + elem = root.find('generations_per_batch') + if elem is not None: + settings.generations_per_batch = int(elem.text) + + # Get keff trigger + elem = root.find('keff_trigger') + if elem is not None: + trigger = elem.findtext('type') + threshold = float(elem.findtext('threshold')) + settings.keff_trigger = {'type': trigger, 'threshold': threshold} + + # Get the source + for elem in root.findall('source'): + settings.source.append(Source.from_xml_element(elem)) + + # Get the output + elem = root.find('output') + if elem is not None: + settings.output = {} + for entry in elem: + key = entry.tag + if key in ('summary', 'tallies'): + value = entry.text == 'true' + else: + value = entry.text + settings.output[key] = value + + # Get the statepoint + elem = root.find('state_point') + if elem is not None: + batches = elem.findtext('batches') + if batches is not None: + settings.statepoint['batches'] = [int(x) for x in batches.split()] + + # Get the sourcepoint + elem = root.find('source_point') + if elem is not None: + for entry in elem: + key = entry.tag + if key in ('separate', 'write', 'overwrite'): + value = entry.text == 'true' + else: + value = [int(x) for x in entry.text.split()] + settings.sourcepoint[key] = value + + # Get confidence intervals + elem = root.find('confidence_intervals') + if elem is not None: + settings.confidence_intervals = elem.text == 'true' + + # Get electron treatment + elem = root.find('electron_treatment') + if elem is not None: + settings.electron_treatment = elem.text + + # Get energy mode + elem = root.find('energy_mode') + if elem is not None: + settings.energy_mode = elem.text + + # Get max order + elem = root.find('max_order') + if elem is not None: + settings.max_order = int(elem.text) + + # Get photon transport + elem = root.find('photon_transport') + if elem is not None: + settings.photon_transport = elem.text == 'true' + + # Get probability tables + elem = root.find('ptables') + if elem is not None: + settings.ptables = elem.text == 'true' + + # Get seed + elem = root.find('seed') + if elem is not None: + settings.seed = int(elem.text) + + # Get survival biasing + elem = root.find('survival_biasing') + if elem is not None: + settings.survival_biasing = elem.text == 'true' + + # Get cutoff + elem = root.find('cutoff') + if elem is not None: + settings.cutoff = {x.tag: float(x.text) for x in elem} + + # Get entropy mesh + elem = root.find('entropy_mesh') + if elem is not None: + settings.entropy_mesh = Mesh.from_xml_element(elem) + + # Get trigger + elem = root.find('trigger') + if elem is not None: + active = elem.find('active') + settings.trigger_active = active.text == 'true' + max_batches = elem.find('max_batches') + if max_batches is not None: + settings.trigger_max_batches = int(max_batches.text) + batch_interval = elem.find('batch_interval') + if batch_interval is not None: + settings.trigger_batch_interval = int(batch_interval.text) + + # Get no reduce + elem = root.find('no_reduce') + if elem is not None: + settings.no_reduce = elem.text == 'true' + + # Get verbosity + elem = root.find('verbosity') + if elem is not None: + settings.verbosity = int(elem.text) + + # Get tabular legendre + elem = root.find('tabular_legendre') + if elem is not None: + enable = elem.findtext('eneable') + settings.tabular_legendre['enable'] = enable == 'true' + num_points = elem.findtext('num_points') + if num_points is not None: + settings.tabular_legendre['num_points'] = int(num_points) + + # Get temperature + elem = root.findtext('temperature_default') + if elem is not None: + settings.temperature['default'] = float(elem) + elem = root.findtext('temperature_tolerance') + if elem is not None: + settings.temperature['tolerance'] = float(elem) + elem = root.findtext('temperature_method') + if elem is not None: + settings.temperature['method'] = elem + elem = root.findtext('temperature_range') + if elem is not None: + settings.temperature['range'] = [float(x) for x in elem.split()] + elem = root.findtext('temperature_multipole') + if elem is not None: + settings.temperature['multipole'] = elem == 'true' + + # Get trace + elem = root.find('trace') + if elem is not None: + settings.trace = [int(x) for x in elem.text.split()] + + # Get track + elem = root.find('track') + if elem is not None: + settings.track = [int(x) for x in elem.text.split()] + + # Get UFS mesh + elem = root.find('ufs_mesh') + if elem is not None: + settings.ufs_mesh = Mesh.from_xml_element(elem) + + # Get resonance scattering + elem = root.find('resonance_scattering') + if elem is not None: + for entry in elem: + key = entry.tag + if key == 'enable': + value = entry.text == 'true' + elif key == 'method': + value = entry.text + elif key == 'energy_min': + value = float(entry.text) + elif key == 'energy_max': + value = float(entry.text) + elif key == 'nuclides': + value = entry.text.split() + settings.resonance_scattering[key] = value + + # TODO: Get volume calculations + + # Get fission neutrons + elem = root.find('create_fission_neutrons') + if elem is not None: + settings.create_fission_neutrons = elem.text == 'true' + + # Get log grid bins + elem = root.find('log_grid_bins') + if elem is not None: + settings.log_grid_bins = int(elem.text) + + # Get dagmc + elem = root.find('dagmc') + if elem is not None: + settings.dagmc = elem.text == 'true' + + return settings diff --git a/openmc/source.py b/openmc/source.py index 6ec882ca6..e278d8089 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -2,8 +2,11 @@ from numbers import Real import sys from xml.etree import ElementTree as ET -from openmc.stats.univariate import Univariate -from openmc.stats.multivariate import UnitSphere, Spatial +from openmc.stats.univariate import (Univariate, Discrete, Uniform, Maxwell, + Watt, Normal, Muir, Tabular) +from openmc.stats.multivariate import (UnitSphere, Spatial, PolarAzimuthal, + Isotropic, Monodirectional, Box, Point, + CartesianIndependent) import openmc.checkvalue as cv @@ -137,3 +140,72 @@ class Source(object): if self.energy is not None: element.append(self.energy.to_xml_element('energy')) return element + + @classmethod + def from_xml_element(cls, elem): + """Generate source from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.Source + Source generated from XML element + + """ + source = cls() + + strength = elem.find('strength') + if strength is not None: + source.strength = float(strength.text) + + particle = elem.find('particle') + if particle is not None: + source.particle = particle.text + + filename = elem.find('file') + if filename is not None: + source.file = filename.text + + space = elem.find('space') + if space is not None: + space_type = space.get('type') + if space_type == 'cartesian': + source.space = CartesianIndependent.from_xml_element(space) + elif space_type == 'box' or space_type == 'fission': + source.space = Box.from_xml_element(space) + elif space_type == 'point': + source.space = Point.from_xml_element(space) + + angle = elem.find('angle') + if angle is not None: + angle_type = angle.get('type') + if angle_type == 'mu-phi': + source.angle = PolarAzimuthal.from_xml_element(angle) + elif angle_type == 'isotropic': + source.angle = Isotropic.from_xml_element(angle) + elif angle_type == 'monodirectional': + source.angle = Monodirectional.from_xml_element(angle) + + energy = elem.find('energy') + if energy is not None: + energy_type = energy.get('type') + if energy_type == 'discrete': + source.energy = Discrete.from_xml_element(energy) + elif energy_type == 'uniform': + source.energy = Uniform.from_xml_element(energy) + elif energy_type == 'maxwell': + source.energy = Maxwell.from_xml_element(energy) + elif energy_type == 'watt': + source.energy = Watt.from_xml_element(energy) + elif energy_type == 'normal': + source.energy = Normal.from_xml_element(energy) + elif energy_type == 'muir': + source.energy = Muir.from_xml_element(energy) + elif energy_type == 'tabular': + source.energy = Tabular.from_xml_element(energy) + + return source diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index ac788c344..61122ada2 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -47,6 +47,11 @@ class UnitSphere(metaclass=ABCMeta): def to_xml_element(self): return '' + @classmethod + @abstractmethod + def from_xml_element(cls, elem): + pass + class PolarAzimuthal(UnitSphere): """Angular distribution represented by polar and azimuthal angles @@ -121,6 +126,29 @@ class PolarAzimuthal(UnitSphere): element.append(self.phi.to_xml_element('phi')) return element + @classmethod + def from_xml_element(cls, elem): + """Generate angular distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.PolarAzimuthal + Angular distribution generated from XML element + + """ + mu_phi = cls() + params = elem.findtext('parameters') + if params is not None: + mu_phi.reference_uvw = [float(x) for x in params.split()] + mu_phi.mu = openmc.stats.Univariate.from_xml_element(elem.find('mu')) + mu_phi.phi = openmc.stats.Univariate.from_xml_element(elem.find('phi')) + return mu_phi + class Isotropic(UnitSphere): """Isotropic angular distribution. @@ -143,6 +171,23 @@ class Isotropic(UnitSphere): element.set("type", "isotropic") return element + @classmethod + def from_xml_element(cls, elem): + """Generate isotropic distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Isotropic + Isotropic distribution generated from XML element + + """ + return cls() + class Monodirectional(UnitSphere): """Monodirectional angular distribution. @@ -178,6 +223,27 @@ class Monodirectional(UnitSphere): element.set("reference_uvw", ' '.join(map(str, self.reference_uvw))) return element + @classmethod + def from_xml_element(cls, elem): + """Generate monodirectional distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Monodirectional + Monodirectional distribution generated from XML element + + """ + monodirectional = cls() + params = elem.findtext('parameters') + if params is not None: + monodirectional.reference_uvw = [float(x) for x in params.split()] + return monodirectional + class Spatial(metaclass=ABCMeta): """Distribution of locations in three-dimensional Euclidean space. @@ -193,6 +259,11 @@ class Spatial(metaclass=ABCMeta): def to_xml_element(self): return '' + @classmethod + @abstractmethod + def from_xml_element(cls, elem): + pass + class CartesianIndependent(Spatial): """Spatial distribution with independent x, y, and z distributions. @@ -270,6 +341,26 @@ class CartesianIndependent(Spatial): element.append(self.z.to_xml_element('z')) return element + @classmethod + def from_xml_element(cls, elem): + """Generate spatial distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.CartesianIndependent + Spatial distribution generated from XML element + + """ + x = openmc.stats.Univariate.from_xml_element(elem.find('x')) + y = openmc.stats.Univariate.from_xml_element(elem.find('y')) + z = openmc.stats.Univariate.from_xml_element(elem.find('z')) + return cls(x, y, z) + class Box(Spatial): """Uniform distribution of coordinates in a rectangular cuboid. @@ -351,6 +442,27 @@ class Box(Spatial): ' '.join(map(str, self.upper_right)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate box distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Box + Box distribution generated from XML element + + """ + only_fissionable = elem.get('type') == 'fission' + params = [float(x) for x in elem.findtext('parameters').split()] + lower_left = params[:len(params)//2] + upper_right = paramx[len(params)//2:] + return cls(lower_left, upper_right, only_fissionable) + class Point(Spatial): """Delta function in three dimensions. @@ -398,3 +510,21 @@ class Point(Spatial): params = ET.SubElement(element, "parameters") params.text = ' '.join(map(str, self.xyz)) return element + + @classmethod + def from_xml_element(cls, elem): + """Generate point distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Point + Point distribution generated from XML element + + """ + xyz = [float(x) for x in elem.findtext('parameters').split()] + return cls(xyz) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index e61f216ef..6a83cc780 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -32,6 +32,11 @@ class Univariate(EqualityMixin, metaclass=ABCMeta): def __len__(self): return 0 + @classmethod + @abstractmethod + def from_xml_element(cls, elem): + pass + class Discrete(Univariate): """Distribution characterized by a probability mass function. @@ -110,6 +115,26 @@ class Discrete(Univariate): return element + @classmethod + def from_xml_element(cls, elem): + """Generate discrete distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Discrete + Discrete distribution generated from XML element + + """ + params = [float(x) for x in elem.findtext('parameters').split()] + x = params[:len(params)//2] + p = params[len(params)//2:] + return cls(x, p) + class Uniform(Univariate): """Distribution with constant probability over a finite interval [a,b] @@ -181,6 +206,26 @@ class Uniform(Univariate): element.set("parameters", '{} {}'.format(self.a, self.b)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate uniform distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Uniform + Uniform distribution generated from XML element + + """ + params = elem.findtext('parameters').split() + a = float(params[0]) + b = float(params[1]) + return cls(a, b) + class Maxwell(Univariate): """Maxwellian distribution in energy. @@ -237,6 +282,24 @@ class Maxwell(Univariate): element.set("parameters", str(self.theta)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate Maxwellian distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Maxwell + Maxwellian distribution generated from XML element + + """ + theta = float(elem.findtext('parameters')) + return cls(theta) + class Watt(Univariate): r"""Watt fission energy spectrum. @@ -308,6 +371,27 @@ class Watt(Univariate): element.set("parameters", '{} {}'.format(self.a, self.b)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate Watt distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Watt + Watt distribution generated from XML element + + """ + params = elem.findtext('parameters').split() + a = float(params[0]) + b = float(params[1]) + return watt(a, b) + + class Normal(Univariate): r"""Normally distributed sampling. @@ -377,6 +461,27 @@ class Normal(Univariate): element.set("parameters", '{} {}'.format(self.mean_value, self.std_dev)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate Normal distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Normal + Normal distribution generated from XML element + + """ + params = elem.findtext('parameters').split() + mean_value = float(params[0]) + std_dev = float(params[1]) + return cls(mean_value, std_dev) + + class Muir(Univariate): """Muir energy spectrum. @@ -465,6 +570,27 @@ class Muir(Univariate): element.set("parameters", '{} {} {}'.format(self._e0, self._m_rat, self._kt)) return element + @classmethod + def from_xml_element(cls, elem): + """Generate Muir distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Muir + Muir distribution generated from XML element + + """ + params = elem.findtext('parameters').split() + e0 = float(params[0]) + m_rat = float(params[1]) + kt = float(params[2]) + return muir(e0, m_rat, kt) + class Tabular(Univariate): """Piecewise continuous probability distribution. @@ -561,6 +687,27 @@ class Tabular(Univariate): return element + @classmethod + def from_xml_element(cls, elem): + """Generate tabular distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Tabular + Tabular distribution generated from XML element + + """ + interpolation = elem.get('interpolation') + params = [float(x) for x in elem.findtext('parameters').split()] + x = params[:len(params)//2] + p = paramx[len(params)//2:] + return cls(x, p, interpolation) + class Legendre(Univariate): r"""Probability density given by a Legendre polynomial expansion @@ -607,6 +754,10 @@ class Legendre(Univariate): def to_xml_element(self, element_name): raise NotImplementedError + @classmethod + def from_xml_element(cls, elem): + raise NotImplementedError + class Mixture(Univariate): """Probability distribution characterized by a mixture of random variables. @@ -660,3 +811,7 @@ class Mixture(Univariate): def to_xml_element(self, element_name): raise NotImplementedError + + @classmethod + def from_xml_element(cls, elem): + raise NotImplementedError From 4d77f6fc943885d87692d5ce6f874af2635c9379 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 Apr 2019 14:00:47 -0400 Subject: [PATCH 40/61] Added a function to do interpolation for scalar input --- openmc/data/function.py | 54 +++++++++++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index f167a3c23..c40e2a699 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -3,6 +3,7 @@ from collections.abc import Iterable, Callable from functools import reduce from itertools import zip_longest from numbers import Real, Integral +from math import exp, log import numpy as np @@ -153,13 +154,11 @@ class Tabulated1D(Function1D): self.y = np.asarray(y) def __call__(self, x): - # Check if input is array or scalar - if isinstance(x, Iterable): - iterable = True - x = np.array(x) - else: - iterable = False - x = np.array([x], dtype=float) + # Check if input is scalar + if not isinstance(x, Iterable): + return self._interpolate_scalar(x) + + x = np.array(x) # Create output array y = np.zeros_like(x) @@ -208,7 +207,46 @@ class Tabulated1D(Function1D): y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0] y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1] - return y if iterable else y[0] + return y + + def _interpolate_scalar(self, x): + if x <= self._x[0]: + return self._y[0] + elif x >= self._x[-1]: + return self._y[-1] + + # Get the index for interpolation + idx = np.searchsorted(self._x, x, side='right') - 1 + + # Loop over interpolation regions + for b, p in zip(self.breakpoints, self.interpolation): + if idx < b - 1: + break + + xi = self._x[idx] # low edge of the corresponding bin + xi1 = self._x[idx + 1] # high edge of the corresponding bin + yi = self._y[idx] + yi1 = self._y[idx + 1] + + if p == 1: + # Histogram + return yi + + elif p == 2: + # Linear-linear + return yi + (x - xi)/(xi1 - xi)*(yi1 - yi) + + elif p == 3: + # Linear-log + return yi + log(x/xi)/log(xi1/xi)*(yi1 - yi) + + elif p == 4: + # Log-linear + return yi*exp((x - xi)/(xi1 - xi)*log(yi1/yi)) + + elif p == 5: + # Log-log + return yi*exp(log(x/xi)/log(xi1/xi)*log(yi1/yi)) def __len__(self): return len(self.x) From 21c693e47326677dd844e6ec27efdc3f86ab38c6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 Apr 2019 17:36:43 -0400 Subject: [PATCH 41/61] update test cross section library for travis --- tools/ci/download-xs.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index ce4d12c1b..38a5bc469 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -3,7 +3,7 @@ set -ex # Download HDF5 data if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget -q -O - https://anl.box.com/shared/static/pzutl4i2717yypv12l78l7fn5nmg6grs.xz | tar -C $HOME -xJ + wget -q -O - https://anl.box.com/shared/static/u1g3n8iai0u1n5f6ev3pg2j3ff941bqa.xz | tar -C $HOME -xJ fi # Download ENDF/B-VII.1 distribution From 489a03c648bb7914b69c2813f996ac37c1f65b77 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 25 Apr 2019 20:07:53 -0500 Subject: [PATCH 42/61] Get value from either attribute or element when reading settings from XML --- openmc/mesh.py | 13 +- openmc/settings.py | 468 +++++++++++++++++++---------------- openmc/source.py | 19 +- openmc/stats/multivariate.py | 13 +- openmc/stats/univariate.py | 17 +- 5 files changed, 289 insertions(+), 241 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index a30b333fb..fb59053ae 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -7,6 +7,7 @@ import numpy as np import openmc.checkvalue as cv import openmc +from openmc._xml import get_text from openmc.mixin import EqualityMixin, IDManagerMixin @@ -261,23 +262,23 @@ class Mesh(IDManagerMixin): Mesh generated from XML element """ - mesh_id = int(elem.get('id')) + mesh_id = int(get_text(elem, 'id')) mesh = cls(mesh_id) - mesh.type = elem.get('type') + mesh.type = get_text(elem, 'type') - dimension = elem.findtext('dimension') + dimension = get_text(elem, 'dimension') if dimension is not None: mesh.dimension = [int(x) for x in dimension.split()] - lower_left = elem.findtext('lower_left') + lower_left = get_text(elem, 'lower_left') if lower_left is not None: mesh.lower_left = [float(x) for x in lower_left.split()] - upper_right = elem.findtext('upper_right') + upper_right = get_text(elem, 'upper_right') if upper_right is not None: mesh.upper_right = [float(x) for x in upper_right.split()] - width = elem.findtext('width') + width = get_text(elem, 'width') if width is not None: mesh.width = [float(x) for x in width.split()] diff --git a/openmc/settings.py b/openmc/settings.py index a7ec08679..39312899b 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -7,7 +7,7 @@ import sys import numpy as np -from openmc._xml import clean_indentation +from openmc._xml import clean_indentation, get_text import openmc.checkvalue as cv from openmc import VolumeCalculation, Source, Mesh @@ -927,6 +927,228 @@ class Settings(object): elem = ET.SubElement(root, "dagmc") elem.text = str(self._dagmc).lower() + def _eigenvalue_from_xml_element(self, root): + elem = root.find('eigenvalue') + if elem is not None: + self._run_mode_from_xml_element(elem) + self._particles_from_xml_element(elem) + self._batches_from_xml_element(elem) + self._inactive_from_xml_element(elem) + self._generations_per_batch_from_xml_element(elem) + + def _run_mode_from_xml_element(self, root): + 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') + if text is not None: + self.particles = int(text) + + def _batches_from_xml_element(self, root): + 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') + if text is not None: + self.inactive = int(text) + + def _generations_per_batch_from_xml_element(self, root): + 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') + if elem is not None: + 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'): + self.source.append(Source.from_xml_element(elem)) + + def _output_from_xml_element(self, root): + elem = root.find('output') + if elem is not None: + self.output = {} + for key in ('summary', 'tallies', 'path'): + value = get_text(elem, key) + if value is not None: + if key in ('summary', 'tallies'): + value = value == 'true' + self.output[key] = value + + def _statepoint_from_xml_element(self, root): + elem = root.find('state_point') + if elem is not None: + text = get_text(elem, 'batches') + if text is not None: + self.statepoint['batches'] = [int(x) for x in text.split()] + + def _sourcepoint_from_xml_element(self, root): + elem = root.find('source_point') + if elem is not None: + for key in ('separate', 'write', 'overwrite', 'batches'): + value = get_text(elem, key) + if value is not None: + if key in ('separate', 'write', 'overwrite'): + value = value == 'true' + else: + value = [int(x) for x in value.split()] + self.sourcepoint[key] = value + + def _confidence_intervals_from_xml_element(self, root): + text = get_text(root, 'confidence_intervals') + if text is not None: + self.confidence_intervals = text == 'true' + + def _electron_treatment_from_xml_element(self, root): + 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') + if text is not None: + self.energy_mode = text + + def _max_order_from_xml_element(self, root): + 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') + if text is not None: + self.photon_transport = text == 'true' + + def _ptables_from_xml_element(self, root): + text = get_text(root, 'ptables') + if text is not None: + self.ptables = text == 'true' + + def _seed_from_xml_element(self, root): + 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') + if text is not None: + self.survival_biasing = text == 'true' + + def _cutoff_from_xml_element(self, root): + 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'): + value = get_text(elem, key) + if value is not None: + self.cutoff[key] = float(value) + + def _entropy_mesh_from_xml_element(self, root): + elem = root.find('entropy_mesh') + if elem is not None: + self.entropy_mesh = Mesh.from_xml_element(elem) + + def _trigger_from_xml_element(self, root): + elem = root.find('trigger') + if elem is not None: + self.trigger_active = get_text(elem, 'active') == 'true' + text = get_text(elem, 'max_batches') + if text is not None: + self.trigger_max_batches = int(text) + 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') + if text is not None: + self.no_reduce = text == 'true' + + def _verbosity_from_xml_element(self, root): + 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') + if elem is not None: + text = get_text(elem, 'enable') + self.tabular_legendre['enable'] = text == 'true' + text = get_text(elem, 'num_points') + if text is not None: + self.tabular_legendre['num_points'] = int(text) + + def _temperature_from_xml_element(self, root): + text = get_text(root, 'temperature_default') + if text is not None: + 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') + if text is not None: + 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') + if text is not None: + self.temperature['multipole'] = text == 'true' + + def _trace_from_xml_element(self, root): + 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') + if text is not None: + self.track = [int(x) for x in text.split()] + + def _ufs_mesh_from_xml_element(self, root): + elem = root.find('ufs_mesh') + if elem is not None: + self.ufs_mesh = Mesh.from_xml_element(elem) + + def _resonance_scattering_from_xml_element(self, root): + elem = root.find('resonance_scattering') + if elem is not None: + 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 == 'true' + elif key in ('energy_min', 'energy_max'): + value = float(value) + 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') + if text is not None: + self.create_fission_neutrons = text == 'true' + + def _log_grid_bins_from_xml_element(self, root): + text = get_text(root, 'log_grid_bins') + if text is not None: + self.log_grid_bins = int(text) + + def _dagmc_from_xml_element(self, root): + text = get_text(root, 'dagmc') + if text is not None: + self.dagmc = text == 'true' + def export_to_xml(self, path='settings.xml'): """Export simulation settings to an XML file. @@ -1005,218 +1227,40 @@ class Settings(object): root = tree.getroot() settings = cls() - - # Get the run mode - elem = root.find('run_mode') - if elem is not None: - settings.run_mode = elem.text - - # Get number of particles - elem = root.find('particles') - if elem is not None: - settings.particles = int(elem.text) - - # Get number of batches - elem = root.find('batches') - if elem is not None: - settings.batches = int(elem.text) - - # Get number of inactive batches - elem = root.find('inactive') - if elem is not None: - settings.inactive = int(elem.text) - - # Get number of generations per batch - elem = root.find('generations_per_batch') - if elem is not None: - settings.generations_per_batch = int(elem.text) - - # Get keff trigger - elem = root.find('keff_trigger') - if elem is not None: - trigger = elem.findtext('type') - threshold = float(elem.findtext('threshold')) - settings.keff_trigger = {'type': trigger, 'threshold': threshold} - - # Get the source - for elem in root.findall('source'): - settings.source.append(Source.from_xml_element(elem)) - - # Get the output - elem = root.find('output') - if elem is not None: - settings.output = {} - for entry in elem: - key = entry.tag - if key in ('summary', 'tallies'): - value = entry.text == 'true' - else: - value = entry.text - settings.output[key] = value - - # Get the statepoint - elem = root.find('state_point') - if elem is not None: - batches = elem.findtext('batches') - if batches is not None: - settings.statepoint['batches'] = [int(x) for x in batches.split()] - - # Get the sourcepoint - elem = root.find('source_point') - if elem is not None: - for entry in elem: - key = entry.tag - if key in ('separate', 'write', 'overwrite'): - value = entry.text == 'true' - else: - value = [int(x) for x in entry.text.split()] - settings.sourcepoint[key] = value - - # Get confidence intervals - elem = root.find('confidence_intervals') - if elem is not None: - settings.confidence_intervals = elem.text == 'true' - - # Get electron treatment - elem = root.find('electron_treatment') - if elem is not None: - settings.electron_treatment = elem.text - - # Get energy mode - elem = root.find('energy_mode') - if elem is not None: - settings.energy_mode = elem.text - - # Get max order - elem = root.find('max_order') - if elem is not None: - settings.max_order = int(elem.text) - - # Get photon transport - elem = root.find('photon_transport') - if elem is not None: - settings.photon_transport = elem.text == 'true' - - # Get probability tables - elem = root.find('ptables') - if elem is not None: - settings.ptables = elem.text == 'true' - - # Get seed - elem = root.find('seed') - if elem is not None: - settings.seed = int(elem.text) - - # Get survival biasing - elem = root.find('survival_biasing') - if elem is not None: - settings.survival_biasing = elem.text == 'true' - - # Get cutoff - elem = root.find('cutoff') - if elem is not None: - settings.cutoff = {x.tag: float(x.text) for x in elem} - - # Get entropy mesh - elem = root.find('entropy_mesh') - if elem is not None: - settings.entropy_mesh = Mesh.from_xml_element(elem) - - # Get trigger - elem = root.find('trigger') - if elem is not None: - active = elem.find('active') - settings.trigger_active = active.text == 'true' - max_batches = elem.find('max_batches') - if max_batches is not None: - settings.trigger_max_batches = int(max_batches.text) - batch_interval = elem.find('batch_interval') - if batch_interval is not None: - settings.trigger_batch_interval = int(batch_interval.text) - - # Get no reduce - elem = root.find('no_reduce') - if elem is not None: - settings.no_reduce = elem.text == 'true' - - # Get verbosity - elem = root.find('verbosity') - if elem is not None: - settings.verbosity = int(elem.text) - - # Get tabular legendre - elem = root.find('tabular_legendre') - if elem is not None: - enable = elem.findtext('eneable') - settings.tabular_legendre['enable'] = enable == 'true' - num_points = elem.findtext('num_points') - if num_points is not None: - settings.tabular_legendre['num_points'] = int(num_points) - - # Get temperature - elem = root.findtext('temperature_default') - if elem is not None: - settings.temperature['default'] = float(elem) - elem = root.findtext('temperature_tolerance') - if elem is not None: - settings.temperature['tolerance'] = float(elem) - elem = root.findtext('temperature_method') - if elem is not None: - settings.temperature['method'] = elem - elem = root.findtext('temperature_range') - if elem is not None: - settings.temperature['range'] = [float(x) for x in elem.split()] - elem = root.findtext('temperature_multipole') - if elem is not None: - settings.temperature['multipole'] = elem == 'true' - - # Get trace - elem = root.find('trace') - if elem is not None: - settings.trace = [int(x) for x in elem.text.split()] - - # Get track - elem = root.find('track') - if elem is not None: - settings.track = [int(x) for x in elem.text.split()] - - # Get UFS mesh - elem = root.find('ufs_mesh') - if elem is not None: - settings.ufs_mesh = Mesh.from_xml_element(elem) - - # Get resonance scattering - elem = root.find('resonance_scattering') - if elem is not None: - for entry in elem: - key = entry.tag - if key == 'enable': - value = entry.text == 'true' - elif key == 'method': - value = entry.text - elif key == 'energy_min': - value = float(entry.text) - elif key == 'energy_max': - value = float(entry.text) - elif key == 'nuclides': - value = entry.text.split() - settings.resonance_scattering[key] = value + settings._eigenvalue_from_xml_element(root) + settings._run_mode_from_xml_element(root) + settings._particles_from_xml_element(root) + settings._batches_from_xml_element(root) + settings._inactive_from_xml_element(root) + settings._generations_per_batch_from_xml_element(root) + settings._keff_trigger_from_xml_element(root) + settings._source_from_xml_element(root) + settings._output_from_xml_element(root) + settings._statepoint_from_xml_element(root) + settings._sourcepoint_from_xml_element(root) + settings._confidence_intervals_from_xml_element(root) + settings._electron_treatment_from_xml_element(root) + settings._energy_mode_from_xml_element(root) + settings._max_order_from_xml_element(root) + settings._photon_transport_from_xml_element(root) + settings._ptables_from_xml_element(root) + settings._seed_from_xml_element(root) + settings._survival_biasing_from_xml_element(root) + settings._cutoff_from_xml_element(root) + settings._entropy_mesh_from_xml_element(root) + settings._trigger_from_xml_element(root) + settings._no_reduce_from_xml_element(root) + settings._verbosity_from_xml_element(root) + settings._tabular_legendre_from_xml_element(root) + settings._temperature_from_xml_element(root) + settings._trace_from_xml_element(root) + settings._track_from_xml_element(root) + settings._ufs_mesh_from_xml_element(root) + settings._resonance_scattering_from_xml_element(root) + settings._create_fission_neutrons_from_xml_element(root) + settings._log_grid_bins_from_xml_element(root) + settings._dagmc_from_xml_element(root) # TODO: Get volume calculations - # Get fission neutrons - elem = root.find('create_fission_neutrons') - if elem is not None: - settings.create_fission_neutrons = elem.text == 'true' - - # Get log grid bins - elem = root.find('log_grid_bins') - if elem is not None: - settings.log_grid_bins = int(elem.text) - - # Get dagmc - elem = root.find('dagmc') - if elem is not None: - settings.dagmc = elem.text == 'true' - return settings diff --git a/openmc/source.py b/openmc/source.py index e278d8089..5bafaf698 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -2,6 +2,7 @@ from numbers import Real import sys from xml.etree import ElementTree as ET +from openmc._xml import get_text from openmc.stats.univariate import (Univariate, Discrete, Uniform, Maxwell, Watt, Normal, Muir, Tabular) from openmc.stats.multivariate import (UnitSphere, Spatial, PolarAzimuthal, @@ -158,21 +159,21 @@ class Source(object): """ source = cls() - strength = elem.find('strength') + strength = get_text(elem, 'strength') if strength is not None: - source.strength = float(strength.text) + source.strength = float(strength) - particle = elem.find('particle') + particle = get_text(elem, 'particle') if particle is not None: - source.particle = particle.text + source.particle = particle - filename = elem.find('file') + filename = get_text(elem, 'file') if filename is not None: - source.file = filename.text + source.file = filename space = elem.find('space') if space is not None: - space_type = space.get('type') + space_type = get_text(space, 'type') if space_type == 'cartesian': source.space = CartesianIndependent.from_xml_element(space) elif space_type == 'box' or space_type == 'fission': @@ -182,7 +183,7 @@ class Source(object): angle = elem.find('angle') if angle is not None: - angle_type = angle.get('type') + angle_type = get_text(angle, 'type') if angle_type == 'mu-phi': source.angle = PolarAzimuthal.from_xml_element(angle) elif angle_type == 'isotropic': @@ -192,7 +193,7 @@ class Source(object): energy = elem.find('energy') if energy is not None: - energy_type = energy.get('type') + energy_type = get_text(energy, 'type') if energy_type == 'discrete': source.energy = Discrete.from_xml_element(energy) elif energy_type == 'uniform': diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 61122ada2..5a836ab63 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -8,6 +8,7 @@ from xml.etree import ElementTree as ET import numpy as np import openmc.checkvalue as cv +from openmc._xml import get_text from openmc.stats.univariate import Univariate, Uniform @@ -142,7 +143,7 @@ class PolarAzimuthal(UnitSphere): """ mu_phi = cls() - params = elem.findtext('parameters') + params = get_text(elem, 'parameters') if params is not None: mu_phi.reference_uvw = [float(x) for x in params.split()] mu_phi.mu = openmc.stats.Univariate.from_xml_element(elem.find('mu')) @@ -239,7 +240,7 @@ class Monodirectional(UnitSphere): """ monodirectional = cls() - params = elem.findtext('parameters') + params = get_text(elem, 'parameters') if params is not None: monodirectional.reference_uvw = [float(x) for x in params.split()] return monodirectional @@ -457,10 +458,10 @@ class Box(Spatial): Box distribution generated from XML element """ - only_fissionable = elem.get('type') == 'fission' - params = [float(x) for x in elem.findtext('parameters').split()] + only_fissionable = get_text(elem, 'type') == 'fission' + params = [float(x) for x in get_text(elem, 'parameters').split()] lower_left = params[:len(params)//2] - upper_right = paramx[len(params)//2:] + upper_right = params[len(params)//2:] return cls(lower_left, upper_right, only_fissionable) @@ -526,5 +527,5 @@ class Point(Spatial): Point distribution generated from XML element """ - xyz = [float(x) for x in elem.findtext('parameters').split()] + xyz = [float(x) for x in get_text(elem, 'parameters').split()] return cls(xyz) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 6a83cc780..3088f8c77 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -7,6 +7,7 @@ from xml.etree import ElementTree as ET import numpy as np import openmc.checkvalue as cv +from openmc._xml import get_text from openmc.mixin import EqualityMixin @@ -130,7 +131,7 @@ class Discrete(Univariate): Discrete distribution generated from XML element """ - params = [float(x) for x in elem.findtext('parameters').split()] + params = [float(x) for x in get_text(elem, 'parameters').split()] x = params[:len(params)//2] p = params[len(params)//2:] return cls(x, p) @@ -221,7 +222,7 @@ class Uniform(Univariate): Uniform distribution generated from XML element """ - params = elem.findtext('parameters').split() + params = get_text(elem, 'parameters').split() a = float(params[0]) b = float(params[1]) return cls(a, b) @@ -297,7 +298,7 @@ class Maxwell(Univariate): Maxwellian distribution generated from XML element """ - theta = float(elem.findtext('parameters')) + theta = float(get_text(elem, 'parameters')) return cls(theta) @@ -386,7 +387,7 @@ class Watt(Univariate): Watt distribution generated from XML element """ - params = elem.findtext('parameters').split() + params = get_text(elem, 'parameters').split() a = float(params[0]) b = float(params[1]) return watt(a, b) @@ -476,7 +477,7 @@ class Normal(Univariate): Normal distribution generated from XML element """ - params = elem.findtext('parameters').split() + params = get_text(elem, 'parameters').split() mean_value = float(params[0]) std_dev = float(params[1]) return cls(mean_value, std_dev) @@ -585,7 +586,7 @@ class Muir(Univariate): Muir distribution generated from XML element """ - params = elem.findtext('parameters').split() + params = get_text(elem, 'parameters').split() e0 = float(params[0]) m_rat = float(params[1]) kt = float(params[2]) @@ -702,8 +703,8 @@ class Tabular(Univariate): Tabular distribution generated from XML element """ - interpolation = elem.get('interpolation') - params = [float(x) for x in elem.findtext('parameters').split()] + interpolation = get_text(elem, 'interpolation') + params = [float(x) for x in get_text(elem, 'parameters').split()] x = params[:len(params)//2] p = paramx[len(params)//2:] return cls(x, p, interpolation) From c0ede1ec8dbd930d6c16bd8c0b72a9b722566f6d Mon Sep 17 00:00:00 2001 From: amandalund Date: Tue, 30 Apr 2019 11:39:17 -0400 Subject: [PATCH 43/61] Add unit tests for settings.from_xml() and a few fixes --- openmc/settings.py | 26 +++++++++----- tests/unit_tests/test_settings.py | 58 +++++++++++++++++++++++++++++-- 2 files changed, 73 insertions(+), 11 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 39312899b..802046c4e 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -174,7 +174,6 @@ class Settings(object): self._source = cv.CheckedList(Source, 'source distributions') self._confidence_intervals = None - self._cross_sections = None self._electron_treatment = None self._photon_transport = None self._ptables = None @@ -993,11 +992,14 @@ class Settings(object): def _sourcepoint_from_xml_element(self, root): elem = root.find('source_point') if elem is not None: - for key in ('separate', 'write', 'overwrite', 'batches'): + for key in ('separate', 'write', 'overwrite_latest', 'batches'): value = get_text(elem, key) if value is not None: - if key in ('separate', 'write', 'overwrite'): + if key in ('separate', 'write'): value = value == 'true' + elif key == 'overwrite_latest': + value = value == 'true' + key = 'overwrite' else: value = [int(x) for x in value.split()] self.sourcepoint[key] = value @@ -1053,9 +1055,12 @@ class Settings(object): self.cutoff[key] = float(value) def _entropy_mesh_from_xml_element(self, root): - elem = root.find('entropy_mesh') - if elem is not None: - self.entropy_mesh = Mesh.from_xml_element(elem) + text = get_text(root, 'entropy_mesh') + if text is not None: + path = "./mesh[@id='{}']".format(int(text)) + elem = root.find(path) + if elem is not None: + self.entropy_mesh = Mesh.from_xml_element(elem) def _trigger_from_xml_element(self, root): elem = root.find('trigger') @@ -1115,9 +1120,12 @@ class Settings(object): self.track = [int(x) for x in text.split()] def _ufs_mesh_from_xml_element(self, root): - elem = root.find('ufs_mesh') - if elem is not None: - self.ufs_mesh = Mesh.from_xml_element(elem) + text = get_text(root, 'ufs_mesh') + if text is not None: + path = "./mesh[@id='{}']".format(int(text)) + elem = root.find(path) + if elem is not None: + self.ufs_mesh = Mesh.from_xml_element(elem) def _resonance_scattering_from_xml_element(self, root): elem = root.find('resonance_scattering') diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index e42f6240f..f99a32a52 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -19,11 +19,12 @@ def test_export_to_xml(run_in_tmpdir): 'write': True, 'overwrite': True} s.statepoint = {'batches': [50, 150, 500, 1000]} s.confidence_intervals = True - s.cross_sections = '/path/to/cross_sections.xml' s.ptables = True s.seed = 17 s.survival_biasing = True - s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy': 1.0e-5} + s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, + 'energy_photon': 1000.0, 'energy_electron': 1.0e-5, + 'energy_positron': 1.0e-5} mesh = openmc.Mesh() mesh.lower_left = (-10., -10., -10.) mesh.upper_right = (10., 10., 10.) @@ -47,6 +48,59 @@ def test_export_to_xml(run_in_tmpdir): upper_right = (10., 10., 10.)) s.create_fission_neutrons = True s.log_grid_bins = 2000 + s.photon_transport = False + s.electron_treatment = 'led' + s.dagmc = False # Make sure exporting XML works s.export_to_xml() + + # Generate settings from XML + s = openmc.Settings.from_xml() + assert s.run_mode == 'fixed source' + assert s.batches == 1000 + assert s.generations_per_batch == 10 + assert s.inactive == 100 + assert s.particles == 1000000 + assert s.keff_trigger == {'type': 'std_dev', 'threshold': 0.001} + assert s.energy_mode == 'continuous-energy' + assert s.max_order == 5 + assert isinstance(s.source[0], openmc.Source) + assert isinstance(s.source[0].space, openmc.stats.Point) + assert s.output == {'summary': True, 'tallies': False, 'path': 'here'} + assert s.verbosity == 7 + assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True, + 'write': True, 'overwrite': True} + assert s.statepoint == {'batches': [50, 150, 500, 1000]} + assert s.confidence_intervals + assert s.ptables + assert s.seed == 17 + assert s.survival_biasing + assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5, + 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, + 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5} + assert isinstance(s.entropy_mesh, openmc.Mesh) + assert s.entropy_mesh.lower_left == [-10., -10., -10.] + assert s.entropy_mesh.upper_right == [10., 10., 10.] + assert s.entropy_mesh.dimension == [5, 5, 5] + assert s.trigger_active + assert s.trigger_max_batches == 10000 + assert s.trigger_batch_interval == 50 + assert not s.no_reduce + assert s.tabular_legendre == {'enable': True, 'num_points': 50} + assert s.temperature == {'default': 293.6, 'method': 'interpolation', + 'multipole': True, 'range': [200., 1000.]} + assert s.trace == [10, 1, 20] + assert s.track == [1, 1, 1, 2, 1, 1] + assert isinstance(s.ufs_mesh, openmc.Mesh) + assert s.ufs_mesh.lower_left == [-10., -10., -10.] + assert s.ufs_mesh.upper_right == [10., 10., 10.] + assert s.ufs_mesh.dimension == [5, 5, 5] + assert s.resonance_scattering == {'enable': True, 'method': 'rvs', + 'energy_min': 1.0, 'energy_max': 1000.0, + 'nuclides': ['U235', 'U238', 'Pu239']} + assert s.create_fission_neutrons + assert s.log_grid_bins == 2000 + assert not s.photon_transport + assert s.electron_treatment == 'led' + assert not s.dagmc From 7d7d9b9f0db6cde385104398b4c65b1e1b211e4d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 30 Apr 2019 23:23:46 -0500 Subject: [PATCH 44/61] Fix for latest pip. --- tools/ci/travis-install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index 2bca50e4b..51e45f64b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -23,7 +23,7 @@ fi python tools/ci/travis-install.py # Install Python API in editable mode -pip install -e .[test,vtk] +pip install --no-use-pep517 -e .[test,vtk] # For uploading to coveralls pip install coveralls From 81674a1c2d5ee2e1c97d02dd1c84e5b936b4ca1a Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 1 May 2019 12:20:46 -0400 Subject: [PATCH 45/61] A couple more fixes for mesh to/from XML --- openmc/mesh.py | 10 +++++++--- openmc/settings.py | 3 ++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index fb59053ae..2489d5faa 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -231,8 +231,9 @@ class Mesh(IDManagerMixin): element.set("id", str(self._id)) element.set("type", self._type) - subelement = ET.SubElement(element, "dimension") - subelement.text = ' '.join(map(str, self._dimension)) + if self._dimension is not None: + subelement = ET.SubElement(element, "dimension") + subelement.text = ' '.join(map(str, self._dimension)) subelement = ET.SubElement(element, "lower_left") subelement.text = ' '.join(map(str, self._lower_left)) @@ -264,7 +265,10 @@ class Mesh(IDManagerMixin): """ mesh_id = int(get_text(elem, 'id')) mesh = cls(mesh_id) - mesh.type = get_text(elem, 'type') + + mesh_type = get_text(elem, 'type') + if mesh_type is not None: + mesh.type = mesh_type dimension = get_text(elem, 'dimension') if dimension is not None: diff --git a/openmc/settings.py b/openmc/settings.py index 802046c4e..063677575 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -551,7 +551,8 @@ class Settings(object): @entropy_mesh.setter def entropy_mesh(self, entropy): cv.check_type('entropy mesh', entropy, Mesh) - cv.check_length('entropy mesh dimension', entropy.dimension, 3) + if entropy.dimension: + cv.check_length('entropy mesh dimension', entropy.dimension, 3) cv.check_length('entropy mesh lower-left corner', entropy.lower_left, 3) cv.check_length('entropy mesh upper-right corner', entropy.upper_right, 3) self._entropy_mesh = entropy From 1086b6cada1a4452dafba8a86e474c927fa97839 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 7 May 2019 06:59:41 -0500 Subject: [PATCH 46/61] Remove --no-use-pep517 on pip install -e --- tools/ci/travis-install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index 51e45f64b..2bca50e4b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -23,7 +23,7 @@ fi python tools/ci/travis-install.py # Install Python API in editable mode -pip install --no-use-pep517 -e .[test,vtk] +pip install -e .[test,vtk] # For uploading to coveralls pip install coveralls From 2a2b7e9b3ccfae7078e3c6aac8de882b97252828 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 May 2019 18:05:04 -0500 Subject: [PATCH 47/61] Fix bug in Plot.__repr__ --- openmc/plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/plots.py b/openmc/plots.py index eb582502e..00d96fd7e 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -439,7 +439,7 @@ class Plot(IDManagerMixin): string += '{: <16}=\t{}\n'.format('\tBasis', self._basis) string += '{: <16}=\t{}\n'.format('\tWidth', self._width) string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin) - string += '{: <16}=\t{}\n'.format('\tPixels', self._origin) + string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels) string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by) string += '{: <16}=\t{}\n'.format('\tBackground', self._background) string += '{: <16}=\t{}\n'.format('\tMask components', From 18b225f14d9ba2b19a565b8c706d36e028300dd1 Mon Sep 17 00:00:00 2001 From: aprilnovak Date: Mon, 29 Apr 2019 15:40:57 -0700 Subject: [PATCH 48/61] Add openmc_cell_get_temperature method. Refs #1223 --- docs/source/capi/index.rst | 11 +++++++++++ include/openmc/capi.h | 1 + openmc/capi/cell.py | 21 +++++++++++++++++++++ src/cell.cpp | 30 ++++++++++++++++++++++++++++++ 4 files changed, 63 insertions(+) diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index 85e5c8ced..274a5240f 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -73,6 +73,17 @@ Functions :return: Return status (negative if an error occurred) :rtype: int +.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) + + Get the temperature of a cell + + :param int32_t index: Index in the cells array + :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature + of the first instance is returned. + :param double* T: temperature of the cell + :return: Return status (negative if an error occurred) + :rtype: int + .. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices) Set the fill for a cell diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 8815d243a..1f6686df1 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -13,6 +13,7 @@ extern "C" { int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n); int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); + int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); int openmc_cell_set_id(int32_t index, int32_t id); int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance); diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 7ad5c75d7..ae536e725 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -24,6 +24,10 @@ _dll.openmc_cell_get_fill.argtypes = [ c_int32, POINTER(c_int), POINTER(POINTER(c_int32)), POINTER(c_int32)] _dll.openmc_cell_get_fill.restype = c_int _dll.openmc_cell_get_fill.errcheck = _error_handler +_dll.openmc_cell_get_temperature.argtypes = [ + c_int32, POINTER(c_int32), POINTER(c_double)] +_dll.openmc_cell_get_temperature.restype = c_int +_dll.openmc_cell_get_temperature.errcheck = _error_handler _dll.openmc_cell_set_fill.argtypes = [ c_int32, c_int, c_int32, POINTER(c_int32)] _dll.openmc_cell_set_fill.restype = c_int @@ -128,6 +132,23 @@ class Cell(_FortranObjectWithID): indices = (c_int32*1)(-1) _dll.openmc_cell_set_fill(self._index, 1, 1, indices) + def get_temperature(self, instance=None): + """Get the temperature of a cell + + Parameters + ---------- + instance: int or None + Which instance of the cell + + """ + + if instance is not None: + instance = c_int32(instance) + + T = c_double() + _dll.openmc_cell_get_temperature(self._index, instance, T) + return T.value + def set_temperature(self, T, instance=None): """Set the temperature of a cell diff --git a/src/cell.cpp b/src/cell.cpp index 59057e084..650a18a8a 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -922,6 +922,36 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) return 0; } +extern "C" int +openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) +{ + if (index < 0 || index >= model::cells.size()) { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + Cell& c {*model::cells[index]}; + + if (c.sqrtkT_.size() < 1) { + strcpy(openmc_err_msg, "Cell temperature has not yet been set."); + return OPENMC_E_UNASSIGNED; + } + + if (instance) { + if (*instance >= 0 && *instance < c.n_instances_) { + double sqrtkT = c.sqrtkT_.size() == 1 ? c.sqrtkT_[0] : c.sqrtkT_[*instance]; + *T = sqrtkT * sqrtkT / K_BOLTZMANN; + } else { + strcpy(openmc_err_msg, "Distribcell instance is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + } else { + *T = c.sqrtkT_[0] * c.sqrtkT_[0] / K_BOLTZMANN; + } + + return 0; +} + //! Return the index in the cells array of a cell with a given ID extern "C" int openmc_get_cell_index(int32_t id, int32_t* index) From 3cb99504fcd352423314765208dbebeecd49b18e Mon Sep 17 00:00:00 2001 From: aprilnovak Date: Mon, 6 May 2019 13:38:07 -0700 Subject: [PATCH 49/61] Unit tests for getting and setting cell temperature. Refs #1223 --- openmc/capi/cell.py | 8 ++++++-- tests/unit_tests/test_capi.py | 8 ++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index ae536e725..959ab08fc 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -143,7 +143,7 @@ class Cell(_FortranObjectWithID): """ if instance is not None: - instance = c_int32(instance) + instance = c_int32(instance) T = c_double() _dll.openmc_cell_get_temperature(self._index, instance, T) @@ -160,7 +160,11 @@ class Cell(_FortranObjectWithID): Which instance of the cell """ - _dll.openmc_cell_set_temperature(self._index, T, c_int32(instance)) + + if instance is not None: + instance = c_int32(instance) + + _dll.openmc_cell_set_temperature(self._index, T, instance) class _CellMapping(Mapping): diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index be88072d9..07fef8516 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -75,6 +75,14 @@ def test_cell(capi_init): assert str(cell) == 'Cell[0]' +def test_cell_temperature(capi_init): + cell = openmc.capi.cells[1] + cell.set_temperature(100.0, 0) + assert cell.get_temperature(0) == 100.0 + cell.set_temperature(200) + assert cell.get_temperature() == 200.0 + + def test_new_cell(capi_init): with pytest.raises(exc.AllocationError): openmc.capi.Cell(1) From 5fbae7cab47f68fd8a1686a8181e336ef22e4033 Mon Sep 17 00:00:00 2001 From: lavistam <20384517+matiaslavista@users.noreply.github.com> Date: Sun, 19 May 2019 21:44:49 -0500 Subject: [PATCH 50/61] Moving special case of keff calculation (n<=3) to cpp side --- openmc/capi/core.py | 15 +++------------ src/eigenvalue.cpp | 12 +++++++++--- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index ad4f77dff..baba0b833 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -230,18 +230,9 @@ def keff(): Mean k-eigenvalue and standard deviation of the mean """ - n = openmc.capi.num_realizations() - if n > 3: - # Use the combined estimator if there are enough realizations - k = (c_double*2)() - _dll.openmc_get_keff(k) - return tuple(k) - else: - # Otherwise, return the tracklength estimator - mean = c_double.in_dll(_dll, 'keff').value - std_dev = c_double.in_dll(_dll, 'keff_std').value \ - if n > 1 else np.inf - return (mean, std_dev) + k = (c_double*2)() + _dll.openmc_get_keff(k) + return tuple(k) def master(): diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 08a4d3f91..0d6aa19b3 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -29,6 +29,7 @@ #include // for sqrt, abs, pow #include // for back_inserter #include +#include //for infinity namespace openmc { @@ -388,10 +389,15 @@ int openmc_get_keff(double* k_combined) k_combined[0] = 0.0; k_combined[1] = 0.0; - // Make sure we have at least four realizations. Notice that at the end, - // there is a N-3 term in a denominator. + //Special case for n <=3. Notice that at the end, + //there is a N-3 term in a denominator. if (simulation::n_realizations <= 3) { - return -1; + k_combined[0] = simulation::keff; + k_combined[1] = simulation::keff_std; + if (simulation::n_realizations <=1){ + k_combined[1] = std::numeric_limits::infinity(); + } + return 0; } // Initialize variables From e15dfb27ede96775ab384a423cc801154ff4f361 Mon Sep 17 00:00:00 2001 From: guillaume Date: Mon, 20 May 2019 14:26:15 -0400 Subject: [PATCH 51/61] Add coverage of cpp source files using cpp-coveralls --- .travis.yml | 3 ++- CMakeLists.txt | 4 ++++ tools/ci/travis-install.py | 3 +++ tools/ci/travis-install.sh | 5 ++++- 4 files changed, 13 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index ade182788..9e51a70eb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -53,4 +53,5 @@ before_script: script: - ./tools/ci/travis-script.sh after_success: - - coveralls + - cpp-coveralls --exclude-pattern "/usr/*" --dump cpp_cov.json + - coveralls --merge=cpp_cov.json diff --git a/CMakeLists.txt b/CMakeLists.txt index 4033d2466..fea014f12 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -95,6 +95,10 @@ if(optimize) list(REMOVE_ITEM cxxflags -O2) list(APPEND cxxflags -O3) endif() +if(coverage) + list(APPEND cxxflags --coverage) + list(APPEND ldflags --coverage) +endif() # Show flags being used message(STATUS "OpenMC C++ flags: ${cxxflags}") diff --git a/tools/ci/travis-install.py b/tools/ci/travis-install.py index 9c12387fe..9a9b06dae 100644 --- a/tools/ci/travis-install.py +++ b/tools/ci/travis-install.py @@ -48,6 +48,9 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False): if dagmc: cmake_cmd.append('-Ddagmc=ON') + # Build in coverage mode for coverage testing + cmake_cmd.append('-Dcoverage=on') + # Build and install cmake_cmd.append('..') print(' '.join(cmake_cmd)) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index 2bca50e4b..4e32f31cc 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -25,5 +25,8 @@ python tools/ci/travis-install.py # Install Python API in editable mode pip install -e .[test,vtk] -# For uploading to coveralls +# For coverage testing of the C++ source files +pip install cpp-coveralls + +# For coverage testing of the Python source files pip install coveralls From 14a95b63731a15ed3ce47a8ad7c63c0fdb576fa4 Mon Sep 17 00:00:00 2001 From: lavistam <20384517+matiaslavista@users.noreply.github.com> Date: Tue, 21 May 2019 13:14:50 -0500 Subject: [PATCH 52/61] Change in Calculating K-Eff from cpp side for cases of batches (n) <= 3 has changed. Therefore filter_distribcell, which ran batches of 1,1,3 and 1, has been updated to reflect the more accurate calcultion. sourcepoint_batches explanation: This test needed to be updated because statepoint files are written for 3 or fewer batches, which will change the k-eff result written in write_eigenvalue_hdf5 for those statepoint files. All other k-effective values are unaffected and the final result of the simulation is the same. --- .../regression_tests/filter_distribcell/case-1/results_true.dat | 2 +- .../regression_tests/filter_distribcell/case-2/results_true.dat | 2 +- .../regression_tests/filter_distribcell/case-3/results_true.dat | 2 +- .../regression_tests/filter_distribcell/case-4/results_true.dat | 2 +- tests/regression_tests/sourcepoint_batch/results_true.dat | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/regression_tests/filter_distribcell/case-1/results_true.dat b/tests/regression_tests/filter_distribcell/case-1/results_true.dat index 2d9835b1a..80036a7c4 100644 --- a/tests/regression_tests/filter_distribcell/case-1/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-1/results_true.dat @@ -1,5 +1,5 @@ k-combined: -0.000000E+00 0.000000E+00 +5.497140E-02 INF tally 1: 1.548980E-02 2.399339E-04 diff --git a/tests/regression_tests/filter_distribcell/case-2/results_true.dat b/tests/regression_tests/filter_distribcell/case-2/results_true.dat index 57e40f496..bbf38b3d1 100644 --- a/tests/regression_tests/filter_distribcell/case-2/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-2/results_true.dat @@ -1,5 +1,5 @@ k-combined: -0.000000E+00 0.000000E+00 +2.149726E-02 INF tally 1: 7.588170E-03 5.758032E-05 diff --git a/tests/regression_tests/filter_distribcell/case-3/results_true.dat b/tests/regression_tests/filter_distribcell/case-3/results_true.dat index 0a4dd9dc7..be197b5e9 100644 --- a/tests/regression_tests/filter_distribcell/case-3/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-3/results_true.dat @@ -1 +1 @@ -11755ecac8355b5e79384f5d72974e618b6f95500c0a5718c01bb340c5d1c8ceafc33de65b9e94b7e9d78f16ec209f8a0fdf6a531eab5430657688d0125db7ef \ No newline at end of file +2ee0162762999f71ad2178936509fd9e054928023a9e0e90c078cede3ecf6583267069486adf15f1b02188333d1bfe1fc41b341bc00aab53feddd1395441f1f8 \ No newline at end of file diff --git a/tests/regression_tests/filter_distribcell/case-4/results_true.dat b/tests/regression_tests/filter_distribcell/case-4/results_true.dat index 078e74bed..43a837e98 100644 --- a/tests/regression_tests/filter_distribcell/case-4/results_true.dat +++ b/tests/regression_tests/filter_distribcell/case-4/results_true.dat @@ -1,5 +1,5 @@ k-combined: -0.000000E+00 0.000000E+00 +1.121246E-01 INF tally 1: 2.265319E-02 5.131668E-04 diff --git a/tests/regression_tests/sourcepoint_batch/results_true.dat b/tests/regression_tests/sourcepoint_batch/results_true.dat index 8b8dd98f2..a9843df53 100644 --- a/tests/regression_tests/sourcepoint_batch/results_true.dat +++ b/tests/regression_tests/sourcepoint_batch/results_true.dat @@ -1,3 +1,3 @@ k-combined: -0.000000E+00 0.000000E+00 +2.976389E-01 3.770725E-03 1.892327E+00 -3.385257E+00 6.702632E-01 From 8ce7a63a045a69908920a33d91df981ccc03c3c8 Mon Sep 17 00:00:00 2001 From: Giud Date: Tue, 21 May 2019 23:05:23 -0400 Subject: [PATCH 53/61] Update .travis.yml --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 9e51a70eb..d0b3de964 100644 --- a/.travis.yml +++ b/.travis.yml @@ -53,5 +53,5 @@ before_script: script: - ./tools/ci/travis-script.sh after_success: - - cpp-coveralls --exclude-pattern "/usr/*" --dump cpp_cov.json + - cpp-coveralls --exclude-pattern "/usr/*" -e "build" --dump cpp_cov.json - coveralls --merge=cpp_cov.json From 55b9db9c446b7082bb11ad8041b10700fc7bc11d Mon Sep 17 00:00:00 2001 From: Giud Date: Wed, 22 May 2019 12:27:46 -0400 Subject: [PATCH 54/61] Update .travis.yml --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index d0b3de964..096459b00 100644 --- a/.travis.yml +++ b/.travis.yml @@ -53,5 +53,5 @@ before_script: script: - ./tools/ci/travis-script.sh after_success: - - cpp-coveralls --exclude-pattern "/usr/*" -e "build" --dump cpp_cov.json + - cpp-coveralls -i src -i include --exclude-pattern "/usr/*" --dump cpp_cov.json - coveralls --merge=cpp_cov.json From 8fd5d76df38bcb232c9789af1a225f6af76162ee Mon Sep 17 00:00:00 2001 From: lavistam <20384517+matiaslavista@users.noreply.github.com> Date: Wed, 22 May 2019 15:17:46 -0500 Subject: [PATCH 55/61] Style --- src/eigenvalue.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 0d6aa19b3..04fdabab1 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -389,12 +389,12 @@ int openmc_get_keff(double* k_combined) k_combined[0] = 0.0; k_combined[1] = 0.0; - //Special case for n <=3. Notice that at the end, - //there is a N-3 term in a denominator. + // Special case for n <=3. Notice that at the end, + // there is a N-3 term in a denominator. if (simulation::n_realizations <= 3) { k_combined[0] = simulation::keff; k_combined[1] = simulation::keff_std; - if (simulation::n_realizations <=1){ + if (simulation::n_realizations <=1) { k_combined[1] = std::numeric_limits::infinity(); } return 0; From 2f09b7ca854ba9a4370687f2bd0521a174c2407f Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 22 May 2019 19:26:15 -0600 Subject: [PATCH 56/61] Fixed some univariate/multivariate from_xml_element issues and updated unit tests --- openmc/settings.py | 28 ++++++------ openmc/source.py | 39 +++------------- openmc/stats/multivariate.py | 26 ++++++++--- openmc/stats/univariate.py | 39 +++++++++------- tests/unit_tests/test_source.py | 8 +++- tests/unit_tests/test_stats.py | 80 +++++++++++++++++++++++++-------- 6 files changed, 130 insertions(+), 90 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 063677575..4fe87d908 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -980,7 +980,7 @@ class Settings(object): value = get_text(elem, key) if value is not None: if key in ('summary', 'tallies'): - value = value == 'true' + value = value in ('true', '1') self.output[key] = value def _statepoint_from_xml_element(self, root): @@ -997,9 +997,9 @@ class Settings(object): value = get_text(elem, key) if value is not None: if key in ('separate', 'write'): - value = value == 'true' + value = value in ('true', '1') elif key == 'overwrite_latest': - value = value == 'true' + value = value in ('true', '1') key = 'overwrite' else: value = [int(x) for x in value.split()] @@ -1008,7 +1008,7 @@ class Settings(object): def _confidence_intervals_from_xml_element(self, root): text = get_text(root, 'confidence_intervals') if text is not None: - self.confidence_intervals = text == 'true' + self.confidence_intervals = text in ('true', '1') def _electron_treatment_from_xml_element(self, root): text = get_text(root, 'electron_treatment') @@ -1028,12 +1028,12 @@ class Settings(object): def _photon_transport_from_xml_element(self, root): text = get_text(root, 'photon_transport') if text is not None: - self.photon_transport = text == 'true' + self.photon_transport = text in ('true', '1') def _ptables_from_xml_element(self, root): text = get_text(root, 'ptables') if text is not None: - self.ptables = text == 'true' + self.ptables = text in ('true', '1') def _seed_from_xml_element(self, root): text = get_text(root, 'seed') @@ -1043,7 +1043,7 @@ class Settings(object): def _survival_biasing_from_xml_element(self, root): text = get_text(root, 'survival_biasing') if text is not None: - self.survival_biasing = text == 'true' + self.survival_biasing = text in ('true', '1') def _cutoff_from_xml_element(self, root): elem = root.find('cutoff') @@ -1066,7 +1066,7 @@ class Settings(object): def _trigger_from_xml_element(self, root): elem = root.find('trigger') if elem is not None: - self.trigger_active = get_text(elem, 'active') == 'true' + 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) @@ -1077,7 +1077,7 @@ class Settings(object): def _no_reduce_from_xml_element(self, root): text = get_text(root, 'no_reduce') if text is not None: - self.no_reduce = text == 'true' + self.no_reduce = text in ('true', '1') def _verbosity_from_xml_element(self, root): text = get_text(root, 'verbosity') @@ -1088,7 +1088,7 @@ class Settings(object): elem = root.find('tabular_legendre') if elem is not None: text = get_text(elem, 'enable') - self.tabular_legendre['enable'] = text == 'true' + 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) @@ -1108,7 +1108,7 @@ class Settings(object): 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 == 'true' + self.temperature['multipole'] = text in ('true', '1') def _trace_from_xml_element(self, root): text = get_text(root, 'trace') @@ -1136,7 +1136,7 @@ class Settings(object): value = get_text(elem, key) if value is not None: if key == 'enable': - value = value == 'true' + value = value in ('true', '1') elif key in ('energy_min', 'energy_max'): value = float(value) elif key == 'nuclides': @@ -1146,7 +1146,7 @@ class Settings(object): def _create_fission_neutrons_from_xml_element(self, root): text = get_text(root, 'create_fission_neutrons') if text is not None: - self.create_fission_neutrons = text == 'true' + self.create_fission_neutrons = text in ('true', '1') def _log_grid_bins_from_xml_element(self, root): text = get_text(root, 'log_grid_bins') @@ -1156,7 +1156,7 @@ class Settings(object): def _dagmc_from_xml_element(self, root): text = get_text(root, 'dagmc') if text is not None: - self.dagmc = text == 'true' + self.dagmc = text in ('true', '1') def export_to_xml(self, path='settings.xml'): """Export simulation settings to an XML file. diff --git a/openmc/source.py b/openmc/source.py index 5bafaf698..88c2f8611 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -3,11 +3,8 @@ import sys from xml.etree import ElementTree as ET from openmc._xml import get_text -from openmc.stats.univariate import (Univariate, Discrete, Uniform, Maxwell, - Watt, Normal, Muir, Tabular) -from openmc.stats.multivariate import (UnitSphere, Spatial, PolarAzimuthal, - Isotropic, Monodirectional, Box, Point, - CartesianIndependent) +from openmc.stats.univariate import Univariate +from openmc.stats.multivariate import UnitSphere, Spatial import openmc.checkvalue as cv @@ -173,40 +170,14 @@ class Source(object): space = elem.find('space') if space is not None: - space_type = get_text(space, 'type') - if space_type == 'cartesian': - source.space = CartesianIndependent.from_xml_element(space) - elif space_type == 'box' or space_type == 'fission': - source.space = Box.from_xml_element(space) - elif space_type == 'point': - source.space = Point.from_xml_element(space) + source.space = Spatial.from_xml_element(space) angle = elem.find('angle') if angle is not None: - angle_type = get_text(angle, 'type') - if angle_type == 'mu-phi': - source.angle = PolarAzimuthal.from_xml_element(angle) - elif angle_type == 'isotropic': - source.angle = Isotropic.from_xml_element(angle) - elif angle_type == 'monodirectional': - source.angle = Monodirectional.from_xml_element(angle) + source.angle = UnitSphere.from_xml_element(angle) energy = elem.find('energy') if energy is not None: - energy_type = get_text(energy, 'type') - if energy_type == 'discrete': - source.energy = Discrete.from_xml_element(energy) - elif energy_type == 'uniform': - source.energy = Uniform.from_xml_element(energy) - elif energy_type == 'maxwell': - source.energy = Maxwell.from_xml_element(energy) - elif energy_type == 'watt': - source.energy = Watt.from_xml_element(energy) - elif energy_type == 'normal': - source.energy = Normal.from_xml_element(energy) - elif energy_type == 'muir': - source.energy = Muir.from_xml_element(energy) - elif energy_type == 'tabular': - source.energy = Tabular.from_xml_element(energy) + source.energy = Univariate.from_xml_element(energy) return source diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 5a836ab63..35afc21dd 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -51,7 +51,13 @@ class UnitSphere(metaclass=ABCMeta): @classmethod @abstractmethod def from_xml_element(cls, elem): - pass + distribution = get_text(elem, 'type') + if distribution == 'mu-phi': + return PolarAzimuthal.from_xml_element(elem) + elif distribution == 'isotropic': + return Isotropic.from_xml_element(elem) + elif distribution == 'monodirectional': + return Monodirectional.from_xml_element(elem) class PolarAzimuthal(UnitSphere): @@ -146,8 +152,8 @@ class PolarAzimuthal(UnitSphere): params = get_text(elem, 'parameters') if params is not None: mu_phi.reference_uvw = [float(x) for x in params.split()] - mu_phi.mu = openmc.stats.Univariate.from_xml_element(elem.find('mu')) - mu_phi.phi = openmc.stats.Univariate.from_xml_element(elem.find('phi')) + mu_phi.mu = Univariate.from_xml_element(elem.find('mu')) + mu_phi.phi = Univariate.from_xml_element(elem.find('phi')) return mu_phi @@ -263,7 +269,13 @@ class Spatial(metaclass=ABCMeta): @classmethod @abstractmethod def from_xml_element(cls, elem): - pass + distribution = get_text(elem, 'type') + if distribution == 'cartesian': + return CartesianIndependent.from_xml_element(elem) + elif distribution == 'box' or distribution == 'fission': + return Box.from_xml_element(elem) + elif distribution == 'point': + return Point.from_xml_element(elem) class CartesianIndependent(Spatial): @@ -357,9 +369,9 @@ class CartesianIndependent(Spatial): Spatial distribution generated from XML element """ - x = openmc.stats.Univariate.from_xml_element(elem.find('x')) - y = openmc.stats.Univariate.from_xml_element(elem.find('y')) - z = openmc.stats.Univariate.from_xml_element(elem.find('z')) + x = Univariate.from_xml_element(elem.find('x')) + y = Univariate.from_xml_element(elem.find('y')) + z = Univariate.from_xml_element(elem.find('z')) return cls(x, y, z) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 3088f8c77..363dd2ee3 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -36,7 +36,25 @@ class Univariate(EqualityMixin, metaclass=ABCMeta): @classmethod @abstractmethod def from_xml_element(cls, elem): - pass + distribution = get_text(elem, 'type') + if distribution == 'discrete': + return Discrete.from_xml_element(elem) + elif distribution == 'uniform': + return Uniform.from_xml_element(elem) + elif distribution == 'maxwell': + return Maxwell.from_xml_element(elem) + elif distribution == 'watt': + return Watt.from_xml_element(elem) + elif distribution == 'normal': + return Normal.from_xml_element(elem) + elif distribution == 'muir': + return Muir.from_xml_element(elem) + elif distribution == 'tabular': + return Tabular.from_xml_element(elem) + elif distribution == 'legendre': + return Legendre.from_xml_element(elem) + elif distribution == 'mixture': + return Mixture.from_xml_element(elem) class Discrete(Univariate): @@ -223,9 +241,7 @@ class Uniform(Univariate): """ params = get_text(elem, 'parameters').split() - a = float(params[0]) - b = float(params[1]) - return cls(a, b) + return cls(*map(float, params)) class Maxwell(Univariate): @@ -388,9 +404,7 @@ class Watt(Univariate): """ params = get_text(elem, 'parameters').split() - a = float(params[0]) - b = float(params[1]) - return watt(a, b) + return cls(*map(float, params)) class Normal(Univariate): @@ -478,9 +492,7 @@ class Normal(Univariate): """ params = get_text(elem, 'parameters').split() - mean_value = float(params[0]) - std_dev = float(params[1]) - return cls(mean_value, std_dev) + return cls(*map(float, params)) class Muir(Univariate): @@ -587,10 +599,7 @@ class Muir(Univariate): """ params = get_text(elem, 'parameters').split() - e0 = float(params[0]) - m_rat = float(params[1]) - kt = float(params[2]) - return muir(e0, m_rat, kt) + return cls(*map(float, params)) class Tabular(Univariate): @@ -706,7 +715,7 @@ class Tabular(Univariate): interpolation = get_text(elem, 'interpolation') params = [float(x) for x in get_text(elem, 'parameters').split()] x = params[:len(params)//2] - p = paramx[len(params)//2:] + p = params[len(params)//2:] return cls(x, p, interpolation) diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 3c963d052..1c70e159d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -11,7 +11,6 @@ def test_source(): assert src.space == space assert src.angle == angle assert src.energy == energy - assert src.strength == 1.0 elem = src.to_xml_element() assert 'strength' in elem.attrib @@ -19,6 +18,13 @@ def test_source(): assert elem.find('angle') is not None assert elem.find('energy') is not None + src = openmc.Source.from_xml_element(elem) + assert isinstance(src.angle, openmc.stats.Isotropic) + assert src.space.xyz == [0.0, 0.0, 0.0] + assert src.energy.x == [1.0e6] + assert src.energy.p == [1.0] + assert src.strength == 1.0 + def test_source_file(): filename = 'source.h5' diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index 553f2410e..e15958478 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -10,10 +10,15 @@ def test_discrete(): x = [0.0, 1.0, 10.0] p = [0.3, 0.2, 0.5] d = openmc.stats.Discrete(x, p) + elem = d.to_xml_element('distribution') + + d = openmc.stats.Discrete.from_xml_element(elem) assert d.x == x assert d.p == p assert len(d) == len(x) - d.to_xml_element('distribution') + + d = openmc.stats.Univariate.from_xml_element(elem) + assert isinstance(d, openmc.stats.Discrete) # Single point d2 = openmc.stats.Discrete(1e6, 1.0) @@ -25,6 +30,9 @@ def test_discrete(): def test_uniform(): a, b = 10.0, 20.0 d = openmc.stats.Uniform(a, b) + elem = d.to_xml_element('distribution') + + d = openmc.stats.Uniform.from_xml_element(elem) assert d.a == a assert d.b == b assert len(d) == 2 @@ -34,35 +42,39 @@ def test_uniform(): assert t.p == [1/(b-a), 1/(b-a)] assert t.interpolation == 'histogram' - d.to_xml_element('distribution') - def test_maxwell(): theta = 1.2895e6 d = openmc.stats.Maxwell(theta) + elem = d.to_xml_element('distribution') + + d = openmc.stats.Maxwell.from_xml_element(elem) assert d.theta == theta assert len(d) == 1 - d.to_xml_element('distribution') def test_watt(): a, b = 0.965e6, 2.29e-6 d = openmc.stats.Watt(a, b) + elem = d.to_xml_element('distribution') + + d = openmc.stats.Watt.from_xml_element(elem) assert d.a == a assert d.b == b assert len(d) == 2 - d.to_xml_element('distribution') def test_tabular(): x = [0.0, 5.0, 7.0] p = [0.1, 0.2, 0.05] d = openmc.stats.Tabular(x, p, 'linear-linear') + elem = d.to_xml_element('distribution') + + d = openmc.stats.Tabular.from_xml_element(elem) assert d.x == x assert d.p == p assert d.interpolation == 'linear-linear' assert len(d) == len(x) - d.to_xml_element('distribution') def test_legendre(): @@ -115,6 +127,15 @@ def test_polar_azimuthal(): assert elem.find('mu') is not None assert elem.find('phi') is not None + d = openmc.stats.PolarAzimuthal.from_xml_element(elem) + assert d.mu.x == [1.] + assert d.mu.p == [1.] + assert d.phi.x == [0.] + assert d.phi.p == [1.] + + d = openmc.stats.UnitSphere.from_xml_element(elem) + assert isinstance(d, openmc.stats.PolarAzimuthal) + def test_isotropic(): d = openmc.stats.Isotropic() @@ -122,24 +143,25 @@ def test_isotropic(): assert elem.tag == 'angle' assert elem.attrib['type'] == 'isotropic' + d = openmc.stats.Isotropic.from_xml_element(elem) + assert isinstance(d, openmc.stats.Isotropic) + def test_monodirectional(): d = openmc.stats.Monodirectional((1., 0., 0.)) - assert d.reference_uvw == pytest.approx((1., 0., 0.)) - elem = d.to_xml_element() assert elem.tag == 'angle' assert elem.attrib['type'] == 'monodirectional' + d = openmc.stats.Monodirectional.from_xml_element(elem) + assert d.reference_uvw == pytest.approx((1., 0., 0.)) + def test_cartesian(): x = openmc.stats.Uniform(-10., 10.) y = openmc.stats.Uniform(-10., 10.) z = openmc.stats.Uniform(0., 20.) d = openmc.stats.CartesianIndependent(x, y, z) - assert d.x == x - assert d.y == y - assert d.z == z elem = d.to_xml_element() assert elem.tag == 'space' @@ -147,55 +169,75 @@ def test_cartesian(): assert elem.find('x') is not None assert elem.find('y') is not None + d = openmc.stats.CartesianIndependent.from_xml_element(elem) + assert d.x == x + assert d.y == y + assert d.z == z + + d = openmc.stats.Spatial.from_xml_element(elem) + assert isinstance(d, openmc.stats.CartesianIndependent) + def test_box(): lower_left = (-10., -10., -10.) upper_right = (10., 10., 10.) d = openmc.stats.Box(lower_left, upper_right) - assert d.lower_left == pytest.approx(lower_left) - assert d.upper_right == pytest.approx(upper_right) - assert not d.only_fissionable elem = d.to_xml_element() assert elem.tag == 'space' assert elem.attrib['type'] == 'box' assert elem.find('parameters') is not None + d = openmc.stats.Box.from_xml_element(elem) + assert d.lower_left == pytest.approx(lower_left) + assert d.upper_right == pytest.approx(upper_right) + assert not d.only_fissionable + # only fissionable parameter d2 = openmc.stats.Box(lower_left, upper_right, True) assert d2.only_fissionable elem = d2.to_xml_element() assert elem.attrib['type'] == 'fission' + d = openmc.stats.Spatial.from_xml_element(elem) + assert isinstance(d, openmc.stats.Box) def test_point(): p = (-4., 2., 10.) d = openmc.stats.Point(p) - assert d.xyz == pytest.approx(p) elem = d.to_xml_element() assert elem.tag == 'space' assert elem.attrib['type'] == 'point' assert elem.find('parameters') is not None + d = openmc.stats.Point.from_xml_element(elem) + assert d.xyz == pytest.approx(p) + def test_normal(): mean = 10.0 std_dev = 2.0 d = openmc.stats.Normal(mean,std_dev) + + elem = d.to_xml_element('distribution') + assert elem.attrib['type'] == 'normal' + + d = openmc.stats.Normal.from_xml_element(elem) assert d.mean_value == pytest.approx(mean) assert d.std_dev == pytest.approx(std_dev) assert len(d) == 2 - elem = d.to_xml_element('distribution') - assert elem.attrib['type'] == 'normal' def test_muir(): mean = 10.0 mass = 5.0 temp = 20000. d = openmc.stats.Muir(mean,mass,temp) + + elem = d.to_xml_element('energy') + assert elem.attrib['type'] == 'muir' + + d = openmc.stats.Muir.from_xml_element(elem) assert d.e0 == pytest.approx(mean) assert d.m_rat == pytest.approx(mass) assert d.kt == pytest.approx(temp) assert len(d) == 3 - elem = d.to_xml_element('energy') - assert elem.attrib['type'] == 'muir' From 827064da85ae549cf2500ab680a2cfa9fbf09aa1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 21 May 2019 21:47:20 -0500 Subject: [PATCH 57/61] Fix bug whereby coherent elastic scattering cross section was zero --- src/thermal.cpp | 2 +- tests/regression_tests/lattice_hex_coincident/results_true.dat | 2 +- tests/regression_tests/salphabeta/results_true.dat | 2 +- tests/regression_tests/triso/results_true.dat | 2 +- tests/regression_tests/void/results_true.dat | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/thermal.cpp b/src/thermal.cpp index 3c2d9ba43..5db578f9d 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -202,7 +202,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, *inelastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1]; // Check for elastic data - if (E < sab.threshold_elastic_) { + if (!sab.elastic_e_in_.empty()) { // Determine whether elastic scattering is given in the coherent or // incoherent approximation. For coherent, the cross section is // represented as P/E whereas for incoherent, it is simply P diff --git a/tests/regression_tests/lattice_hex_coincident/results_true.dat b/tests/regression_tests/lattice_hex_coincident/results_true.dat index ba047f674..8bf1c539e 100644 --- a/tests/regression_tests/lattice_hex_coincident/results_true.dat +++ b/tests/regression_tests/lattice_hex_coincident/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.741370E+00 1.384609E-03 +1.744283E+00 1.734930E-03 diff --git a/tests/regression_tests/salphabeta/results_true.dat b/tests/regression_tests/salphabeta/results_true.dat index e1121bce5..fb9d2d77b 100644 --- a/tests/regression_tests/salphabeta/results_true.dat +++ b/tests/regression_tests/salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.403447E-01 2.461538E-02 +8.757212E-01 5.686265E-02 diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index eb06a771c..0c5186e45 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.707485E+00 9.795497E-02 +1.701412E+00 3.180877E-02 diff --git a/tests/regression_tests/void/results_true.dat b/tests/regression_tests/void/results_true.dat index 48be2778a..c0e8184b6 100644 --- a/tests/regression_tests/void/results_true.dat +++ b/tests/regression_tests/void/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.062505E+00 2.674375E-02 +9.612556E-01 1.990135E-02 From e4d7417f981ef25b6caaff4c713f307c375f82e5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 21 May 2019 22:03:23 -0500 Subject: [PATCH 58/61] Fix default constant on Watt spectrum when no source is specified --- src/settings.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/settings.cpp b/src/settings.cpp index 4c041a080..32032ed2a 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -400,7 +400,7 @@ void read_settings_xml() SourceDistribution source { UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})}, UPtrAngle{new Isotropic()}, - UPtrDist{new Watt(0.988, 2.249e-6)} + UPtrDist{new Watt(0.988e6, 2.249e-6)} }; model::external_sources.push_back(std::move(source)); } From 7b4d9d268d08df1a07873d794f4675aecfc75a73 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 06:19:59 -0500 Subject: [PATCH 59/61] Change calculation of thermal elastic/inelastic xs --- src/thermal.cpp | 4 ++-- .../asymmetric_lattice/results_true.dat | 2 +- tests/regression_tests/cmfd_feed_2g/results_true.dat | 12 ++++++------ 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/thermal.cpp b/src/thermal.cpp index 5db578f9d..40bb19e89 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -199,7 +199,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, // Calculate S(a,b) inelastic scattering cross section auto& xs = sab.inelastic_sigma_; - *inelastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1]; + *inelastic = xs[i_grid] + f * (xs[i_grid + 1] - xs[i_grid]); // Check for elastic data if (!sab.elastic_e_in_.empty()) { @@ -231,7 +231,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, // Calculate S(a,b) elastic scattering cross section auto& xs = sab.elastic_P_; - *elastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1]; + *elastic = xs[i_grid] + f*(xs[i_grid + 1] - xs[i_grid]); } } else { // No elastic data diff --git a/tests/regression_tests/asymmetric_lattice/results_true.dat b/tests/regression_tests/asymmetric_lattice/results_true.dat index ecc30a056..43f59aa5d 100644 --- a/tests/regression_tests/asymmetric_lattice/results_true.dat +++ b/tests/regression_tests/asymmetric_lattice/results_true.dat @@ -1 +1 @@ -4b75e203d06d0fc1b4c4dfcb8c180d6f3df8fa2bc44e9775b59bbfd8f7a3785f956f9a9f301526f69c0f8a963ce2e553e0c62c7f6c696656ce1d47b415af6076 \ No newline at end of file +4401f503237c94e9d9cfc9f60e0269d5ae5bb67be3225e18c5510ed08616482964e2962a06268751f66a455fac3ddd5faf91555638dfb56fcd09eee60219edff \ No newline at end of file diff --git a/tests/regression_tests/cmfd_feed_2g/results_true.dat b/tests/regression_tests/cmfd_feed_2g/results_true.dat index 8161b2645..4b3dc3c96 100644 --- a/tests/regression_tests/cmfd_feed_2g/results_true.dat +++ b/tests/regression_tests/cmfd_feed_2g/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.038883E+00 1.017026E-02 +1.038883E+00 1.017030E-02 tally 1: 1.167304E+02 1.362680E+03 @@ -49,7 +49,7 @@ tally 3: 0.000000E+00 2.033842E-02 4.375294E-05 -4.296338E+00 +4.296339E+00 9.280716E-01 3.493776E+00 6.157094E-01 @@ -74,7 +74,7 @@ tally 3: 9.717439E+01 4.724254E+02 8.435691E-01 -3.666151E-02 +3.666152E-02 6.192881E+01 1.924016E+02 0.000000E+00 @@ -82,7 +82,7 @@ tally 3: 1.796382E-02 3.190854E-05 4.343238E+00 -9.514039E-01 +9.514040E-01 3.504683E+00 6.157615E-01 0.000000E+00 @@ -106,7 +106,7 @@ tally 3: 9.885292E+01 4.888720E+02 9.509199E-01 -4.670951E-02 +4.670952E-02 tally 4: 0.000000E+00 0.000000E+00 @@ -383,7 +383,7 @@ cmfd balance 1.08402E-03 1.09178E-03 5.45977E-04 -4.45554E-04 +4.45555E-04 4.01147E-04 3.71025E-04 3.57715E-04 From 8231751b888a7d427282ba501cd884d4e4ead80f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 May 2019 08:24:39 -0500 Subject: [PATCH 60/61] Some updates to setup.py and MANIFEST.in --- MANIFEST.in | 10 ++++++++++ setup.py | 14 ++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/MANIFEST.in b/MANIFEST.in index be82928a1..5715bfb44 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,7 +1,10 @@ include CMakeLists.txt include LICENSE +include CODE_OF_CONDUCT.md +include CONTRIBUTING.md include schemas.xml include pyproject.toml +include pytest.ini include openmc/data/reconstruct.pyx include docs/source/_templates/layout.html include docs/sphinxext/LICENSE @@ -33,4 +36,11 @@ recursive-include tests *.dat recursive-include tests *.h5 recursive-include tests *.py recursive-include tests *.xml +recursive-include vendor CMakeLists.txt +recursive-include vendor *.cmake.in +recursive-include vendor *.cc +recursive-include vendor *.cpp +recursive-include vendor *.hh +recursive-include vendor *.hpp prune docs/build +prune docs/source/pythonapi/generated/ diff --git a/setup.py b/setup.py index 38ee57f71..712244efb 100755 --- a/setup.py +++ b/setup.py @@ -39,7 +39,13 @@ kwargs = { 'author': 'The OpenMC Development Team', 'author_email': 'openmc-dev@googlegroups.com', 'description': 'OpenMC', - 'url': 'https://github.com/openmc-dev/openmc', + 'url': 'https://openmc.org', + 'download_url': 'https://github.com/openmc-dev/openmc/releases', + 'project_urls': { + 'Issue Tracker': 'https://github.com/openmc-dev/openmc/issues', + 'Documentation': 'https://openmc.readthedocs.io', + 'Source Code': 'https://github.com/openmc-dev/openmc', + }, 'classifiers': [ 'Development Status :: 4 - Beta', 'Intended Audience :: Developers', @@ -48,6 +54,7 @@ kwargs = { 'License :: OSI Approved :: MIT License', 'Natural Language :: English', 'Topic :: Scientific/Engineering' + 'Programming Language :: C++', 'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3.4', 'Programming Language :: Python :: 3.5', @@ -55,13 +62,12 @@ kwargs = { 'Programming Language :: Python :: 3.7', ], - # Required dependencies + # Dependencies + 'python_requires': '>=3.4', 'install_requires': [ 'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib', 'pandas', 'lxml', 'uncertainties' ], - - # Optional dependencies 'extras_require': { 'test': ['pytest', 'pytest-cov'], 'vtk': ['vtk'], From d1e051c8570bf8cff82e9142cdb87c7c9b0ec4e4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 May 2019 15:44:02 -0500 Subject: [PATCH 61/61] Update test results for Ubuntu 16.04 --- docs/source/methods/neutron_physics.rst | 2 +- .../cmfd_feed_ng/results_true.dat | 6 +++--- .../lattice_hex_coincident/inputs_true.dat | 6 +++--- .../lattice_hex_coincident/results_true.dat | 2 +- .../lattice_hex_coincident/test.py | 15 ++++++++------- tests/regression_tests/salphabeta/inputs_true.dat | 2 +- .../regression_tests/salphabeta/results_true.dat | 2 +- tests/regression_tests/salphabeta/test.py | 2 +- tests/regression_tests/triso/results_true.dat | 2 +- 9 files changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst index 3451ae83f..465929ba0 100644 --- a/docs/source/methods/neutron_physics.rst +++ b/docs/source/methods/neutron_physics.rst @@ -86,7 +86,7 @@ Elastic Scattering ------------------ Note that the multi-group mode makes no distinction between elastic or -inelastic scattering reactions. The spceific multi-group scattering +inelastic scattering reactions. The specific multi-group scattering implementation is discussed in the :ref:`multi-group-scatter` section. Elastic scattering refers to the process by which a neutron scatters off a diff --git a/tests/regression_tests/cmfd_feed_ng/results_true.dat b/tests/regression_tests/cmfd_feed_ng/results_true.dat index 8308e8634..707411459 100644 --- a/tests/regression_tests/cmfd_feed_ng/results_true.dat +++ b/tests/regression_tests/cmfd_feed_ng/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.029540E+00 1.765356E-02 +1.029540E+00 1.765357E-02 tally 1: 1.144958E+02 1.311468E+03 @@ -93,7 +93,7 @@ tally 3: 0.000000E+00 6.939496E+01 3.014077E+02 -6.294738E-01 +6.294737E-01 2.532146E-02 4.991526E+01 1.566280E+02 @@ -129,7 +129,7 @@ tally 3: 0.000000E+00 6.885265E+01 2.965191E+02 -6.537576E-01 +6.537575E-01 2.783733E-02 4.958341E+01 1.546656E+02 diff --git a/tests/regression_tests/lattice_hex_coincident/inputs_true.dat b/tests/regression_tests/lattice_hex_coincident/inputs_true.dat index 8d73f6374..0388a4f55 100644 --- a/tests/regression_tests/lattice_hex_coincident/inputs_true.dat +++ b/tests/regression_tests/lattice_hex_coincident/inputs_true.dat @@ -39,7 +39,7 @@ - + @@ -67,8 +67,8 @@ eigenvalue 1000 - 10 - 5 + 5 + 2 -0.9899494936611666 -0.9899494936611666 0.0 0.9899494936611666 0.9899494936611666 10.0 diff --git a/tests/regression_tests/lattice_hex_coincident/results_true.dat b/tests/regression_tests/lattice_hex_coincident/results_true.dat index 8bf1c539e..1f95edaf8 100644 --- a/tests/regression_tests/lattice_hex_coincident/results_true.dat +++ b/tests/regression_tests/lattice_hex_coincident/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.744283E+00 1.734930E-03 +1.910374E+00 2.685991E-02 diff --git a/tests/regression_tests/lattice_hex_coincident/test.py b/tests/regression_tests/lattice_hex_coincident/test.py index 65786e2ec..eaea74b19 100644 --- a/tests/regression_tests/lattice_hex_coincident/test.py +++ b/tests/regression_tests/lattice_hex_coincident/test.py @@ -1,4 +1,5 @@ -import numpy as np +from math import sqrt + import openmc from tests.testing_harness import PyAPITestHarness @@ -13,7 +14,7 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness): materials.append(fuel_mat) matrix = openmc.Material() - matrix.set_density('atom/b-cm', 8.7742E-02) + matrix.set_density('atom/b-cm', 1.7742E-02) matrix.add_element('C', 1.0, 'ao') matrix.add_s_alpha_beta('c_Graphite') materials.append(matrix) @@ -94,7 +95,7 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness): coolant_univ.add_cells(coolant_channel) half_width = assembly_pitch # cm - edge_length = (2./np.sqrt(3.0)) * half_width + edge_length = (2./sqrt(3.0)) * half_width inf_mat = openmc.Cell() inf_mat.fill = matrix @@ -132,19 +133,19 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness): settings.run_mode = 'eigenvalue' source = openmc.Source() - corner_dist = np.sqrt(2) * pin_rad + corner_dist = sqrt(2) * pin_rad ll = [-corner_dist, -corner_dist, 0.0] ur = [corner_dist, corner_dist, 10.0] source.space = openmc.stats.Box(ll, ur) source.strength = 1.0 settings.source = source settings.output = {'summary' : False} - settings.batches = 10 - settings.inactive = 5 + settings.batches = 5 + settings.inactive = 2 settings.particles = 1000 settings.seed = 22 settings.export_to_xml() def test_lattice_hex_coincident_surf(): - harness = HexLatticeCoincidentTestHarness('statepoint.10.h5') + harness = HexLatticeCoincidentTestHarness('statepoint.5.h5') harness.main() diff --git a/tests/regression_tests/salphabeta/inputs_true.dat b/tests/regression_tests/salphabeta/inputs_true.dat index ede96d376..29d8065cf 100644 --- a/tests/regression_tests/salphabeta/inputs_true.dat +++ b/tests/regression_tests/salphabeta/inputs_true.dat @@ -49,7 +49,7 @@ eigenvalue - 1000 + 400 5 0 diff --git a/tests/regression_tests/salphabeta/results_true.dat b/tests/regression_tests/salphabeta/results_true.dat index fb9d2d77b..708e6e726 100644 --- a/tests/regression_tests/salphabeta/results_true.dat +++ b/tests/regression_tests/salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.757212E-01 5.686265E-02 +8.474822E-01 1.767966E-02 diff --git a/tests/regression_tests/salphabeta/test.py b/tests/regression_tests/salphabeta/test.py index e69ada5ba..ab5305126 100644 --- a/tests/regression_tests/salphabeta/test.py +++ b/tests/regression_tests/salphabeta/test.py @@ -66,7 +66,7 @@ def make_model(): # Settings model.settings.batches = 5 model.settings.inactive = 0 - model.settings.particles = 1000 + model.settings.particles = 400 model.settings.source = openmc.Source(space=openmc.stats.Box( [-4, -4, -4], [4, 4, 4])) diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index 0c5186e45..c8e683234 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.701412E+00 3.180877E-02 +1.701412E+00 3.180881E-02