mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #1595 from paulromano/depletion-photon-fix
Fix for depletion with photon transport turned on
This commit is contained in:
Amanda Lund
GitHub
commit
ba2563f868
15 changed files with 257 additions and 174 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
<https://doi.org/10.1016/j.pnucene.2020.103418>`_", *Prog. Nucl. Energy*,
|
||||
**126**, 103418 (2020).
|
||||
|
||||
- Stefano Segantin, Raffaella Testoni, and Massimo Zucchetti, "`ARC reactor --
|
||||
Neutron irradiation analysis <https://doi.org/10.1016/j.fusengdes.2020.111792>`_",
|
||||
*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
|
||||
<https://doi.org/10.1088/1742-6596/1493/1/012007>`_", *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
|
||||
<https://doi.org/10.1016/j.net.2019.08.008>`_", Nucl. Eng. Technol., **52**,
|
||||
<https://doi.org/10.1016/j.net.2019.08.008>`_", *Nucl. Eng. Technol.*, **52**,
|
||||
469-476 (2020).
|
||||
|
||||
- J. Rolando Granada, J. Ignacio Marquez Damian, and Christian Helman, "`Studies
|
||||
|
|
@ -162,6 +184,10 @@ Miscellaneous
|
|||
<https://doi.org/10.13140/RG.2.2.26088.01281>`_," 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
|
||||
<https://doi.org/10.1080/00295639.2020.1743578>`_", *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
|
||||
<https://doi.org/10.1016/j.anucene.2020.107514>`_", *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 <https://doi.org/10.1080/00295450.2019.1571828>`_,"
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@
|
|||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
#include <gsl/gsl>
|
||||
#include <hdf5.h>
|
||||
|
||||
#include "openmc/constants.h"
|
||||
|
|
@ -34,8 +35,9 @@ public:
|
|||
std::vector<double> energy;
|
||||
};
|
||||
|
||||
// Constructors
|
||||
Nuclide(hid_t group, const std::vector<double>& temperature, int i_nuclide);
|
||||
// Constructors/destructors
|
||||
Nuclide(hid_t group, const std::vector<double>& 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<double> kTs_; //!< temperatures in eV (k*T)
|
||||
|
|
|
|||
|
|
@ -4,9 +4,11 @@
|
|||
#include "openmc/endf.h"
|
||||
#include "openmc/particle.h"
|
||||
|
||||
#include <gsl/gsl>
|
||||
#include <hdf5.h>
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
#include <memory> // for unique_ptr
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
#include <utility> // 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<double, 1> energy_;
|
||||
|
|
@ -117,7 +120,7 @@ namespace data {
|
|||
extern xt::xtensor<double, 1> compton_profile_pz; //! Compton profile momentum grid
|
||||
|
||||
//! Photon interaction data for each element
|
||||
extern std::vector<PhotonInteraction> elements;
|
||||
extern std::vector<std::unique_ptr<PhotonInteraction>> elements;
|
||||
extern std::unordered_map<std::string, int> element_map;
|
||||
|
||||
} // namespace data
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
||||
|
|
|
|||
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
|
|||
|
|
@ -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<std::vector<double>>& 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<Nuclide>(
|
||||
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<int>(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<int>(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<std::vector<double>>& nuc_temps,
|
|||
mat->finalize();
|
||||
} // materials
|
||||
|
||||
|
||||
// Set up logarithmic grid for nuclides
|
||||
for (auto& nuc : data::nuclides) {
|
||||
nuc->init_grid();
|
||||
}
|
||||
int neutron = static_cast<int>(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<int>(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<int>(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);
|
||||
|
|
|
|||
|
|
@ -475,7 +475,7 @@ void Material::collision_stopping_power(double* s_col, bool positron)
|
|||
std::vector<double> 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<std::string>& name,
|
|||
if (n != nuclide_.size()) {
|
||||
nuclide_.resize(n);
|
||||
atom_density_ = xt::zeros<double>({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<std::string>& 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
|
||||
|
|
|
|||
113
src/nuclide.cpp
113
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<double>& temperature, int i_nuclide)
|
||||
: i_nuclide_{i_nuclide}
|
||||
Nuclide::Nuclide(hid_t group, const std::vector<double>& 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<double>& 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<int64_t>(i_nuclide_ + 1), *p.current_seed());
|
||||
//replaces, the seed is set with index_ + 1 instead of index_
|
||||
double r = future_prn(static_cast<int64_t>(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<double> temperature;
|
||||
int i_nuclide = data::nuclides.size();
|
||||
data::nuclides.push_back(std::make_unique<Nuclide>(
|
||||
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<double> temperature{temps, temps + n};
|
||||
data::nuclides.push_back(std::make_unique<Nuclide>(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<PhotonInteraction>(group));
|
||||
|
||||
close_group(group);
|
||||
file_close(file_id);
|
||||
}
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -27,7 +27,7 @@ namespace data {
|
|||
|
||||
xt::xtensor<double, 1> compton_profile_pz;
|
||||
|
||||
std::vector<PhotonInteraction> elements;
|
||||
std::vector<std::unique_ptr<PhotonInteraction>> elements;
|
||||
std::unordered_map<std::string, int> element_map;
|
||||
|
||||
} // namespace data
|
||||
|
|
@ -36,11 +36,14 @@ std::unordered_map<std::string, int> 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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -34,6 +34,7 @@
|
|||
#endif
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <string>
|
||||
|
||||
|
||||
|
|
@ -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<int>(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<int>(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<int>(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<int>(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<int>(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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue