From ce54673f1d3afc4e725ed34e2040f3e47b9e2956 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 16 Mar 2019 22:01:41 -0400 Subject: [PATCH 01/43] 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 47bdbf1d96..444136b317 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/43] 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 902d242321..0fbdbc0c9e 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 79345b90e0..86bf4390ac 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/43] 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 d9307abe20..dd9ecf211c 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 eec3a2b66b..9529a31f21 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 df0976d94a..0c7ebcbda3 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 86bf4390ac..20000c4a9a 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 d2405d899e..222822230d 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 682db35405..d26725db22 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 7e235fa458..d574ed109f 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/43] 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 2ced8acf0a..6ab4af84cd 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/43] 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 444136b317..9d453929c8 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/43] 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 9d453929c8..ea0e0cb924 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/43] 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 ea0e0cb924..9c046c525c 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/43] 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 9c046c525c..6c66efc012 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/43] 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 6c66efc012..2722d03247 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/43] 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 2722d03247..663c6a2b7d 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/43] 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 663c6a2b7d..73c48527d1 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/43] 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 d574ed109f..9fd31925f0 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/43] 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 73c48527d1..f8400b736d 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/43] 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 f8400b736d..4e266bc9d7 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/43] 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 4e266bc9d7..b2edebacd2 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/43] 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 9fd31925f0..820b514017 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/43] 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 b2edebacd2..57f68510c6 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/43] 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 57f68510c6..78b48593ee 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/43] 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 820b514017..f249bc1215 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/43] 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 f249bc1215..b0bf8f4483 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/43] 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 78b48593ee..fd6a8fd177 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 b0bf8f4483..8fb85db5db 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/43] fix photon analog --- src/physics.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/physics.cpp b/src/physics.cpp index 1cfd316580..258cea9248 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/43] 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 dd9ecf211c..a6397b2a7a 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 8ec53de447..fcaf209a07 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 4e0066b20e..42889a981b 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 67af6dc421..018677a547 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/43] 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 8fb85db5db..3a759496ba 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/43] 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 018677a547..e8611dbd72 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 d6ccff0f1b..2cb0945ab3 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 3a759496ba..ea046ed042 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/43] 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 ea046ed042..b2e0761dd7 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/43] 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 fb40b0a6b8..974b9323cf 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 edf50f189586fe43c4bbab193a9284f0a072cffe Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Apr 2019 21:04:41 -0500 Subject: [PATCH 28/43] Don't build redundant reactions when calling IncidentNeutron.from_hdf5 --- openmc/data/neutron.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index cd1e406047..7b340f0b0c 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -543,20 +543,6 @@ class IncidentNeutron(EqualityMixin): tgroup = group['total_nu'] rx.derived_products.append(Product.from_hdf5(tgroup)) - # Build redundant reactions. Start from the highest MT number because - # high MTs never depend on lower MTs. - for mt_sum in sorted(SUM_RULES, reverse=True): - if mt_sum not in data: - rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data] - if len(rxs) > 0: - data.reactions[mt_sum] = rx = Reaction(mt_sum) - rx.redundant = True - if rx.mt == 18 and 'total_nu' in group: - tgroup = group['total_nu'] - rx.derived_products.append(Product.from_hdf5(tgroup)) - for T in data.temperatures: - rx.xs[T] = Sum([rx_i.xs[T] for rx_i in rxs]) - # Read unresolved resonance probability tables if 'urr' in group: urr_group = group['urr'] From 63527fe9087982b3b34bbd2e76ec19e666552392 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 6 Apr 2019 08:02:07 -0400 Subject: [PATCH 29/43] 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 b2e0761dd7..e7cde6bd70 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 30/43] 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 974b9323cf..aa01eac9b4 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 31/43] 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 e7cde6bd70..578aa004c4 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 32/43] 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 107380eb84..f0e6dde015 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 e8611dbd72..9b4a1a9c5f 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 33/43] 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 a6397b2a7a..49f980cc29 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 fcaf209a07..62923c11f1 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 f0e6dde015..107380eb84 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 42889a981b..cbc6217caa 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 9b4a1a9c5f..cd29e97c26 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 34/43] 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 35c25e5480..db7635d0a5 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 2216b37daf..5d8d18c219 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 de0d7e5de5..74d777fbd6 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 cd29e97c26..0eb8c494e7 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 c1a203e48b..40a20f2321 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 578aa004c4..74cf9dcd44 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 35/43] 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 81e4f4bcd7..957bfa3a3c 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 6ab4af84cd..f005b50606 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 36/43] 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 141838e49a..cd1d8c05cc 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 6eb35341a593fe878175809725d932307f07eee9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 12 Apr 2019 06:56:33 -0500 Subject: [PATCH 37/43] Fix bug re: first nuclide with zero density --- src/material.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/material.cpp b/src/material.cpp index de0d7e5de5..7d12578449 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -356,8 +356,8 @@ void Material::finalize() void Material::normalize_density() { - bool percent_in_atom = (atom_density_(0) > 0.0); - bool density_in_atom = (density_ > 0.0); + bool percent_in_atom = (atom_density_(0) >= 0.0); + bool density_in_atom = (density_ >= 0.0); for (int i = 0; i < nuclide_.size(); ++i) { // determine atomic weight ratio From acc69ea6cb6891076fcc9424b4f1f7a0e087d1a3 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 12 Apr 2019 20:38:36 -0400 Subject: [PATCH 38/43] 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 1d43c89137..f82b4a6a78 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 7df7e5dbcf..a7725576c8 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 47215fff30..c6a1fc3908 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 b4a6d796d56f5a9202ad24ee7bd1efeef7e75eab Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Apr 2019 07:25:29 -0500 Subject: [PATCH 39/43] Add two papers to publication list --- docs/source/publications.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/docs/source/publications.rst b/docs/source/publications.rst index e1c0f6de00..d554fe5a14 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -139,6 +139,18 @@ Miscellaneous `_," *Nucl. Sci. Tech.*, **30**: 63 (2019). +- M. Sajjad, Muhammad Rizwan Ali, M. Naveed Ashraf, Rustam Khan, Tasneem Fatima, + "`KANUPP Reactor Core Model and its Validation + `_," International Conference on + Power Generation Systems and Renewable Energy Technologies, Islamabad, + Pakistan, Sep. 10-12 (2018). + +- Muhammad Waqas Tariq, Muhammad Sohail, and Sikander Majid Mirza, "`Calculation + of Neutronic Parameters using OpenMC for Potential Dispersed Fuels of MNSR + `_," International Conference on + Power Generation Systems and Renewable Energy Technologies, Islamabad, + Pakistan, Sep. 10-12 (2018). + - Amanda L. Lund and Paul K. Romano, "`Implementation and Validation of Photon Transport in OpenMC `_", Argonne National Laboratory, Technical Report ANL/MCS-TM-381 (2018). From 435b21eecb92d759f212954942f158404f0b1cea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 17 Apr 2019 10:40:11 -0500 Subject: [PATCH 40/43] 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 cd1d8c05cc..1de84721ee 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 41/43] 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 db7635d0a5..a60fd06a82 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 1de84721ee..da576b5c8f 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 40a20f2321..614e38767b 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 4d77f6fc943885d87692d5ce6f874af2635c9379 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 19 Apr 2019 14:00:47 -0400 Subject: [PATCH 42/43] 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 f167a3c236..c40e2a699c 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 43/43] 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 ce4d12c1b2..38a5bc4692 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