crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

clearing up const and making use of generic base class methods

Browse source
This commit is contained in:
Adam G Nelson 2018-06-16 09:25:49 -04:00
parent 2b8c2075d0
commit 863f91b1e7
6 changed files with 232 additions and 387 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/api.F90
View file

@ -129,7 +129,7 @@ contains
index_ufs_mesh = -1
keff = ONE
legendre_to_tabular = .true.
legendre_to_tabular_points = 33
legendre_to_tabular_points = C_NONE
n_batch_interval = 1
n_lost_particles = 0
n_particles = 0

505
src/scattdata.cpp
View file

@ -41,6 +41,125 @@ ScattData::generic_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
void
ScattData::generic_combine(const 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)
{
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.));
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattData* that = those_scatts[i];
// Build the dense matrix for that object
double_3dvec 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++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
}
i_gout++;
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
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];
}
}
}
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.
for (int gin = 0; gin < groups; gin++) {
// Find the minimum and maximum group boundaries
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
// treat the case of all values being 0
if (gmin_ > gmax_) {
gmin_ = gin;
gmax_ = gin;
}
// Store the group bounds
in_gmin[gin] = gmin_;
in_gmax[gin] = gmax_;
// Store the data in the compressed format
sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
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];
i_gout++;
}
}
}
//==============================================================================
void
ScattData::sample_energy(int gin, int& gout, int& i_gout)
{
@ -60,7 +179,8 @@ ScattData::sample_energy(int gin, int& gout, int& i_gout)
//==============================================================================
double
ScattData::get_xs(const int xstype, int gin, int* gout, double* mu)
ScattData::get_xs(const int xstype, const int gin, const int* gout,
const double* mu)
{
// Set the outgoing group offset index as needed
int i_gout = 0;
@ -214,7 +334,7 @@ ScattDataLegendre::update_max_val()
//==============================================================================
double
ScattDataLegendre::calc_f(int gin, int gout, double mu)
ScattDataLegendre::calc_f(const int gin, const int gout, const double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -230,7 +350,7 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
//==============================================================================
void
ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
ScattDataLegendre::sample(const int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -261,8 +381,8 @@ ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
void
ScattDataLegendre::combine(std::vector<ScattData*>& those_scatts,
double_1dvec& scalars)
ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = 0;
@ -277,120 +397,18 @@ ScattDataLegendre::combine(std::vector<ScattData*>& those_scatts,
}
max_order++; // Add one since this is a Legendre
// Get the groups as a shorthand
int groups = dynamic_cast<ScattDataLegendre*>(those_scatts[0])->energy.size();
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.));
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattDataLegendre* that = dynamic_cast<ScattDataLegendre*>(those_scatts[i]);
// Build the dense matrix for that object
double_3dvec 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++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
}
i_gout++;
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
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];
}
}
}
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.
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
for (int gin = 0; gin < groups; gin++) {
// Find the minimum and maximum group boundaries
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
// treat the case of all values being 0
if (gmin_ > gmax_) {
gmin_ = gin;
gmax_ = gin;
}
// Store the group bounds
in_gmin[gin] = gmin_;
in_gmax[gin] = gmax_;
// Store the data in the compressed format
sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
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];
i_gout++;
}
}
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::generic_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -399,7 +417,7 @@ ScattDataLegendre::combine(std::vector<ScattData*>& those_scatts,
//==============================================================================
double_3dvec
ScattDataLegendre::get_matrix(int max_order)
ScattDataLegendre::get_matrix(const int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
@ -504,7 +522,7 @@ ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
double
ScattDataHistogram::calc_f(int gin, int gout, double mu)
ScattDataHistogram::calc_f(const int gin, const int gout, const double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -528,7 +546,7 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu)
//==============================================================================
void
ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
ScattDataHistogram::sample(const int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -563,7 +581,7 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataHistogram::get_matrix(int max_order)
ScattDataHistogram::get_matrix(const int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
@ -587,8 +605,8 @@ ScattDataHistogram::get_matrix(int max_order)
//==============================================================================
void
ScattDataHistogram::combine(std::vector<ScattData*>& those_scatts,
double_1dvec& scalars)
ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order;
@ -605,120 +623,18 @@ ScattDataHistogram::combine(std::vector<ScattData*>& those_scatts,
}
}
// Get the groups as a shorthand
int groups = dynamic_cast<ScattDataHistogram*>(those_scatts[0])->energy.size();
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.));
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattDataHistogram* that = dynamic_cast<ScattDataHistogram*>(those_scatts[i]);
// Build the dense matrix for that object
double_3dvec 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++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
}
i_gout++;
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
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];
}
}
}
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.
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
for (int gin = 0; gin < groups; gin++) {
// Find the minimum and maximum group boundaries
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
// treat the case of all values being 0
if (gmin_ > gmax_) {
gmin_ = gin;
gmax_ = gin;
}
// Store the group bounds
in_gmin[gin] = gmin_;
in_gmax[gin] = gmax_;
// Store the data in the compressed format
sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
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];
i_gout++;
}
}
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::generic_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -818,7 +734,7 @@ ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
double
ScattDataTabular::calc_f(int gin, int gout, double mu)
ScattDataTabular::calc_f(const int gin, const int gout, const double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -843,7 +759,7 @@ ScattDataTabular::calc_f(int gin, int gout, double mu)
//==============================================================================
void
ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
ScattDataTabular::sample(const int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -891,7 +807,7 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataTabular::get_matrix(int max_order)
ScattDataTabular::get_matrix(const int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
@ -915,8 +831,8 @@ ScattDataTabular::get_matrix(int max_order)
//==============================================================================
void
ScattDataTabular::combine(std::vector<ScattData*>& those_scatts,
double_1dvec& scalars)
ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order;
@ -933,120 +849,18 @@ ScattDataTabular::combine(std::vector<ScattData*>& those_scatts,
}
}
// Get the groups as a shorthand
int groups = dynamic_cast<ScattDataTabular*>(those_scatts[0])->energy.size();
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.));
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
// To combine from nuclidic data we need to use the final relationship
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// Developed as follows:
// mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'}
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'})
// mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
// sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
// nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member
// variables
for (int i = 0; i < those_scatts.size(); i++) {
ScattDataTabular* that = dynamic_cast<ScattDataTabular*>(those_scatts[i]);
// Build the dense matrix for that object
double_3dvec 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++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
}
i_gout++;
}
}
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
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];
}
}
}
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.
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
for (int gin = 0; gin < groups; gin++) {
// Find the minimum and maximum group boundaries
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
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.) {
non_zero = true;
break;
}
}
if (non_zero) break;
}
// treat the case of all values being 0
if (gmin_ > gmax_) {
gmin_ = gin;
gmax_ = gin;
}
// Store the group bounds
in_gmin[gin] = gmin_;
in_gmax[gin] = gmax_;
// Store the data in the compressed format
sparse_scatter[gin].resize(gmax_ - gmin_ + 1);
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];
i_gout++;
}
}
// The rest of the steps do not depend on the type of angular representation
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::generic_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -1060,6 +874,17 @@ void
convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu)
{
// 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) {
n_mu = 2;
} else {
n_mu = DEFAULT_NMU;
}
}
tab.generic_init(n_mu, leg.gmin, leg.gmax, leg.energy, leg.mult);
tab.scattxs = leg.scattxs;

