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/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). diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 2ced8acf0a..f005b50606 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -256,9 +256,12 @@ 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 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 | diff --git a/include/openmc/constants.h b/include/openmc/constants.h index d9307abe20..49f980cc29 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_NONE {-1}; // ============================================================================ // CROSS SECTION RELATED CONSTANTS @@ -127,6 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2}; // Reaction types // TODO: Convert to enum +constexpr int REACTION_NONE {0}; constexpr int TOTAL_XS {1}; constexpr int ELASTIC {2}; constexpr int N_NONELASTIC {3}; @@ -230,7 +232,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}; @@ -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}; @@ -366,6 +369,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/include/openmc/particle.h b/include/openmc/particle.h index 8ec53de447..62923c11f1 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 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/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/include/openmc/string_utils.h b/include/openmc/string_utils.h index 35c25e5480..a60fd06a82 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(const std::string& name); + void to_lower(std::string& str); int word_count(std::string const& str); 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/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) 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'] diff --git a/openmc/data/photon.py b/openmc/data/photon.py index e038ca05d2..da576b5c8f 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -3,12 +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 openmc.mixin import EqualityMixin import openmc.checkvalue as cv @@ -19,6 +21,14 @@ 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 +# classical electron radius in cm +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', 'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3', @@ -33,6 +43,7 @@ _REACTION_NAME = { 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'), @@ -150,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): @@ -379,6 +391,49 @@ 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 + + Parameters + ---------- + shell : str + The subshell to compute + + Returns + ------- + float + Energy of fluorescent photons + + """ + + 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] + 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 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 the photoelectron hole + e_row += self.energy_fluorescence(primary) + + # Expected fluorescent photon energy + e += e_row * prob + + self._e_fluorescence[shell] = e + return e class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -499,9 +554,14 @@ 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 [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)]) + # Compton profiles n_shell = ace.nxs[5] if n_shell != 0: @@ -631,6 +691,9 @@ class IncidentPhoton(EqualityMixin): # Add bremsstrahlung DCS data data._add_bremsstrahlung() + # Add heating cross sections + data._compute_heating() + return data @classmethod @@ -751,7 +814,7 @@ class IncidentPhoton(EqualityMixin): designators = [] for mt, rx in self.reactions.items(): name, key = _REACTION_NAME[mt] - if mt in [502, 504, 515, 517, 522]: + if mt in (502, 504, 515, 517, 522, 525): sub_group = group.create_group(key) elif mt >= 534 and mt <= 572: # Subshell @@ -858,6 +921,80 @@ class IncidentPhoton(EqualityMixin): self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy'] self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number]) + 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: + + .. 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}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV + + \overline{E}_{PE} &= E(\text{fluorescent photons}) + + 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 in self: + 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): + k = E / MASS_ELECTRON_EV + 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.0 - mu)) + return Eout * dsigma_dmu(mu, E) + + def eout_average(E): + 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) + 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: + # 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]) + self.reactions[525] = heat_rx class PhotonReaction(EqualityMixin): """Photon-induced reaction @@ -972,14 +1109,21 @@ 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 == 525: + # Get heating factors in [eV per collision] + xs *= EV_PER_MEV + else: + nonzero = (xs != 0.0) + xs[nonzero] = np.exp(xs[nonzero]) rx.xs = Tabulated1D(energy, xs, [n], [5]) # Get form factors for incoherent/coherent scattering 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..b5cb2c47a7 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}; @@ -356,8 +355,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 diff --git a/src/output.cpp b/src/output.cpp index ee74a27e15..bbc2588ac4 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -617,6 +617,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/particle.cpp b/src/particle.cpp index 4e0066b20e..cbc6217caa 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_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 // 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/photon.cpp b/src/photon.cpp index 8a0ddf9aa9..929822bfcf 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -96,6 +96,15 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) read_dataset(rgroup, "xs", photoelectric_total_); close_group(rgroup); + // Read heating + 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"); std::vector designators; @@ -280,6 +289,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, diff --git a/src/physics.cpp b/src/physics.cpp index 67af6dc421..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" @@ -61,12 +62,16 @@ 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 { + } else if (p->type_ == Particle::Type::photon) { msg << " " << reaction_name(p->event_mt_) << " with " << - data::elements[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."; } write_message(msg, 1); } @@ -209,7 +214,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]}; @@ -226,6 +230,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 +273,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 +325,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 +351,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 +369,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 +395,7 @@ void sample_positron_reaction(Particle* p) p->E_ = 0.0; p->alive_ = false; + p->event_ = EVENT_ABSORB; } int sample_nuclide(const Particle* p) @@ -432,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 for tally purpose + p->event_nuclide_ = mat->nuclide_[i]; + + return i_element; + } } // If we made it here, no element was sampled 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/string_utils.cpp b/src/string_utils.cpp index c1a203e48b..614e38767b 100644 --- a/src/string_utils.cpp +++ b/src/string_utils.cpp @@ -26,6 +26,12 @@ char* strtrim(char* c_str) } +std::string to_element(const std::string& 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.cpp b/src/tallies/tally.cpp index 84ec223ed9..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" @@ -109,6 +110,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 +211,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; @@ -771,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 8d51aa7f28..a69ab8ca5c 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -7,10 +7,12 @@ #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" #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" @@ -1149,6 +1151,156 @@ score_general_ce(Particle* p, int i_tally, int start_index, break; + case SCORE_HEATING: + score = 0.; + if (p->type_ == Particle::Type::neutron) { + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // 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 + score = p->wgt_last_ + p->wgt_absorb_; + } else { + score = p->wgt_last_; + } + 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]); + } + } + 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) { + 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) { + 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_[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) { + 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) { + // 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]; + 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 { + // 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_; + 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; + 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; + } + } + } + } + } + break; + default: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need @@ -1979,12 +2131,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 { 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): 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