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

Saving state

Browse source
This commit is contained in:
Adam G Nelson 2018-08-31 18:18:01 -04:00
parent 6a739abe52
commit 4e92988433
10 changed files with 536 additions and 467 deletions
Download patch file Download diff file Expand all files Collapse all files

5
include/openmc/constants.h
View file

@ -16,11 +16,6 @@ typedef std::vector<double> double_1dvec;
typedef std::vector<std::vector<double> > double_2dvec;
typedef std::vector<std::vector<std::vector<double> > > double_3dvec;
typedef std::vector<std::vector<std::vector<std::vector<double> > > > double_4dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > double_5dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > > double_6dvec;
typedef std::vector<int> int_1dvec;
typedef std::vector<std::vector<int> > int_2dvec;
typedef std::vector<std::vector<std::vector<int> > > int_3dvec;
// ============================================================================
// VERSIONING NUMBERS

27
include/openmc/hdf5_interface.h
View file

@ -46,37 +46,36 @@ hid_t file_open(const std::string& filename, char mode, bool parallel=false);
void write_string(hid_t group_id, const char* name, const std::string& buffer,
bool indep);
void
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 1>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 2>& 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,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 2>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<int> > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 3>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<double> > > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 3>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 4>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 5>& result,
bool must_have = false);
std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);

26
include/openmc/mgxs.h
View file