60
src/scattdata.h
View file

@ -16,35 +16,42 @@
namespace openmc {
//==============================================================================
// SCATTDATA contains all the data needed to describe the scattering energy and
// angular distribution data
//==============================================================================
// temporary declaations so we can name our friend functions
class ScattDataLegendre;
class ScattDataTabular;
//==============================================================================
// SCATTDATA contains all the data needed to describe the scattering energy and
// angular distribution data
//==============================================================================
class ScattData {
protected:
void generic_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_energy, double_2dvec& in_mult);
void generic_combine(const 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);
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
public:
double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
virtual double calc_f(int gin, int gout, double mu) = 0;
virtual void sample(int gin, int& gout, double& mu, double& wgt) = 0;
virtual double calc_f(const int gin, const int gout, const double mu) = 0;
virtual void sample(const int gin, int& gout, double& mu, double& wgt) = 0;
virtual void init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_mult, double_3dvec& coeffs) = 0;
void sample_energy(int gin, int& gout, int& i_gout);
double get_xs(const int xstype, int gin, int* gout, double* mu);
void generic_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_energy, double_2dvec& in_mult);
virtual void combine(std::vector<ScattData*>& those_scatts,
double_1dvec& scalars) = 0;
double get_xs(const int xstype, const int gin, const int* gout,
const double* mu);
virtual void combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars) = 0;
virtual int get_order() = 0;
virtual double_3dvec get_matrix(int max_order) = 0;
virtual double_3dvec get_matrix(const int max_order) = 0;
};
//==============================================================================
@ -61,11 +68,12 @@ class ScattDataLegendre: public ScattData {
void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
void update_max_val();
double calc_f(int gin, int gout, double mu);
void sample(int gin, int& gout, double& mu, double& wgt);
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
double calc_f(const int gin, const int gout, const double mu);
void sample(const int gin, int& gout, double& mu, double& wgt);
void combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
int get_order() {return dist[0][0].size() - 1;};
double_3dvec get_matrix(int max_order);
double_3dvec get_matrix(const int max_order);
};
//==============================================================================
@ -81,11 +89,12 @@ class ScattDataHistogram: public ScattData {
public:
void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
double calc_f(int gin, int gout, double mu);
void sample(int gin, int& gout, double& mu, double& wgt);
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
double calc_f(const int gin, const int gout, const double mu);
void sample(const int gin, int& gout, double& mu, double& wgt);
void combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
int get_order() {return dist[0][0].size();};
double_3dvec get_matrix(int max_order);
double_3dvec get_matrix(const int max_order);
};
//==============================================================================
@ -103,11 +112,12 @@ class ScattDataTabular: public ScattData {
public:
void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
double calc_f(int gin, int gout, double mu);
void sample(int gin, int& gout, double& mu, double& wgt);
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
double calc_f(const int gin, const int gout, const double mu);
void sample(const int gin, int& gout, double& mu, double& wgt);
void combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
int get_order() {return dist[0][0].size();};
double_3dvec get_matrix(int max_order);
double_3dvec get_matrix(const int max_order);
};
//==============================================================================

