crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fixed issue in MG-Mode which surfaced when running problems which have microscopic cross sections defined and multiple fissile isotopes in a material

Browse source
This commit is contained in:
Adam Nelson 2020-10-15 13:46:05 -05:00
parent a70eda04ff
commit 77015c0b14
8 changed files with 120 additions and 109 deletions
Download patch file Download diff file Expand all files Collapse all files

4
include/openmc/constants.h
View file

@ -289,9 +289,9 @@ enum class MgxsType {
ABSORPTION,
INVERSE_VELOCITY,
DECAY_RATE,
NU_SCATTER,
SCATTER,
SCATTER_MULT,
SCATTER_FMU_MULT,
NU_SCATTER_FMU,
SCATTER_FMU,
FISSION,
KAPPA_FISSION,

2
include/openmc/mgxs.h
View file

@ -117,7 +117,7 @@ class Mgxs {
int num_group, int num_delay);
//! \brief Constructor that initializes and populates all data to build a
//! macroscopic cross section from microscopic cross section.
//! macroscopic cross section from microscopic cross sections.
//!
//! @param in_name Name of the object.
//! @param mat_kTs temperatures (in units of eV) that data is needed.

3
include/openmc/scattdata.h
View file

@ -33,7 +33,8 @@ class ScattData {
//! \brief Combines microscopic ScattDatas into a macroscopic one.
void
base_combine(size_t max_order, const std::vector<ScattData*>& those_scatts,
base_combine(size_t max_order, size_t order_dim,
const std::vector<ScattData*>& those_scatts,
const std::vector<double>& scalars, xt::xtensor<int, 1>& in_gmin,
xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter);

6
src/material.cpp
View file

@ -211,10 +211,12 @@ Material::Material(pugi::xml_node node)
for (int i = 0; i < n; ++i) {
const auto& name {names[i]};
// Check that this nuclide is listed in the cross_sections.xml file
// Check that this nuclide is listed in the nuclear data library
// (cross_sections.xml for CE and the MGXS HDF5 for MG)
LibraryKey key {Library::Type::neutron, name};
if (data::library_map.find(key) == data::library_map.end()) {
fatal_error("Could not find nuclide " + name + " in cross_sections.xml.");
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

60
src/mgxs.cpp
View file

@ -318,7 +318,8 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
// Get the minimum data needed to initialize:
// Dont need awr, but lets just initialize it anyways
double in_awr = -1.;
// start with the assumption it is not fissionable
// start with the assumption it is not fissionable and set
// the fissionable status if we learn differently
bool in_fissionable = false;
for (int m = 0; m < micros.size(); m++) {
if (micros[m]->fissionable) in_fissionable = true;
@ -377,25 +378,43 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
} // end switch
} // end microscopic temperature loop
// We are about to loop through each of the microscopic objects
// and incorporate the contribution of each microscopic data at
// one of the two temperature interpolants to this macroscopic quantity.
// If we are doing nearest temperature interpolation, then we don't need
// to do the 2nd temperature
int num_interp_points = 2;
if (settings::temperature_method == TemperatureMethod::NEAREST) num_interp_points = 1;
std::vector<double> interp(micros.size());
std::vector<int> temp_indices(micros.size());
for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
for (int m = 0; m < micros.size(); m++) {
interp[m] = (1. - micro_t_interp[m]) * atom_densities[m];
temp_indices[m] = micro_t[m] + interp_point;
micro_t_interp[m] = 1. - micro_t_interp[m];
// Now combine the microscopic data at each relevant temperature
// We will do this by treating the multiple temperatures of a nuclide as
// a different nuclide. Mathematically this just means the temperature
// interpolant is included in the number density.
// These interpolants are contained within interp.
std::vector<double> interpolant; // the interpolant for the Mgxs
std::vector<int> temp_indices; // the temperature index for each Mgxs
std::vector<Mgxs*> mgxs_to_combine; // The Mgxs to combine
// Now go through and build the above vectors so that we can use them to
// combine the data. We will step through each microscopic data and
// add in its lower and upper temperature points
for (int m = 0; m < micros.size(); m++) {
if (settings::temperature_method == TemperatureMethod::NEAREST) {
// Nearest interpolation only has one temperature point per isotope
// and so we dont need to include a temperature interpolant in
// the interpolant vector
interpolant.push_back(atom_densities[m]);
temp_indices.push_back(micro_t[m]);
mgxs_to_combine.push_back(micros[m]);
} else {
// This will be an interpolation between two points so get both these
// points
// Start with the low point
interpolant.push_back((1. - micro_t_interp[m]) * atom_densities[m]);
temp_indices.push_back(micro_t[m]);
mgxs_to_combine.push_back(micros[m]);
// The higher point
interpolant.push_back((micro_t_interp[m]) * atom_densities[m]);
temp_indices.push_back(micro_t[m] + 1);
mgxs_to_combine.push_back(micros[m]);
}
}
combine(micros, interp, temp_indices, t);
} // end loop to sum all micros across the temperatures
// And finally, combine the data
combine(mgxs_to_combine, interpolant, temp_indices, t);
} // end temperature (t) loop
}
//==============================================================================
@ -407,9 +426,6 @@ Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scala
// Build the vector of pointers to the xs objects within micros
std::vector<XsData*> those_xs(micros.size());
for (int i = 0; i < micros.size(); i++) {
if (!xs[this_t].equiv(micros[i]->xs[micro_ts[i]])) {
fatal_error("Cannot combine the Mgxs objects!");
}
those_xs[i] = &(micros[i]->xs[micro_ts[i]]);
}
@ -448,9 +464,9 @@ Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
case MgxsType::KAPPA_FISSION:
val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
break;
case MgxsType::NU_SCATTER:
case MgxsType::SCATTER:
case MgxsType::SCATTER_MULT:
case MgxsType::SCATTER_FMU_MULT:
case MgxsType::NU_SCATTER_FMU:
case MgxsType::SCATTER_FMU:
val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
break;

121
src/scattdata.cpp
View file

@ -5,6 +5,7 @@
#include <cmath>
#include "xtensor/xbuilder.hpp"
#include "xtensor/xview.hpp"
#include "openmc/constants.h"
#include "openmc/error.h"
@ -36,6 +37,13 @@ ScattData::base_init(int order, const xt::xtensor<int, 1>& in_gmin,
energy[gin] = in_energy[gin];
mult[gin] = in_mult[gin];
// Make sure the multiplicity does not have 0s
for (int go = 0; go < mult[gin].size(); go++) {
if (mult[gin][go] == 0.) {
mult[gin][go] = 1.;
}
}
// Make sure the energy is normalized
double norm = std::accumulate(energy[gin].begin(), energy[gin].end(), 0.);
@ -54,77 +62,45 @@ ScattData::base_init(int order, const xt::xtensor<int, 1>& in_gmin,
//==============================================================================
void
ScattData::base_combine(size_t max_order,
const std::vector<ScattData*>& those_scatts, const std::vector<double>& scalars,
xt::xtensor<int, 1>& in_gmin, xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
ScattData::base_combine(size_t max_order, size_t order_dim,
const std::vector<ScattData*>& those_scatts,
const std::vector<double>& scalars, xt::xtensor<int, 1>& in_gmin,
xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter)
{
size_t groups = those_scatts[0] -> energy.size();
// Now allocate and zero our storage spaces
xt::xtensor<double, 3> this_matrix({groups, groups, max_order}, 0.);
xt::xtensor<double, 2> mult_numer({groups, groups}, 0.);
xt::xtensor<double, 2> mult_denom({groups, groups}, 0.);
// TODO: Need to review this:
if (this->scattxs.size() > 0) {
this_matrix = this->get_matrix(max_order);
}
xt::xtensor<double, 3> this_nuscatt_matrix({groups, groups, order_dim}, 0.);
xt::xtensor<double, 2> this_nuscatt_P0({groups, groups}, 0.);
xt::xtensor<double, 2> this_scatt_P0({groups, groups}, 0.);
xt::xtensor<double, 2> this_mult({groups, groups}, 1.);
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattData* that = those_scatts[i];
// Build the dense matrix for that object
xt::xtensor<double, 3> that_matrix = that->get_matrix(max_order);
// Now add that to this for the scattering and multiplicity
// Now add that to this for the nu-scatter matrix
this_nuscatt_matrix += scalars[i] * that_matrix;
// Do the same with the P0 matrices
for (int gin = 0; gin < groups; gin++) {
// Only spend time adding that's gmin to gmax data since the rest will
// be zeros
int i_gout = 0;
for (int gout = that->gmin(gin); gout <= that->gmax(gin); gout++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix(gin, gout, l) += scalars[i] * that_matrix(gin, gout, l);
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs(gin) * that->energy[gin][i_gout];
mult_numer(gin, gout) += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom(gin, gout) += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom(gin, gout) += scalars[i];
}
i_gout++;
for (int go = 0; go < groups; go++) {
this_nuscatt_P0(gin, go) +=
scalars[i] * that->get_xs(MgxsType::NU_SCATTER, gin, &go, nullptr);
this_scatt_P0(gin, go) +=
scalars[i] *
that->get_xs(MgxsType::SCATTER, gin, &go, nullptr);
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
xt::xtensor<double, 2> this_mult({groups, groups}, 1.);
// TODO: Need to check this too
for (int gin = 0; gin < groups; gin++) {
for (int gout = 0; gout < groups; gout++) {
if (std::abs(mult_denom(gin, gout)) > 0.0) {
this_mult(gin, gout) = mult_numer(gin, gout) / mult_denom(gin, gout);
} else {
if (mult_numer(gin, gout) == 0.0) {
this_mult(gin, gout) = 1.0;
}
}
}
}
// Now we have the dense nuscatt and scatt, we can easily compute the
// multiplicity matrix by dividing the two and fixing any nans
this_mult = xt::nan_to_num(this_nuscatt_P0 / this_scatt_P0);
// We have the data, now we need to convert to a jagged array and then use
// the initialize function to store it on the object.
@ -133,8 +109,8 @@ ScattData::base_combine(size_t max_order,
int gmin_;
for (gmin_ = 0; gmin_ < groups; gmin_++) {
bool non_zero = false;
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmin_, l) != 0.) {
for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) {
if (this_nuscatt_matrix(gin, gmin_, l) != 0.) {
non_zero = true;
break;
}
@ -144,8 +120,8 @@ ScattData::base_combine(size_t max_order,
int gmax_;
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
bool non_zero = false;
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmax_, l) != 0.) {
for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) {
if (this_nuscatt_matrix(gin, gmax_, l) != 0.) {
non_zero = true;
break;
}
@ -168,9 +144,9 @@ ScattData::base_combine(size_t max_order,
sparse_mult[gin].resize(gmax_ - gmin_ + 1);
int i_gout = 0;
for (int gout = gmin_; gout <= gmax_; gout++) {
sparse_scatter[gin][i_gout].resize(this_matrix.shape()[2]);
for (int l = 0; l < this_matrix.shape()[2]; l++) {
sparse_scatter[gin][i_gout][l] = this_matrix(gin, gout, l);
sparse_scatter[gin][i_gout].resize(this_nuscatt_matrix.shape()[2]);
for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) {
sparse_scatter[gin][i_gout][l] = this_nuscatt_matrix(gin, gout, l);
}
sparse_mult[gin][i_gout] = this_mult(gin, gout);
i_gout++;
@ -178,6 +154,7 @@ ScattData::base_combine(size_t max_order,
}
}
//==============================================================================
void
@ -212,10 +189,10 @@ ScattData::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu)
double val = scattxs[gin];
switch(xstype) {
case MgxsType::SCATTER:
case MgxsType::NU_SCATTER:
if (gout != nullptr) val *= energy[gin][i_gout];
break;
case MgxsType::SCATTER_MULT:
case MgxsType::SCATTER:
if (gout != nullptr) {
val *= energy[gin][i_gout] / mult[gin][i_gout];
} else {
@ -223,7 +200,7 @@ ScattData::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu)
energy[gin].begin(), 0.0);
}
break;
case MgxsType::SCATTER_FMU_MULT:
case MgxsType::NU_SCATTER_FMU:
if ((gout != nullptr) && (mu != nullptr)) {
val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu);
} else {
@ -407,7 +384,6 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
size_t that_order = that->get_order();
if (that_order > max_order) max_order = that_order;
}
max_order++; // Add one since this is a Legendre
size_t groups = those_scatts[0] -> energy.size();
@ -419,8 +395,9 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
size_t order_dim = max_order + 1;
ScattData::base_combine(max_order, order_dim, those_scatts, scalars, in_gmin,
in_gmax, sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -636,8 +613,9 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
size_t order_dim = max_order;
ScattData::base_combine(max_order, order_dim, those_scatts, scalars, in_gmin,
in_gmax, sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -854,8 +832,9 @@ ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
size_t order_dim = max_order;
ScattData::base_combine(max_order, order_dim, those_scatts, scalars, in_gmin,
in_gmax, sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);

16
src/tallies/tally_scoring.cpp
View file

@ -1566,17 +1566,17 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index,
score = p.wgt_last_ * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(
MgxsType::SCATTER_FMU_MULT, p.g_last_, &p.g_, &p.mu_, nullptr)
MgxsType::SCATTER_FMU, p.g_last_, &p.g_, &p.mu_, nullptr)
/ macro_xs.get_xs(
MgxsType::SCATTER_FMU_MULT, p.g_last_, &p.g_, &p.mu_, nullptr);
MgxsType::SCATTER_FMU, p.g_last_, &p.g_, &p.mu_, nullptr);
}
} else {
if (i_nuclide >= 0) {
score = atom_density * flux * nuc_xs.get_xs(
MgxsType::SCATTER_MULT, p_g, nullptr, &p.mu_, nullptr);
MgxsType::SCATTER, p_g, nullptr, &p.mu_, nullptr);
} else {
score = flux * macro_xs.get_xs(
MgxsType::SCATTER_MULT, p_g, nullptr, &p.mu_, nullptr);
MgxsType::SCATTER, p_g, nullptr, &p.mu_, nullptr);
}
}
break;
@ -1595,16 +1595,16 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index,
// adjust the score by the actual probability for that nuclide.
if (i_nuclide >= 0) {
score *= atom_density
* nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last_, &p.g_,
* nuc_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last_, &p.g_,
&p.mu_, nullptr)
/ macro_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last_, &p.g_,
/ macro_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last_, &p.g_,
&p.mu_, nullptr);
}
} else {
if (i_nuclide >= 0) {
score = atom_density * flux * nuc_xs.get_xs(MgxsType::SCATTER, p_g);
score = atom_density * flux * nuc_xs.get_xs(MgxsType::NU_SCATTER, p_g);
} else {
score = flux * macro_xs.get_xs(MgxsType::SCATTER, p_g);
score = flux * macro_xs.get_xs(MgxsType::NU_SCATTER, p_g);
}
}
break;

