From 77015c0b1490e96db397f09f9ba2507cf377be8d Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 15 Oct 2020 13:46:05 -0500 Subject: [PATCH] Fixed issue in MG-Mode which surfaced when running problems which have microscopic cross sections defined and multiple fissile isotopes in a material --- include/openmc/constants.h | 4 +- include/openmc/mgxs.h | 2 +- include/openmc/scattdata.h | 3 +- src/material.cpp | 6 +- src/mgxs.cpp | 60 ++++++++++------- src/scattdata.cpp | 121 ++++++++++++++-------------------- src/tallies/tally_scoring.cpp | 16 ++--- src/xsdata.cpp | 17 ++++- 8 files changed, 120 insertions(+), 109 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 13b62ee87..e75d8c03c 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -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, diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index f9c4722e9..ae7a54ad3 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -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. diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index 75fd7a2fd..aa91d5688 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -33,7 +33,8 @@ class ScattData { //! \brief Combines microscopic ScattDatas into a macroscopic one. void - base_combine(size_t max_order, const std::vector& those_scatts, + base_combine(size_t max_order, size_t order_dim, + const std::vector& those_scatts, const std::vector& scalars, xt::xtensor& in_gmin, xt::xtensor& in_gmax, double_2dvec& sparse_mult, double_3dvec& sparse_scatter); diff --git a/src/material.cpp b/src/material.cpp index 1d4e18108..dd3f7ffd7 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -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 diff --git a/src/mgxs.cpp b/src/mgxs.cpp index d5c144658..5b61e5ded 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -318,7 +318,8 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector& 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& 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 interp(micros.size()); - std::vector 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 interpolant; // the interpolant for the Mgxs + std::vector temp_indices; // the temperature index for each Mgxs + std::vector 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& micros, const std::vector& scala // Build the vector of pointers to the xs objects within micros std::vector 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; diff --git a/src/scattdata.cpp b/src/scattdata.cpp index e3c97512c..dca906cb1 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -5,6 +5,7 @@ #include #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& 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& in_gmin, //============================================================================== void -ScattData::base_combine(size_t max_order, - const std::vector& those_scatts, const std::vector& scalars, - xt::xtensor& in_gmin, xt::xtensor& in_gmax, double_2dvec& sparse_mult, +ScattData::base_combine(size_t max_order, size_t order_dim, + const std::vector& those_scatts, + const std::vector& scalars, xt::xtensor& in_gmin, + xt::xtensor& 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 this_matrix({groups, groups, max_order}, 0.); - xt::xtensor mult_numer({groups, groups}, 0.); - xt::xtensor mult_denom({groups, groups}, 0.); - // TODO: Need to review this: - if (this->scattxs.size() > 0) { - this_matrix = this->get_matrix(max_order); - } + xt::xtensor this_nuscatt_matrix({groups, groups, order_dim}, 0.); + xt::xtensor this_nuscatt_P0({groups, groups}, 0.); + xt::xtensor this_scatt_P0({groups, groups}, 0.); + xt::xtensor 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 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 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& 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& 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& 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& 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); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 396ab2f82..894cdecea 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -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; diff --git a/src/xsdata.cpp b/src/xsdata.cpp index d4698e2a0..ff360c052 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -529,12 +529,25 @@ XsData::combine(const std::vector& 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