2
src/settings.F90
View file

@ -36,7 +36,7 @@ module settings
logical :: legendre_to_tabular = .true.
! Number of points to use in the Legendre to tabular conversion
integer(C_INT) :: legendre_to_tabular_points = 33
integer(C_INT) :: legendre_to_tabular_points = C_NONE
! ============================================================================
! SIMULATION VARIABLES

27
src/xsdata.cpp
View file

@ -6,8 +6,9 @@ namespace openmc {
// XsData class methods
//==============================================================================
XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi)
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)
{
// check to make sure scatter format is OK before we allocate
if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
@ -69,9 +70,10 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
//==============================================================================
void
XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic)
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)
{
// Reconstruct the dimension information so it doesn't need to be passed
int n_pol = total.size();
@ -132,8 +134,9 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
//==============================================================================
void
XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int delayed_groups, bool is_isotropic)
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)
{
// Get the fission and kappa_fission data xs; these are optional
@ -532,9 +535,10 @@ XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
//==============================================================================
void
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)
XsData::_scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
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)
{
if (!object_exists(xsdata_grp, "scatter_data")) {
fatal_error("Must provide scatter_data group!");
@ -671,7 +675,8 @@ XsData::_scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
//==============================================================================
void
XsData::combine(std::vector<XsData*> those_xs, double_1dvec& scalars)
XsData::combine(const std::vector<XsData*> those_xs,
const double_1dvec& scalars)
{
// Combine the non-scattering data
for (int i = 0; i < those_xs.size(); i++) {

23
src/xsdata.h
View file

@ -24,9 +24,10 @@ namespace openmc {
class XsData {
private:
void _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);
void _scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
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);
public:
@ -55,12 +56,16 @@ class XsData {
std::vector<std::vector<ScattData*> > scatter;
XsData() = default;
XsData(int num_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi);
void from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic);
void combine(std::vector<XsData*> those_xs, double_1dvec& scalars);
XsData(const int num_groups, const int num_delayed_groups,
const bool fissionable, const int scatter_format, const int n_pol,
const int n_azi);
void 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);
void combine(const std::vector<XsData*> those_xs,
const double_1dvec& scalars);
bool equiv(const XsData& that);
};