@ -7,6 +7,8 @@
#include <string>
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xsdata.h"
@ -35,7 +37,7 @@ struct CacheData {
class Mgxs {
private:
double_1dvec kTs; // temperature in eV (k * T)
xt::xtensor<double, 1> kTs; // temperature in eV (k * T)
int scatter_format; // flag for if this is legendre, histogram, or tabular
int num_delayed_groups; // number of delayed neutron groups
int num_groups; // number of energy groups
@ -44,8 +46,8 @@ class Mgxs {
bool is_isotropic; // used to skip search for angle indices if isotropic
int n_pol;
int n_azi;
double_1dvec polar;
double_1dvec azimuthal;
std::vector<double> polar;
std::vector<double> azimuthal;
//! \brief Initializes the Mgxs object metadata
//!
@ -62,10 +64,10 @@ class Mgxs {
//! @param in_polar Polar angle grid.
//! @param in_azimuthal Azimuthal angle grid.
void
init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs,
init(const std::string& in_name, double in_awr, const std::vector<double>& in_kTs,
bool in_fissionable, int in_scatter_format, int in_num_groups,
int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal);
const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal);
//! \brief Initializes the Mgxs object metadata from the HDF5 file
//!
@ -80,8 +82,8 @@ class Mgxs {
//! @param method Method of choosing nearest temperatures.
void
metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim,
int in_num_delayed_groups, const std::vector<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim,
int& method);
//! \brief Performs the actual act of combining the microscopic data for a
@ -93,8 +95,8 @@ class Mgxs {
//! corresponds to the temperature of interest.
//! @param this_t The temperature index of the macroscopic object.
void
combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t);
combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scalars,
const std::vector<int>& micro_ts, int this_t);
//! \brief Checks to see if this and that are able to be combined
//!
@ -128,7 +130,7 @@ class Mgxs {
//! provides the number of points to use in the tabular representation.
//! @param method Method of choosing nearest temperatures.
Mgxs(hid_t xs_id, int energy_groups,
int delayed_groups, const double_1dvec& temperature, double tolerance,
int delayed_groups, const std::vector<double>& temperature, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points, int& method);
@ -141,8 +143,8 @@ class Mgxs {
//! @param atom_densities Atom densities of those microscopic quantities.
//! @param tolerance Tolerance of temperature selection method.
//! @param method Method of choosing nearest temperatures.
Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method);
//! \brief Provides a cross section value given certain parameters

52
include/openmc/scattdata.h
View file

@ -6,6 +6,8 @@
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/constants.h"
namespace openmc {
@ -25,23 +27,25 @@ class ScattData {
protected:
//! \brief Initializes the attributes of the base class.
void
base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_energy, const double_2dvec& in_mult);
base_init(int order, const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult);
//! \brief Combines microscopic ScattDatas into a macroscopic one.
void
base_combine(int max_order, const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& sparse_mult, double_3dvec& sparse_scatter);
const double_1dvec& scalars, xt::xtensor<int, 1>& in_gmin,
xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter);
public:
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
int_1dvec gmin; // minimum outgoing group
int_1dvec gmax; // maximum outgoing group
double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
xt::xtensor<double, 1> gmin; // minimum outgoing group
xt::xtensor<double, 1> gmax; // maximum outgoing group
xt::xtensor<double, 1> scattxs; // Isotropic Sigma_{s,g_{in}}
//! \brief Calculates the value of normalized f(mu).
//!
@ -72,7 +76,7 @@ class ScattData {
//! @param in_mult Input sparse multiplicity matrix
//! @param coeffs Input sparse scattering matrix
virtual void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs) = 0;
//! \brief Combines the microscopic data.
@ -97,7 +101,7 @@ class ScattData {
//! @param max_order If Legendre this is the maximum value of "n" in "Pn"
//! requested; ignored otherwise.
//! @return The dense scattering matrix.
virtual double_3dvec
virtual xt::xtensor<double, 3>
get_matrix(int max_order) = 0;
//! \brief Samples the outgoing energy from the ScattData info.
@ -142,7 +146,7 @@ class ScattDataLegendre: public ScattData {
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
@ -163,7 +167,7 @@ class ScattDataLegendre: public ScattData {
int
get_order() {return dist[0][0].size() - 1;};
double_3dvec
xt::xtensor<double, 3>
get_matrix(int max_order);
};
@ -176,14 +180,14 @@ class ScattDataHistogram: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the spacing between the mu bin points
double_3dvec fmu; // The angular distribution histogram
xt::xtensor<double, 1> mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the mu spacing
double_3dvec fmu; // The angular distribution histogram
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
@ -199,7 +203,7 @@ class ScattDataHistogram: public ScattData {
int
get_order() {return dist[0][0].size();};
double_3dvec
xt::xtensor<double, 3>
get_matrix(int max_order);
};
@ -212,9 +216,9 @@ class ScattDataTabular: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu grid points
double dmu; // Quick storage of the spacing between the mu points
double_3dvec fmu; // The angular distribution function
xt::xtensor<double, 1> mu; // Angle distribution mu grid points
double dmu; // Quick storage of the mu spacing
double_3dvec fmu; // The angular distribution function
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
@ -225,7 +229,7 @@ class ScattDataTabular: public ScattData {
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
@ -241,7 +245,7 @@ class ScattDataTabular: public ScattData {
int
get_order() {return dist[0][0].size();};
double_3dvec get_matrix(int max_order);
xt::xtensor<double, 3> get_matrix(int max_order);
};
//==============================================================================

24
include/openmc/xsdata.h
View file

@ -7,6 +7,8 @@
#include <memory>
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/hdf5_interface.h"
#include "openmc/scattdata.h"
@ -35,26 +37,26 @@ class XsData {
// The following quantities have the following dimensions:
// [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;
xt::xtensor<double, 2> total;
xt::xtensor<double, 2> absorption;
xt::xtensor<double, 2> nu_fission;
xt::xtensor<double, 2> prompt_nu_fission;
xt::xtensor<double, 2> kappa_fission;
xt::xtensor<double, 2> fission;
xt::xtensor<double, 2> inverse_velocity;
// decay_rate has the following dimensions:
// [angle][delayed group]
double_2dvec decay_rate;
xt::xtensor<double, 2> decay_rate;
// delayed_nu_fission has the following dimensions:
// [angle][incoming group][delayed group]
double_3dvec delayed_nu_fission;
xt::xtensor<double, 3> delayed_nu_fission;
// chi_prompt has the following dimensions:
// [angle][incoming group][outgoing group]
double_3dvec chi_prompt;
xt::xtensor<double, 3> chi_prompt;
// chi_delayed has the following dimensions:
// [angle][incoming group][outgoing group][delayed group]
double_4dvec chi_delayed;
xt::xtensor<double, 4> chi_delayed;
// scatter has the following dimensions: [angle]
std::vector<std::shared_ptr<ScattData> > scatter;

192
src/hdf5_interface.cpp
View file

@ -5,6 +5,9 @@
#include <sstream>
#include <string>
#include "xtensor/xtensor.hpp"
#include "xtensor/xarray.hpp"
#include "hdf5.h"
#include "hdf5_hl.h"
#ifdef OPENMC_MPI
@ -474,6 +477,106 @@ read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
if (name) H5Dclose(dset);
}
// //*****************************************************************************
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 1>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// read_double(obj_id, name, result.data(), true);
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 2>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<double> temp;
// read_double(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 3>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<double> temp;
// read_double(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 4>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<double> temp;
// read_double(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 5>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<double> temp;
// read_double(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 2>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<int> temp;
// read_int(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
// void
// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 3>& result,
// bool must_have)
// {
// if (object_exists(obj_id, name)) {
// xt::xarray<int> temp;
// read_int(obj_id, name, temp.data(), true);
// result = xt::adapt(temp, result.shape());
// } else if (must_have) {
// fatal_error(std::string("Must provide " + std::string(name) + "!"));
// }
// }
//*****************************************************************************
void
read_double(hid_t obj_id, const char* name, double* buffer, bool indep)
@ -545,19 +648,38 @@ read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result, bool must_have)
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 1>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim1 = result.shape()[0];
double temp_arr[dim1];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
result(i) = 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, xt::xtensor<double, 2>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
double temp_arr[dim1 * dim2];
read_double(obj_id, name, temp_arr, 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++];
result(i, j) = temp_arr[temp_idx++];
}
}
} else if (must_have) {
@ -567,19 +689,19 @@ 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)
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 2>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
int temp_arr[dim1 * dim2];
read_int(obj_id, name, temp_arr, 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++];
result(i, j) = temp_arr[temp_idx++];
}
}
} else if (must_have) {
@ -589,14 +711,13 @@ 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<std::vector<double> > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 3>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
int dim3 = result.shape()[2];
double temp_arr[dim1 * dim2 * dim3];
read_double(obj_id, name, temp_arr, true);
@ -604,7 +725,7 @@ read_nd_vector(hid_t obj_id, const char* name,
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
result(i, j, k) = temp_arr[temp_idx++];
}
}
}
@ -614,14 +735,13 @@ 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<std::vector<int> > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<int, 3>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
int dim3 = result.shape()[2];
int temp_arr[dim1 * dim2 * dim3];
read_int(obj_id, name, temp_arr, true);
@ -629,7 +749,7 @@ read_nd_vector(hid_t obj_id, const char* name,
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
result(i, j, k) = temp_arr[temp_idx++];
}
}
}
@ -639,15 +759,14 @@ 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<std::vector<std::vector<double> > > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 4>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
int dim3 = result.shape()[2];
int dim4 = result.shape()[3];
double temp_arr[dim1 * dim2 * dim3 * dim4];
read_double(obj_id, name, temp_arr, true);
@ -656,7 +775,7 @@ read_nd_vector(hid_t obj_id, const char* name,
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
result[i][j][k][l] = temp_arr[temp_idx++];
result(i, j, k, l) = temp_arr[temp_idx++];
}
}
}
@ -667,16 +786,15 @@ 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<std::vector<std::vector<std::vector<double> > > > >& result,
read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<double, 5>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
int dim5 = result[0][0][0][0].size();
int dim1 = result.shape()[0];
int dim2 = result.shape()[1];
int dim3 = result.shape()[2];
int dim4 = result.shape()[3];
int dim5 = result.shape()[4];
double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5];
read_double(obj_id, name, temp_arr, true);
@ -686,7 +804,7 @@ read_nd_vector(hid_t obj_id, const char* name,
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
for (int m = 0; m < dim5; m++) {
result[i][j][k][l][m] = temp_arr[temp_idx++];
result(i, j, k, l, m) = temp_arr[temp_idx++];
}
}
}

172
src/mgxs.cpp
View file

@ -3,12 +3,17 @@
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <valarray>
#include <sstream>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xmath.hpp"
#include "xtensor/xsort.hpp"
#include "xtensor/xadapt.hpp"
#include "xtensor/xview.hpp"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
@ -28,14 +33,15 @@ std::vector<Mgxs> macro_xs;
void
Mgxs::init(const std::string& in_name, double in_awr,
const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format,
const std::vector<double>& in_kTs, bool in_fissionable, int in_scatter_format,
int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal)
const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal)
{
// Set the metadata
name = in_name;
awr = in_awr;
kTs = in_kTs;
//TODO: Remove adapt when in_KTs is an xtensor
kTs = xt::adapt(in_kTs);
fissionable = in_fissionable;
scatter_format = in_scatter_format;
num_groups = in_num_groups;
@ -61,8 +67,8 @@ Mgxs::init(const std::string& in_name, double in_awr,
void
Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method)
int in_num_delayed_groups, const std::vector<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim, int& method)
{
// get name
char char_name[MAX_WORD_LEN];
@ -87,7 +93,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
dset_names[i] = new char[151];
}
get_datasets(kT_group, dset_names);
double_1dvec available_temps(num_temps);
xt::xarray<double> available_temps(num_temps);
for (int i = 0; i < num_temps; i++) {
read_double(kT_group, dset_names[i], &available_temps[i], true);
@ -110,24 +116,20 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
switch(method) {
case TEMPERATURE_NEAREST:
// Find the minimum difference
for (int i = 0; i < temperature.size(); i++) {
std::valarray<double> temp_diff(available_temps.data(),
available_temps.size());
temp_diff = std::abs(temp_diff - temperature[i]);
int i_closest = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// Determine actual temperatures to read
for (const auto& T : temperature) {
auto i_closest = xt::argmin(xt::abs(available_temps - T))[0];
double temp_actual = available_temps[i_closest];
if (std::abs(temp_actual - temperature[i]) < tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(temp_actual)) == temps_to_read.end()) {
if (std::fabs(temp_actual - T) < tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual))
== temps_to_read.end()) {
temps_to_read.push_back(std::round(temp_actual));
} else {
fatal_error("MGXS Library does not contain cross section for " +
in_name + " at or near " +
std::to_string(std::round(temperature[i])) + " K.");
}
} else {
std::stringstream msg;
msg << "MGXS library does not contain cross sections for "
<< in_name << " at or near " << std::round(T) << " K.";
fatal_error(msg);
}
}
break;
@ -160,7 +162,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
// Get the library's temperatures
int n_temperature = temps_to_read.size();
double_1dvec in_kTs(n_temperature);
std::vector<double> in_kTs(n_temperature);
for (int i = 0; i < n_temperature; i++) {
std::string temp_str(std::to_string(temps_to_read[i]) + "K");
@ -254,12 +256,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
}
// Set the angular bins to use equally-spaced bins
double_1dvec in_polar(in_n_pol);
std::vector<double> in_polar(in_n_pol);
double dangle = PI / in_n_pol;
for (int p = 0; p < in_n_pol; p++) {
in_polar[p] = (p + 0.5) * dangle;
}
double_1dvec in_azimuthal(in_n_azi);
std::vector<double> in_azimuthal(in_n_azi);
dangle = 2. * PI / in_n_azi;
for (int a = 0; a < in_n_azi; a++) {
in_azimuthal[a] = (a + 0.5) * dangle - PI;
@ -274,12 +276,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
//==============================================================================
Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
const double_1dvec& temperature, double tolerance, int max_order,
const std::vector<double>& temperature, double tolerance, int max_order,
bool legendre_to_tabular, int legendre_to_tabular_points, int& method)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
std::vector<int> temps_to_read;
metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
tolerance, temps_to_read, order_data, method);
@ -310,8 +312,8 @@ Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
//==============================================================================
Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method)
{
// Get the minimum data needed to initialize:
@ -328,8 +330,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
int in_num_groups = micros[0]->num_groups;
int in_num_delayed_groups = micros[0]->num_delayed_groups;
bool in_is_isotropic = micros[0]->is_isotropic;
double_1dvec in_polar = micros[0]->polar;
double_1dvec in_azimuthal = micros[0]->azimuthal;
std::vector<double> in_polar = micros[0]->polar;
std::vector<double> 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_is_isotropic, in_polar,
@ -345,33 +347,27 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
// Create the list of temperature indices and interpolation factors for
// each microscopic data at the material temperature
int_1dvec micro_t(micros.size(), 0);
double_1dvec micro_t_interp(micros.size(), 0.);
std::vector<int> micro_t(micros.size(), 0);
std::vector<double> micro_t_interp(micros.size(), 0.);
for (int m = 0; m < micros.size(); m++) {
switch(method) {
case TEMPERATURE_NEAREST:
{
// Find the nearest temperature
std::valarray<double> 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]];
micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0];
auto 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 " +
name + " at or near " +
std::to_string(std::round(temp_desired / K_BOLTZMANN))
+ " K.");
std::stringstream msg;
msg << "MGXS Library does not contain cross section for " << name
<< " at or near " << std::round(temp_desired / K_BOLTZMANN) << "K.";
fatal_error(msg);
}
}
break;
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++) {
for (int k = 0; k < micros[m]->kTs.shape()[0] - 1; k++) {
if ((micros[m]->kTs[k] <= temp_desired) &&
(temp_desired < micros[m]->kTs[k + 1])) {
micro_t[m] = k;
@ -394,8 +390,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
int num_interp_points = 2;
if (method == TEMPERATURE_NEAREST) num_interp_points = 1;
for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
double_1dvec interp(micros.size());
double_1dvec temp_indices(micros.size());
std::vector<double> interp(micros.size());
std::vector<double> temp_indices(micros.size());
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;
@ -409,8 +405,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
//==============================================================================
void
Mgxs::combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t)
Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scalars,
const std::vector<int>& micro_ts, int this_t)
{
// Build the vector of pointers to the xs objects within micros
std::vector<XsData*> those_xs(micros.size());
@ -441,19 +437,19 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs_t->total[a][gin];
val = xs_t->total(a, gin);
break;
case MG_GET_XS_NU_FISSION:
val = fissionable ? xs_t->nu_fission[a][gin] : 0.;
val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
break;
case MG_GET_XS_ABSORPTION:
val = xs_t->absorption[a][gin];
val = xs_t->absorption(a, gin);;
break;
case MG_GET_XS_FISSION:
val = fissionable ? xs_t->fission[a][gin] : 0.;
val = fissionable ? xs_t->fission(a, gin) : 0.;
break;
case MG_GET_XS_KAPPA_FISSION:
val = fissionable ? xs_t->kappa_fission[a][gin] : 0.;
val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
@ -462,16 +458,16 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_PROMPT_NU_FISSION:
val = fissionable ? xs_t->prompt_nu_fission[a][gin] : 0.;
val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
break;
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs_t->delayed_nu_fission[a][gin][*dg];
val = xs_t->delayed_nu_fission(a, gin, *dg);
} else {
val = 0.;
for (auto& num : xs_t->delayed_nu_fission[a][gin]) {
val += num;
for (int d = 0; d < xs_t->delayed_nu_fission.shape()[2]; d++) {
val += xs_t->delayed_nu_fission(a, gin, d);
}
}
} else {
@ -481,12 +477,12 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs_t->chi_prompt[a][gin][*gout];
val = xs_t->chi_prompt(a, gin, *gout);
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs_t->chi_prompt[a][gin]) {
val += num;
for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) {
val += xs_t->chi_prompt(a, gin, g);
}
}
} else {
@ -497,21 +493,21 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs_t->chi_delayed[a][gin][*gout][*dg];
val = xs_t->chi_delayed(a, gin, *gout, *dg);
} else {
val = xs_t->chi_delayed[a][gin][*gout][0];
val = xs_t->chi_delayed(a, gin, *gout, 0);
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
val += xs_t->chi_delayed[a][gin][i][*dg];
for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
val += xs_t->delayed_nu_fission(a, gin, g, *dg);
}
} else {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
for (auto& num : xs_t->chi_delayed[a][gin][i]) {
val += num;
for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) {
val += xs_t->delayed_nu_fission(a, gin, g, d);
}
}
}
@ -521,13 +517,13 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
}
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs_t->inverse_velocity[a][gin];
val = xs_t->inverse_velocity(a, gin);
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs_t->decay_rate[a][*dg + 1];
val = xs_t->decay_rate(a, *dg + 1);
} else {
val = xs_t->decay_rate[a][0];
val = xs_t->decay_rate(a, 0);
}
break;
default:
@ -548,11 +544,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
int tid = 0;
#endif
XsData* xs_t = &xs[cache[tid].t];
double nu_fission = xs_t->nu_fission[cache[tid].a][gin];
double nu_fission = xs_t->nu_fission(cache[tid].a, gin);
// Find the probability of having a prompt neutron
double prob_prompt =
xs_t->prompt_nu_fission[cache[tid].a][gin];
xs_t->prompt_nu_fission(cache[tid].a, gin);
// sample random numbers
double xi_pd = prn() * nu_fission;
@ -568,10 +564,10 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_prompt[cache[tid].a][gin][gout];
xs_t->chi_prompt(cache[tid].a, gin, gout);
while (prob_gout < xi_gout) {
gout++;
prob_gout += xs_t->chi_prompt[cache[tid].a][gin][gout];
prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout);
}
} else {
@ -582,7 +578,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs_t->delayed_nu_fission[cache[tid].a][gin][dg];
xs_t->delayed_nu_fission(cache[tid].a, gin, dg);
}
// adjust dg in case of round-off error
@ -591,11 +587,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
xs_t->chi_delayed(cache[tid].a, gin, gout, dg);
while (prob_gout < xi_gout) {
gout++;
prob_gout +=
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
xs_t->chi_delayed(cache[tid].a, gin, gout, dg);
}
}
}
@ -630,10 +626,10 @@ Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3],
set_temperature_index(sqrtkT);
set_angle_index(uvw);
XsData* xs_t = &xs[cache[tid].t];
total_xs = xs_t->total[cache[tid].a][gin];
abs_xs = xs_t->absorption[cache[tid].a][gin];
total_xs = xs_t->total(cache[tid].a, gin);
abs_xs = xs_t->absorption(cache[tid].a, gin);
nu_fiss_xs = fissionable ? xs_t->nu_fission[cache[tid].a][gin] : 0.;
nu_fiss_xs = fissionable ? xs_t->nu_fission(cache[tid].a, gin) : 0.;
}
//==============================================================================
@ -662,17 +658,7 @@ Mgxs::set_temperature_index(double sqrtkT)
int tid = 0;
#endif
if (sqrtkT != cache[tid].sqrtkT) {
double kT = sqrtkT * sqrtkT;
// initialize vector for storage of the differences
std::valarray<double> temp_diff(kTs.data(), kTs.size());
// Find the minimum difference of kT and kTs
temp_diff = std::abs(temp_diff - kT);
cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// store this temperature as the last one used
cache[tid].t = xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
cache[tid].sqrtkT = sqrtkT;
}
}

137
src/scattdata.cpp
View file

@ -4,6 +4,8 @@
#include <numeric>
#include <cmath>
#include "xtensor/xbuilder.hpp"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -16,8 +18,8 @@ namespace openmc {
//==============================================================================
void
ScattData::base_init(int order, const int_1dvec& in_gmin,
const int_1dvec& in_gmax, const double_2dvec& in_energy,
ScattData::base_init(int order, const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult)
{
int groups = in_energy.size();
@ -53,16 +55,15 @@ ScattData::base_init(int order, const int_1dvec& in_gmin,
void
ScattData::base_combine(int max_order,
const std::vector<ScattData*>& those_scatts, const double_1dvec& scalars,
int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult,
xt::xtensor<int, 1>& in_gmin, xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter)
{
int groups = those_scatts[0] -> energy.size();
// Now allocate and zero our storage spaces
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.));
xt::xtensor<double, 3> this_matrix({groups, groups, max_order}, 0.);
xt::xtensor<double, 2> mult_numer({groups, groups}, 0.);
xt::xtensor<double, 2> mult_denom({groups, groups}, 0.);
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
@ -80,26 +81,26 @@ ScattData::base_combine(int max_order,
ScattData* that = those_scatts[i];
// Build the dense matrix for that object
double_3dvec that_matrix = that->get_matrix(max_order);
xt::xtensor<double, 3> that_matrix = that->get_matrix(max_order);
// Now add that to this for the scattering and multiplicity
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++) {
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];
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;
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];
mult_denom(gin, gout) += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
mult_denom(gin, gout) += scalars[i];
}
i_gout++;
}
@ -107,16 +108,14 @@ ScattData::base_combine(int max_order,
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
xt::xtensor<double, 2> this_mult = xt::ones<double>({groups, groups});
for (int gin = 0; gin < groups; gin++) {
for (int gout = 0; gout < groups; gout++) {
if (mult_denom[gin][gout] > 0.) {
this_mult[gin][gout] = mult_numer[gin][gout] / mult_denom[gin][gout];
if (mult_denom(gin, gout) > 0.) {
this_mult(gin, gout) = mult_numer(gin, gout) / mult_denom(gin, gout);
}
}
}
mult_numer.clear();
mult_denom.clear();
// 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.
@ -125,8 +124,8 @@ ScattData::base_combine(int max_order,
int gmin_;
for (gmin_ = 0; gmin_ < groups; gmin_++) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmin_].size(); l++) {
if (this_matrix[gin][gmin_][l] != 0.) {
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmin_, l) != 0.) {
non_zero = true;
break;
}
@ -136,8 +135,8 @@ ScattData::base_combine(int max_order,
int gmax_;
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmax_].size(); l++) {
if (this_matrix[gin][gmax_][l] != 0.) {
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmax_, l) != 0.) {
non_zero = true;
break;
}
@ -160,8 +159,11 @@ ScattData::base_combine(int 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] = this_matrix[gin][gout];
sparse_mult[gin][i_gout] = this_mult[gin][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_mult[gin][i_gout] = this_mult(gin, gout);
i_gout++;
}
}
@ -241,8 +243,9 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
//==============================================================================
void
ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataLegendre::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -252,10 +255,9 @@ ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
// Get the scattering cross section value by summing the un-normalized P0
// coefficient in the variable matrix over all outgoing groups.
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
scattxs[gin] += matrix[gin][i_gout][0];
}
@ -401,8 +403,8 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups});
xt::xtensor<int, 1> in_gmax({groups});
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -418,20 +420,19 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
double_3dvec
xt::xtensor<double, 3>
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(groups,
double_1dvec(order_dim, 0.)));
xt::xtensor<double, 3> matrix({groups, groups, order_dim}, 0.);
for (int gin = 0; gin < groups; gin++) {
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] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
dist[gin][i_gout][l];
}
}
@ -444,8 +445,9 @@ ScattDataLegendre::get_matrix(int max_order)
//==============================================================================
void
ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataHistogram::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -455,9 +457,8 @@ ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
// Get the scattering cross section value by summing the distribution
// over all the histogram bins in angle and outgoing energy groups
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
scattxs[gin] += std::accumulate(matrix[gin][i_gout].begin(),
matrix[gin][i_gout].end(), 0.);
@ -484,12 +485,8 @@ ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Build the angular distribution mu values
mu = double_1dvec(order);
mu = xt::linspace(-1., 1., order + 1);
dmu = 2. / order;
mu[0] = -1.;
for (int imu = 1; imu < order; imu++) {
mu[imu] = -1. + imu * dmu;
}
// Calculate f(mu) and integrate it so we can avoid rejection sampling
fmu.resize(groups);
@ -581,21 +578,20 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
xt::xtensor<double, 3>
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(groups,
double_1dvec(order_dim, 0.)));
xt::xtensor<double, 3> matrix = xt::zeros<double>({groups, groups, order_dim});
for (int gin = 0; gin < groups; gin++) {
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] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
@ -624,8 +620,8 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups});
xt::xtensor<int, 1> in_gmax({groups});
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -633,7 +629,7 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
// 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);
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -644,8 +640,9 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
void
ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataTabular::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -654,19 +651,13 @@ ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
double_3dvec matrix = coeffs;
// Build the angular distribution mu values
mu = double_1dvec(order);
mu = xt::linspace(-1., 1., order);
dmu = 2. / (order - 1);
mu[0] = -1.;
for (int imu = 1; imu < order - 1; imu++) {
mu[imu] = -1. + imu * dmu;
}
mu[order - 1] = 1.;
// Get the scattering cross section value by integrating the distribution
// over all mu points and then combining over all outgoing groups
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
for (int imu = 1; imu < order; imu++) {
scattxs[gin] += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
@ -743,7 +734,7 @@ ScattDataTabular::calc_f(int gin, int gout, double mu)
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
imu = this->mu.shape()[0] - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
@ -764,7 +755,7 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
sample_energy(gin, gout, i_gout);
// Determine the outgoing cosine bin
int NP = this->mu.size();
int NP = this->mu.shape()[0];
double xi = prn();
double c_k = dist[gin][i_gout][0];
@ -804,21 +795,20 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
xt::xtensor<double, 3>
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(groups,
double_1dvec(order_dim, 0.)));
xt::xtensor<double, 3> matrix({groups, groups, order_dim}, 0.);
for (int gin = 0; gin < groups; gin++) {
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] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
@ -847,8 +837,8 @@ ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
int groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups});
xt::xtensor<int, 1> in_gmax({groups});
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -856,7 +846,7 @@ ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
// 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);
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -885,13 +875,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
tab.scattxs = leg.scattxs;
// Build mu and dmu
tab.mu = double_1dvec(n_mu);
tab.mu = xt::linspace(-1., 1., n_mu);
tab.dmu = 2. / (n_mu - 1);
tab.mu[0] = -1.;
for (int imu = 1; imu < n_mu - 1; imu++) {
tab.mu[imu] = -1. + imu * tab.dmu;
}
tab.mu[n_mu - 1] = 1.;
// Calculate f(mu) and integrate it so we can avoid rejection sampling
int groups = tab.energy.size();

