mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Implemented sample_scatter in C++
This commit is contained in:
Adam G Nelson
parent
59c4bc3fc4
commit
6c8de73e3e
8 changed files with 156 additions and 69 deletions
|
|
@ -10,7 +10,6 @@
|
|||
|
||||
namespace openmc {
|
||||
|
||||
typedef std::array<double, 3> dir_arr;
|
||||
typedef std::vector<double> double_1dvec;
|
||||
typedef std::vector<std::vector<double> > double_2dvec;
|
||||
typedef std::vector<std::vector<std::vector<double> > > double_3dvec;
|
||||
|
|
|
|||
72
src/mgxs.cpp
72
src/mgxs.cpp
|
|
@ -24,6 +24,13 @@ void Mgxs::init(const std::string& in_name, const double in_awr,
|
|||
azimuthal = in_azimuthal;
|
||||
n_pol = polar.size();
|
||||
n_azi = azimuthal.size();
|
||||
index_temp = 0;
|
||||
last_sqrtkT = 0.;
|
||||
index_pol = 0;
|
||||
index_azi = 0;
|
||||
last_uvw[0] = 1.;
|
||||
last_uvw[1] = 0.;
|
||||
last_uvw[2] = 0.;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -383,7 +390,8 @@ void Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
|
|||
}
|
||||
|
||||
|
||||
double Mgxs::get_xs(const char* xstype, int gin, int* gout, double* mu, int* dg)
|
||||
double Mgxs::get_xs(const char* xstype, const int gin, int* gout, double* mu,
|
||||
int* dg)
|
||||
{
|
||||
// This method assumes that the temperature and angle indices are set
|
||||
double val;
|
||||
|
|
@ -469,7 +477,8 @@ double Mgxs::get_xs(const char* xstype, int gin, int* gout, double* mu, int* dg)
|
|||
}
|
||||
|
||||
|
||||
void Mgxs::sample_fission_energy(int gin, double nu_fission, int& dg, int& gout)
|
||||
void Mgxs::sample_fission_energy(const int gin, const double nu_fission,
|
||||
int& dg, int& gout)
|
||||
{
|
||||
// This method assumes that the temperature and angle indices are set
|
||||
// Find the probability of having a prompt neutron
|
||||
|
|
@ -525,8 +534,7 @@ void Mgxs::sample_fission_energy(int gin, double nu_fission, int& dg, int& gout)
|
|||
}
|
||||
|
||||
|
||||
void Mgxs::sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
|
||||
double& wgt)
|
||||
void Mgxs::sample_scatter(const int gin, int& gout, double& mu, double& wgt)
|
||||
{
|
||||
// This method assumes that the temperature and angle indices are set
|
||||
// Sample the data
|
||||
|
|
@ -534,8 +542,8 @@ void Mgxs::sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
|
|||
}
|
||||
|
||||
|
||||
void Mgxs::calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
|
||||
double& abs_xs, double& nu_fiss_xs)
|
||||
void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
|
||||
double& total_xs, double& abs_xs, double& nu_fiss_xs)
|
||||
{
|
||||
// Set our indices
|
||||
set_temperature_index(sqrtkT);
|
||||
|
|
@ -543,10 +551,14 @@ void Mgxs::calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
|
|||
total_xs = xs[index_temp].total[index_pol][index_azi][gin];
|
||||
abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin];
|
||||
|
||||
// nu-fission is made up of the prompt and all the delayed nu_fission data
|
||||
nu_fiss_xs = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
|
||||
for (auto& val : xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin]) {
|
||||
nu_fiss_xs += val;
|
||||
if (fissionable) {
|
||||
// nu-fission is made up of the prompt and all the delayed nu_fission data
|
||||
nu_fiss_xs = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
|
||||
for (auto& val : xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin]) {
|
||||
nu_fiss_xs += val;
|
||||
}
|
||||
} else {
|
||||
nu_fiss_xs = 0.;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -568,7 +580,7 @@ bool Mgxs::equiv(const Mgxs& that)
|
|||
}
|
||||
|
||||
|
||||
inline void Mgxs::set_temperature_index(double sqrtkT)
|
||||
inline void Mgxs::set_temperature_index(const double sqrtkT)
|
||||
{
|
||||
// See if we need to find the new index
|
||||
if (sqrtkT != last_sqrtkT) {
|
||||
|
|
@ -588,27 +600,30 @@ inline void Mgxs::set_temperature_index(double sqrtkT)
|
|||
}
|
||||
|
||||
|
||||
inline void Mgxs::set_angle_index(dir_arr& uvw)
|
||||
inline void Mgxs::set_angle_index(const double uvw[3])
|
||||
{
|
||||
// See if we need to find the new index
|
||||
if (uvw != last_uvw) {
|
||||
if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) ||
|
||||
(uvw[2] != last_uvw[2])) {
|
||||
// convert uvw to polar and azimuthal angles
|
||||
double my_pol = std::acos(uvw[2]);
|
||||
double my_azi = std::atan2(uvw[1], uvw[0]);
|
||||
|
||||
// Find the location, assuming equal-bin angles
|
||||
double delta_angle = PI / n_pol;
|
||||
index_pol = std::floor(my_pol / delta_angle + 1.);
|
||||
delta_angle = PI / n_azi;
|
||||
index_azi = std::floor((my_azi + PI) / delta_angle + 1.);
|
||||
index_pol = std::floor(my_pol / delta_angle);
|
||||
delta_angle = 2. * PI / n_azi;
|
||||
index_azi = std::floor((my_azi + PI) / delta_angle);
|
||||
|
||||
// store this direction as the last one used
|
||||
last_uvw = uvw;
|
||||
last_uvw[0] = uvw[0];
|
||||
last_uvw[1] = uvw[1];
|
||||
last_uvw[2] = uvw[2];
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Mgxs data loading methods
|
||||
// Mgxs data loading interface methods
|
||||
//==============================================================================
|
||||
|
||||
void add_mgxs(hid_t file_id, char* name, int energy_groups,
|
||||
|
|
@ -680,4 +695,25 @@ void create_macro_xs(char* mat_name, const int n_nuclides,
|
|||
macro_xs.push_back(macro);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Mgxs tracking/transport/tallying interface methods
|
||||
//==============================================================================
|
||||
|
||||
void calculate_xs(const int i_mat, const int gin, const double sqrtkT,
|
||||
const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
|
||||
{
|
||||
macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
|
||||
nu_fiss_xs);
|
||||
}
|
||||
|
||||
void scatter(const int i_mat, const int gin, int& gout, double& mu,
|
||||
double& wgt, double uvw[3])
|
||||
{
|
||||
int gout_c = gout - 1;
|
||||
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
|
||||
gout = gout_c + 1;
|
||||
|
||||
rotate_angle_c(uvw, mu, nullptr);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
25
src/mgxs.h
25
src/mgxs.h
|
|
@ -44,7 +44,7 @@ class Mgxs {
|
|||
int index_azi;
|
||||
double_1dvec polar;
|
||||
double_1dvec azimuthal;
|
||||
dir_arr last_uvw;
|
||||
double last_uvw[3];
|
||||
void _metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
|
||||
const int in_num_delayed_groups, double_1dvec& temperature,
|
||||
int& method, const double tolerance, int_1dvec& temps_to_read,
|
||||
|
|
@ -66,16 +66,16 @@ class Mgxs {
|
|||
double_1dvec& temperature, int& method, double tolerance,
|
||||
int max_order, bool legendre_to_tabular,
|
||||
int legendre_to_tabular_points);
|
||||
double get_xs(const char* xstype, int gin, int* gout, double* mu,
|
||||
double get_xs(const char* xstype, const int gin, int* gout, double* mu,
|
||||
int* dg);
|
||||
void sample_fission_energy(int gin, double nu_fission, int& dg, int& gout);
|
||||
void sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu,
|
||||
double& wgt);
|
||||
void calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs,
|
||||
double& abs_xs, double& nu_fiss_xs);
|
||||
void sample_fission_energy(const int gin, const double nu_fission,
|
||||
int& dg, int& gout);
|
||||
void sample_scatter(const int gin, int& gout, double& mu, double& wgt);
|
||||
void calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
|
||||
double& total_xs, double& abs_xs, double& nu_fiss_xs);
|
||||
bool equiv(const Mgxs& that);
|
||||
inline void set_temperature_index(double sqrtkT);
|
||||
inline void set_angle_index(dir_arr& uvw);
|
||||
inline void set_temperature_index(const double sqrtkT);
|
||||
inline void set_angle_index(const double uvw[3]);
|
||||
};
|
||||
|
||||
extern "C" void add_mgxs(hid_t file_id, char* name, int energy_groups,
|
||||
|
|
@ -89,6 +89,13 @@ extern "C" void create_macro_xs(char* mat_name, const int n_nuclides,
|
|||
const int i_nuclides[], const int n_temps, const double temps[],
|
||||
const double atom_densities[], int& method, const double tolerance);
|
||||
|
||||
extern "C" void calculate_xs(const int i_mat, const int gin,
|
||||
const double sqrtkT, const double uvw[3], double& total_xs,
|
||||
double& abs_xs, double& nu_fiss_xs);
|
||||
|
||||
extern "C" void scatter(const int i_mat, const int gin, int& gout, double& mu,
|
||||
double& wgt, double uvw[3]);
|
||||
|
||||
|
||||
// Storage for the MGXS data
|
||||
std::vector<Mgxs> nuclides_MG;
|
||||
|
|
|
|||
|
|
@ -166,10 +166,10 @@ module nuclide_header
|
|||
!===============================================================================
|
||||
|
||||
type MaterialMacroXS
|
||||
real(8) :: total ! macroscopic total xs
|
||||
real(8) :: absorption ! macroscopic absorption xs
|
||||
real(8) :: fission ! macroscopic fission xs
|
||||
real(8) :: nu_fission ! macroscopic production xs
|
||||
real(C_DOUBLE) :: total ! macroscopic total xs
|
||||
real(C_DOUBLE) :: absorption ! macroscopic absorption xs
|
||||
real(C_DOUBLE) :: fission ! macroscopic fission xs
|
||||
real(C_DOUBLE) :: nu_fission ! macroscopic production xs
|
||||
end type MaterialMacroXS
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -22,6 +22,20 @@ module physics_mg
|
|||
|
||||
implicit none
|
||||
|
||||
interface
|
||||
subroutine scatter_c(i_mat, gin, gout, mu, wgt, uvw) &
|
||||
bind(C, name='scatter')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: i_mat
|
||||
integer(C_INT), value, intent(in) :: gin
|
||||
integer(C_INT), intent(inout) :: gout
|
||||
real(C_DOUBLE), intent(inout) :: mu
|
||||
real(C_DOUBLE), intent(inout) :: wgt
|
||||
real(C_DOUBLE), intent(inout) :: uvw(1:3)
|
||||
end subroutine scatter_c
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -143,16 +157,18 @@ contains
|
|||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
|
||||
p % last_g, p % g, &
|
||||
p % mu, p % wgt)
|
||||
! call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, &
|
||||
! p % last_g, p % g, p % mu, p % wgt)
|
||||
|
||||
! Convert change in angle (mu) to new direction
|
||||
! p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
|
||||
|
||||
call scatter_c(p % material, p % last_g, p % g, p % mu, p % wgt, &
|
||||
p % coord(1) % uvw)
|
||||
|
||||
! Update energy value for downstream compatability (in tallying)
|
||||
p % E = energy_bin_avg(p % g)
|
||||
|
||||
! Convert change in angle (mu) to new direction
|
||||
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu)
|
||||
|
||||
! Set event component
|
||||
p % event = EVENT_SCATTER
|
||||
|
||||
|
|
|
|||
|
|
@ -44,7 +44,8 @@ void ScattData::sample_energy(int gin, int& gout, int& i_gout)
|
|||
{
|
||||
// Sample the outgoing group
|
||||
double xi = prn();
|
||||
i_gout = 0; //TODO: + 1?
|
||||
|
||||
i_gout = 0;
|
||||
gout = gmin[gin];
|
||||
double prob = energy[gin][i_gout];
|
||||
while((prob < xi) && (gout < gmax[gin])) {
|
||||
|
|
@ -247,7 +248,7 @@ void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
|
|||
}
|
||||
|
||||
|
||||
void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
|
||||
void ScattDataLegendre::combine(std::vector<ScattData*>& those_scatts,
|
||||
double_1dvec& scalars)
|
||||
{
|
||||
// Find the max order in the data set and make sure we can combine the sets
|
||||
|
|
@ -294,6 +295,7 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
|
|||
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++) {
|
||||
|
|
@ -301,13 +303,14 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
|
|||
}
|
||||
|
||||
// Incorporate that's contribution to the multiplicity matrix data
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
|
||||
mult_numer[gin][gout] += scalars[i] * nuscatt;
|
||||
if (that->mult[gin][gout] > 0.) {
|
||||
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][gout];
|
||||
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++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -336,14 +339,14 @@ void ScattDataLegendre::combine(std::vector<ScattData*> those_scatts,
|
|||
for (gmin_ = 0; gmin_ < groups; gmin_++) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmin_].begin(),
|
||||
this_matrix[gin][gmin_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
int gmax_;
|
||||
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmax_].begin(),
|
||||
this_matrix[gin][gmax_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
|
||||
|
|
@ -425,6 +428,7 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
|
|||
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
|
||||
double norm = std::accumulate(matrix[gin][i_gout].begin(),
|
||||
matrix[gin][i_gout].end(), 0.);
|
||||
in_energy[gin][i_gout] = norm;
|
||||
if (norm != 0.) {
|
||||
for (auto& n : matrix[gin][i_gout]) n /= norm;
|
||||
}
|
||||
|
|
@ -440,7 +444,7 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
|
|||
dmu = 2. / order;
|
||||
mu[0] = -1.;
|
||||
for (int imu = 1; imu < order; imu++) {
|
||||
mu[imu] = -1. + (imu - 1) * dmu;
|
||||
mu[imu] = -1. + imu * dmu;
|
||||
}
|
||||
|
||||
// Calculate f(mu) and integrate it so we can avoid rejection sampling
|
||||
|
|
@ -507,7 +511,7 @@ void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
|
|||
|
||||
int imu;
|
||||
if (xi < dist[gin][i_gout][0]) {
|
||||
imu = 1;
|
||||
imu = 0;
|
||||
} else {
|
||||
// TODO lower_bound? + 1?
|
||||
imu = std::upper_bound(dist[gin][i_gout].begin(),
|
||||
|
|
@ -551,7 +555,7 @@ double_3dvec ScattDataHistogram::get_matrix(int max_order)
|
|||
}
|
||||
|
||||
|
||||
void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
|
||||
void ScattDataHistogram::combine(std::vector<ScattData*>& those_scatts,
|
||||
double_1dvec& scalars)
|
||||
{
|
||||
// Find the max order in the data set and make sure we can combine the sets
|
||||
|
|
@ -600,6 +604,7 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
|
|||
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++) {
|
||||
|
|
@ -607,13 +612,14 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
|
|||
}
|
||||
|
||||
// Incorporate that's contribution to the multiplicity matrix data
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
|
||||
mult_numer[gin][gout] += scalars[i] * nuscatt;
|
||||
if (that->mult[gin][gout] > 0.) {
|
||||
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][gout];
|
||||
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++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -642,14 +648,14 @@ void ScattDataHistogram::combine(std::vector<ScattData*> those_scatts,
|
|||
for (gmin_ = 0; gmin_ < groups; gmin_++) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmin_].begin(),
|
||||
this_matrix[gin][gmin_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
int gmax_;
|
||||
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmax_].begin(),
|
||||
this_matrix[gin][gmax_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
|
||||
|
|
@ -695,7 +701,7 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
|
|||
dmu = 2. / (order - 1);
|
||||
mu[0] = -1.;
|
||||
for (int imu = 1; imu < order - 1; imu++) {
|
||||
mu[imu] = -1. + (imu - 1) * dmu;
|
||||
mu[imu] = -1. + imu * dmu;
|
||||
}
|
||||
mu[order - 1] = 1.;
|
||||
|
||||
|
|
@ -724,9 +730,7 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
|
|||
norm += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
|
||||
matrix[gin][i_gout][imu]);
|
||||
}
|
||||
if (norm != 0.) {
|
||||
for (auto& n : matrix[gin][i_gout]) n /= norm;
|
||||
}
|
||||
in_energy[gin][i_gout] = norm;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -806,14 +810,14 @@ void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
|
|||
|
||||
double c_k = dist[gin][i_gout][0];
|
||||
int k;
|
||||
for (k = 0; k < NP - 2; k++) {
|
||||
for (k = 0; k < NP - 1; k++) {
|
||||
double c_k1 = dist[gin][i_gout][k + 1];
|
||||
if (xi < c_k1) break;
|
||||
c_k = c_k1;
|
||||
}
|
||||
|
||||
// Check to make sure k is <= NP - 1
|
||||
k = std::min(k, NP - 1);
|
||||
k = std::min(k, NP - 2);
|
||||
|
||||
// Find the pdf values we want
|
||||
double p0 = fmu[gin][i_gout][k];
|
||||
|
|
@ -861,7 +865,7 @@ double_3dvec ScattDataTabular::get_matrix(int max_order)
|
|||
return matrix;
|
||||
}
|
||||
|
||||
void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
|
||||
void ScattDataTabular::combine(std::vector<ScattData*>& those_scatts,
|
||||
double_1dvec& scalars)
|
||||
{
|
||||
// Find the max order in the data set and make sure we can combine the sets
|
||||
|
|
@ -910,6 +914,7 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
|
|||
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++) {
|
||||
|
|
@ -917,13 +922,14 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
|
|||
}
|
||||
|
||||
// Incorporate that's contribution to the multiplicity matrix data
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][gout];
|
||||
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
|
||||
mult_numer[gin][gout] += scalars[i] * nuscatt;
|
||||
if (that->mult[gin][gout] > 0.) {
|
||||
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][gout];
|
||||
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++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -952,14 +958,14 @@ void ScattDataTabular::combine(std::vector<ScattData*> those_scatts,
|
|||
for (gmin_ = 0; gmin_ < groups; gmin_++) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmin_].begin(),
|
||||
this_matrix[gin][gmin_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
int gmax_;
|
||||
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
|
||||
bool non_zero = std::all_of(this_matrix[gin][gmax_].begin(),
|
||||
this_matrix[gin][gmax_].end(),
|
||||
[](double val){return val > 0.;});
|
||||
[](double val){return val != 0.;});
|
||||
if (non_zero) break;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -41,7 +41,7 @@ class ScattData {
|
|||
double get_xs(const char* xstype, int gin, int* gout, double* mu);
|
||||
void generic_init(int order, int_1dvec in_gmin, int_1dvec in_gmax,
|
||||
double_2dvec in_energy, double_2dvec in_mult);
|
||||
virtual void combine(std::vector<ScattData*> those_scatts,
|
||||
virtual void combine(std::vector<ScattData*>& those_scatts,
|
||||
double_1dvec& scalars) = 0;
|
||||
virtual int get_order() = 0;
|
||||
virtual double_3dvec get_matrix(int max_order) = 0;
|
||||
|
|
@ -59,7 +59,7 @@ class ScattDataLegendre: public ScattData {
|
|||
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);
|
||||
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
|
||||
int get_order() {return dist[0][0].size() - 1;};
|
||||
double_3dvec get_matrix(int max_order);
|
||||
};
|
||||
|
|
@ -74,7 +74,7 @@ class ScattDataHistogram: public ScattData {
|
|||
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);
|
||||
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
|
||||
int get_order() {return dist[0][0].size();};
|
||||
double_3dvec get_matrix(int max_order);
|
||||
};
|
||||
|
|
@ -91,7 +91,7 @@ class ScattDataTabular: public ScattData {
|
|||
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);
|
||||
void combine(std::vector<ScattData*>& those_scatts, double_1dvec& scalars);
|
||||
int get_order() {return dist[0][0].size();};
|
||||
double_3dvec get_matrix(int max_order);
|
||||
};
|
||||
|
|
|
|||
|
|
@ -1,5 +1,7 @@
|
|||
module tracking
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants
|
||||
use error, only: warning, write_message
|
||||
use geometry_header, only: cells
|
||||
|
|
@ -27,6 +29,21 @@ module tracking
|
|||
|
||||
implicit none
|
||||
|
||||
interface
|
||||
subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
|
||||
nu_fiss_xs) bind(C, name='calculate_xs')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: i_mat
|
||||
integer(C_INT), value, intent(in) :: gin
|
||||
real(C_DOUBLE), value, intent(in) :: sqrtkT
|
||||
real(C_DOUBLE), intent(in) :: uvw(1:3)
|
||||
real(C_DOUBLE), intent(inout) :: total_xs
|
||||
real(C_DOUBLE), intent(inout) :: abs_xs
|
||||
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
|
||||
end subroutine calculate_xs_c
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -112,8 +129,14 @@ contains
|
|||
end if
|
||||
else
|
||||
! Get the MG data
|
||||
!!TODO: Remove Fortran call - needed until I'm done replacing Fortran
|
||||
!!with C++ code because it sets index_temp
|
||||
call macro_xs(p % material) % obj % calculate_xs(p % g, p % sqrtkT, &
|
||||
p % coord(p % n_coord) % uvw, material_xs)
|
||||
call calculate_xs_c(p % material, p % g, p % sqrtkT, &
|
||||
p % coord(p % n_coord) % uvw, material_xs % total, &
|
||||
material_xs % absorption, material_xs % nu_fission)
|
||||
|
||||
|
||||
! Finally, update the particle group while we have already checked
|
||||
! for if multi-group
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue