From ce919a0ad85b9ee0023ded1224bb05686edb1e94 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sat, 9 Jun 2018 08:53:40 -0400 Subject: [PATCH] extended to being able to combine microscopic data into macroscopic. Seems to work, next step is to actually incorporate the C++ data into the transport and tallying process. That will be the true test --- src/input_xml.F90 | 3 +- src/material_header.F90 | 2 +- src/mgxs.cpp | 84 ++++++++++++++++------------ src/mgxs.h | 12 ++-- src/mgxs_data.F90 | 47 ++++++++++++++++ src/scattdata.cpp | 119 +++++++++++++++++----------------------- src/scattdata.h | 19 +++---- src/xsdata.cpp | 67 +++++++++++----------- 8 files changed, 199 insertions(+), 154 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8d4c8de817..507bdf8a9c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,7 +19,7 @@ module input_xml use material_header use mesh_header use message_passing - use mgxs_data, only: create_macro_xs, read_mgxs, read_mgxs2 + use mgxs_data, only: create_macro_xs, read_mgxs, read_mgxs2, create_macro_xs2 use mgxs_header use nuclide_header use output, only: title, header, print_plot @@ -86,6 +86,7 @@ contains call read_mgxs() call read_mgxs2() call create_macro_xs() + call create_macro_xs2() end if call time_read_xs % stop() end if diff --git a/src/material_header.F90 b/src/material_header.F90 index 038e529bee..7fbfe553b3 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -34,7 +34,7 @@ module material_header integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm - real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm + real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm real(8) :: density_gpcc ! total density in g/cm^3 ! To improve performance of tallying, we store an array (direct address diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 9c11da8ac8..35ac42f171 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -266,11 +266,14 @@ void Mgxs::from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups, close_group(xsdata_grp); } // end temperature loop + + // Make sure the scattering format is updated to the final case + scatter_format = final_scatter_format; } void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, - std::vector& micros, double_1dvec& atom_densities, + std::vector& micros, double_1dvec& atom_densities, int& method, double tolerance) { // Get the minimum data needed to initialize: @@ -279,21 +282,25 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, // start with the assumption it is not fissionable bool in_fissionable = false; for (int m = 0; m < micros.size(); m++) { - if (micros[m].fissionable) in_fissionable = true; + if (micros[m]->fissionable) in_fissionable = true; } // Force all of the following data to be the same; these will be verified // to be true later - int in_scatter_format = micros[0].scatter_format; - int in_num_groups = micros[0].num_groups; - int in_num_delayed_groups = micros[0].num_delayed_groups; - double_1dvec in_polar = micros[0].polar; - double_1dvec in_azimuthal = micros[0].azimuthal; + int in_scatter_format = micros[0]->scatter_format; + int in_num_groups = micros[0]->num_groups; + int in_num_delayed_groups = micros[0]->num_delayed_groups; + double_1dvec in_polar = micros[0]->polar; + double_1dvec in_azimuthal = micros[0]->azimuthal; - init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format, in_num_groups, - in_num_delayed_groups, in_polar, in_azimuthal); + init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format, + in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal); // Create the xs data for each temperature for (int t = 0; t < mat_kTs.size(); t++) { + xs[t]= XsData(in_num_groups, in_num_delayed_groups, in_fissionable, + in_scatter_format, in_polar.size(), in_azimuthal.size()); + + // Find the right temperature index to use double temp_desired = mat_kTs[t]; // Create the list of temperature indices and interpolation factors for @@ -305,13 +312,13 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, case TEMPERATURE_NEAREST: { // Find the nearest temperature - std::valarray temp_diff(micros[m].kTs.data(), - micros[m].kTs.size()); + std::valarray temp_diff(micros[m]->kTs.data(), + micros[m]->kTs.size()); temp_diff = std::abs(temp_diff - temp_desired); micro_t[m] = std::min_element(std::begin(temp_diff), std::end(temp_diff)) - std::begin(temp_diff); - double temp_actual = micros[m].kTs[micro_t[m]]; + double temp_actual = micros[m]->kTs[micro_t[m]]; if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) { fatal_error("MGXS Library does not contain cross section for " + @@ -324,13 +331,13 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, case TEMPERATURE_INTERPOLATION: // Get a list of bounding temperatures for each actual temperature // present in the model - for (int k = 0; k < micros[m].kTs.size() - 1; k++) { - if ((micros[m].kTs[k] <= temp_desired) && - (temp_desired < micros[m].kTs[k + 1])) { + for (int k = 0; k < micros[m]->kTs.size() - 1; k++) { + if ((micros[m]->kTs[k] <= temp_desired) && + (temp_desired < micros[m]->kTs[k + 1])) { micro_t[m] = k; if (k == 0) { - micro_t_interp[m] = (temp_desired - micros[m].kTs[k]) / - (micros[m].kTs[k + 1] - micros[m].kTs[k]); + micro_t_interp[m] = (temp_desired - micros[m]->kTs[k]) / + (micros[m]->kTs[k + 1] - micros[m]->kTs[k]); } else { micro_t_interp[m] = 1.; } @@ -353,23 +360,23 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs, interp[m] = (1. - micro_t_interp[m]) * atom_densities[m]; temp_indices[m] = micro_t[m] + interp_point; } + combine(micros, interp, micro_t, t); } // end loop to sum all micros across the temperatures } // end temperature (t) loop - } -void Mgxs::combine(std::vector& micros, double_1dvec& scalars, +void Mgxs::combine(std::vector& micros, double_1dvec& scalars, int_1dvec& micro_ts, int this_t) { // 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]])) { + 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]]); + those_xs[i] = &(micros[i]->xs[micro_ts[i]]); } xs[this_t].combine(those_xs, scalars); @@ -646,24 +653,31 @@ bool query_fissionable(const int n_nuclides, const int i_nuclides[]) } -void create_macro_xs(int n_materials, double_2dvec& mat_kTs, - std::vector& mat_names, - double_1dvec& atom_densities, int& method, - double tolerance) +void create_macro_xs(char* mat_name, const int n_nuclides, + const int i_nuclides[], const int n_temps, const double temps[], + const double atom_densities[], int& method, const double tolerance) { - // TODO mat_kTs needs to be converted from Fortran - // it is currently an array of type(VectorReal), a wrapper should convert to - // the vector. - macro_xs.resize(n_materials); + Mgxs macro; + if (n_temps > 0) { + // // Convert temps to a vector + double_1dvec temperature; + temperature.assign(temps, temps + n_temps); - for (int m = 0; m < n_materials; m++) - { - if (mat_kTs[m].size() > 0) { - macro_xs[m].build_macro(mat_names[m], mat_kTs[m], nuclides_MG, - atom_densities, method, tolerance); + // Convert atom_densities to a vector + double_1dvec atom_densities_vec; + atom_densities_vec.assign(atom_densities, atom_densities + n_nuclides); + + // Build array of pointers to nuclides_MG's Mgxs objects needed for this + // material + std::vector mgxs_ptr(n_nuclides); + for (int n = 0; n < n_nuclides; n++) { + mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1]; } + + macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec, + method, tolerance); } + macro_xs.push_back(macro); } - } // namespace openmc \ No newline at end of file diff --git a/src/mgxs.h b/src/mgxs.h index 8c13cf7cb4..90ed520d9a 100644 --- a/src/mgxs.h +++ b/src/mgxs.h @@ -58,9 +58,9 @@ class Mgxs { const int in_num_delayed_groups, const double_1dvec& in_polar, const double_1dvec& in_azimuthal); void build_macro(const std::string& in_name, double_1dvec& mat_kTs, - std::vector& micros, double_1dvec& atom_densities, + std::vector& micros, double_1dvec& atom_densities, int& method, double tolerance); - void combine(std::vector& micros, double_1dvec& scalars, + void combine(std::vector& micros, double_1dvec& scalars, int_1dvec& micro_ts, int this_t); void from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups, double_1dvec& temperature, int& method, double tolerance, @@ -82,10 +82,12 @@ extern "C" void add_mgxs(hid_t file_id, char* name, int energy_groups, int delayed_groups, int n_temps, double temps[], int& method, double tolerance, int max_order, bool legendre_to_tabular, int legendre_to_tabular_points); + extern "C" bool query_fissionable(const int n_nuclides, const int i_nuclides[]); -void create_macro_xs(int n_materials, double_2dvec& mat_kTs, - std::vector& mat_names, double_1dvec& atom_densities, - int& method, double tolerance); + +extern "C" void create_macro_xs(char* mat_name, const int n_nuclides, + const int i_nuclides[], const int n_temps, const double temps[], + const double atom_densities[], int& method, const double tolerance); // Storage for the MGXS data diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 6be463a91e..a61bb65a2b 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -45,6 +45,20 @@ module mgxs_data integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) logical(C_BOOL) :: result end function query_fissionable_c + + subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, & + atom_densities, method, tolerance) bind(C, name='create_macro_xs') + use ISO_C_BINDING + implicit none + character(kind=C_CHAR),intent(in) :: name(*) + integer(C_INT), value, intent(in) :: n_nuclides + integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) + integer(C_INT), value, intent(in) :: n_temps + real(C_DOUBLE), intent(in) :: temps(1:n_temps) + real(C_DOUBLE), intent(in) :: atom_densities(1:n_nuclides) + integer(C_INT), intent(inout) :: method + real(C_DOUBLE), value, intent(in) :: tolerance + end subroutine create_macro_xs_c end interface contains @@ -316,6 +330,39 @@ contains end subroutine create_macro_xs + + subroutine create_macro_xs2() + integer :: i_mat ! index in materials array + type(Material), pointer :: mat ! current material + type(VectorReal), allocatable :: kTs(:) + character(MAX_WORD_LEN) :: name ! name of material + + ! Get temperatures to read for each material + call get_mat_kTs(kTs) + + ! Force all nuclides in a material to be the same representation. + ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs). + ! At the same time, we will find the scattering type, as that will dictate + ! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials)) + do i_mat = 1, n_materials + + ! Get the material + mat => materials(i_mat) + + name = trim(mat % name) // C_NULL_CHAR + + ! Do not read materials which we do not actually use in the problem to + ! reduce storage + if (allocated(kTs(i_mat) % data)) then + call create_macro_xs_c(name, mat % n_nuclides, mat % nuclide, & + kTs(i_mat) % size(), kTs(i_mat) % data, mat % atom_density, & + temperature_method, temperature_tolerance) + end if + end do + + end subroutine create_macro_xs2 + + !=============================================================================== ! GET_MAT_kTs returns a list of temperatures (in eV) that each ! material appears at in the model. diff --git a/src/scattdata.cpp b/src/scattdata.cpp index cb72f42298..dc4d169a7c 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -109,11 +109,11 @@ double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu) // ScattDataLegendre methods //============================================================================== -void ScattDataLegendre::init(int_1dvec in_gmin, int_1dvec in_gmax, - double_2dvec in_mult, double_3dvec coeffs) +void ScattDataLegendre::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); - int order = coeffs[0].size(); + int order = coeffs[0][0].size(); // make a copy of coeffs that we can use to both extract data and normalize double_3dvec matrix = coeffs; @@ -250,13 +250,12 @@ void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) void ScattDataLegendre::combine(std::vector those_scatts, double_1dvec& scalars) { - int groups = energy.size(); - // Find the maximum order in the data set - int max_order = get_order(); + // Find the max order in the data set and make sure we can combine the sets + int max_order = 0; for (int i = 0; i < those_scatts.size(); i++) { // Lets also make sure these items are combineable ScattDataLegendre* that = dynamic_cast(those_scatts[i]); - if (!equiv(*that)) { + if (!that) { fatal_error("Cannot combine the ScattData objects!"); } int that_order = that->get_order(); @@ -264,6 +263,9 @@ void ScattDataLegendre::combine(std::vector those_scatts, } max_order++; // Add one since this is a Legendre + // Get the groups as a shorthand + int groups = dynamic_cast(those_scatts[0])->energy.size(); + // Now allocate and zero our storage spaces double_3dvec this_matrix = get_matrix(max_order); double_2dvec mult_numer(groups, double_1dvec(groups, 0.)); @@ -369,23 +371,16 @@ void ScattDataLegendre::combine(std::vector those_scatts, } -bool ScattDataLegendre::equiv(const ScattDataLegendre& that) -{ - // ensure that the number of groups match - return (this->energy.size() == that.energy.size()); -} - - double_3dvec ScattDataLegendre::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); int order_dim = max_order + 1; - double_3dvec matrix = double_3dvec(groups, double_2dvec(order_dim, + double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, double_1dvec(order_dim, 0.))); for (int gin = 0; gin < groups; gin++) { - for (int i_gout = 0; i_gout < energy[0].size(); i_gout++) { + for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * @@ -400,11 +395,11 @@ double_3dvec ScattDataLegendre::get_matrix(int max_order) // ScattDataHistogram methods //============================================================================== -void ScattDataHistogram::init(int_1dvec in_gmin, int_1dvec in_gmax, - double_2dvec in_mult, double_3dvec coeffs) +void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); - int order = coeffs[0].size(); + int order = coeffs[0][0].size(); // make a copy of coeffs that we can use to both extract data and normalize double_3dvec matrix = coeffs; @@ -452,7 +447,7 @@ void ScattDataHistogram::init(int_1dvec in_gmin, int_1dvec in_gmax, for (int gin = 0; gin < groups; gin++) { int num_groups = gmax[gin] - gmin[gin] + 1; fmu[gin].resize(num_groups); - for (int i_gout = 0; i_gout < num_groups; i_gout) { + for (int i_gout = 0; i_gout < num_groups; i_gout++) { fmu[gin][i_gout].resize(order); // The variable matrix contains f(mu); so directly assign it fmu[gin][i_gout] = matrix[gin][i_gout]; @@ -533,29 +528,17 @@ void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) } -bool ScattDataHistogram::equiv(const ScattDataHistogram& that) -{ - bool match = false; - if (this->energy.size() == that.energy.size() && - this->dmu == that.dmu && - std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) { - match = true; - } - return match; -} - - double_3dvec ScattDataHistogram::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations int order_dim = get_order(); - double_3dvec matrix = double_3dvec(groups, double_2dvec(order_dim, + double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, double_1dvec(order_dim, 0.))); for (int gin = 0; gin < groups; gin++) { - for (int i_gout = 0; i_gout < energy[0].size(); i_gout++) { + for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * @@ -570,19 +553,23 @@ double_3dvec ScattDataHistogram::get_matrix(int max_order) void ScattDataHistogram::combine(std::vector those_scatts, double_1dvec& scalars) { - int groups = energy.size(); - // Find the maximum order in the data set - int max_order = get_order(); + // Find the max order in the data set and make sure we can combine the sets + int max_order; for (int i = 0; i < those_scatts.size(); i++) { // Lets also make sure these items are combineable ScattDataHistogram* that = dynamic_cast(those_scatts[i]); - if (!equiv(*that)) { + if (!that) { + fatal_error("Cannot combine the ScattData objects!"); + } + if (i == 0) { + max_order = that->get_order(); + } else if (max_order != that->get_order()) { fatal_error("Cannot combine the ScattData objects!"); } - int that_order = that->get_order(); - if (that_order > max_order) max_order = that_order; } - max_order++; // Add one since this is a Legendre + + // Get the groups as a shorthand + int groups = dynamic_cast(those_scatts[0])->energy.size(); // Now allocate and zero our storage spaces double_3dvec this_matrix = get_matrix(max_order); @@ -692,11 +679,11 @@ void ScattDataHistogram::combine(std::vector those_scatts, // ScattDataTabular methods //============================================================================== -void ScattDataTabular::init(int_1dvec in_gmin, int_1dvec in_gmax, - double_2dvec in_mult, double_3dvec coeffs) +void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) { int groups = coeffs.size(); - int order = coeffs[0].size(); + int order = coeffs[0][0].size(); // make a copy of coeffs that we can use to both extract data and normalize double_3dvec matrix = coeffs; @@ -742,15 +729,14 @@ void ScattDataTabular::init(int_1dvec in_gmin, int_1dvec in_gmax, } // Initialize the base class attributes - ScattData::generic_init(order, in_gmin, in_gmax, in_energy, - in_mult); + ScattData::generic_init(order, in_gmin, in_gmax, in_energy, in_mult); // Calculate f(mu) and integrate it so we can avoid rejection sampling fmu.resize(groups); for (int gin = 0; gin < groups; gin++) { int num_groups = gmax[gin] - gmin[gin] + 1; fmu[gin].resize(num_groups); - for (int i_gout = 0; i_gout < num_groups; i_gout) { + for (int i_gout = 0; i_gout < num_groups; i_gout++) { fmu[gin][i_gout].resize(order); // The variable matrix contains f(mu); so directly assign it fmu[gin][i_gout] = matrix[gin][i_gout]; @@ -852,29 +838,17 @@ void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) } -bool ScattDataTabular::equiv(const ScattDataTabular& that) -{ - bool match = false; - if (this->energy.size() == that.energy.size() && - this->dmu == that.dmu && - std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) { - match = true; - } - return match; -} - - double_3dvec ScattDataTabular::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations int order_dim = get_order(); - double_3dvec matrix = double_3dvec(groups, double_2dvec(order_dim, + double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, double_1dvec(order_dim, 0.))); for (int gin = 0; gin < groups; gin++) { - for (int i_gout = 0; i_gout < energy[0].size(); i_gout++) { + for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * @@ -888,22 +862,27 @@ double_3dvec ScattDataTabular::get_matrix(int max_order) void ScattDataTabular::combine(std::vector those_scatts, double_1dvec& scalars) { - int groups = energy.size(); - // Find the maximum order in the data set - int max_order = get_order(); + // Find the max order in the data set and make sure we can combine the sets + int max_order; for (int i = 0; i < those_scatts.size(); i++) { // Lets also make sure these items are combineable ScattDataTabular* that = dynamic_cast(those_scatts[i]); - if (!equiv(*that)) { + if (!that) { + fatal_error("Cannot combine the ScattData objects!"); + } + if (i == 0) { + max_order = that->get_order(); + } else if (max_order != that->get_order()) { fatal_error("Cannot combine the ScattData objects!"); } - int that_order = that->get_order(); - if (that_order > max_order) max_order = that_order; } - max_order++; // Add one since this is a Legendre + + // Get the groups as a shorthand + int groups = dynamic_cast(those_scatts[0])->energy.size(); // Now allocate and zero our storage spaces - double_3dvec this_matrix = get_matrix(max_order); + double_3dvec this_matrix = double_3dvec(groups, double_2dvec(groups, + double_1dvec(max_order, 0.))); double_2dvec mult_numer(groups, double_1dvec(groups, 0.)); double_2dvec mult_denom(groups, double_1dvec(groups, 0.)); @@ -999,6 +978,7 @@ void ScattDataTabular::combine(std::vector those_scatts, for (int gout = gmin_; gout <= gmax_; gout++) { sparse_scatter[gin][i_gout] = this_matrix[gin][gout]; sparse_mult[gin][i_gout] = this_mult[gin][gout]; + i_gout++; } } @@ -1011,6 +991,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, int n_mu) { tab.generic_init(n_mu, leg.gmin, leg.gmax, leg.energy, leg.mult); + tab.scattxs = leg.scattxs; // Build mu and dmu tab.mu = double_1dvec(n_mu); diff --git a/src/scattdata.h b/src/scattdata.h index 787a8f005d..6e24a50bc9 100644 --- a/src/scattdata.h +++ b/src/scattdata.h @@ -35,8 +35,8 @@ class ScattData { double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}} virtual double calc_f(int gin, int gout, double mu) = 0; virtual void sample(int gin, int& gout, double& mu, double& wgt) = 0; - virtual void init(int_1dvec in_gmin, int_1dvec in_gmax, - double_2dvec in_mult, double_3dvec coeffs) = 0; + virtual void init(int_1dvec& in_gmin, int_1dvec& in_gmax, + double_2dvec& in_mult, double_3dvec& coeffs) = 0; void sample_energy(int gin, int& gout, int& i_gout); double get_xs(const char* xstype, int gin, int* gout, double* mu); void generic_init(int order, int_1dvec in_gmin, int_1dvec in_gmax, @@ -54,12 +54,11 @@ class ScattDataLegendre: public ScattData { friend void convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, int n_mu); public: - void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, - double_3dvec coeffs); + void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult, + double_3dvec& coeffs); void update_max_val(); double calc_f(int gin, int gout, double mu); void sample(int gin, int& gout, double& mu, double& wgt); - bool equiv(const ScattDataLegendre& that); void combine(std::vector those_scatts, double_1dvec& scalars); int get_order() {return dist[0][0].size() - 1;}; double_3dvec get_matrix(int max_order); @@ -71,12 +70,11 @@ class ScattDataHistogram: public ScattData { double dmu; double_3dvec fmu; public: - void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, - double_3dvec coeffs); + void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult, + double_3dvec& coeffs); double calc_f(int gin, int gout, double mu); void sample(int gin, int& gout, double& mu, double& wgt); void combine(std::vector those_scatts, double_1dvec& scalars); - bool equiv(const ScattDataHistogram& that); int get_order() {return dist[0][0].size();}; double_3dvec get_matrix(int max_order); }; @@ -89,12 +87,11 @@ class ScattDataTabular: public ScattData { friend void convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, int n_mu); public: - void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, - double_3dvec coeffs); + void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult, + double_3dvec& coeffs); double calc_f(int gin, int gout, double mu); void sample(int gin, int& gout, double& mu, double& wgt); void combine(std::vector those_scatts, double_1dvec& scalars); - bool equiv(const ScattDataTabular& that); int get_order() {return dist[0][0].size();}; double_3dvec get_matrix(int max_order); }; diff --git a/src/xsdata.cpp b/src/xsdata.cpp index 42e32baa37..49b07453df 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -666,26 +666,27 @@ void XsData::combine(std::vector those_xs, double_1dvec& scalars) absorption[p][a][gin] += scalar * that->absorption[p][a][gin]; inverse_velocity[p][a][gin] += scalar * that->inverse_velocity[p][a][gin]; + if (that->prompt_nu_fission.size() > 0) { + prompt_nu_fission[p][a][gin] += + scalar * that->prompt_nu_fission[p][a][gin]; + kappa_fission[p][a][gin] += + scalar * that->kappa_fission[p][a][gin]; + fission[p][a][gin] += + scalar * that->fission[p][a][gin]; - prompt_nu_fission[p][a][gin] += - scalar * that->prompt_nu_fission[p][a][gin]; - kappa_fission[p][a][gin] += - scalar * that->kappa_fission[p][a][gin]; - fission[p][a][gin] += - scalar * that->fission[p][a][gin]; + for (int dg = 0; dg < delayed_nu_fission[p][a][gin].size(); dg++) { + delayed_nu_fission[p][a][gin][dg] += + scalar * that->delayed_nu_fission[p][a][gin][dg]; + } - for (int dg = 0; dg < delayed_nu_fission[p][a][gin].size(); dg++) { - delayed_nu_fission[p][a][gin][dg] += - scalar * that->delayed_nu_fission[p][a][gin][dg]; - } + for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { + chi_prompt[p][a][gin][gout] += + scalar * that->chi_prompt[p][a][gin][gout]; - for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { - chi_prompt[p][a][gin][gout] += - scalar * that->chi_prompt[p][a][gin][gout]; - - for (int dg = 0; dg < chi_delayed[p][a][gin][gout].size(); dg++) { - chi_delayed[p][a][gin][gout][dg] += - scalar * that->chi_delayed[p][a][gin][gout][dg]; + for (int dg = 0; dg < chi_delayed[p][a][gin][gout].size(); dg++) { + chi_delayed[p][a][gin][gout][dg] += + scalar * that->chi_delayed[p][a][gin][gout][dg]; + } } } } @@ -695,23 +696,25 @@ void XsData::combine(std::vector those_xs, double_1dvec& scalars) } // Normalize chi - for (int gin = 0; gin < chi_prompt[p][a].size(); gin++) { - double norm = std::accumulate(chi_prompt[p][a][gin].begin(), - chi_prompt[p][a][gin].end(), 0.); - if (norm > 0.) { - for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { - chi_prompt[p][a][gin][gout] /= norm; - } - } - - for (int dg = 0; dg < chi_delayed[p][a][gin][0].size(); dg++) { - norm = 0.; - for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) { - norm += chi_delayed[p][a][gin][gout][dg]; - } + if (chi_prompt.size() > 0) { + for (int gin = 0; gin < chi_prompt[p][a].size(); gin++) { + double norm = std::accumulate(chi_prompt[p][a][gin].begin(), + chi_prompt[p][a][gin].end(), 0.); if (norm > 0.) { + for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { + chi_prompt[p][a][gin][gout] /= norm; + } + } + + for (int dg = 0; dg < chi_delayed[p][a][gin][0].size(); dg++) { + norm = 0.; for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) { - chi_delayed[p][a][gin][gout][dg] /= norm; + norm += chi_delayed[p][a][gin][gout][dg]; + } + if (norm > 0.) { + for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) { + chi_delayed[p][a][gin][gout][dg] /= norm; + } } } }