diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index 63e7028aa2..e77c70b1c5 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -313,11 +313,15 @@ Functions :return: Return status (negative if an error occurs) :rtype: int -.. c:function:: int openmc_load_nuclide(char name[]) +.. c:function:: int openmc_load_nuclide(const char* name, const double* temps, int n) Load data for a nuclide from the HDF5 data library. - :param char[] name: Name of the nuclide. + :param name: Name of the nuclide. + :type name: const char* + :param temps: Temperatures in [K] to load data at + :type temps: const double* + :param int n: Number of temperatures :return: Return status (negative if an error occurs) :rtype: int @@ -373,7 +377,7 @@ Functions :return: Return status (negative if an error occurs) :rtype: int -.. c:function:: int openmc_material_set_densities(int32_t index, int n, const char** name, const double density*) +.. c:function:: int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density) :param int32_t index: Index in the materials array :param int n: Length of name/density diff --git a/docs/source/publications.rst b/docs/source/publications.rst index a6746741e3..89df0d29f8 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -121,6 +121,9 @@ Coupling and Multi-physics Geometry and Visualization -------------------------- +- Patrick C. Shriwise, Xiaokang Zhang, and Andrew Davis, "DAG-OpenMC: CAD-Based + Geometry in OpenMC", *Trans. Am. Nucl. Soc.*, **122**, 395-398 (2020). + - Sterling Harper, Paul Romano, Benoit Forget, and Kord Smith, "Efficient dynamic threadsafe neighbor lists for Monte Carlo ray tracing," *Proc. M&C*, 918-926, Portland, Oregon, Aug. 25-29 (2019). @@ -146,10 +149,29 @@ Geometry and Visualization Miscellaneous ------------- -- Ned Xoubi, Sharif Abu Darda, Abdelfattah Y. Soliman, and Tareq Abulfaraj, +- Jiankai Yu, Qiudong Wang, Ding She, and Benoit Forget, "Modelling of the + HTR-PM Pebble-bed Reactor using OpenMC", *Trans. Am. Nucl. Soc.*, **122**, + 643-646 (2020). + +- Sharif Abu Darda, Abdelfattah Y. Soliman, Mohammed S. Aljohani, and Ned Xoubi, + "`Technical feasibility study of BAEC TRIGA reactor (BTRR) as a neutron source + for BNCT using OpenMC Monte Carlo code + `_", *Prog. Nucl. Energy*, + **126**, 103418 (2020). + +- Stefano Segantin, Raffaella Testoni, and Massimo Zucchetti, "`ARC reactor -- + Neutron irradiation analysis `_", + *Fus. Eng. Design*, **159**, 111792 (2020). + +- Muhammad Ilham, Helen Raflis, and Zaki Suud, "`Full Core Optimization of Small + Modular Gas-Cooled Fast Reactors Using OpenMC Program Code + `_", *J. Phys.: Conf. Series*, + **1493**, 012007 (2020). + +- Ned Xoubi, Sharif Abu Darda, Abdelfattah Y. Soliman, and Tareq Abulfaraj, "`An investigative study of enrichment reduction impact on the neutron flux in the in-core flux-trap facility of MTR research reactors - `_", Nucl. Eng. Technol., **52**, + `_", *Nucl. Eng. Technol.*, **52**, 469-476 (2020). - J. Rolando Granada, J. Ignacio Marquez Damian, and Christian Helman, "`Studies @@ -162,6 +184,10 @@ Miscellaneous `_," M.S. Thesis, KTH Royal Institute of Technology (2019). +- Ilham Variansyah, Benjamin R. Betzler, and William R. Martin, + "α-weighted transition rate matrix method", *Proc. M&C*, 1368-1377, Portland, + Oregon, Aug. 25-29 (2019). + - Shikhar Kumar, Benoit Forget, and Kord Smith, "Analysis of fission source convergence for a 3-D SMR core using functional expansion tallies," *Proc. M&C*, 937-947, Portland, Oregon, Aug. 25-29 (2019). @@ -249,6 +275,16 @@ Miscellaneous Multigroup Cross Section Generation ----------------------------------- +- Ilham Variansyah, Benjamin R. Betzler, and William R. Martin, "`Multigroup + Constant Calculation with Static α-Eigenvalue Monte Carlo for + Time-Dependent Neutron Transport Simulation + `_", *Nucl. Sci. Eng.*, 2020. + +- Chenghui Wan, Tianliang Hu, and Liangzhi Cao, "`Multi-physics numerical + analysis of the fuel-addition transients in the liquid-fuel molten salt reactor + `_", *Ann. Nucl. Energy*, + **144**, 107514 (2020). + - William Boyd, Adam Nelson, Paul K. Romano, Samuel Shaner, Benoit Forget, and Kord Smith, "`Multigroup Cross-Section Generation with the OpenMC Monte Carlo Particle Transport Code `_," diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 0d2432d290..6d82929d17 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -55,7 +55,7 @@ extern "C" { bool openmc_is_statepoint_batch(); int openmc_legendre_filter_get_order(int32_t index, int* order); int openmc_legendre_filter_set_order(int32_t index, int order); - int openmc_load_nuclide(const char* name); + int openmc_load_nuclide(const char* name, const double* temps, int n); int openmc_material_add_nuclide(int32_t index, const char name[], double density); int openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n); int openmc_material_get_id(int32_t index, int32_t* id); diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 597dbb11b7..60bac81d2d 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -9,6 +9,7 @@ #include #include +#include #include #include "openmc/constants.h" @@ -34,8 +35,9 @@ public: std::vector energy; }; - // Constructors - Nuclide(hid_t group, const std::vector& temperature, int i_nuclide); + // Constructors/destructors + Nuclide(hid_t group, const std::vector& temperature); + ~Nuclide(); //! Initialize logarithmic grid for energy searches void init_grid(); @@ -62,7 +64,7 @@ public: int A_; //!< Mass number int metastable_; //!< Metastable state double awr_; //!< Atomic weight ratio - int i_nuclide_; //!< Index in the nuclides array + gsl::index index_; //!< Index in the nuclides array // Temperature dependent cross section data std::vector kTs_; //!< temperatures in eV (k*T) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index b652a4d8b7..0b3688bbae 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -4,9 +4,11 @@ #include "openmc/endf.h" #include "openmc/particle.h" +#include #include #include "xtensor/xtensor.hpp" +#include // for unique_ptr #include #include #include // for pair @@ -38,8 +40,9 @@ public: class PhotonInteraction { public: - // Constructors - PhotonInteraction(hid_t group, int i_element); + // Constructors/destructor + PhotonInteraction(hid_t group); + ~PhotonInteraction(); // Methods void calculate_xs(Particle& p) const; @@ -57,7 +60,7 @@ public: // Data members std::string name_; //!< Name of element, e.g. "Zr" int Z_; //!< Atomic number - int i_element_; //!< Index in global elements vector + gsl::index index_; //!< Index in global elements vector // Microscopic cross sections xt::xtensor energy_; @@ -117,7 +120,7 @@ namespace data { extern xt::xtensor compton_profile_pz; //! Compton profile momentum grid //! Photon interaction data for each element -extern std::vector elements; +extern std::vector> elements; extern std::unordered_map element_map; } // namespace data diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index a3032ee010..32e2c67303 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -57,6 +57,9 @@ void allocate_banks(); //! Determine number of particles to transport per process void calculate_work(); +//! Initialize nuclear data before a simulation +void initialize_data(); + //! Initialize a batch void initialize_batch(); diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 27e253d306..040a28cbfc 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -965,6 +965,8 @@ class Chain: new_nuclide = Nuclide(previous.name) new_nuclide.half_life = previous.half_life new_nuclide.decay_energy = previous.decay_energy + if hasattr(previous, '_fpy'): + new_nuclide._fpy = previous._fpy new_decay = [] for mode in previous.decay_modes: diff --git a/openmc/lib/nuclide.py b/openmc/lib/nuclide.py index f9a4c1fadd..571dcc6551 100644 --- a/openmc/lib/nuclide.py +++ b/openmc/lib/nuclide.py @@ -1,5 +1,5 @@ from collections.abc import Mapping -from ctypes import c_int, c_char_p, POINTER, c_size_t +from ctypes import c_int, c_double, c_char_p, POINTER, c_size_t from weakref import WeakValueDictionary from ..exceptions import DataError, AllocationError @@ -14,7 +14,7 @@ __all__ = ['Nuclide', 'nuclides', 'load_nuclide'] _dll.openmc_get_nuclide_index.argtypes = [c_char_p, POINTER(c_int)] _dll.openmc_get_nuclide_index.restype = c_int _dll.openmc_get_nuclide_index.errcheck = _error_handler -_dll.openmc_load_nuclide.argtypes = [c_char_p] +_dll.openmc_load_nuclide.argtypes = [c_char_p, POINTER(c_double), c_int] _dll.openmc_load_nuclide.restype = c_int _dll.openmc_load_nuclide.errcheck = _error_handler _dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)] @@ -32,7 +32,7 @@ def load_nuclide(name): Name of the nuclide, e.g. 'U235' """ - _dll.openmc_load_nuclide(name.encode()) + _dll.openmc_load_nuclide(name.encode(), None, 0) class Nuclide(_FortranObject): diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 1d36b5cff4..6267dc5f9b 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -1,5 +1,6 @@ #include "openmc/cross_sections.h" +#include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/container_util.h" #ifdef DAGMC @@ -207,78 +208,11 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // If we've already read this nuclide, skip it if (already_read.find(name) != already_read.end()) continue; - LibraryKey key {Library::Type::neutron, name}; - int idx = data::library_map[key]; - std::string& filename = data::libraries[idx].path_; + const auto& temps = nuc_temps[i_nuc]; + int err = openmc_load_nuclide(name.c_str(), temps.data(), temps.size()); + if (err < 0) throw std::runtime_error{openmc_err_msg}; - write_message("Reading " + name + " from " + filename, 6); - - // Open file and make sure version is sufficient - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read nuclide data from HDF5 - hid_t group = open_group(file_id, name.c_str()); - int i_nuclide = data::nuclides.size(); - data::nuclides.push_back(std::make_unique( - group, nuc_temps[i_nuc], i_nuclide)); - - close_group(group); - file_close(file_id); - - // Determine if minimum/maximum energy for this nuclide is greater/less - // than the previous - if (data::nuclides[i_nuclide]->grid_.size() >= 1) { - int neutron = static_cast(Particle::Type::neutron); - data::energy_min[neutron] = std::max(data::energy_min[neutron], - data::nuclides[i_nuclide]->grid_[0].energy.front()); - data::energy_max[neutron] = std::min(data::energy_max[neutron], - data::nuclides[i_nuclide]->grid_[0].energy.back()); - } - - // Add name and alias to dictionary already_read.insert(name); - - // Check if elemental data has been read, if needed - 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 - LibraryKey key {Library::Type::photon, element}; - int idx = data::library_map[key]; - std::string& filename = data::libraries[idx].path_; - write_message("Reading " + element + " from " + filename, 6); - - // Open file and make sure version is sufficient - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read element data from HDF5 - hid_t group = open_group(file_id, element.c_str()); - data::elements.emplace_back(group, data::elements.size()); - - // Determine if minimum/maximum energy for this element is greater/less than - // the previous - const auto& elem {data::elements.back()}; - if (elem.energy_.size() >= 1) { - int photon = static_cast(Particle::Type::photon); - int n = elem.energy_.size(); - data::energy_min[photon] = std::max(data::energy_min[photon], - std::exp(elem.energy_(1))); - data::energy_max[photon] = std::min(data::energy_max[photon], - std::exp(elem.energy_(n - 1))); - } - - close_group(group); - file_close(file_id); - - // Add element to set - already_read.insert(element); - } - } - - // Read multipole file into the appropriate entry on the nuclides array - if (settings::temperature_multipole) read_multipole_data(i_nuclide); } } @@ -316,49 +250,11 @@ read_ce_cross_sections(const std::vector>& nuc_temps, mat->finalize(); } // materials - - // Set up logarithmic grid for nuclides - for (auto& nuc : data::nuclides) { - nuc->init_grid(); - } - int neutron = static_cast(Particle::Type::neutron); - simulation::log_spacing = std::log(data::energy_max[neutron] / - data::energy_min[neutron]) / settings::n_log_bins; - if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) { - // Determine if minimum/maximum energy for bremsstrahlung is greater/less - // than the current minimum/maximum - if (data::ttb_e_grid.size() >= 1) { - int photon = static_cast(Particle::Type::photon); - int n_e = data::ttb_e_grid.size(); - data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); - data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); - } - // Take logarithm of energies since they are log-log interpolated data::ttb_e_grid = xt::log(data::ttb_e_grid); } - // Show which nuclide results in lowest energy for neutron transport - for (const auto& nuc : data::nuclides) { - // If a nuclide is present in a material that's not used in the model, its - // grid has not been allocated - if (nuc->grid_.size() > 0) { - double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(Particle::Type::neutron); - if (max_E == data::energy_max[neutron]) { - write_message("Maximum neutron transport energy: " + - std::to_string(data::energy_max[neutron]) + " eV for " + - nuc->name_, 7); - if (mpi::master && data::energy_max[neutron] < 20.0e6) { - warning("Maximum neutron energy is below 20 MeV. This may bias " - " the results."); - } - break; - } - } - } - // Show minimum/maximum temperature write_message("Minimum neutron data temperature: " + std::to_string(data::temperature_min) + " K", 4); diff --git a/src/material.cpp b/src/material.cpp index 6ee8749351..ea5571df8d 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -475,7 +475,7 @@ void Material::collision_stopping_power(double* s_col, bool positron) std::vector e_b_sq; for (int i = 0; i < element_.size(); ++i) { - const auto& elm = data::elements[element_[i]]; + const auto& elm = *data::elements[element_[i]]; double awr = data::nuclides[nuclide_[i]]->awr_; // Get atomic density of nuclide given atom/weight percent @@ -589,7 +589,7 @@ void Material::init_bremsstrahlung() // Bragg's additivity rule. for (int i = 0; i < n; ++i) { // Get pointer to current element - const auto& elm = data::elements[element_[i]]; + const auto& elm = *data::elements[element_[i]]; double awr = data::nuclides[nuclide_[i]]->awr_; // Get atomic density and mass density of nuclide given atom/weight percent @@ -838,7 +838,7 @@ void Material::calculate_photon_xs(Particle& p) const // Calculate microscopic cross section for this nuclide const auto& micro {p.photon_xs_[i_element]}; if (p.E_ != micro.last_E) { - data::elements[i_element].calculate_xs(p); + data::elements[i_element]->calculate_xs(p); } // ======================================================================== @@ -936,13 +936,14 @@ void Material::set_densities(const std::vector& name, if (n != nuclide_.size()) { nuclide_.resize(n); atom_density_ = xt::zeros({n}); + if (settings::photon_transport) element_.resize(n); } double sum_density = 0.0; for (gsl::index i = 0; i < n; ++i) { const auto& nuc {name[i]}; if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) { - int err = openmc_load_nuclide(nuc.c_str()); + int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0); if (err < 0) throw std::runtime_error{openmc_err_msg}; } @@ -950,11 +951,21 @@ void Material::set_densities(const std::vector& name, Expects(density[i] > 0.0); atom_density_(i) = density[i]; sum_density += density[i]; + + if (settings::photon_transport) { + auto element_name = to_element(nuc); + element_[i] = data::element_map.at(element_name); + } } // Set total density to the sum of the vector this->set_density(sum_density, "atom/b-cm"); + // Generate material bremsstrahlung data for electrons and positrons + if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) { + this->init_bremsstrahlung(); + } + // Assign S(a,b) tables this->init_thermal(); } @@ -1042,13 +1053,19 @@ void Material::add_nuclide(const std::string& name, double density) } // If nuclide wasn't found, extend nuclide/density arrays - int err = openmc_load_nuclide(name.c_str()); + int err = openmc_load_nuclide(name.c_str(), nullptr, 0); if (err < 0) throw std::runtime_error{openmc_err_msg}; // Append new nuclide/density int i_nuc = data::nuclide_map[name]; nuclide_.push_back(i_nuc); + // Append new element if photon transport is on + if (settings::photon_transport) { + int i_elem = data::element_map[to_element(name)]; + element_.push_back(i_elem); + } + auto n = nuclide_.size(); // Create copy of atom_density_ array with one extra entry diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 22fd568554..4e5d41b4be 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -46,11 +46,14 @@ int Nuclide::XS_FISSION {2}; int Nuclide::XS_NU_FISSION {3}; int Nuclide::XS_PHOTON_PROD {4}; -Nuclide::Nuclide(hid_t group, const std::vector& temperature, int i_nuclide) - : i_nuclide_{i_nuclide} +Nuclide::Nuclide(hid_t group, const std::vector& temperature) { + // Set index of nuclide in global vector + index_ = data::nuclides.size(); + // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); + data::nuclide_map[name_] = index_; read_attribute(group, "Z", Z_); read_attribute(group, "A", A_); @@ -276,6 +279,11 @@ Nuclide::Nuclide(hid_t group, const std::vector& temperature, int i_nucl this->create_derived(prompt_photons.get(), delayed_photons.get()); } +Nuclide::~Nuclide() +{ + data::nuclide_map.erase(name_); +} + void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons) { for (const auto& grid : grid_) { @@ -489,7 +497,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const void Nuclide::calculate_elastic_xs(Particle& p) const { // Get temperature index, grid index, and interpolation factor - auto& micro {p.neutron_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[index_]}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -525,7 +533,7 @@ double Nuclide::elastic_xs_0K(double E) const void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p) { - auto& micro {p.neutron_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[index_]}; // Initialize cached cross sections to zero micro.elastic = CACHE_INVALID; @@ -732,7 +740,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) { - auto& micro {p.neutron_xs_[i_nuclide_]}; + auto& micro {p.neutron_xs_[index_]}; // Set flag that S(a,b) treatment should be used for scattering micro.index_sab = i_sab; @@ -761,7 +769,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const { - auto& micro = p.neutron_xs_[i_nuclide_]; + auto& micro = p.neutron_xs_[index_]; micro.use_ptable = true; // Create a shorthand for the URR data @@ -779,8 +787,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const // therefore preserving correlation of temperature in probability tables. p.stream_ = STREAM_URR_PTABLE; //TODO: to maintain the same random number stream as the Fortran code this - //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_ - double r = future_prn(static_cast(i_nuclide_ + 1), *p.current_seed()); + //replaces, the seed is set with index_ + 1 instead of index_ + double r = future_prn(static_cast(index_ + 1), *p.current_seed()); p.stream_ = STREAM_TRACKING; int i_low = 0; @@ -919,42 +927,67 @@ nuclides_size() // C API //============================================================================== -extern "C" int openmc_load_nuclide(const char* name) +extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n) { - if (data::nuclide_map.find(name) == data::nuclide_map.end()) { - const auto& it = data::library_map.find({Library::Type::neutron, name}); - if (it != data::library_map.end()) { - // Get filename for library containing nuclide - int idx = it->second; - std::string& filename = data::libraries[idx].path_; - write_message("Reading " + std::string{name} + " from " + filename, 6); - - // Open file and make sure version is sufficient - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read nuclide data from HDF5 - hid_t group = open_group(file_id, name); - std::vector temperature; - int i_nuclide = data::nuclides.size(); - data::nuclides.push_back(std::make_unique( - group, temperature, i_nuclide)); - - close_group(group); - file_close(file_id); - - // Add entry to nuclide dictionary - data::nuclide_map[name] = i_nuclide; - - // Initialize nuclide grid - data::nuclides.back()->init_grid(); - - // Read multipole file into the appropriate entry on the nuclides array - if (settings::temperature_multipole) read_multipole_data(i_nuclide); - } else { + if (data::nuclide_map.find(name) == data::nuclide_map.end() || + data::nuclide_map.at(name) >= data::elements.size()) { + LibraryKey key {Library::Type::neutron, name}; + const auto& it = data::library_map.find(key); + if (it == data::library_map.end()) { set_errmsg("Nuclide '" + std::string{name} + "' is not present in library."); return OPENMC_E_DATA; } + + // Get filename for library containing nuclide + int idx = it->second; + const auto& filename = data::libraries[idx].path_; + write_message("Reading " + std::string{name} + " from " + filename, 6); + + // Open file and make sure version is sufficient + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read nuclide data from HDF5 + hid_t group = open_group(file_id, name); + std::vector temperature{temps, temps + n}; + data::nuclides.push_back(std::make_unique(group, temperature)); + + close_group(group); + file_close(file_id); + + // Read multipole file into the appropriate entry on the nuclides array + int i_nuclide = data::nuclide_map.at(name); + if (settings::temperature_multipole) read_multipole_data(i_nuclide); + + // Read elemental data, if necessary + if (settings::photon_transport) { + auto element = to_element(name); + if (data::element_map.find(element) == data::element_map.end() || + data::element_map.at(element) >= data::elements.size()) { + // Read photon interaction data from HDF5 photon library + LibraryKey key {Library::Type::photon, element}; + const auto& it = data::library_map.find(key); + if (it == data::library_map.end()) { + set_errmsg("Element '" + std::string{element} + "' is not present in library."); + return OPENMC_E_DATA; + } + + int idx = it->second; + const auto& filename = data::libraries[idx].path_; + write_message("Reading " + element + " from " + filename, 6); + + // Open file and make sure version is sufficient + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read element data from HDF5 + hid_t group = open_group(file_id, element.c_str()); + data::elements.push_back(std::make_unique(group)); + + close_group(group); + file_close(file_id); + } + } } return 0; } diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 3347030baa..b0f739bd06 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -74,6 +74,9 @@ void run_particle_restart() // Set verbosity high settings::verbosity = 10; + // Initialize nuclear data (energy limits, log grid, etc.) + initialize_data(); + // Initialize the particle to be tracked Particle p; diff --git a/src/photon.cpp b/src/photon.cpp index 7af8687dfb..fa142609af 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -27,7 +27,7 @@ namespace data { xt::xtensor compton_profile_pz; -std::vector elements; +std::vector> elements; std::unordered_map element_map; } // namespace data @@ -36,11 +36,14 @@ std::unordered_map element_map; // PhotonInteraction implementation //============================================================================== -PhotonInteraction::PhotonInteraction(hid_t group, int i_element) - : i_element_{i_element} +PhotonInteraction::PhotonInteraction(hid_t group) { + // Set index of element in global vector + index_ = data::elements.size(); + // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); + data::element_map[name_] = index_; // Get atomic number read_attribute(group, "Z", Z_); @@ -294,6 +297,11 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) heating_ = xt::where(heating_ > 0.0, xt::log(heating_), -500.0); } +PhotonInteraction::~PhotonInteraction() +{ + data::element_map.erase(name_); +} + void PhotonInteraction::compton_scatter(double alpha, bool doppler, double* alpha_out, double* mu, int* i_shell, uint64_t* seed) const { @@ -464,7 +472,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const // calculate interpolation factor double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); - auto& xs {p.photon_xs_[i_element_]}; + auto& xs {p.photon_xs_[index_]}; xs.index_grid = i_grid; xs.interp_factor = f; diff --git a/src/physics.cpp b/src/physics.cpp index b307d9487e..1000c3c4b0 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -258,7 +258,7 @@ void sample_photon_reaction(Particle& p) // Sample element within material int i_element = sample_element(p); const auto& micro {p.photon_xs_[i_element]}; - const auto& element {data::elements[i_element]}; + const auto& element {*data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent double alpha = p.E_/MASS_ELECTRON_EV; diff --git a/src/simulation.cpp b/src/simulation.cpp index d3e0f48c71..1b39c13296 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -34,6 +34,7 @@ #endif #include +#include #include @@ -68,6 +69,11 @@ int openmc_simulation_init() // Skip if simulation has already been initialized if (simulation::initialized) return 0; + // Initialize nuclear data (energy limits, log grid) + if (settings::run_CE) { + initialize_data(); + } + // Determine how much work each process should do calculate_work(); @@ -555,6 +561,76 @@ void calculate_work() } } +void initialize_data() +{ + // Determine minimum/maximum energy for incident neutron/photon data + data::energy_max = {INFTY, INFTY}; + data::energy_min = {0.0, 0.0}; + for (const auto& nuc : data::nuclides) { + if (nuc->grid_.size() >= 1) { + int neutron = static_cast(Particle::Type::neutron); + data::energy_min[neutron] = std::max(data::energy_min[neutron], + nuc->grid_[0].energy.front()); + data::energy_max[neutron] = std::min(data::energy_max[neutron], + nuc->grid_[0].energy.back()); + } + } + + if (settings::photon_transport) { + for (const auto& elem : data::elements) { + if (elem->energy_.size() >= 1) { + int photon = static_cast(Particle::Type::photon); + int n = elem->energy_.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(elem->energy_(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(elem->energy_(n - 1))); + } + } + + if (settings::electron_treatment == ElectronTreatment::TTB) { + // Determine if minimum/maximum energy for bremsstrahlung is greater/less + // than the current minimum/maximum + if (data::ttb_e_grid.size() >= 1) { + int photon = static_cast(Particle::Type::photon); + int n_e = data::ttb_e_grid.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(data::ttb_e_grid(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(data::ttb_e_grid(n_e - 1))); + } + } + } + + // Show which nuclide results in lowest energy for neutron transport + for (const auto& nuc : data::nuclides) { + // If a nuclide is present in a material that's not used in the model, its + // grid has not been allocated + if (nuc->grid_.size() > 0) { + double max_E = nuc->grid_[0].energy.back(); + int neutron = static_cast(Particle::Type::neutron); + if (max_E == data::energy_max[neutron]) { + write_message("Maximum neutron transport energy: " + + std::to_string(data::energy_max[neutron]) + " eV for " + + nuc->name_, 7); + if (mpi::master && data::energy_max[neutron] < 20.0e6) { + warning("Maximum neutron energy is below 20 MeV. This may bias " + "the results."); + } + break; + } + } + } + + // Set up logarithmic grid for nuclides + for (auto& nuc : data::nuclides) { + nuc->init_grid(); + } + int neutron = static_cast(Particle::Type::neutron); + simulation::log_spacing = std::log(data::energy_max[neutron] / + data::energy_min[neutron]) / settings::n_log_bins; +} + #ifdef OPENMC_MPI void broadcast_results() { // Broadcast tally results so that each process has access to results