crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 06:25:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5
OpenMC/src/material.cpp
Gavin Ridley aa4de82258
remove gsl-lite dependency (#3225) 
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
2025-02-20 01:03:20 +00:00

1596 lines
49 KiB
C++
Raw Blame History

#include "openmc/material.h"
#include <algorithm> // for min, max, sort, fill
#include <cassert>
#include <cmath>
#include <iterator>
#include <sstream>
#include <string>
#include <unordered_set>
#include "xtensor/xbuilder.hpp"
#include "xtensor/xoperation.hpp"
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/container_util.h"
#include "openmc/cross_sections.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/thermal.h"
#include "openmc/xml_interface.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace model {
std::unordered_map<int32_t, int32_t> material_map;
vector<unique_ptr<Material>> materials;
} // namespace model
//==============================================================================
// Material implementation
//==============================================================================
Material::Material(pugi::xml_node node)
{
index_ = model::materials.size(); // Avoids warning about narrowing
if (check_for_node(node, "id")) {
this->set_id(std::stoi(get_node_value(node, "id")));
} else {
fatal_error("Must specify id of material in materials XML file.");
}
if (check_for_node(node, "name")) {
name_ = get_node_value(node, "name");
}
if (check_for_node(node, "cfg")) {
auto cfg = get_node_value(node, "cfg");
write_message(
5, "NCrystal config string for material #{}: '{}'", this->id(), cfg);
ncrystal_mat_ = NCrystalMat(cfg);
}
if (check_for_node(node, "depletable")) {
depletable_ = get_node_value_bool(node, "depletable");
}
bool sum_density {false};
pugi::xml_node density_node = node.child("density");
std::string units;
if (density_node) {
units = get_node_value(density_node, "units");
if (units == "sum") {
sum_density = true;
} else if (units == "macro") {
if (check_for_node(density_node, "value")) {
density_ = std::stod(get_node_value(density_node, "value"));
} else {
density_ = 1.0;
}
} else {
double val = std::stod(get_node_value(density_node, "value"));
if (val <= 0.0) {
fatal_error("Need to specify a positive density on material " +
std::to_string(id_) + ".");
}
if (units == "g/cc" || units == "g/cm3") {
density_ = -val;
} else if (units == "kg/m3") {
density_ = -1.0e-3 * val;
} else if (units == "atom/b-cm") {
density_ = val;
} else if (units == "atom/cc" || units == "atom/cm3") {
density_ = 1.0e-24 * val;
} else {
fatal_error("Unknown units '" + units + "' specified on material " +
std::to_string(id_) + ".");
}
}
} else {
fatal_error("Must specify <density> element in material " +
std::to_string(id_) + ".");
}
if (node.child("element")) {
fatal_error(
"Unable to add an element to material " + std::to_string(id_) +
" since the element option has been removed from the xml input. "
"Elements can only be added via the Python API, which will expand "
"elements into their natural nuclides.");
}
// =======================================================================
// READ AND PARSE <nuclide> TAGS
// Check to ensure material has at least one nuclide
if (!check_for_node(node, "nuclide") &&
!check_for_node(node, "macroscopic")) {
fatal_error("No macroscopic data or nuclides specified on material " +
std::to_string(id_));
}
// Create list of macroscopic x/s based on those specified, just treat
// them as nuclides. This is all really a facade so the user thinks they
// are entering in macroscopic data but the code treats them the same
// as nuclides internally.
// Get pointer list of XML <macroscopic>
auto node_macros = node.children("macroscopic");
int num_macros = std::distance(node_macros.begin(), node_macros.end());
vector<std::string> names;
vector<double> densities;
if (settings::run_CE && num_macros > 0) {
fatal_error("Macroscopic can not be used in continuous-energy mode.");
} else if (num_macros > 1) {
fatal_error("Only one macroscopic object permitted per material, " +
std::to_string(id_));
} else if (num_macros == 1) {
pugi::xml_node node_nuc = *node_macros.begin();
// Check for empty name on nuclide
if (!check_for_node(node_nuc, "name")) {
fatal_error("No name specified on macroscopic data in material " +
std::to_string(id_));
}
// store nuclide name
std::string name = get_node_value(node_nuc, "name", false, true);
names.push_back(name);
// Set density for macroscopic data
if (units == "macro") {
densities.push_back(density_);
} else {
fatal_error("Units can only be macro for macroscopic data " + name);
}
} else {
// Create list of nuclides based on those specified
for (auto node_nuc : node.children("nuclide")) {
// Check for empty name on nuclide
if (!check_for_node(node_nuc, "name")) {
fatal_error(
"No name specified on nuclide in material " + std::to_string(id_));
}
// store nuclide name
std::string name = get_node_value(node_nuc, "name", false, true);
names.push_back(name);
// Check if no atom/weight percents were specified or if both atom and
// weight percents were specified
if (units == "macro") {
densities.push_back(density_);
} else {
bool has_ao = check_for_node(node_nuc, "ao");
bool has_wo = check_for_node(node_nuc, "wo");
if (!has_ao && !has_wo) {
fatal_error(
"No atom or weight percent specified for nuclide: " + name);
} else if (has_ao && has_wo) {
fatal_error("Cannot specify both atom and weight percents for a "
"nuclide: " +
name);
}
// Copy atom/weight percents
if (has_ao) {
densities.push_back(std::stod(get_node_value(node_nuc, "ao")));
} else {
densities.push_back(-std::stod(get_node_value(node_nuc, "wo")));
}
}
}
}
// =======================================================================
// READ AND PARSE <isotropic> element
vector<std::string> iso_lab;
if (check_for_node(node, "isotropic")) {
iso_lab = get_node_array<std::string>(node, "isotropic");
}
// ========================================================================
// COPY NUCLIDES TO ARRAYS IN MATERIAL
// allocate arrays in Material object
auto n = names.size();
nuclide_.reserve(n);
atom_density_ = xt::empty<double>({n});
if (settings::photon_transport)
element_.reserve(n);
for (int i = 0; i < n; ++i) {
const auto& name {names[i]};
// Check that this nuclide is listed in the nuclear data library
// (cross_sections.xml for CE and the MGXS HDF5 for MG)
if (settings::run_mode != RunMode::PLOTTING) {
LibraryKey key {Library::Type::neutron, name};
if (data::library_map.find(key) == data::library_map.end()) {
fatal_error("Could not find nuclide " + name +
" in the "
"nuclear data library.");
}
}
// If this nuclide hasn't been encountered yet, we need to add its name
// and alias to the nuclide_dict
if (data::nuclide_map.find(name) == data::nuclide_map.end()) {
int index = data::nuclide_map.size();
data::nuclide_map[name] = index;
nuclide_.push_back(index);
} else {
nuclide_.push_back(data::nuclide_map[name]);
}
// 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) {
std::string element = to_element(name);
// Make sure photon cross section data is available
if (settings::run_mode != RunMode::PLOTTING) {
LibraryKey key {Library::Type::photon, element};
if (data::library_map.find(key) == data::library_map.end()) {
fatal_error(
"Could not find element " + element + " in cross_sections.xml.");
}
}
if (data::element_map.find(element) == data::element_map.end()) {
int index = data::element_map.size();
data::element_map[element] = index;
element_.push_back(index);
} else {
element_.push_back(data::element_map[element]);
}
}
// Copy atom/weight percent
atom_density_(i) = densities[i];
}
if (settings::run_CE) {
// By default, isotropic-in-lab is not used
if (iso_lab.size() > 0) {
p0_.resize(n);
// Apply isotropic-in-lab treatment to specified nuclides
for (int j = 0; j < n; ++j) {
for (const auto& nuc : iso_lab) {
if (names[j] == nuc) {
p0_[j] = true;
break;
}
}
}
}
}
// Check to make sure either all atom percents or all weight percents are
// given
if (!(xt::all(atom_density_ >= 0.0) || xt::all(atom_density_ <= 0.0))) {
fatal_error(
"Cannot mix atom and weight percents in material " + std::to_string(id_));
}
// Determine density if it is a sum value
if (sum_density)
density_ = xt::sum(atom_density_)();
if (check_for_node(node, "temperature")) {
temperature_ = std::stod(get_node_value(node, "temperature"));
}
if (check_for_node(node, "volume")) {
volume_ = std::stod(get_node_value(node, "volume"));
}
// =======================================================================
// READ AND PARSE <sab> TAG FOR THERMAL SCATTERING DATA
if (settings::run_CE) {
// Loop over <sab> elements
vector<std::string> sab_names;
for (auto node_sab : node.children("sab")) {
// Determine name of thermal scattering table
if (!check_for_node(node_sab, "name")) {
fatal_error("Need to specify <name> for thermal scattering table.");
}
std::string name = get_node_value(node_sab, "name");
sab_names.push_back(name);
// Read the fraction of nuclei affected by this thermal scattering table
double fraction = 1.0;
if (check_for_node(node_sab, "fraction")) {
fraction = std::stod(get_node_value(node_sab, "fraction"));
}
// Check that the thermal scattering table is listed in the
// cross_sections.xml file
if (settings::run_mode != RunMode::PLOTTING) {
LibraryKey key {Library::Type::thermal, name};
if (data::library_map.find(key) == data::library_map.end()) {
fatal_error("Could not find thermal scattering data " + name +
" in cross_sections.xml file.");
}
}
// Determine index of thermal scattering data in global
// data::thermal_scatt array
int index_table;
if (data::thermal_scatt_map.find(name) == data::thermal_scatt_map.end()) {
index_table = data::thermal_scatt_map.size();
data::thermal_scatt_map[name] = index_table;
} else {
index_table = data::thermal_scatt_map[name];
}
// Add entry to thermal tables vector. For now, we put the nuclide index
// as zero since we don't know which nuclides the table is being applied
// to yet (this is assigned in init_thermal)
thermal_tables_.push_back({index_table, 0, fraction});
}
}
}
Material::~Material()
{
model::material_map.erase(id_);
}
Material& Material::clone()
{
std::unique_ptr<Material> mat = std::make_unique<Material>();
// set all other parameters to whatever the calling Material has
mat->name_ = name_;
mat->nuclide_ = nuclide_;
mat->element_ = element_;
mat->ncrystal_mat_ = ncrystal_mat_;
mat->atom_density_ = atom_density_;
mat->density_ = density_;
mat->density_gpcc_ = density_gpcc_;
mat->volume_ = volume_;
mat->fissionable() = fissionable_;
mat->depletable() = depletable_;
mat->p0_ = p0_;
mat->mat_nuclide_index_ = mat_nuclide_index_;
mat->thermal_tables_ = thermal_tables_;
mat->temperature_ = temperature_;
if (ttb_)
mat->ttb_ = std::make_unique<Bremsstrahlung>(*ttb_);
mat->index_ = model::materials.size();
mat->set_id(C_NONE);
model::materials.push_back(std::move(mat));
return *model::materials.back();
}
void Material::finalize()
{
// Set fissionable if any nuclide is fissionable
if (settings::run_CE) {
for (const auto& i_nuc : nuclide_) {
if (data::nuclides[i_nuc]->fissionable_) {
fissionable_ = true;
break;
}
}
// Generate material bremsstrahlung data for electrons and positrons
if (settings::photon_transport &&
settings::electron_treatment == ElectronTreatment::TTB) {
this->init_bremsstrahlung();
}
// Assign thermal scattering tables
this->init_thermal();
}
// Normalize density
this->normalize_density();
}
void Material::normalize_density()
{
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
int i_nuc = nuclide_[i];
double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_
: data::mg.nuclides_[i_nuc].awr;
// if given weight percent, convert all values so that they are divided
// by awr. thus, when a sum is done over the values, it's actually
// sum(w/awr)
if (!percent_in_atom)
atom_density_(i) = -atom_density_(i) / awr;
}
// determine normalized atom percents. if given atom percents, this is
// straightforward. if given weight percents, the value is w/awr and is
// divided by sum(w/awr)
atom_density_ /= xt::sum(atom_density_)();
// Change density in g/cm^3 to atom/b-cm. Since all values are now in
// atom percent, the sum needs to be re-evaluated as 1/sum(x*awr)
if (!density_in_atom) {
double sum_percent = 0.0;
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_
: data::mg.nuclides_[i_nuc].awr;
sum_percent += atom_density_(i) * awr;
}
sum_percent = 1.0 / sum_percent;
density_ = -density_ * N_AVOGADRO / MASS_NEUTRON * sum_percent;
}
// Calculate nuclide atom densities
atom_density_ *= density_;
// Calculate density in g/cm^3.
density_gpcc_ = 0.0;
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0;
density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
}
}
void Material::init_thermal()
{
vector<ThermalTable> tables;
std::unordered_set<int> already_checked;
for (const auto& table : thermal_tables_) {
// Make sure each S(a,b) table only gets checked once
if (already_checked.find(table.index_table) != already_checked.end()) {
continue;
}
already_checked.insert(table.index_table);
// In order to know which nuclide the S(a,b) table applies to, we need
// to search through the list of nuclides for one which has a matching
// name
bool found = false;
for (int j = 0; j < nuclide_.size(); ++j) {
const auto& name {data::nuclides[nuclide_[j]]->name_};
if (contains(data::thermal_scatt[table.index_table]->nuclides_, name)) {
tables.push_back({table.index_table, j, table.fraction});
found = true;
}
}
// Check to make sure thermal scattering table matched a nuclide
if (!found) {
fatal_error("Thermal scattering table " +
data::thermal_scatt[table.index_table]->name_ +
" did not match any nuclide on material " +
std::to_string(id_));
}
}
// Make sure each nuclide only appears in one table.
for (int j = 0; j < tables.size(); ++j) {
for (int k = j + 1; k < tables.size(); ++k) {
if (tables[j].index_nuclide == tables[k].index_nuclide) {
int index = nuclide_[tables[j].index_nuclide];
auto name = data::nuclides[index]->name_;
fatal_error(
name + " in material " + std::to_string(id_) +
" was found "
"in multiple thermal scattering tables. Each nuclide can appear in "
"only one table per material.");
}
}
}
// If there are multiple S(a,b) tables, we need to make sure that the
// entries in i_sab_nuclides are sorted or else they won't be applied
// correctly in the cross_section module.
std::sort(tables.begin(), tables.end(), [](ThermalTable a, ThermalTable b) {
return a.index_nuclide < b.index_nuclide;
});
// Update the list of thermal tables
thermal_tables_ = tables;
}
void Material::collision_stopping_power(double* s_col, bool positron)
{
// Average electron number and average atomic weight
double electron_density = 0.0;
double mass_density = 0.0;
// Log of the mean excitation energy of the material
double log_I = 0.0;
// Effective number of conduction electrons in the material
double n_conduction = 0.0;
// Oscillator strength and square of the binding energy for each oscillator
// in material
vector<double> f;
vector<double> e_b_sq;
for (int i = 0; i < element_.size(); ++i) {
const auto& elm = *data::elements[element_[i]];
double awr = data::nuclides[nuclide_[i]]->awr_;
// Get atomic density of nuclide given atom/weight percent
double atom_density =
(atom_density_[0] > 0.0) ? atom_density_[i] : -atom_density_[i] / awr;
electron_density += atom_density * elm.Z_;
mass_density += atom_density * awr * MASS_NEUTRON;
log_I += atom_density * elm.Z_ * std::log(elm.I_);
for (int j = 0; j < elm.n_electrons_.size(); ++j) {
if (elm.n_electrons_[j] < 0) {
n_conduction -= elm.n_electrons_[j] * atom_density;
continue;
}
e_b_sq.push_back(elm.ionization_energy_[j] * elm.ionization_energy_[j]);
f.push_back(elm.n_electrons_[j] * atom_density);
}
}
log_I /= electron_density;
n_conduction /= electron_density;
for (auto& f_i : f)
f_i /= electron_density;
// Get density in g/cm^3 if it is given in atom/b-cm
double density = (density_ < 0.0) ? -density_ : mass_density / N_AVOGADRO;
// Calculate the square of the plasma energy
double e_p_sq =
PLANCK_C * PLANCK_C * PLANCK_C * N_AVOGADRO * electron_density * density /
(2.0 * PI * PI * FINE_STRUCTURE * MASS_ELECTRON_EV * mass_density);
// Get the Sternheimer adjustment factor
double rho =
sternheimer_adjustment(f, e_b_sq, e_p_sq, n_conduction, log_I, 1.0e-6, 100);
// Classical electron radius in cm
constexpr double CM_PER_ANGSTROM {1.0e-8};
constexpr double r_e =
CM_PER_ANGSTROM * PLANCK_C / (2.0 * PI * FINE_STRUCTURE * MASS_ELECTRON_EV);
// Constant in expression for collision stopping power
constexpr double BARN_PER_CM_SQ {1.0e24};
double c =
BARN_PER_CM_SQ * 2.0 * PI * r_e * r_e * MASS_ELECTRON_EV * electron_density;
// Loop over incident charged particle energies
for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
double E = data::ttb_e_grid(i);
// Get the density effect correction
double delta =
density_effect(f, e_b_sq, e_p_sq, n_conduction, rho, E, 1.0e-6, 100);
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) /
((E + MASS_ELECTRON_EV) * (E + MASS_ELECTRON_EV));
double tau = E / MASS_ELECTRON_EV;
double F;
if (positron) {
double t = tau + 2.0;
F = std::log(4.0) - (beta_sq / 12.0) * (23.0 + 14.0 / t + 10.0 / (t * t) +
4.0 / (t * t * t));
} else {
F = (1.0 - beta_sq) *
(1.0 + tau * tau / 8.0 - (2.0 * tau + 1.0) * std::log(2.0));
}
// Calculate the collision stopping power for this energy
s_col[i] =
c / beta_sq *
(2.0 * (std::log(E) - log_I) + std::log(1.0 + tau / 2.0) + F - delta);
}
}
void Material::init_bremsstrahlung()
{
// Create new object
ttb_ = make_unique<Bremsstrahlung>();
// Get the size of the energy grids
auto n_k = data::ttb_k_grid.size();
auto n_e = data::ttb_e_grid.size();
// Determine number of elements
int n = element_.size();
for (int particle = 0; particle < 2; ++particle) {
// Loop over logic twice, once for electron, once for positron
BremsstrahlungData* ttb =
(particle == 0) ? &ttb_->electron : &ttb_->positron;
bool positron = (particle == 1);
// Allocate arrays for TTB data
ttb->pdf = xt::zeros<double>({n_e, n_e});
ttb->cdf = xt::zeros<double>({n_e, n_e});
ttb->yield = xt::empty<double>({n_e});
// Allocate temporary arrays
xt::xtensor<double, 1> stopping_power_collision({n_e}, 0.0);
xt::xtensor<double, 1> stopping_power_radiative({n_e}, 0.0);
xt::xtensor<double, 2> dcs({n_e, n_k}, 0.0);
double Z_eq_sq = 0.0;
double sum_density = 0.0;
// Get the collision stopping power of the material
this->collision_stopping_power(stopping_power_collision.data(), positron);
// Calculate the molecular DCS and the molecular radiative stopping power
// using Bragg's additivity rule.
for (int i = 0; i < n; ++i) {
// Get pointer to current element
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
double atom_density =
(atom_density_[0] > 0.0) ? atom_density_[i] : -atom_density_[i] / awr;
// Calculate the "equivalent" atomic number Zeq of the material
Z_eq_sq += atom_density * elm.Z_ * elm.Z_;
sum_density += atom_density;
// Accumulate material DCS
dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_;
// Accumulate material radiative stopping power
stopping_power_radiative += atom_density * elm.stopping_power_radiative_;
}
Z_eq_sq /= sum_density;
// Calculate the positron DCS and radiative stopping power. These are
// obtained by multiplying the electron DCS and radiative stopping powers by
// a factor r, which is a numerical approximation of the ratio of the
// radiative stopping powers for positrons and electrons. Source: F. Salvat,
// J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
// Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
// Issy-les-Moulineaux, France (2011).
if (positron) {
for (int i = 0; i < n_e; ++i) {
double t = std::log(
1.0 + 1.0e6 * data::ttb_e_grid(i) / (Z_eq_sq * MASS_ELECTRON_EV));
double r =
1.0 -
std::exp(-1.2359e-1 * t + 6.1274e-2 * std::pow(t, 2) -
3.1516e-2 * std::pow(t, 3) + 7.7446e-3 * std::pow(t, 4) -
1.0595e-3 * std::pow(t, 5) + 7.0568e-5 * std::pow(t, 6) -
1.808e-6 * std::pow(t, 7));
stopping_power_radiative(i) *= r;
auto dcs_i = xt::view(dcs, i, xt::all());
dcs_i *= r;
}
}
// Total material stopping power
xt::xtensor<double, 1> stopping_power =
stopping_power_collision + stopping_power_radiative;
// Loop over photon energies
xt::xtensor<double, 1> f({n_e}, 0.0);
xt::xtensor<double, 1> z({n_e}, 0.0);
for (int i = 0; i < n_e - 1; ++i) {
double w = data::ttb_e_grid(i);
// Loop over incident particle energies
for (int j = i; j < n_e; ++j) {
double e = data::ttb_e_grid(j);
// Reduced photon energy
double k = w / e;
// Find the lower bounding index of the reduced photon energy
int i_k = lower_bound_index(
data::ttb_k_grid.cbegin(), data::ttb_k_grid.cend(), k);
// Get the interpolation bounds
double k_l = data::ttb_k_grid(i_k);
double k_r = data::ttb_k_grid(i_k + 1);
double x_l = dcs(j, i_k);
double x_r = dcs(j, i_k + 1);
// Find the value of the DCS using linear interpolation in reduced
// photon energy k
double x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l);
// Square of the ratio of the speed of light to the velocity of the
// charged particle
double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) /
((e + MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
// Compute the integrand of the PDF
f(j) = x / (beta_sq * stopping_power(j) * w);
}
// Number of points to integrate
int n = n_e - i;
// Integrate the PDF using cubic spline integration over the incident
// particle energy
if (n > 2) {
spline(n, &data::ttb_e_grid(i), &f(i), &z(i));
double c = 0.0;
for (int j = i; j < n_e - 1; ++j) {
c += spline_integrate(n, &data::ttb_e_grid(i), &f(i), &z(i),
data::ttb_e_grid(j), data::ttb_e_grid(j + 1));
ttb->pdf(j + 1, i) = c;
}
// Integrate the last two points using trapezoidal rule in log-log space
} else {
double e_l = std::log(data::ttb_e_grid(i));
double e_r = std::log(data::ttb_e_grid(i + 1));
double x_l = std::log(f(i));
double x_r = std::log(f(i + 1));
ttb->pdf(i + 1, i) =
0.5 * (e_r - e_l) * (std::exp(e_l + x_l) + std::exp(e_r + x_r));
}
}
// Loop over incident particle energies
for (int j = 1; j < n_e; ++j) {
// Set last element of PDF to small non-zero value to enable log-log
// interpolation
ttb->pdf(j, j) = std::exp(-500.0);
// Loop over photon energies
double c = 0.0;
for (int i = 0; i < j; ++i) {
// Integrate the CDF from the PDF using the trapezoidal rule in log-log
// space
double w_l = std::log(data::ttb_e_grid(i));
double w_r = std::log(data::ttb_e_grid(i + 1));
double x_l = std::log(ttb->pdf(j, i));
double x_r = std::log(ttb->pdf(j, i + 1));
c += 0.5 * (w_r - w_l) * (std::exp(w_l + x_l) + std::exp(w_r + x_r));
ttb->cdf(j, i + 1) = c;
}
// Set photon number yield
ttb->yield(j) = c;
}
// Use logarithm of number yield since it is log-log interpolated
ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0);
}
}
void Material::init_nuclide_index()
{
int n = settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
mat_nuclide_index_.resize(n);
std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
for (int i = 0; i < nuclide_.size(); ++i) {
mat_nuclide_index_[nuclide_[i]] = i;
}
}
void Material::calculate_xs(Particle& p) const
{
// Set all material macroscopic cross sections to zero
p.macro_xs().total = 0.0;
p.macro_xs().absorption = 0.0;
p.macro_xs().fission = 0.0;
p.macro_xs().nu_fission = 0.0;
if (p.type() == ParticleType::neutron) {
this->calculate_neutron_xs(p);
} else if (p.type() == ParticleType::photon) {
this->calculate_photon_xs(p);
}
}
void Material::calculate_neutron_xs(Particle& p) const
{
// Find energy index on energy grid
int neutron = static_cast<int>(ParticleType::neutron);
int i_grid =
std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;
// Determine if this material has S(a,b) tables
bool check_sab = (thermal_tables_.size() > 0);
// Initialize position in i_sab_nuclides
int j = 0;
// Calculate NCrystal cross section
double ncrystal_xs = -1.0;
if (ncrystal_mat_ && p.E() < NCRYSTAL_MAX_ENERGY) {
ncrystal_xs = ncrystal_mat_.xs(p);
}
// Add contribution from each nuclide in material
for (int i = 0; i < nuclide_.size(); ++i) {
// ======================================================================
// CHECK FOR S(A,B) TABLE
int i_sab = C_NONE;
double sab_frac = 0.0;
// Check if this nuclide matches one of the S(a,b) tables specified.
// This relies on thermal_tables_ being sorted by .index_nuclide
if (check_sab) {
const auto& sab {thermal_tables_[j]};
if (i == sab.index_nuclide) {
// Get index in sab_tables
i_sab = sab.index_table;
sab_frac = sab.fraction;
// If particle energy is greater than the highest energy for the
// S(a,b) table, then don't use the S(a,b) table
if (p.E() > data::thermal_scatt[i_sab]->energy_max_)
i_sab = C_NONE;
// Increment position in thermal_tables_
++j;
// Don't check for S(a,b) tables if there are no more left
if (j == thermal_tables_.size())
check_sab = false;
}
}
// ======================================================================
// CALCULATE MICROSCOPIC CROSS SECTION
// Get nuclide index
int i_nuclide = nuclide_[i];
// Update microscopic cross section for this nuclide
p.update_neutron_xs(i_nuclide, i_grid, i_sab, sab_frac, ncrystal_xs);
auto& micro = p.neutron_xs(i_nuclide);
// ======================================================================
// ADD TO MACROSCOPIC CROSS SECTION
// Copy atom density of nuclide in material
double atom_density = atom_density_(i);
// Add contributions to cross sections
p.macro_xs().total += atom_density * micro.total;
p.macro_xs().absorption += atom_density * micro.absorption;
p.macro_xs().fission += atom_density * micro.fission;
p.macro_xs().nu_fission += atom_density * micro.nu_fission;
}
}
void Material::calculate_photon_xs(Particle& p) const
{
p.macro_xs().coherent = 0.0;
p.macro_xs().incoherent = 0.0;
p.macro_xs().photoelectric = 0.0;
p.macro_xs().pair_production = 0.0;
// Add contribution from each nuclide in material
for (int i = 0; i < nuclide_.size(); ++i) {
// ========================================================================
// CALCULATE MICROSCOPIC CROSS SECTION
// Determine microscopic cross sections for this nuclide
int i_element = element_[i];
// 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);
}
// ========================================================================
// ADD TO MACROSCOPIC CROSS SECTION
// Copy atom density of nuclide in material
double atom_density = atom_density_(i);
// Add contributions to material macroscopic cross sections
p.macro_xs().total += atom_density * micro.total;
p.macro_xs().coherent += atom_density * micro.coherent;
p.macro_xs().incoherent += atom_density * micro.incoherent;
p.macro_xs().photoelectric += atom_density * micro.photoelectric;
p.macro_xs().pair_production += atom_density * micro.pair_production;
}
}
void Material::set_id(int32_t id)
{
assert(id >= 0 || id == C_NONE);
// Clear entry in material map if an ID was already assigned before
if (id_ != C_NONE) {
model::material_map.erase(id_);
id_ = C_NONE;
}
// Make sure no other material has same ID
if (model::material_map.find(id) != model::material_map.end()) {
throw std::runtime_error {
"Two materials have the same ID: " + std::to_string(id)};
}
// If no ID specified, auto-assign next ID in sequence
if (id == C_NONE) {
id = 0;
for (const auto& m : model::materials) {
id = std::max(id, m->id_);
}
++id;
}
// Update ID and entry in material map
id_ = id;
model::material_map[id] = index_;
}
void Material::set_density(double density, const std::string& units)
{
assert(density >= 0.0);
if (nuclide_.empty()) {
throw std::runtime_error {"No nuclides exist in material yet."};
}
if (units == "atom/b-cm") {
// Set total density based on value provided
density_ = density;
// Determine normalized atom percents
double sum_percent = xt::sum(atom_density_)();
atom_density_ /= sum_percent;
// Recalculate nuclide atom densities based on given density
atom_density_ *= density;
// Calculate density in g/cm^3.
density_gpcc_ = 0.0;
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
double awr = data::nuclides[i_nuc]->awr_;
density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
}
} else if (units == "g/cm3" || units == "g/cc") {
// Determine factor by which to change densities
double previous_density_gpcc = density_gpcc_;
double f = density / previous_density_gpcc;
// Update densities
density_gpcc_ = density;
density_ *= f;
atom_density_ *= f;
} else {
throw std::invalid_argument {
"Invalid units '" + std::string(units.data()) + "' specified."};
}
}
void Material::set_densities(
const vector<std::string>& name, const vector<double>& density)
{
auto n = name.size();
assert(n > 0);
assert(n == density.size());
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 (int64_t 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(), nullptr, 0);
if (err < 0)
throw std::runtime_error {openmc_err_msg};
}
nuclide_[i] = data::nuclide_map.at(nuc);
assert(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();
}
double Material::volume() const
{
if (volume_ < 0.0) {
throw std::runtime_error {
"Volume for material with ID=" + std::to_string(id_) + " not set."};
}
return volume_;
}
double Material::temperature() const
{
// If material doesn't have an assigned temperature, use global default
return temperature_ >= 0 ? temperature_ : settings::temperature_default;
}
void Material::to_hdf5(hid_t group) const
{
hid_t material_group = create_group(group, "material " + std::to_string(id_));
write_attribute(material_group, "depletable", static_cast<int>(depletable()));
if (volume_ > 0.0) {
write_attribute(material_group, "volume", volume_);
}
if (temperature_ > 0.0) {
write_attribute(material_group, "temperature", temperature_);
}
write_dataset(material_group, "name", name_);
write_dataset(material_group, "atom_density", density_);
// Copy nuclide/macro name for each nuclide to vector
vector<std::string> nuc_names;
vector<std::string> macro_names;
vector<double> nuc_densities;
if (settings::run_CE) {
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
nuc_names.push_back(data::nuclides[i_nuc]->name_);
nuc_densities.push_back(atom_density_(i));
}
} else {
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
if (data::mg.nuclides_[i_nuc].awr != MACROSCOPIC_AWR) {
nuc_names.push_back(data::mg.nuclides_[i_nuc].name);
nuc_densities.push_back(atom_density_(i));
} else {
macro_names.push_back(data::mg.nuclides_[i_nuc].name);
}
}
}
// Write vector to 'nuclides'
if (!nuc_names.empty()) {
write_dataset(material_group, "nuclides", nuc_names);
write_dataset(material_group, "nuclide_densities", nuc_densities);
}
// Write vector to 'macroscopics'
if (!macro_names.empty()) {
write_dataset(material_group, "macroscopics", macro_names);
}
if (!thermal_tables_.empty()) {
vector<std::string> sab_names;
for (const auto& table : thermal_tables_) {
sab_names.push_back(data::thermal_scatt[table.index_table]->name_);
}
write_dataset(material_group, "sab_names", sab_names);
}
close_group(material_group);
}
void Material::export_properties_hdf5(hid_t group) const
{
hid_t material_group = create_group(group, "material " + std::to_string(id_));
write_attribute(material_group, "atom_density", density_);
write_attribute(material_group, "mass_density", density_gpcc_);
close_group(material_group);
}
void Material::import_properties_hdf5(hid_t group)
{
hid_t material_group = open_group(group, "material " + std::to_string(id_));
double density;
read_attribute(material_group, "atom_density", density);
this->set_density(density, "atom/b-cm");
close_group(material_group);
}
void Material::add_nuclide(const std::string& name, double density)
{
// Check if nuclide is already in material
for (int i = 0; i < nuclide_.size(); ++i) {
int i_nuc = nuclide_[i];
if (data::nuclides[i_nuc]->name_ == name) {
double awr = data::nuclides[i_nuc]->awr_;
density_ += density - atom_density_(i);
density_gpcc_ +=
(density - atom_density_(i)) * awr * MASS_NEUTRON / N_AVOGADRO;
atom_density_(i) = density;
return;
}
}
// If nuclide wasn't found, extend nuclide/density arrays
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
xt::xtensor<double, 1> atom_density = xt::zeros<double>({n});
xt::view(atom_density, xt::range(0, n - 1)) = atom_density_;
atom_density(n - 1) = density;
atom_density_ = atom_density;
density_ += density;
density_gpcc_ +=
density * data::nuclides[i_nuc]->awr_ * MASS_NEUTRON / N_AVOGADRO;
}
//==============================================================================
// Non-method functions
//==============================================================================
double sternheimer_adjustment(const vector<double>& f,
const vector<double>& e_b_sq, double e_p_sq, double n_conduction,
double log_I, double tol, int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Calculate the Sternheimer adjustment factor using Newton's method
double rho = 2.0;
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double rho_0 = rho;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
// Square of resonance energy of a bound-shell oscillator
double e_r_sq = e_b_sq[i] * rho * rho + 2.0 / 3.0 * f[i] * e_p_sq;
g += f[i] * std::log(e_r_sq);
gp += e_b_sq[i] * f[i] * rho / e_r_sq;
}
// Include conduction electrons
if (n_conduction > 0.0) {
g += n_conduction * std::log(n_conduction * e_p_sq);
}
// Set the next guess: rho_n+1 = rho_n - g(rho_n)/g'(rho_n)
rho -= (g - 2.0 * log_I) / (2.0 * gp);
// If the initial guess is too large, rho can be negative
if (rho < 0.0)
rho = rho_0 / 2.0;
// Check for convergence
if (std::abs(rho - rho_0) / rho_0 < tol)
break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded.");
rho = 1.0e-6;
}
return rho;
}
double density_effect(const vector<double>& f, const vector<double>& e_b_sq,
double e_p_sq, double n_conduction, double rho, double E, double tol,
int max_iter)
{
// Get the total number of oscillators
int n = f.size();
// Square of the ratio of the speed of light to the velocity of the charged
// particle
double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) /
((E + MASS_ELECTRON_EV) * (E + MASS_ELECTRON_EV));
// For nonmetals, delta = 0 for beta < beta_0, where beta_0 is obtained by
// setting the frequency w = 0.
double beta_0_sq = 0.0;
if (n_conduction == 0.0) {
for (int i = 0; i < n; ++i) {
beta_0_sq += f[i] * e_p_sq / (e_b_sq[i] * rho * rho);
}
beta_0_sq = 1.0 / (1.0 + beta_0_sq);
}
double delta = 0.0;
if (beta_sq < beta_0_sq)
return delta;
// Compute the square of the frequency w^2 using Newton's method, with the
// initial guess of w^2 equal to beta^2 * gamma^2
double w_sq = E / MASS_ELECTRON_EV * (E / MASS_ELECTRON_EV + 2);
int iter;
for (iter = 0; iter < max_iter; ++iter) {
double w_sq_0 = w_sq;
// Function to find the root of and its derivative
double g = 0.0;
double gp = 0.0;
for (int i = 0; i < n; ++i) {
double c = e_b_sq[i] * rho * rho / e_p_sq + w_sq;
g += f[i] / c;
gp -= f[i] / (c * c);
}
// Include conduction electrons
g += n_conduction / w_sq;
gp -= n_conduction / (w_sq * w_sq);
// Set the next guess: w_n+1 = w_n - g(w_n)/g'(w_n)
w_sq -= (g + 1.0 - 1.0 / beta_sq) / gp;
// If the initial guess is too large, w can be negative
if (w_sq < 0.0)
w_sq = w_sq_0 / 2.0;
// Check for convergence
if (std::abs(w_sq - w_sq_0) / w_sq_0 < tol)
break;
}
// Did not converge
if (iter >= max_iter) {
warning("Maximum Newton-Raphson iterations exceeded: setting density "
"effect correction to zero.");
return delta;
}
// Solve for the density effect correction
for (int i = 0; i < n; ++i) {
double l_sq = e_b_sq[i] * rho * rho / e_p_sq + 2.0 / 3.0 * f[i];
delta += f[i] * std::log((l_sq + w_sq) / l_sq);
}
// Include conduction electrons
if (n_conduction > 0.0) {
delta += n_conduction * std::log((n_conduction + w_sq) / n_conduction);
}
return delta - w_sq * (1.0 - beta_sq);
}
void read_materials_xml()
{
write_message("Reading materials XML file...", 5);
pugi::xml_document doc;
// Check if materials.xml exists
std::string filename = settings::path_input + "materials.xml";
if (!file_exists(filename)) {
fatal_error("Material XML file '" + filename + "' does not exist!");
}
// Parse materials.xml file and get root element
doc.load_file(filename.c_str());
// Loop over XML material elements and populate the array.
pugi::xml_node root = doc.document_element();
read_materials_xml(root);
}
void read_materials_xml(pugi::xml_node root)
{
for (pugi::xml_node material_node : root.children("material")) {
model::materials.push_back(make_unique<Material>(material_node));
}
model::materials.shrink_to_fit();
}
void free_memory_material()
{
model::materials.clear();
model::material_map.clear();
}
//==============================================================================
// C API
//==============================================================================
extern "C" int openmc_get_material_index(int32_t id, int32_t* index)
{
auto it = model::material_map.find(id);
if (it == model::material_map.end()) {
set_errmsg("No material exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
} else {
*index = it->second;
return 0;
}
}
extern "C" int openmc_material_add_nuclide(
int32_t index, const char* name, double density)
{
int err = 0;
if (index >= 0 && index < model::materials.size()) {
try {
model::materials[index]->add_nuclide(name, density);
} catch (const std::runtime_error& e) {
return OPENMC_E_DATA;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return err;
}
extern "C" int openmc_material_get_densities(
int32_t index, const int** nuclides, const double** densities, int* n)
{
if (index >= 0 && index < model::materials.size()) {
auto& mat = model::materials[index];
if (!mat->nuclides().empty()) {
*nuclides = mat->nuclides().data();
*densities = mat->densities().data();
*n = mat->nuclides().size();
return 0;
} else {
set_errmsg("Material atom density array has not been allocated.");
return OPENMC_E_ALLOCATE;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_get_density(int32_t index, double* density)
{
if (index >= 0 && index < model::materials.size()) {
auto& mat = model::materials[index];
*density = mat->density_gpcc();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_get_fissionable(int32_t index, bool* fissionable)
{
if (index >= 0 && index < model::materials.size()) {
*fissionable = model::materials[index]->fissionable();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_get_id(int32_t index, int32_t* id)
{
if (index >= 0 && index < model::materials.size()) {
*id = model::materials[index]->id();
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_get_temperature(
int32_t index, double* temperature)
{
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*temperature = model::materials[index]->temperature();
return 0;
}
extern "C" int openmc_material_get_volume(int32_t index, double* volume)
{
if (index >= 0 && index < model::materials.size()) {
try {
*volume = model::materials[index]->volume();
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
return 0;
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_set_density(
int32_t index, double density, const char* units)
{
if (index >= 0 && index < model::materials.size()) {
try {
model::materials[index]->set_density(density, units);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int openmc_material_set_densities(
int32_t index, int n, const char** name, const double* density)
{
if (index >= 0 && index < model::materials.size()) {
try {
model::materials[index]->set_densities(
{name, name + n}, {density, density + n});
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int openmc_material_set_id(int32_t index, int32_t id)
{
if (index >= 0 && index < model::materials.size()) {
try {
model::materials.at(index)->set_id(id);
} catch (const std::exception& e) {
set_errmsg(e.what());
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int openmc_material_get_name(int32_t index, const char** name)
{
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*name = model::materials[index]->name().data();
return 0;
}
extern "C" int openmc_material_set_name(int32_t index, const char* name)
{
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::materials[index]->set_name(name);
return 0;
}
extern "C" int openmc_material_set_volume(int32_t index, double volume)
{
if (index >= 0 && index < model::materials.size()) {
auto& m {model::materials[index]};
if (volume >= 0.0) {
m->volume_ = volume;
return 0;
} else {
set_errmsg("Volume must be non-negative");
return OPENMC_E_INVALID_ARGUMENT;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int openmc_material_get_depletable(int32_t index, bool* depletable)
{
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*depletable = model::materials[index]->depletable();
return 0;
}
extern "C" int openmc_material_set_depletable(int32_t index, bool depletable)
{
if (index < 0 || index >= model::materials.size()) {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::materials[index]->depletable() = depletable;
return 0;
}
extern "C" int openmc_extend_materials(
int32_t n, int32_t* index_start, int32_t* index_end)
{
if (index_start)
*index_start = model::materials.size();
if (index_end)
*index_end = model::materials.size() + n - 1;
for (int32_t i = 0; i < n; i++) {
model::materials.push_back(make_unique<Material>());
}
return 0;
}
extern "C" size_t n_materials()
{
return model::materials.size();
}
} // namespace openmc
Reference in a new issue View git blame Copy permalink