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

Merge pull request #1694 from nelsonag/nuclidic_mg

Browse source 
Improvement to microscopic MG mode calculations
This commit is contained in:
Paul Romano 2020-10-19 08:33:22 -05:00 • committed by GitHub
commit 29dfc0ae5f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 136 additions and 116 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);

14
openmc/mgxs/mgxs.py
View file

@ -2872,7 +2872,7 @@ class DiffusionCoefficient(TransportXS):
\phi \rangle}{\langle \phi \rangle} \\
D = \frac{1}{3 \sigma_{tr}}
\end{aligned}
To incorporate the effect of scattering multiplication in the above
relation, the `nu` parameter can be set to `True`.
@ -3057,7 +3057,7 @@ class DiffusionCoefficient(TransportXS):
filter_bins=[('P1',)],
squeeze=True)
p1_tally._scores = ['scatter-1']
total = self.tallies['total'] / self.tallies['flux (tracklength)']
trans_corr = p1_tally / self.tallies['flux (analog)']
transport = (total - trans_corr)
@ -3066,7 +3066,7 @@ class DiffusionCoefficient(TransportXS):
else:
dif_coef = self.rxn_rate_tally
flux_tally = condensed_xs.tallies['flux (tracklength)']
condensed_xs._tallies = OrderedDict()
condensed_xs._tallies[self._rxn_type] = dif_coef
@ -5633,8 +5633,8 @@ class Chi(MGXS):
"""
# Store whether or not the number density should be removed for microscopic
# values of this data; since this chi data is normalized to 1.0, the
# data should not be divided by the number density
# values of this data; since this chi data is normalized to 1.0, the
# data should not be divided by the number density
_divide_by_density = False
def __init__(self, domain=None, domain_type=None, groups=None,
@ -6554,7 +6554,7 @@ class MeshSurfaceMGXS(MGXS):
surfaces = df.pop(col_key)
df.insert(len(self.domain.dimension), col_key, surfaces)
if len(self.domain.dimension) == 1:
df.sort_values(by=[(mesh_str, 'x'), (mesh_str, 'surf')]
df.sort_values(by=[(mesh_str, 'x'), (mesh_str, 'surf')]
+ columns, inplace=True)
elif len(self.domain.dimension) == 2:
df.sort_values(by=[(mesh_str, 'x'), (mesh_str, 'y'),
@ -6642,7 +6642,7 @@ class Current(MeshSurfaceMGXS):
is None unless the multi-group cross section has been computed.
num_subdomains : int
The number of subdomains is equal to the number of mesh surfaces times
two to account for both the incoming and outgoing current from the
two to account for both the incoming and outgoing current from the
mesh cell surfaces.
num_nuclides : int
Unused in MeshSurfaceMGXS

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 interpolant.
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;

130
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,22 @@ 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
unsigned long int num_converted = 0;
for (int go = 0; go < mult[gin].size(); go++) {
if (mult[gin][go] == 0.) {
num_converted += 1;
mult[gin][go] = 1.;
}
}
if (num_converted > 0) {
// Raise a warning to the user if we did have to do the conversion
std::string msg = std::to_string(num_converted) +
"entries in the Multiplicity Matrix were changed from 0 to 1";
warning(msg);
}
// Make sure the energy is normalized
double norm = std::accumulate(energy[gin].begin(), energy[gin].end(), 0.);
@ -54,77 +71,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 +118,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 +129,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 +153,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 +163,7 @@ ScattData::base_combine(size_t max_order,
}
}
//==============================================================================
void
@ -212,10 +198,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 +209,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 +393,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 +404,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 +622,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 +841,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