Replaced all polar and azimuthal indices with just a 1d angle index. Also fixed mg_max_order

This commit is contained in:
Adam G Nelson 2018-06-16 13:37:50 -04:00
parent 863f91b1e7
commit 35c62affe6
7 changed files with 431 additions and 477 deletions

View file

@ -481,6 +481,28 @@ read_nd_vector(hid_t obj_id, const char* name,
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
std::vector<int> temp_arr = std::vector<int>(dim1 * dim2);
read_int(obj_id, name, &temp_arr[0], true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
result[i][j] = temp_arr[temp_idx++];
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,

View file

@ -65,6 +65,10 @@ read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,

View file

@ -34,11 +34,10 @@ Mgxs::init(const std::string& in_name, const double in_awr,
for (int thread = 0; thread < n_threads; thread++) {
cache[thread].sqrtkT = 0.;
cache[thread].t = 0;
cache[thread].p = 0;
cache[thread].a = 0;
cache[thread].uvw[0] = 1.;
cache[thread].uvw[1] = 0.;
cache[thread].uvw[2] = 0.;
cache[thread].u = 0.;
cache[thread].v = 0.;
cache[thread].w = 0.;
}
}
@ -48,7 +47,7 @@ void
Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
const int in_num_delayed_groups, double_1dvec& temperature, int& method,
const double tolerance, int_1dvec& temps_to_read, int& order_dim,
bool& is_isotropic, const int n_threads)
const int n_threads)
{
// get name
char char_name[MAX_WORD_LEN];
@ -179,7 +178,7 @@ Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
fatal_error("Invalid scatter_shape option!");
}
}
//TODO: do i even need this flag? - it should be easy to self-determine
bool in_fissionable = false;
if (attribute_exists(xs_id, "fissionable")) {
int int_fiss;
@ -264,9 +263,8 @@ Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
bool is_isotropic;
_metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
method, tolerance, temps_to_read, order_data, is_isotropic, n_threads);
method, tolerance, temps_to_read, order_data, n_threads);
// Set number of energy and delayed groups
int final_scatter_format = scatter_format;
@ -284,7 +282,7 @@ Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
final_scatter_format, order_data, max_order,
legendre_to_tabular_points, is_isotropic);
legendre_to_tabular_points, is_isotropic, n_pol, n_azi);
close_group(xsdata_grp);
} // end temperature loop
@ -419,45 +417,41 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
val = xs[cache[tid].t].total[cache[tid].a][gin];
break;
case MG_GET_XS_ABSORPTION:
val = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs[cache[tid].t].inverse_velocity[cache[tid].p][cache[tid].a][gin];
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][*dg + 1];
case MG_GET_XS_NU_FISSION:
if (fissionable) {
val = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
} else {
val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][0];
val = 0.;
}
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs[cache[tid].t].scatter[cache[tid].p]
[cache[tid].a]->get_xs(xstype, gin, gout, mu);
case MG_GET_XS_ABSORPTION:
val = xs[cache[tid].t].absorption[cache[tid].a][gin];
break;
case MG_GET_XS_FISSION:
if (fissionable) {
val = xs[cache[tid].t].fission[cache[tid].p][cache[tid].a][gin];
val = xs[cache[tid].t].fission[cache[tid].a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_KAPPA_FISSION:
if (fissionable) {
val = xs[cache[tid].t].kappa_fission[cache[tid].p][cache[tid].a][gin];
val = xs[cache[tid].t].kappa_fission[cache[tid].a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs[cache[tid].t].scatter[cache[tid].a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_PROMPT_NU_FISSION:
if (fissionable) {
val = xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
val = xs[cache[tid].t].prompt_nu_fission[cache[tid].a][gin];
} else {
val = 0.;
}
@ -465,10 +459,10 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][*dg];
val = xs[cache[tid].t].delayed_nu_fission[cache[tid].a][gin][*dg];
} else {
val = 0.;
for (auto& num : xs[cache[tid].t].delayed_nu_fission[cache[tid].p]
for (auto& num : xs[cache[tid].t].delayed_nu_fission
[cache[tid].a][gin]) {
val += num;
}
@ -477,21 +471,14 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
val = 0.;
}
break;
case MG_GET_XS_NU_FISSION:
if (fissionable) {
val = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][*gout];
val = xs[cache[tid].t].chi_prompt[cache[tid].a][gin][*gout];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin]) {
for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].a][gin]) {
val += num;
}
}
@ -503,22 +490,22 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][*dg];
val = xs[cache[tid].t].chi_delayed[cache[tid].a][gin][*gout][*dg];
} else {
val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][0];
val = xs[cache[tid].t].chi_delayed[cache[tid].a][gin][*gout][0];
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
for (int i = 0; i < xs[cache[tid].t].chi_delayed
[cache[tid].a][gin].size(); i++) {
val += xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][i][*dg];
val += xs[cache[tid].t].chi_delayed[cache[tid].a][gin][i][*dg];
}
} else {
val = 0.;
for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
for (int i = 0; i < xs[cache[tid].t].chi_delayed
[cache[tid].a][gin].size(); i++) {
for (auto& num : xs[cache[tid].t].chi_delayed[cache[tid].p]
for (auto& num : xs[cache[tid].t].chi_delayed
[cache[tid].a][gin][i]) {
val += num;
}
@ -529,6 +516,16 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
val = 0.;
}
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs[cache[tid].t].inverse_velocity[cache[tid].a][gin];
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs[cache[tid].t].decay_rate[cache[tid].a][*dg + 1];
} else {
val = xs[cache[tid].t].decay_rate[cache[tid].a][0];
}
break;
default:
val = 0.;
}
@ -541,11 +538,11 @@ void
Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
{
// This method assumes that the temperature and angle indices are set
double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
// Find the probability of having a prompt neutron
double prob_prompt =
xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
xs[cache[tid].t].prompt_nu_fission[cache[tid].a][gin];
// sample random numbers
double xi_pd = prn() * nu_fission;
@ -561,10 +558,10 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
xs[cache[tid].t].chi_prompt[cache[tid].a][gin][gout];
while (prob_gout < xi_gout) {
gout++;
prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].a][gin][gout];
}
} else {
@ -575,7 +572,7 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][dg];
xs[cache[tid].t].delayed_nu_fission[cache[tid].a][gin][dg];
}
// adjust dg in case of round-off error
@ -584,11 +581,11 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
xs[cache[tid].t].chi_delayed[cache[tid].a][gin][gout][dg];
while (prob_gout < xi_gout) {
gout++;
prob_gout +=
xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
xs[cache[tid].t].chi_delayed[cache[tid].a][gin][gout][dg];
}
}
}
@ -601,7 +598,7 @@ Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu,
{
// This method assumes that the temperature and angle indices are set
// Sample the data
xs[cache[tid].t].scatter[cache[tid].p][cache[tid].a]->sample(gin, gout, mu, wgt);
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt);
}
//==============================================================================
@ -613,11 +610,11 @@ Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT,
// Set our indices
set_temperature_index(tid, sqrtkT);
set_angle_index(tid, uvw);
total_xs = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
abs_xs = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
total_xs = xs[cache[tid].t].total[cache[tid].a][gin];
abs_xs = xs[cache[tid].t].absorption[cache[tid].a][gin];
if (fissionable) {
nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].a][gin];
} else {
nu_fiss_xs = 0.;
}
@ -670,22 +667,25 @@ void
Mgxs::set_angle_index(const int tid, const double uvw[3])
{
// See if we need to find the new index
if ((uvw[0] != cache[tid].uvw[0]) || (uvw[1] != cache[tid].uvw[1]) ||
(uvw[2] != cache[tid].uvw[2])) {
if (!is_isotropic &&
((uvw[0] != cache[tid].u) || (uvw[1] != cache[tid].v) ||
(uvw[2] != cache[tid].w))) {
// convert uvw to polar and azimuthal angles
double my_pol = std::acos(uvw[2]);
double my_azi = std::atan2(uvw[1], uvw[0]);
// Find the location, assuming equal-bin angles
double delta_angle = PI / n_pol;
cache[tid].p = std::floor(my_pol / delta_angle);
int p = std::floor(my_pol / delta_angle);
delta_angle = 2. * PI / n_azi;
cache[tid].a = std::floor((my_azi + PI) / delta_angle);
int a = std::floor((my_azi + PI) / delta_angle);
cache[tid].a = n_azi * p + a;
// store this direction as the last one used
cache[tid].uvw[0] = uvw[0];
cache[tid].uvw[1] = uvw[1];
cache[tid].uvw[2] = uvw[2];
cache[tid].u = uvw[0];
cache[tid].v = uvw[1];
cache[tid].w = uvw[2];
}
}

View file

@ -29,9 +29,11 @@ namespace openmc {
struct CacheData {
double sqrtkT; // last temperature corresponding to t
int t; // temperature index
int p; // polar angle index
int a; // azimuthal angle index
double uvw[3]; // last angle that corresponds to p and a
int a; // angle index
// last angle that corresponds to p and a
double u;
double v;
double w;
};
//==============================================================================
@ -45,6 +47,8 @@ class Mgxs {
int num_delayed_groups; // number of delayed neutron groups
int num_groups; // number of energy groups
std::vector<XsData> xs; // Cross section data
// MGXS Incoming Flux Angular grid information
bool is_isotropic; // used to skip search for angle indices if isotropic
int n_pol;
int n_azi;
double_1dvec polar;
@ -52,7 +56,7 @@ class Mgxs {
void _metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
const int in_num_delayed_groups, double_1dvec& temperature,
int& method, const double tolerance, int_1dvec& temps_to_read,
int& order_dim, bool& is_isotropic, const int n_threads);
int& order_dim, const int n_threads);
bool equiv(const Mgxs& that);
public:

View file

@ -876,9 +876,9 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
{
// See if the user wants us to figure out how many points to use
if (n_mu == C_NONE) {
// then we will use 2 pts if its P0 or P1 (super fast), or the default if
// a higher order
if (leg.get_order() <= 1) {
// then we will use 2 pts if its P0, or the default if a higher order
// TODO use an error minimization algorithm that also picks n_mu
if (leg.get_order() == 0) {
n_mu = 2;
} else {
n_mu = DEFAULT_NMU;

View file

@ -10,59 +10,49 @@ XsData::XsData(const int energy_groups, const int num_delayed_groups,
const bool fissionable, const int scatter_format,
const int n_pol, const int n_azi)
{
int n_ang = n_pol * n_azi;
// check to make sure scatter format is OK before we allocate
if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
scatter_format != ANGLE_LEGENDRE) {
fatal_error("Invalid scatter_format!");
}
// allocate all [temperature][phi][theta][in group] quantities
total = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
absorption = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
inverse_velocity = double_3dvec(n_pol,
double_2dvec(n_azi, double_1dvec(energy_groups, 0.)));
total = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
absorption = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
inverse_velocity = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
if (fissionable) {
fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
nu_fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
prompt_nu_fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
kappa_fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
prompt_nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
kappa_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.));
}
// allocate decay_rate; [temperature][phi][theta][delayed group]
decay_rate = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(num_delayed_groups, 0.)));
decay_rate = double_2dvec(n_ang, double_1dvec(num_delayed_groups, 0.));
if (fissionable) {
// allocate delayed_nu_fission; [temperature][phi][theta][in group][delay group]
delayed_nu_fission = double_4dvec(n_pol, double_3dvec(n_azi,
double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
delayed_nu_fission = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(num_delayed_groups, 0.)));
// chi_prompt; [temperature][phi][theta][in group][delayed group]
chi_prompt = double_4dvec(n_pol, double_3dvec(n_azi,
double_2dvec(energy_groups, double_1dvec(energy_groups, 0.))));
chi_prompt = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups, 0.)));
// chi_delayed; [temperature][phi][theta][in group][out group][delay group]
chi_delayed = double_5dvec(n_pol, double_4dvec(n_azi,
double_3dvec(energy_groups, double_2dvec(energy_groups,
double_1dvec(num_delayed_groups, 0.)))));
chi_delayed = double_4dvec(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
}
scatter.resize(n_pol);
for (int p = 0; p < n_pol; p++) {
scatter[p].resize(n_azi);
for (int a = 0; a < n_azi; a++) {
if (scatter_format == ANGLE_HISTOGRAM) {
scatter[p][a] = new ScattDataHistogram;
} else if (scatter_format == ANGLE_TABULAR) {
scatter[p][a] = new ScattDataTabular;
} else if (scatter_format == ANGLE_LEGENDRE) {
scatter[p][a] = new ScattDataLegendre;
}
scatter.resize(n_ang);
for (int a = 0; a < n_ang; a++) {
if (scatter_format == ANGLE_HISTOGRAM) {
scatter[a] = new ScattDataHistogram;
} else if (scatter_format == ANGLE_TABULAR) {
scatter[a] = new ScattDataTabular;
} else if (scatter_format == ANGLE_LEGENDRE) {
scatter[a] = new ScattDataLegendre;
}
}
}
@ -73,13 +63,13 @@ void
XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
const int scatter_format, const int final_scatter_format,
const int order_data, const int max_order,
const int legendre_to_tabular_points, const bool is_isotropic)
const int legendre_to_tabular_points, const bool is_isotropic,
const int n_pol, const int n_azi)
{
// Reconstruct the dimension information so it doesn't need to be passed
int n_pol = total.size();
int n_azi = total[0].size();
int energy_groups = total[0][0].size();
int delayed_groups = decay_rate[0][0].size();
int n_ang = n_pol * n_azi;
int energy_groups = total[0].size();
int delayed_groups = decay_rate[0].size();
// Set the fissionable-specific data
if (fissionable) {
@ -98,11 +88,9 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
// Check absorption to ensure it is not 0 since it is often the
// denominator in tally methods
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (absorption[p][a][gin] == 0.) absorption[p][a][gin] = 1.e-10;
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (absorption[a][gin] == 0.) absorption[a][gin] = 1.e-10;
}
}
@ -110,23 +98,17 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
if (object_exists(xsdata_grp, "total")) {
read_nd_vector(xsdata_grp, "total", total);
} else {
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
total[p][a][gin] = absorption[p][a][gin] +
scatter[p][a]->scattxs[gin];
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
total[a][gin] = absorption[a][gin] + scatter[a]->scattxs[gin];
}
}
}
// Check total to ensure it is not 0 since it is often the denominator in
// tally methods
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (total[p][a][gin] == 0.) total[p][a][gin] = 1.e-10;
}
// Fix if total is 0, since it is in the denominator when tallying
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
if (total[a][gin] == 0.) total[a][gin] = 1.e-10;
}
}
}
@ -138,15 +120,14 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
const int n_azi, const int energy_groups, const int delayed_groups,
const bool is_isotropic)
{
int n_ang = n_pol * n_azi;
// Get the fission and kappa_fission data xs; these are optional
read_nd_vector(xsdata_grp, "fission", fission);
read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);
// Set/get beta
double_4dvec temp_beta =
double_4dvec(n_pol, double_3dvec(n_azi,
double_2dvec(energy_groups, double_1dvec(delayed_groups, 0.))));
double_3dvec temp_beta =double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(delayed_groups, 0.)));
if (object_exists(xsdata_grp, "beta")) {
hid_t xsdata = open_dataset(xsdata_grp, "beta");
int ndims = dataset_ndims(xsdata);
@ -160,14 +141,12 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
// Broadcast to all incoming groups
int temp_idx = 0;
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set the first group index and copy the rest
temp_beta[p][a][0][dg] = temp_arr[temp_idx++];
for (int gin = 1; gin < energy_groups; gin++) {
temp_beta[p][a][gin] = temp_beta[p][a][0];
}
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set the first group index and copy the rest
temp_beta[a][0][dg] = temp_arr[temp_idx++];
for (int gin = 1; gin < energy_groups; gin++) {
temp_beta[a][gin] = temp_beta[a][0];
}
}
}
@ -181,41 +160,38 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
// If chi is provided, set chi-prompt and chi-delayed
if (object_exists(xsdata_grp, "chi")) {
double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups)));
double_2dvec temp_arr = double_2dvec(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi", temp_arr);
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
// First set the first group
for (int a = 0; a < n_ang; a++) {
// First set the first group
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][0][gout] = temp_arr[a][gout];
}
// Now normalize this data
double chi_sum = std::accumulate(chi_prompt[a][0].begin(),
chi_prompt[a][0].end(),
0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi for a group that is <= 0!");
}
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][0][gout] /= chi_sum;
}
// And extend to the remaining incoming groups
for (int gin = 1; gin < energy_groups; gin++) {
chi_prompt[a][gin] = chi_prompt[a][0];
}
// Finally set chi-delayed equal to chi-prompt
// Set chi-delayed to chi-prompt
for(int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][0][gout] = temp_arr[p][a][gout];
}
// Now normalize this data
double chi_sum = std::accumulate(chi_prompt[p][a][0].begin(),
chi_prompt[p][a][0].end(),
0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi for a group that is <= 0!");
}
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][0][gout] /= chi_sum;
}
// And extend to the remaining incoming groups
for (int gin = 1; gin < energy_groups; gin++) {
chi_prompt[p][a][gin] = chi_prompt[p][a][0];
}
// Finally set chi-delayed equal to chi-prompt
// Set chi-delayed to chi-prompt
for(int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[p][a][gin][gout][dg] =
chi_prompt[p][a][gin][gout];
}
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] =
chi_prompt[a][gin][gout];
}
}
}
@ -234,21 +210,18 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission);
// set delayed-nu-fission and correct prompt-nu-fission with beta
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[p][a][gin][dg] =
temp_beta[p][a][gin][dg] *
prompt_nu_fission[p][a][gin];
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] * prompt_nu_fission[a][gin];
}
// Correct the prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(),
temp_beta[p][a][gin].end(), 0.);
prompt_nu_fission[p][a][gin] *= (1. - beta_sum);
}
// Correct the prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
}
}
}
@ -258,47 +231,45 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
read_nd_vector(xsdata_grp, "nu-fission", chi_prompt);
// Normalize the chi info so the CDF is 1.
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(),
chi_prompt[p][a][gin].end(), 0.);
// Set the vector nu-fission from the matrix nu-fission
prompt_nu_fission[p][a][gin] = chi_sum;
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
// Set the vector nu-fission from the matrix nu-fission
prompt_nu_fission[a][gin] = chi_sum;
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi for a group that is <= 0!");
}
}
// set chi-delayed to chi-prompt
for (int gin = 0; gin < energy_groups; gin++) {
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[p][a][gin][gout][dg] =
chi_prompt[p][a][gin][gout];
}
chi_prompt[a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi for a group that is <= 0!");
}
}
// set chi-delayed to chi-prompt
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] =
chi_prompt[a][gin][gout];
}
}
}
// Set the delayed-nu-fission and correct prompt-nu-fission with beta
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[p][a][gin][dg] =
temp_beta[p][a][gin][dg] *
prompt_nu_fission[p][a][gin];
}
// Set the delayed-nu-fission and correct prompt-nu-fission with beta
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] *
prompt_nu_fission[a][gin];
}
// Correct prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(),
temp_beta[p][a][gin].end(), 0.);
prompt_nu_fission[p][a][gin] *= (1. - beta_sum);
}
// Correct prompt-nu-fission using the delayed neutron fraction
if (delayed_groups > 0) {
double beta_sum = std::accumulate(temp_beta[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
}
}
}
@ -311,27 +282,24 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
// If chi-prompt is provided, set chi-prompt
if (object_exists(xsdata_grp, "chi-prompt")) {
double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups)));
double_2dvec temp_arr = double_2dvec(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi-prompt", temp_arr);
for (int a = 0; a < n_azi; a++) {
for (int p = 0; p < n_pol; p++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][gin][gout] = temp_arr[p][a][gout];
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] = temp_arr[a][gout];
}
// Normalize chi so its CDF goes to 1
double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(),
chi_prompt[p][a][gin].end(), 0.);
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0.!");
// Normalize chi so its CDF goes to 1
double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] /= chi_sum;
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0.!");
}
}
}
@ -346,26 +314,22 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
if (ndims == 3) {
// chi-delayed is a [in group] vector
double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups)));
double_2dvec temp_arr = double_2dvec(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "chi-delayed", temp_arr);
for (int a = 0; a < n_azi; a++) {
for (int p = 0; p < n_pol; p++) {
// normalize the chi CDF to 1
double chi_sum = std::accumulate(temp_arr[p][a].begin(),
temp_arr[p][a].end(), 0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
for (int a = 0; a < n_ang; a++) {
// normalize the chi CDF to 1
double chi_sum = std::accumulate(temp_arr[a].begin(),
temp_arr[a].end(), 0.);
if (chi_sum <= 0.) {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
// set chi-delayed
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[p][a][gin][gout][dg] =
temp_arr[p][a][gout] / chi_sum;
}
// set chi-delayed
for (int gin = 0; gin < energy_groups; gin++) {
for (int gout = 0; gout < energy_groups; gout++) {
for (int dg = 0; dg < delayed_groups; dg++) {
chi_delayed[a][gin][gout][dg] = temp_arr[a][gout] / chi_sum;
}
}
}
@ -375,22 +339,20 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
read_nd_vector(xsdata_grp, "chi-delayed", chi_delayed);
// Normalize the chi info so the CDF is 1.
for (int a = 0; a < n_azi; a++) {
for (int p = 0; p < n_pol; p++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_delayed[p][a][gin][gout][dg];
}
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_delayed[a][gin][gout][dg];
}
if (chi_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[p][a][gin][gout][dg] /= chi_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
if (chi_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[a][gin][gout][dg] /= chi_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
}
}
@ -414,30 +376,27 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
} else if (ndims == 4) {
// prompt nu fission is a matrix,
// so set prompt_nu_fiss & chi_prompt
double_4dvec temp_arr = double_4dvec(n_pol, double_3dvec(n_azi,
double_2dvec(energy_groups, double_1dvec(energy_groups))));
double_3dvec temp_arr = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups)));
read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_arr);
// The prompt_nu_fission vector from the matrix form
for (int a = 0; a < n_azi; a++) {
for (int p = 0; p < n_pol; p++) {
for (int gin = 0; gin < energy_groups; gin++) {
double prompt_sum = std::accumulate(temp_arr[p][a][gin].begin(),
temp_arr[p][a][gin].end(), 0.);
prompt_nu_fission[p][a][gin] = prompt_sum;
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
double prompt_sum = std::accumulate(temp_arr[a][gin].begin(),
temp_arr[a][gin].end(), 0.);
prompt_nu_fission[a][gin] = prompt_sum;
}
// The chi_prompt data is just the normalized fission matrix
for (int gin= 0; gin < energy_groups; gin++) {
if (prompt_nu_fission[p][a][gin] > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[p][a][gin][gout] =
temp_arr[p][a][gin][gout] /
prompt_nu_fission[p][a][gin];
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0!");
// The chi_prompt data is just the normalized fission matrix
for (int gin= 0; gin < energy_groups; gin++) {
if (prompt_nu_fission[a][gin] > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] =
temp_arr[a][gin][gout] / prompt_nu_fission[a][gin];
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0!");
}
}
}
@ -455,23 +414,20 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
if (is_isotropic) ndims += 2;
if (ndims == 3) {
// delayed-nu-fission is a [in group] vector
if (temp_beta[0][0][0][0] == 0.) {
// delayed-nu-fission is an [in group] vector
if (temp_beta[0][0][0] == 0.) {
fatal_error("cannot set delayed-nu-fission with a 1D array if "
"beta is not provided");
}
double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups)));
double_2dvec temp_arr = double_2dvec(n_ang, double_1dvec(energy_groups));
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set delayed-nu-fission using beta
delayed_nu_fission[p][a][gin][dg] =
temp_beta[p][a][gin][dg] * temp_arr[p][a][gin];
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
for (int dg = 0; dg < delayed_groups; dg++) {
// Set delayed-nu-fission using beta
delayed_nu_fission[a][gin][dg] =
temp_beta[a][gin][dg] * temp_arr[a][gin];
}
}
}
@ -482,32 +438,28 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
} else if (ndims == 5) {
// This will contain delayed-nu-fision and chi-delayed data
double_5dvec temp_arr = double_5dvec(n_pol, double_4dvec(n_azi,
double_3dvec(energy_groups, double_2dvec(energy_groups,
double_1dvec(delayed_groups)))));
double_4dvec temp_arr = double_4dvec(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(delayed_groups))));
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
// Set the 4D delayed-nu-fission matrix and 5D chi-delayed matrix
// from the 5D delayed-nu-fission matrix
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double gout_sum = 0.;
// Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix
// from the 4D delayed-nu-fission matrix
for (int a = 0; a < n_ang; a++) {
for (int dg = 0; dg < delayed_groups; dg++) {
for (int gin = 0; gin < energy_groups; gin++) {
double gout_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
gout_sum += temp_arr[a][gin][gout][dg];
chi_delayed[a][gin][gout][dg] = temp_arr[a][gin][gout][dg];
}
delayed_nu_fission[a][gin][dg] = gout_sum;
// Normalize chi-delayed
if (gout_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
gout_sum += temp_arr[p][a][gin][gout][dg];
chi_delayed[p][a][gin][gout][dg] =
temp_arr[p][a][gin][gout][dg];
}
delayed_nu_fission[p][a][gin][dg] = gout_sum;
// Normalize chi-delayed
if (gout_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[p][a][gin][gout][dg] /= gout_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
chi_delayed[a][gin][gout][dg] /= gout_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
}
}
@ -520,14 +472,12 @@ XsData::_fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
}
// Combine prompt_nu_fission and delayed_nu_fission into nu_fission
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
nu_fission[p][a][gin] =
std::accumulate(delayed_nu_fission[p][a][gin].begin(),
delayed_nu_fission[p][a][gin].end(),
prompt_nu_fission[p][a][gin]);
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
nu_fission[a][gin] =
std::accumulate(delayed_nu_fission[a][gin].begin(),
delayed_nu_fission[a][gin].end(),
prompt_nu_fission[a][gin]);
}
}
}
@ -540,37 +490,32 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
const int final_scatter_format, const int order_data, const int max_order,
const int legendre_to_tabular_points)
{
int n_ang = n_pol * n_azi;
if (!object_exists(xsdata_grp, "scatter_data")) {
fatal_error("Must provide scatter_data group!");
}
hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");
// Get the outgoing group boundary indices
int_3dvec gmin = int_3dvec(n_pol, int_2dvec(n_azi,
int_1dvec(energy_groups)));
int_2dvec gmin = int_2dvec(n_ang, int_1dvec(energy_groups));
read_nd_vector(scatt_grp, "g_min", gmin, true);
int_3dvec gmax = int_3dvec(n_pol, int_2dvec(n_azi,
int_1dvec(energy_groups)));
int_2dvec gmax = int_2dvec(n_ang, int_1dvec(energy_groups));
read_nd_vector(scatt_grp, "g_max", gmax, true);
// Make gmin and gmax start from 0 vice 1 as they do in the library
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
gmin[p][a][gin] -= 1;
gmax[p][a][gin] -= 1;
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
gmin[a][gin] -= 1;
gmax[a][gin] -= 1;
}
}
// Now use this info to find the length of a vector to hold the flattened
// data.
int length = 0;
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
length += order_data * (gmax[p][a][gin] - gmin[p][a][gin] + 1);
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
length += order_data * (gmax[a][gin] - gmin[a][gin] + 1);
}
}
double_1dvec temp_arr = double_1dvec(length);
@ -587,55 +532,48 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
// convert the flattened temp_arr to a jagged array for passing to
// scatt data
double_5dvec input_scatt =
double_5dvec(n_pol, double_4dvec(n_azi, double_3dvec(energy_groups)));
double_4dvec input_scatt =
double_4dvec(n_ang, double_3dvec(energy_groups));
int temp_idx = 0;
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
input_scatt[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
for (int i_gout = 0; i_gout < input_scatt[p][a][gin].size(); i_gout++) {
input_scatt[p][a][gin][i_gout].resize(order_dim);
for (int l = 0; l < order_dim; l++) {
input_scatt[p][a][gin][i_gout][l] = temp_arr[temp_idx++];
}
// Adjust index for the orders we didnt take
temp_idx += (order_data - order_dim);
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
input_scatt[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
input_scatt[a][gin][i_gout].resize(order_dim);
for (int l = 0; l < order_dim; l++) {
input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++];
}
// Adjust index for the orders we didnt take
temp_idx += (order_data - order_dim);
}
}
}
temp_arr.clear();
// Get multiplication matrix
double_4dvec temp_mult = double_4dvec(n_pol, double_3dvec(n_azi,
double_2dvec(energy_groups)));
double_3dvec temp_mult = double_3dvec(n_ang, double_2dvec(energy_groups));
if (object_exists(scatt_grp, "multiplicity_matrix")) {
temp_arr.resize(length / order_data);
read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);
// convert the flat temp_arr to a jagged array for passing to scatt data
int temp_idx = 0;
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) {
temp_mult[p][a][gin][i_gout] = temp_arr[temp_idx++];
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
}
}
}
} else {
// Use a default: multiplicities are 1.0.
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) {
temp_mult[p][a][gin][i_gout] = 1.;
}
for (int a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1);
for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
temp_mult[a][gin][i_gout] = 1.;
}
}
}
@ -646,28 +584,22 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
// Finally, convert the Legendre data to tabular, if needed
if (scatter_format == ANGLE_LEGENDRE &&
final_scatter_format == ANGLE_TABULAR) {
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
ScattDataLegendre legendre_scatt;
legendre_scatt.init(gmin[p][a], gmax[p][a], temp_mult[p][a],
input_scatt[p][a]);
for (int a = 0; a < n_ang; a++) {
ScattDataLegendre legendre_scatt;
legendre_scatt.init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
// Now create a tabular version of legendre_scatt
convert_legendre_to_tabular(legendre_scatt,
*static_cast<ScattDataTabular*>(scatter[p][a]),
legendre_to_tabular_points);
// Now create a tabular version of legendre_scatt
convert_legendre_to_tabular(legendre_scatt,
*static_cast<ScattDataTabular*>(scatter[a]),
legendre_to_tabular_points);
scatter_format = final_scatter_format;
}
scatter_format = final_scatter_format;
}
} else {
// We are sticking with the current representation
// Initialize the ScattData object with this data
for (int p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
scatter[p][a]->init(gmin[p][a], gmax[p][a], temp_mult[p][a],
input_scatt[p][a]);
}
for (int a = 0; a < n_ang; a++) {
scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
}
}
}
@ -675,7 +607,7 @@ XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
//==============================================================================
void
XsData::combine(const std::vector<XsData*> those_xs,
XsData::combine(const std::vector<XsData*>& those_xs,
const double_1dvec& scalars)
{
// Combine the non-scattering data
@ -683,64 +615,60 @@ XsData::combine(const std::vector<XsData*> those_xs,
XsData* that = those_xs[i];
if (!equiv(*that)) fatal_error("Cannot combine the XsData objects!");
double scalar = scalars[i];
for (int p = 0; p < total.size(); p++) {
for (int a = 0; a < total[p].size(); a++) {
for (int gin = 0; gin < total[p][a].size(); gin++) {
total[p][a][gin] += scalar * that->total[p][a][gin];
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) {
nu_fission[p][a][gin] +=
scalar * that->nu_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 a = 0; a < total.size(); a++) {
for (int gin = 0; gin < total[a].size(); gin++) {
total[a][gin] += scalar * that->total[a][gin];
absorption[a][gin] += scalar * that->absorption[a][gin];
if (i == 0) {
inverse_velocity[a][gin] = that->inverse_velocity[a][gin];
}
if (that->prompt_nu_fission.size() > 0) {
nu_fission[a][gin] += scalar * that->nu_fission[a][gin];
prompt_nu_fission[a][gin] +=
scalar * that->prompt_nu_fission[a][gin];
kappa_fission[a][gin] += scalar * that->kappa_fission[a][gin];
fission[a][gin] += scalar * that->fission[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[a][gin].size(); dg++) {
delayed_nu_fission[a][gin][dg] +=
scalar * that->delayed_nu_fission[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[a][gin].size(); gout++) {
chi_prompt[a][gin][gout] +=
scalar * that->chi_prompt[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[a][gin][gout].size(); dg++) {
chi_delayed[a][gin][gout][dg] +=
scalar * that->chi_delayed[a][gin][gout][dg];
}
}
}
}
for (int dg = 0; dg < decay_rate[p][a].size(); dg++) {
decay_rate[p][a][dg] += scalar * that->decay_rate[p][a][dg];
}
for (int dg = 0; dg < decay_rate[a].size(); dg++) {
decay_rate[a][dg] += scalar * that->decay_rate[a][dg];
}
// Normalize chi
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;
}
// Normalize chi
if (chi_prompt.size() > 0) {
for (int gin = 0; gin < chi_prompt[a].size(); gin++) {
double norm = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
if (norm > 0.) {
for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) {
chi_prompt[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 (norm > 0.) {
for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) {
chi_delayed[p][a][gin][gout][dg] /= norm;
}
for (int dg = 0; dg < chi_delayed[a][gin][0].size(); dg++) {
norm = 0.;
for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
norm += chi_delayed[a][gin][gout][dg];
}
if (norm > 0.) {
for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) {
chi_delayed[a][gin][gout][dg] /= norm;
}
}
}
@ -750,17 +678,15 @@ XsData::combine(const std::vector<XsData*> those_xs,
}
// Allow the ScattData object to combine itself
for (int p = 0; p < total.size(); p++) {
for (int a = 0; a < total[p].size(); a++) {
// Build vector of the scattering objects to incorporate
std::vector<ScattData*> those_scatts(those_xs.size());
for (int i = 0; i < those_xs.size(); i++) {
those_scatts[i] = those_xs[i]->scatter[p][a];
}
// Now combine these guys
scatter[p][a]->combine(those_scatts, scalars);
for (int a = 0; a < total.size(); a++) {
// Build vector of the scattering objects to incorporate
std::vector<ScattData*> those_scatts(those_xs.size());
for (int i = 0; i < those_xs.size(); i++) {
those_scatts[i] = those_xs[i]->scatter[a];
}
// Now combine these guys
scatter[a]->combine(those_scatts, scalars);
}
}
@ -770,12 +696,8 @@ bool
XsData::equiv(const XsData& that)
{
bool match = false;
// check n_pol (total.size()), n_azi (total[0].size()), and
// groups (total[0][0].size())
// This assumes correct initializatino of the remaining cross sections
if ((total.size() == that.total.size()) &&
(total[0].size() == that.total[0].size()) &&
(total[0][0].size() == that.total[0][0].size())) {
if ((absorption.size() == that.absorption.size()) &&
(absorption[0].size() == that.absorption[0].size())) {
match = true;
}
return match;

View file

@ -28,32 +28,34 @@ class XsData {
const int n_azi, const int energy_groups, int scatter_format,
const int final_scatter_format, const int order_data,
const int max_order, const int legendre_to_tabular_points);
void _fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int delayed_groups, bool is_isotropic);
void _fissionable_from_hdf5(const hid_t xsdata_grp, const int n_pol,
const int n_azi, const int energy_groups, const int delayed_groups,
const bool is_isotropic);
public:
// The following quantities have the following dimensions:
// [phi][theta][incoming group]
double_3dvec total;
double_3dvec absorption;
double_3dvec nu_fission;
double_3dvec prompt_nu_fission;
double_3dvec kappa_fission;
double_3dvec fission;
double_3dvec inverse_velocity;
// [angle][incoming group]
double_2dvec total;
double_2dvec absorption;
double_2dvec nu_fission;
double_2dvec prompt_nu_fission;
double_2dvec kappa_fission;
double_2dvec fission;
double_2dvec inverse_velocity;
// decay_rate has the following dimensions:
// [phi][theta][delayed group]
double_3dvec decay_rate;
// [angle][delayed group]
double_2dvec decay_rate;
// delayed_nu_fission has the following dimensions:
// [phi][theta][incoming group][delayed group]
double_4dvec delayed_nu_fission;
// [angle][incoming group][delayed group]
double_3dvec delayed_nu_fission;
// chi_prompt has the following dimensions:
// [phi][theta][incoming group][outgoing group]
double_4dvec chi_prompt;
// [angle][incoming group][outgoing group]
double_3dvec chi_prompt;
// chi_delayed has the following dimensions:
// [phi][theta][incoming group][outgoing group][delayed group]
double_5dvec chi_delayed;
// scatter has the following dimensions: [phi][theta]
std::vector<std::vector<ScattData*> > scatter;
// [angle][incoming group][outgoing group][delayed group]
double_4dvec chi_delayed;
// scatter has the following dimensions: [angle]
std::vector<ScattData*> scatter;
XsData() = default;
XsData(const int num_groups, const int num_delayed_groups,
@ -63,8 +65,8 @@ class XsData {
const int scatter_format, const int final_scatter_format,
const int order_data, const int max_order,
const int legendre_to_tabular_points,
const bool is_isotropic);
void combine(const std::vector<XsData*> those_xs,
const bool is_isotropic, const int n_pol, const int n_azi);
void combine(const std::vector<XsData*>& those_xs,
const double_1dvec& scalars);
bool equiv(const XsData& that);
};