17
src/xsdata.cpp
View file

@ -529,12 +529,25 @@ XsData::combine(const std::vector<XsData*>& those_xs,
kappa_fission += scalar * that->kappa_fission;
fission += scalar * that->fission;
delayed_nu_fission += scalar * that->delayed_nu_fission;
chi_prompt += scalar * that->chi_prompt;
chi_delayed += scalar * that->chi_delayed;
chi_prompt += scalar *
xt::view(xt::sum(that->prompt_nu_fission, {1}),
xt::all(), xt::newaxis(), xt::newaxis()) *
that->chi_prompt;
chi_delayed += scalar *
xt::view(xt::sum(that->delayed_nu_fission, {2}),
xt::all(), xt::all(), xt::newaxis(), xt::newaxis()) *
that->chi_delayed;
}
decay_rate += scalar * that->decay_rate;
}
// Ensure the chi_prompt and chi_delayed are normalized to 1 for each
// azimuthal angle and delayed group (for chi_delayed)
chi_prompt /=
xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
chi_delayed /= xt::view(xt::sum(chi_delayed, {3}), xt::all(), xt::all(),
xt::all(), xt::newaxis());
// Allow the ScattData object to combine itself
for (size_t a = 0; a < total.shape()[0]; a++) {
// Build vector of the scattering objects to incorporate