359
src/xsdata.cpp
View file

@ -5,6 +5,9 @@
#include <algorithm>
#include <numeric>
#include "xtensor/xbuilder.hpp"
#include "xtensor/xview.hpp"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -20,7 +23,7 @@ namespace openmc {
XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi)
{
int n_ang = n_pol * n_azi;
size_t 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 &&
@ -28,31 +31,33 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
fatal_error("Invalid scatter_format!");
}
// allocate all [temperature][phi][theta][in group] quantities
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.));
std::vector<size_t> shape = {n_ang, energy_groups};
total = xt::zeros<double>(shape);
absorption = xt::zeros<double>(shape);
inverse_velocity = xt::zeros<double>(shape);
if (fissionable) {
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.));
fission = xt::zeros<double>(shape);
nu_fission = xt::zeros<double>(shape);
prompt_nu_fission = xt::zeros<double>(shape);
kappa_fission = xt::zeros<double>(shape);
}
// allocate decay_rate; [temperature][phi][theta][delayed group]
decay_rate = double_2dvec(n_ang, double_1dvec(num_delayed_groups, 0.));
shape[1] = num_delayed_groups;
decay_rate = xt::zeros<double>(shape);
if (fissionable) {
shape = {n_ang, energy_groups, num_delayed_groups};
// allocate delayed_nu_fission; [temperature][phi][theta][in group][delay group]
delayed_nu_fission = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(num_delayed_groups, 0.)));
delayed_nu_fission = xt::zeros<double>(shape);
// chi_prompt; [temperature][phi][theta][in group][delayed group]
chi_prompt = double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups, 0.)));
// chi_prompt; [temperature][phi][theta][in group][out group]
shape = {n_ang, energy_groups, energy_groups};
chi_prompt = xt::zeros<double>(shape);
// chi_delayed; [temperature][phi][theta][in group][out group][delay group]
chi_delayed = double_4dvec(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.))));
shape = {n_ang, energy_groups, energy_groups, num_delayed_groups};
chi_delayed = xt::zeros<double>(shape);
}
@ -79,8 +84,8 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
{
// Reconstruct the dimension information so it doesn't need to be passed
int n_ang = n_pol * n_azi;
int energy_groups = total[0].size();
int delayed_groups = decay_rate[0].size();
int energy_groups = total.shape()[1];
int delayed_groups = decay_rate.shape()[1];
// Set the fissionable-specific data
if (fissionable) {
@ -100,7 +105,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
// denominator in tally methods
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;
if (absorption(a, gin) == 0.) absorption(a, gin) = 1.e-10;
}
}
@ -110,7 +115,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
} else {
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];
total(a, gin) = absorption(a, gin) + scatter[a]->scattxs[gin];
}
}
}
@ -118,7 +123,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
// 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;
if (total(a, gin) == 0.) total(a, gin) = 1.e-10;
}
}
}
@ -129,14 +134,13 @@ void
XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int delayed_groups, bool is_isotropic)
{
int n_ang = n_pol * n_azi;
size_t 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_3dvec temp_beta =double_3dvec(n_ang, double_2dvec(energy_groups,
double_1dvec(delayed_groups, 0.)));
xt::xtensor<double, 3> temp_beta({n_ang, energy_groups, delayed_groups}, 0.);
if (object_exists(xsdata_grp, "beta")) {
hid_t xsdata = open_dataset(xsdata_grp, "beta");
int ndims = dataset_ndims(xsdata);
@ -146,7 +150,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
if (ndims == 3) {
// Beta is input as [delayed group]
double_1dvec temp_arr(n_pol * n_azi * delayed_groups);
std::vector<double> temp_arr({n_pol * n_azi * delayed_groups});
read_nd_vector(xsdata_grp, "beta", temp_arr);
// Broadcast to all incoming groups
@ -154,9 +158,9 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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++];
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];
temp_beta(a, gin, dg) = temp_beta(a, 0, dg);
}
}
}
@ -170,29 +174,33 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// If chi is provided, set chi-prompt and chi-delayed
if (object_exists(xsdata_grp, "chi")) {
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
xt::xtensor<double, 2> temp_arr ({n_ang, energy_groups});
read_nd_vector(xsdata_grp, "chi", temp_arr);
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];
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.);
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_prompt(a, 0, gout);
}
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;
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];
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt(a, gin, gout) = chi_prompt(a, 0, gout);
}
}
// Finally set chi-delayed equal to chi-prompt
@ -200,8 +208,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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];
chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout);
}
}
}
@ -224,15 +231,18 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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];
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[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
double beta_sum = 0.;
for (int gin = 0; gin < energy_groups; gin++) {
beta_sum += temp_beta(a, gin);
}
prompt_nu_fission(a, gin) *= (1. - beta_sum);
}
}
}
@ -244,14 +254,17 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// Normalize the chi info so the CDF is 1.
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.);
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_prompt(a, gin, gout);
}
// Set the vector nu-fission from the matrix nu-fission
prompt_nu_fission[a][gin] = chi_sum;
prompt_nu_fission(a, gin) = chi_sum;
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] /= chi_sum;
chi_prompt(a, gin, gout) /= chi_sum;
}
} else {
fatal_error("Encountered chi for a group that is <= 0!");
@ -262,8 +275,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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];
chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout);
}
}
}
@ -271,16 +283,17 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// 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];
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[a][gin].begin(),
temp_beta[a][gin].end(), 0.);
prompt_nu_fission[a][gin] *= (1. - beta_sum);
double beta_sum = 0.;
for (int dg = 0; dg < delayed_groups; dg++) {
beta_sum += temp_beta(a, gin, dg);
}
prompt_nu_fission(a, gin) *= (1. - beta_sum);
}
}
}
@ -293,21 +306,24 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// If chi-prompt is provided, set chi-prompt
if (object_exists(xsdata_grp, "chi-prompt")) {
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
xt::xtensor<double, 2> temp_arr({n_ang, energy_groups});
read_nd_vector(xsdata_grp, "chi-prompt", temp_arr);
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];
chi_prompt(a, gin, gout) = temp_arr(a, gout);
}
// Normalize chi so its CDF goes to 1
double chi_sum = std::accumulate(chi_prompt[a][gin].begin(),
chi_prompt[a][gin].end(), 0.);
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += chi_prompt(a, gin, gout);
}
if (chi_sum >= 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_prompt[a][gin][gout] /= chi_sum;
chi_prompt(a, gin, gout) /= chi_sum;
}
} else {
fatal_error("Encountered chi-prompt for a group that is <= 0.!");
@ -326,13 +342,16 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
if (ndims == 3) {
// chi-delayed is a [in group] vector
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
xt::xtensor<double, 2> temp_arr({n_ang, energy_groups});
read_nd_vector(xsdata_grp, "chi-delayed", temp_arr);
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.);
double chi_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
chi_sum += temp_arr(a, gout);
}
if (chi_sum <= 0.) {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
}
@ -341,7 +360,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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;
chi_delayed(a, gin, gout, dg) = temp_arr(a, gout) / chi_sum;
}
}
}
@ -356,12 +375,12 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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];
chi_sum += chi_delayed(a, gin, gout, dg);
}
if (chi_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[a][gin][gout][dg] /= chi_sum;
chi_delayed(a, gin, gout, dg) /= chi_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
@ -384,29 +403,30 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
if (ndims == 3) {
// prompt-nu-fission is a [in group] vector
read_nd_vector(xsdata_grp, "prompt-nu-fission",
prompt_nu_fission);
read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission);
} else if (ndims == 4) {
// prompt nu fission is a matrix,
// so set prompt_nu_fiss & chi_prompt
double_3dvec temp_arr(n_ang, double_2dvec(energy_groups,
double_1dvec(energy_groups)));
xt::xtensor<double, 3> temp_arr({n_ang, energy_groups, 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_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;
double prompt_sum = 0.;
for (int gout = 0; gout < energy_groups; gout++) {
prompt_sum += temp_arr(a, gin, gout);
}
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[a][gin] > 0.) {
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];
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!");
@ -429,19 +449,19 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
if (ndims == 3) {
// delayed-nu-fission is an [in group] vector
if (temp_beta[0][0][0] == 0.) {
if (temp_beta(0, 0, 0) == 0.) {
fatal_error("cannot set delayed-nu-fission with a 1D array if "
"beta is not provided");
}
double_2dvec temp_arr(n_ang, double_1dvec(energy_groups));
xt::xtensor<double, 2> temp_arr({n_ang, energy_groups});
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
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];
delayed_nu_fission(a, gin, dg) =
temp_beta(a, gin, dg) * temp_arr(a, gin);
}
}
}
@ -451,9 +471,9 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
delayed_nu_fission);
} else if (ndims == 5) {
// This will contain delayed-nu-fision and chi-delayed data
double_4dvec temp_arr(n_ang, double_3dvec(energy_groups,
double_2dvec(energy_groups, double_1dvec(delayed_groups))));
// This will contain delayed-nu-fission and chi-delayed data
xt::xtensor<double, 4> temp_arr({n_ang, energy_groups, energy_groups,
delayed_groups});
read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr);
// Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix
@ -463,14 +483,14 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
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];
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;
delayed_nu_fission(a, gin, dg) = gout_sum;
// Normalize chi-delayed
if (gout_sum > 0.) {
for (int gout = 0; gout < energy_groups; gout++) {
chi_delayed[a][gin][gout][dg] /= gout_sum;
chi_delayed(a, gin, gout, dg) /= gout_sum;
}
} else {
fatal_error("Encountered chi-delayed for a group that is <= 0!");
@ -488,10 +508,10 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// Combine prompt_nu_fission and delayed_nu_fission into nu_fission
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]);
nu_fission(a, gin) = prompt_nu_fission(a, gin);
for (int dg = 0; dg < delayed_groups; dg++) {
nu_fission(a, gin) += delayed_nu_fission(a, gin, dg);
}
}
}
}
@ -503,94 +523,92 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int scatter_format, int final_scatter_format,
int order_data, int max_order, int legendre_to_tabular_points)
{
int n_ang = n_pol * n_azi;
size_t 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_2dvec gmin(n_ang, int_1dvec(energy_groups));
xt::xtensor<int, 2> gmin({n_ang, energy_groups});
read_nd_vector(scatt_grp, "g_min", gmin, true);
int_2dvec gmax(n_ang, int_1dvec(energy_groups));
xt::xtensor<int, 2> gmax({n_ang, 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 a = 0; a < n_ang; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
gmin[a][gin] -= 1;
gmax[a][gin] -= 1;
}
}
gmin -= 1;
gmax -= 1;
// Now use this info to find the length of a vector to hold the flattened
// data.
int length = 0;
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);
length += order_data * (gmax(a, gin) - gmin(a, gin) + 1);
}
}
double_1dvec temp_arr(length);
read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
// Compare the number of orders given with the max order of the problem;
// strip off the superfluous orders if needed
int order_dim;
if (scatter_format == ANGLE_LEGENDRE) {
order_dim = std::min(order_data - 1, max_order) + 1;
} else {
order_dim = order_data;
}
// convert the flattened temp_arr to a jagged array for passing to
// scatt data
double_4dvec input_scatt(n_ang, double_3dvec(energy_groups));
//temp_arr scope
{
std::vector<double> temp_arr(length);
read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
int temp_idx = 0;
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++];
// Compare the number of orders given with the max order of the problem;
// strip off the superfluous orders if needed
int order_dim;
if (scatter_format == ANGLE_LEGENDRE) {
order_dim = std::min(order_data - 1, max_order) + 1;
} else {
order_dim = order_data;
}
// convert the flattened temp_arr to a jagged array for passing to
// scatt data
int temp_idx = 0;
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);
}
// Adjust index for the orders we didnt take
temp_idx += (order_data - order_dim);
}
}
}
temp_arr.clear();
// Get multiplication matrix
double_3dvec temp_mult(n_ang, double_2dvec(energy_groups));
if (object_exists(scatt_grp, "multiplicity_matrix")) {
temp_arr.resize(length / order_data);
std::vector<double> temp_arr(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 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);
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++];
}
}
}
temp_arr.clear();
} else {
// Use a default: multiplicities are 1.0.
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);
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.;
}
}
}
}
temp_arr.clear();
close_group(scatt_grp);
// Finally, convert the Legendre data to tabular, if needed
@ -598,7 +616,10 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
final_scatter_format == ANGLE_TABULAR) {
for (int a = 0; a < n_ang; a++) {
ScattDataLegendre legendre_scatt;
legendre_scatt.init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
xt::xtensor<int, 1> in_gmax = xt::view(gmax, a, xt::all());
legendre_scatt.init(in_gmin, in_gmax,
temp_mult[a], input_scatt[a]);
// Now create a tabular version of legendre_scatt
convert_legendre_to_tabular(legendre_scatt,
@ -611,7 +632,9 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
// We are sticking with the current representation
// Initialize the ScattData object with this data
for (int a = 0; a < n_ang; a++) {
scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]);
scatter[a]->init(xt::view(gmin, a, xt::all()),
xt::view(gmax, a, xt::all()),
temp_mult[a], input_scatt[a]);
}
}
}
@ -627,70 +650,25 @@ 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 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[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[a][gin].size(); gout++) {
chi_prompt[a][gin][gout] +=
scalar * that->chi_prompt[a][gin][gout];
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[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[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[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;
}
}
}
}
}
total += scalar * that->total;
absorption += scalar * that->absorption;
if (i == 0) {
inverse_velocity = that->inverse_velocity;
}
if (that->prompt_nu_fission.shape()[0] > 0) {
nu_fission += scalar * that->nu_fission;
prompt_nu_fission += scalar * that->prompt_nu_fission;
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;
}
decay_rate += scalar * that->decay_rate;
}
// Allow the ScattData object to combine itself
for (int a = 0; a < total.size(); a++) {
for (int a = 0; a < total.shape()[0]; 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++) {
@ -707,8 +685,7 @@ XsData::combine(const std::vector<XsData*>& those_xs,
bool
XsData::equiv(const XsData& that)
{
return ((absorption.size() == that.absorption.size()) &&
(absorption[0].size() == that.absorption[0].size()));
return (absorption.shape() == that.absorption.shape());
}
} //namespace openmc

9
tests/regression_tests/mgxs_library_ce_to_mg/test.py
View file

@ -71,10 +71,11 @@ class MGXSTestHarness(PyAPITestHarness):
openmc.run(openmc_exec=config['exe'])
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
pass
# super()._cleanup()
# f = 'mgxs.h5'
# if os.path.exists(f):
# os.remove(f)
def test_mgxs_library_ce_to_mg():