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

It works! Replaced all mgxs_header functionality with the C++ version and a mgxs_interface module to act as the go-between the C++ and Fortran

Browse source
This commit is contained in:
Adam G Nelson 2018-06-13 20:22:51 -04:00
parent 6c8de73e3e
commit b1c73918a8
29 changed files with 1111 additions and 674 deletions
Download patch file Download diff file Expand all files Collapse all files

5
CMakeLists.txt
View file

@ -363,7 +363,7 @@ set(LIBOPENMC_FORTRAN_SRC
src/mesh_header.F90
src/message_passing.F90
src/mgxs_data.F90
src/mgxs_header.F90
src/mgxs_interface.F90
src/multipole_header.F90
src/nuclide_header.F90
src/output.F90
@ -380,7 +380,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/reaction_header.F90
src/relaxng
src/sab_header.F90
src/scattdata_header.F90
src/secondary_correlated.F90
src/secondary_kalbach.F90
src/secondary_nbody.F90
@ -439,11 +438,13 @@ set(LIBOPENMC_CXX_SRC
src/math_functions.cpp
src/message_passing.cpp
src/mgxs.cpp
src/mgxs_interface.cpp
src/plot.cpp
src/random_lcg.cpp
src/scattdata.cpp
src/simulation.cpp
src/state_point.cpp
src/string_functions.cpp
src/surface.cpp
src/xml_interface.cpp
src/xsdata.cpp

2
src/api.F90
View file

@ -309,7 +309,6 @@ contains
subroutine free_memory()
use cmfd_header
use mgxs_header
use plot_header
use sab_header
use settings
@ -327,7 +326,6 @@ contains
call free_memory_simulation()
call free_memory_nuclide()
call free_memory_settings()
call free_memory_mgxs()
call free_memory_sab()
call free_memory_source()
call free_memory_mesh()

4
src/cmfd_input.F90
View file

@ -3,8 +3,8 @@ module cmfd_input
use, intrinsic :: ISO_C_BINDING
use cmfd_header
use mesh_header, only: mesh_dict
use mgxs_header, only: energy_bins
use mesh_header, only: mesh_dict
use mgxs_interface, only: energy_bins, num_energy_groups
use tally
use tally_header
use timer_header

19
src/constants.F90
View file

@ -402,6 +402,25 @@ module constants
DIFF_NUCLIDE_DENSITY = 2, &
DIFF_TEMPERATURE = 3
! Mgxs::get_xs enumerated types
integer(C_INT), parameter :: &
MG_GET_XS_TOTAL = 0, &
MG_GET_XS_ABSORPTION = 1, &
MG_GET_XS_INVERSE_VELOCITY = 2, &
MG_GET_XS_DECAY_RATE = 3, &
MG_GET_XS_SCATTER = 4, &
MG_GET_XS_SCATTER_MULT = 5, &
MG_GET_XS_SCATTER_FMU_MULT = 6, &
MG_GET_XS_SCATTER_FMU = 7, &
MG_GET_XS_FISSION = 8, &
MG_GET_XS_KAPPA_FISSION = 9, &
MG_GET_XS_PROMPT_NU_FISSION = 10, &
MG_GET_XS_DELAYED_NU_FISSION = 11, &
MG_GET_XS_NU_FISSION = 12, &
MG_GET_XS_CHI_PROMPT = 13, &
MG_GET_XS_CHI_DELAYED = 14
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS

17
src/constants.h
View file

@ -63,6 +63,23 @@ constexpr double PI {3.1415926535898};
const double SQRT_PI {std::sqrt(PI)};
// Mgxs::get_xs enumerated types
constexpr int MG_GET_XS_TOTAL {0};
constexpr int MG_GET_XS_ABSORPTION {1};
constexpr int MG_GET_XS_INVERSE_VELOCITY {2};
constexpr int MG_GET_XS_DECAY_RATE {3};
constexpr int MG_GET_XS_SCATTER {4};
constexpr int MG_GET_XS_SCATTER_MULT {5};
constexpr int MG_GET_XS_SCATTER_FMU_MULT {6};
constexpr int MG_GET_XS_SCATTER_FMU {7};
constexpr int MG_GET_XS_FISSION {8};
constexpr int MG_GET_XS_KAPPA_FISSION {9};
constexpr int MG_GET_XS_PROMPT_NU_FISSION {10};
constexpr int MG_GET_XS_DELAYED_NU_FISSION {11};
constexpr int MG_GET_XS_NU_FISSION {12};
constexpr int MG_GET_XS_CHI_PROMPT {13};
constexpr int MG_GET_XS_CHI_DELAYED {14};
} // namespace openmc

10
src/input_xml.F90
View file

@ -19,8 +19,8 @@ module input_xml
use material_header
use mesh_header
use message_passing
use mgxs_data, only: create_macro_xs, read_mgxs, read_mgxs2, create_macro_xs2
use mgxs_header
use mgxs_data, only: create_macro_xs, read_mgxs
use mgxs_interface
use nuclide_header
use output, only: title, header, print_plot
use plot_header
@ -84,9 +84,7 @@ contains
else
! Create material macroscopic data for MGXS
call read_mgxs()
call read_mgxs2()
call create_macro_xs()
call create_macro_xs2()
end if
call time_read_xs % stop()
end if
@ -3856,7 +3854,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
awr = get_awr_c(mat % nuclide(j))
end if
! if given weight percent, convert all values so that they are divided
@ -3881,7 +3879,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
awr = get_awr_c(mat % nuclide(j))
end if
x = mat % atom_density(j)
sum_percent = sum_percent + x*awr

265
src/mgxs.cpp
View file

@ -2,6 +2,11 @@
namespace openmc {
// Storage for the MGXS data
std::vector<Mgxs> nuclides_MG;
std::vector<Mgxs> macro_xs;
//==============================================================================
// Mgxs base-class methods
//==============================================================================
@ -390,112 +395,149 @@ void Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
}
double Mgxs::get_xs(const char* xstype, const int gin, int* gout, double* mu,
double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
int* dg)
{
// This method assumes that the temperature and angle indices are set
double val;
if (std::strcmp(xstype, "total")) {
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs[index_temp].total[index_pol][index_azi][gin];
} else if (std::strcmp(xstype, "absorption")) {
break;
case MG_GET_XS_ABSORPTION:
val = xs[index_temp].absorption[index_pol][index_azi][gin];
} else if (std::strcmp(xstype, "inverse-velocity")) {
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs[index_temp].inverse_velocity[index_pol][index_azi][gin];
} else if (std::strcmp(xstype, "decay rate")) {
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs[index_temp].decay_rate[index_pol][index_azi][*dg + 1];
} else {
val = xs[index_temp].decay_rate[index_pol][index_azi][0];
}
} else if ((std::strcmp(xstype, "scatter")) ||
(std::strcmp(xstype, "scatter/mult")) ||
(std::strcmp(xstype, "scatter*f_mu/mult")) ||
(std::strcmp(xstype, "scatter*f_mu"))) {
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs[index_temp].scatter[index_pol]
[index_azi]->get_xs(xstype, gin, gout, mu);
} else if (fissionable && std::strcmp(xstype, "fission")) {
val = xs[index_temp].fission[index_pol][index_azi][gin];
} else if (fissionable && std::strcmp(xstype, "kappa-fission")) {
val = xs[index_temp].kappa_fission[index_pol][index_azi][gin];
} else if (fissionable && std::strcmp(xstype, "prompt-nu-fission")) {
val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
} else if (fissionable && std::strcmp(xstype, "delayed-nu-fission")) {
if (dg != nullptr) {
val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg];
break;
case MG_GET_XS_FISSION:
if (fissionable) {
val = xs[index_temp].fission[index_pol][index_azi][gin];
} else {
val = 0.;
for (auto& num : xs[index_temp].delayed_nu_fission[index_pol]
[index_azi][gin]) {
val += num;
}
}
} else if (fissionable && std::strcmp(xstype, "nu-fission")) {
val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
for (auto& num : xs[index_temp].delayed_nu_fission[index_pol]
[index_azi][gin]) {
val += num;
}
} else if (fissionable && std::strcmp(xstype, "chi-prompt")) {
if (gout != nullptr) {
val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout];
break;
case MG_GET_XS_KAPPA_FISSION:
if (fissionable) {
val = xs[index_temp].kappa_fission[index_pol][index_azi][gin];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) {
val += num;
}
}
} else if (fissionable && std::strcmp(xstype, "chi-delayed")) {
if (gout != nullptr) {
break;
case MG_GET_XS_PROMPT_NU_FISSION:
if (fissionable) {
val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg];
val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg];
} else {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0];
val = 0.;
for (auto& num : xs[index_temp].delayed_nu_fission[index_pol]
[index_azi][gin]) {
val += num;
}
}
} else {
if (dg != nullptr) {
val = 0.;
}
break;
case MG_GET_XS_NU_FISSION:
if (fissionable) {
val = xs[index_temp].nu_fission[index_pol][index_azi][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
val += xs[index_temp].chi_delayed[index_pol][index_azi][gin][i][*dg];
for (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) {
val += num;
}
}
} else {
val = 0.;
}
break;
case MG_GET_XS_CHI_DELAYED:
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg];
} else {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0];
}
} else {
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
for (auto& num : xs[index_temp].chi_delayed[index_pol]
[index_azi][gin][i]) {
val += num;
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
val += xs[index_temp].chi_delayed[index_pol][index_azi][gin][i][*dg];
}
} else {
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
for (auto& num : xs[index_temp].chi_delayed[index_pol]
[index_azi][gin][i]) {
val += num;
}
}
}
}
} else {
val = 0.;
}
} else {
break;
default:
val = 0.;
}
return val;
}
void Mgxs::sample_fission_energy(const int gin, const double nu_fission,
int& dg, int& gout)
void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
{
// This method assumes that the temperature and angle indices are set
double nu_fission = xs[index_temp].nu_fission[index_pol][index_azi][gin];
// Find the probability of having a prompt neutron
double prob_prompt =
xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin] /
nu_fission;
xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
// sample random numbers
double xi_pd = prn();
double xi_pd = prn() * nu_fission;
double xi_gout = prn();
// Select whether the neutron is prompt or delayed
if (xi_pd <= prob_prompt) {
// the neutron is prompt
// set the delayed group for the particle to be 0, indicating prompt
dg = 0;
// set the delayed group for the particle to be -1, indicating prompt
dg = -1;
// sample the outgoing energy group
gout = 0;
@ -514,12 +556,11 @@ void Mgxs::sample_fission_energy(const int gin, const double nu_fission,
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg] /
nu_fission;
xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg];
}
// adjust dg in case of round-off error
dg = std::min(dg, num_delayed_groups);
dg = std::min(dg, num_delayed_groups - 1);
// sample the outgoing energy group
gout = 0;
@ -552,11 +593,7 @@ void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin];
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;
}
nu_fiss_xs = xs[index_temp].nu_fission[index_pol][index_azi][gin];
} else {
nu_fiss_xs = 0.;
}
@ -580,7 +617,7 @@ bool Mgxs::equiv(const Mgxs& that)
}
inline void Mgxs::set_temperature_index(const double sqrtkT)
void Mgxs::set_temperature_index(const double sqrtkT)
{
// See if we need to find the new index
if (sqrtkT != last_sqrtkT) {
@ -600,7 +637,7 @@ inline void Mgxs::set_temperature_index(const double sqrtkT)
}
inline void Mgxs::set_angle_index(const double uvw[3])
void Mgxs::set_angle_index(const double uvw[3])
{
// See if we need to find the new index
if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) ||
@ -622,98 +659,4 @@ inline void Mgxs::set_angle_index(const double uvw[3])
}
}
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
void add_mgxs(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points)
{
//!! mgxs_data.F90 will be modified to just create the list of names
//!! in the order needed
// Convert temps to a vector for the from_hdf5 function
double_1dvec temperature;
temperature.assign(temps, temps + n_temps);
// TODO: C++ replacement for write_message
// write_message("Loading " + std::string(names[i]) + " data...", 6);
// Check to make sure cross section set exists in the library
hid_t xs_grp;
if (object_exists(file_id, name)) {
xs_grp = open_group(file_id, name);
} else {
fatal_error("Data for " + std::string(name) + " does not exist in "
+ "provided MGXS Library");
}
Mgxs mg;
mg.from_hdf5(xs_grp, energy_groups, delayed_groups,
temperature, method, tolerance, max_order, legendre_to_tabular,
legendre_to_tabular_points);
nuclides_MG.push_back(mg);
}
bool query_fissionable(const int n_nuclides, const int i_nuclides[])
{
bool result = false;
for (int n = 0; n < n_nuclides; n++) {
if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
}
return result;
}
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)
{
Mgxs macro;
if (n_temps > 0) {
// // Convert temps to a vector
double_1dvec temperature;
temperature.assign(temps, temps + n_temps);
// Convert atom_densities to a vector
double_1dvec atom_densities_vec;
atom_densities_vec.assign(atom_densities, atom_densities + n_nuclides);
// Build array of pointers to nuclides_MG's Mgxs objects needed for this
// material
std::vector<Mgxs*> mgxs_ptr(n_nuclides);
for (int n = 0; n < n_nuclides; n++) {
mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
}
macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
method, tolerance);
}
macro_xs.push_back(macro);
}
//==============================================================================
// 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

47
src/mgxs.h
View file

@ -29,29 +29,31 @@ namespace openmc {
class Mgxs {
private:
std::string name; // name of dataset, e.g., UO2
double awr; // atomic weight ratio
double_1dvec 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
int index_temp; // cache of temperature index
double last_sqrtkT; // cache of the temperature corresponding to index_temp
std::vector<XsData> xs; // Cross section data
int n_pol;
int n_azi;
int index_pol; // cache fof the angle indices
int index_azi;
double_1dvec polar;
double_1dvec azimuthal;
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,
int& order_dim, bool& is_isotropic);
bool equiv(const Mgxs& that);
public:
std::string name; // name of dataset, e.g., UO2
double awr; // atomic weight ratio
bool fissionable; // Is this fissionable
// TODO: The following attributes be private when Fortran is fully replaced
int index_pol; // cache for the angle indices
int index_azi;
double last_uvw[3]; // cache of the angle corresponding to the above indices
int index_temp; // cache of temperature index
void init(const std::string& in_name, const double in_awr,
const double_1dvec& in_kTs, const bool in_fissionable,
const int in_scatter_format, const int in_num_groups,
@ -66,40 +68,15 @@ 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, const int gin, int* gout, double* mu,
double get_xs(const int xstype, const int gin, int* gout, double* mu,
int* dg);
void sample_fission_energy(const int gin, const double nu_fission,
int& dg, int& gout);
void sample_fission_energy(const int gin, 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(const double sqrtkT);
inline void set_angle_index(const double uvw[3]);
void set_temperature_index(const double sqrtkT);
void set_angle_index(const double uvw[3]);
};
extern "C" void add_mgxs(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points);
extern "C" bool query_fissionable(const int n_nuclides, const int i_nuclides[]);
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;
std::vector<Mgxs> macro_xs;
} // namespace openmc
#endif // MGXS_H

414
src/mgxs_data.F90
View file

@ -9,7 +9,7 @@ module mgxs_data
use geometry_header, only: get_temperatures, cells
use hdf5_interface
use material_header, only: Material, materials, n_materials
use mgxs_header
use mgxs_interface
use nuclide_header, only: n_nuclides
use set_header, only: SetChar
use settings
@ -17,50 +17,6 @@ module mgxs_data
use string, only: to_lower
implicit none
interface
subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
n_temps, temps, method, tolerance, max_order, legendre_to_tabular, &
legendre_to_tabular_points) bind(C, name='add_mgxs')
use ISO_C_BINDING
import HID_T
implicit none
integer(HID_T), value, intent(in) :: file_id
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: energy_groups
integer(C_INT), value, intent(in) :: delayed_groups
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), value, intent(in) :: max_order
logical(C_BOOL),value, intent(in) :: legendre_to_tabular
integer(C_INT), value, intent(in) :: legendre_to_tabular_points
end subroutine add_mgxs_c
function query_fissionable_c(n_nuclides, i_nuclides) result(result) &
bind(C, name='query_fissionable')
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
logical(C_BOOL) :: result
end function query_fissionable_c
subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
atom_densities, method, tolerance) bind(C, name='create_macro_xs')
use ISO_C_BINDING
implicit none
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
real(C_DOUBLE), intent(in) :: atom_densities(1:n_nuclides)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
end subroutine create_macro_xs_c
end interface
contains
!===============================================================================
@ -68,150 +24,150 @@ contains
! nuclides and sab_tables arrays
!===============================================================================
! subroutine read_mgxs()
! integer :: i ! index in materials array
! integer :: j ! index over nuclides in material
! integer :: i_nuclide ! index in nuclides array
! character(20) :: name ! name of library to load
! integer :: representation ! Data representation
! character(MAX_LINE_LEN) :: temp_str
! type(Material), pointer :: mat
! type(SetChar) :: already_read
! integer(HID_T) :: file_id
! integer(HID_T) :: xsdata_group
! logical :: file_exists
! type(VectorReal), allocatable :: temps(:)
! character(MAX_WORD_LEN) :: word
! integer, allocatable :: array(:)
! ! Check if MGXS Library exists
! inquire(FILE=path_cross_sections, EXIST=file_exists)
! if (.not. file_exists) then
! ! Could not find MGXS Library file
! call fatal_error("Cross sections HDF5 file '" &
! // trim(path_cross_sections) // "' does not exist!")
! end if
! call write_message("Loading cross section data...", 5)
! ! Get temperatures
! call get_temperatures(temps)
! ! Open file for reading
! file_id = file_open(path_cross_sections, 'r', parallel=.true.)
! ! Read filetype
! call read_attribute(word, file_id, "filetype")
! if (word /= 'mgxs') then
! call fatal_error("Provided MGXS Library is not a MGXS Library file.")
! end if
! ! Read revision number for the MGXS Library file and make sure it matches
! ! with the current version
! call read_attribute(array, file_id, "version")
! if (any(array /= VERSION_MGXS_LIBRARY)) then
! call fatal_error("MGXS Library file version does not match current &
! &version supported by OpenMC.")
! end if
! ! allocate arrays for MGXS storage and cross section cache
! allocate(nuclides_MG(n_nuclides))
! ! ==========================================================================
! ! READ ALL MGXS CROSS SECTION TABLES
! ! Loop over all files
! MATERIAL_LOOP: do i = 1, n_materials
! mat => materials(i)
! NUCLIDE_LOOP: do j = 1, mat % n_nuclides
! name = mat % names(j)
! if (.not. already_read % contains(name)) then
! i_nuclide = mat % nuclide(j)
! call write_message("Loading " // trim(name) // " data...", 6)
! ! Check to make sure cross section set exists in the library
! if (object_exists(file_id, trim(name))) then
! xsdata_group = open_group(file_id, trim(name))
! else
! call fatal_error("Data for '" // trim(name) // "' does not exist in "&
! &// trim(path_cross_sections))
! end if
! ! First find out the data representation
! if (attribute_exists(xsdata_group, "representation")) then
! call read_attribute(temp_str, xsdata_group, "representation")
! if (trim(temp_str) == 'isotropic') then
! representation = MGXS_ISOTROPIC
! else if (trim(temp_str) == 'angle') then
! representation = MGXS_ANGLE
! else
! call fatal_error("Invalid Data Representation!")
! end if
! else
! ! Default to isotropic representation
! representation = MGXS_ISOTROPIC
! end if
! ! Now allocate accordingly
! select case(representation)
! case(MGXS_ISOTROPIC)
! allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj)
! case(MGXS_ANGLE)
! allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj)
! end select
! ! Now read in the data specific to the type we just declared
! call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
! num_energy_groups, num_delayed_groups, temps(i_nuclide), &
! temperature_method, temperature_tolerance, max_order, &
! legendre_to_tabular, legendre_to_tabular_points)
! ! Add name to dictionary
! call already_read % add(name)
! call close_group(xsdata_group)
! end if
! end do NUCLIDE_LOOP
! end do MATERIAL_LOOP
! ! Avoid some valgrind leak errors
! call already_read % clear()
! ! Loop around material
! MATERIAL_LOOP3: do i = 1, n_materials
! ! Get material
! mat => materials(i)
! ! Loop around nuclides in material
! NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! ! Is this fissionable?
! if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
! mat % fissionable = .true.
! end if
! if (mat % fissionable) then
! exit NUCLIDE_LOOP2
! end if
! end do NUCLIDE_LOOP2
! end do MATERIAL_LOOP3
! call file_close(file_id)
! end subroutine read_mgxs
subroutine read_mgxs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
character(MAX_LINE_LEN) :: temp_str
type(Material), pointer :: mat
type(SetChar) :: already_read
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
type(VectorReal), allocatable :: temps(:)
character(MAX_WORD_LEN) :: word
integer, allocatable :: array(:)
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find MGXS Library file
call fatal_error("Cross sections HDF5 file '" &
// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Loading cross section data...", 5)
! Get temperatures
call get_temperatures(temps)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
! Read filetype
call read_attribute(word, file_id, "filetype")
if (word /= 'mgxs') then
call fatal_error("Provided MGXS Library is not a MGXS Library file.")
end if
! Read revision number for the MGXS Library file and make sure it matches
! with the current version
call read_attribute(array, file_id, "version")
if (any(array /= VERSION_MGXS_LIBRARY)) then
call fatal_error("MGXS Library file version does not match current &
&version supported by OpenMC.")
end if
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides))
! ==========================================================================
! READ ALL MGXS CROSS SECTION TABLES
! Loop over all files
MATERIAL_LOOP: do i = 1, n_materials
mat => materials(i)
NUCLIDE_LOOP: do j = 1, mat % n_nuclides
name = mat % names(j)
if (.not. already_read % contains(name)) then
i_nuclide = mat % nuclide(j)
call write_message("Loading " // trim(name) // " data...", 6)
! Check to make sure cross section set exists in the library
if (object_exists(file_id, trim(name))) then
xsdata_group = open_group(file_id, trim(name))
else
call fatal_error("Data for '" // trim(name) // "' does not exist in "&
&// trim(path_cross_sections))
end if
! First find out the data representation
if (attribute_exists(xsdata_group, "representation")) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
else if (trim(temp_str) == 'angle') then
representation = MGXS_ANGLE
else
call fatal_error("Invalid Data Representation!")
end if
else
! Default to isotropic representation
representation = MGXS_ISOTROPIC
end if
! Now allocate accordingly
select case(representation)
case(MGXS_ISOTROPIC)
allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj)
case(MGXS_ANGLE)
allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj)
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
num_energy_groups, num_delayed_groups, temps(i_nuclide), &
temperature_method, temperature_tolerance, max_order, &
legendre_to_tabular, legendre_to_tabular_points)
! Add name to dictionary
call already_read % add(name)
call close_group(xsdata_group)
end if
end do NUCLIDE_LOOP
end do MATERIAL_LOOP
! Avoid some valgrind leak errors
call already_read % clear()
! Loop around material
MATERIAL_LOOP3: do i = 1, n_materials
! Get material
mat => materials(i)
! Loop around nuclides in material
NUCLIDE_LOOP2: do j = 1, mat % n_nuclides
! Is this fissionable?
if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then
mat % fissionable = .true.
end if
if (mat % fissionable) then
exit NUCLIDE_LOOP2
end if
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
call file_close(file_id)
end subroutine read_mgxs
subroutine read_mgxs2()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_nuclide ! index in nuclides array
@ -286,55 +242,55 @@ contains
! Avoid some valgrind leak errors
call already_read % clear()
end subroutine read_mgxs2
end subroutine read_mgxs
!===============================================================================
! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data
!===============================================================================
! subroutine create_macro_xs()
! integer :: i_mat ! index in materials array
! type(Material), pointer :: mat ! current material
! type(VectorReal), allocatable :: kTs(:)
! allocate(macro_xs(n_materials))
! ! Get temperatures to read for each material
! call get_mat_kTs(kTs)
! ! Force all nuclides in a material to be the same representation.
! ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs).
! ! At the same time, we will find the scattering type, as that will dictate
! ! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials))
! do i_mat = 1, n_materials
! ! Get the material
! mat => materials(i_mat)
! ! Get the scattering type for the first nuclide
! select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
! type is (MgxsIso)
! allocate(MgxsIso :: macro_xs(i_mat) % obj)
! type is (MgxsAngle)
! allocate(MgxsAngle :: macro_xs(i_mat) % obj)
! end select
! ! Do not read materials which we do not actually use in the problem to
! ! reduce storage
! if (allocated(kTs(i_mat) % data)) then
! call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
! num_energy_groups, num_delayed_groups, max_order, &
! temperature_tolerance, temperature_method)
! end if
! end do
! end subroutine create_macro_xs
subroutine create_macro_xs()
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
allocate(macro_xs(n_materials))
! Get temperatures to read for each material
call get_mat_kTs(kTs)
! Force all nuclides in a material to be the same representation.
! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs).
! At the same time, we will find the scattering type, as that will dictate
! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials))
do i_mat = 1, n_materials
! Get the material
mat => materials(i_mat)
! Get the scattering type for the first nuclide
select type(nuc => nuclides_MG(mat % nuclide(1)) % obj)
type is (MgxsIso)
allocate(MgxsIso :: macro_xs(i_mat) % obj)
type is (MgxsAngle)
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
! Do not read materials which we do not actually use in the problem to
! reduce storage
if (allocated(kTs(i_mat) % data)) then
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
num_energy_groups, num_delayed_groups, max_order, &
temperature_tolerance, temperature_method)
end if
end do
end subroutine create_macro_xs
subroutine create_macro_xs2()
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
character(MAX_WORD_LEN) :: name ! name of material
! Get temperatures to read for each material
@ -360,7 +316,7 @@ contains
end if
end do
end subroutine create_macro_xs2
end subroutine create_macro_xs
!===============================================================================

1
src/mgxs_header.F90
View file

@ -11,7 +11,6 @@ module mgxs_header
use math, only: evaluate_legendre
use nuclide_header, only: MaterialMacroXS
use random_lcg, only: prn
use scattdata_header
use string
use stl_vector, only: VectorInt, VectorReal

201
src/mgxs_interface.F90 Normal file
View file

@ -0,0 +1,201 @@
module mgxs_interface
use, intrinsic :: ISO_C_BINDING
use hdf5_interface
implicit none
interface
subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
n_temps, temps, method, tolerance, max_order, legendre_to_tabular, &
legendre_to_tabular_points) bind(C)
use ISO_C_BINDING
import HID_T
implicit none
integer(HID_T), value, intent(in) :: file_id
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: energy_groups
integer(C_INT), value, intent(in) :: delayed_groups
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), value, intent(in) :: max_order
logical(C_BOOL),value, intent(in) :: legendre_to_tabular
integer(C_INT), value, intent(in) :: legendre_to_tabular_points
end subroutine add_mgxs_c
function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
logical(C_BOOL) :: result
end function query_fissionable_c
subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
atom_densities, method, tolerance) bind(C)
use ISO_C_BINDING
implicit none
character(kind=C_CHAR),intent(in) :: name(*)
integer(C_INT), value, intent(in) :: n_nuclides
integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides)
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
real(C_DOUBLE), intent(in) :: atom_densities(1:n_nuclides)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
end subroutine create_macro_xs_c
subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
nu_fiss_xs) bind(C)
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
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
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 sample_scatter_c
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
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) :: dg
integer(C_INT), intent(inout) :: gout
end subroutine sample_fission_energy_c
subroutine get_name_c(index, name_len, name) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: name_len
character(kind=C_CHAR), intent(inout) :: name(name_len)
end subroutine get_name_c
function get_awr_c(index) result(awr) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE) :: awr
end function get_awr_c
function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_nuclide_xs_c
function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_macro_xs_c
subroutine set_nuclide_angle_index_c(index, uvw, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
integer(C_INT), intent(inout) :: last_pol
integer(C_INT), intent(inout) :: last_azi
real(C_DOUBLE), intent(inout) :: last_uvw(1:3)
end subroutine set_nuclide_angle_index_c
subroutine reset_nuclide_angle_index_c(index, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_pol
integer(C_INT), value, intent(in) :: last_azi
real(C_DOUBLE), intent(in) :: last_uvw(1:3)
end subroutine reset_nuclide_angle_index_c
subroutine set_macro_angle_index_c(index, uvw, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
integer(C_INT), intent(inout) :: last_pol
integer(C_INT), intent(inout) :: last_azi
real(C_DOUBLE), intent(inout) :: last_uvw(1:3)
end subroutine set_macro_angle_index_c
subroutine reset_macro_angle_index_c(index, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_pol
integer(C_INT), value, intent(in) :: last_azi
real(C_DOUBLE), intent(in) :: last_uvw(1:3)
end subroutine reset_macro_angle_index_c
function set_nuclide_temperature_index_c(index, sqrtkT) result(last_temp) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: sqrtkT
integer(C_INT) :: last_temp
end function set_nuclide_temperature_index_c
subroutine reset_nuclide_temperature_index_c(index, last_temp) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_temp
end subroutine reset_nuclide_temperature_index_c
end interface
! Number of energy groups
integer(C_INT) :: num_energy_groups
! Number of delayed groups
integer(C_INT) :: num_delayed_groups
! Energy group structure with decreasing energy
real(8), allocatable :: energy_bins(:)
! Midpoint of the energy group structure
real(8), allocatable :: energy_bin_avg(:)
! Energy group structure with increasing energy
real(8), allocatable :: rev_energy_bins(:)
end module mgxs_interface

247
src/mgxs_interface.cpp Normal file
View file

@ -0,0 +1,247 @@
#include "mgxs_interface.h"
namespace openmc {
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points)
{
//!! mgxs_data.F90 will be modified to just create the list of names
//!! in the order needed
// Convert temps to a vector for the from_hdf5 function
double_1dvec temperature;
temperature.assign(temps, temps + n_temps);
// TODO: C++ replacement for write_message
// write_message("Loading " + std::string(names[i]) + " data...", 6);
// Check to make sure cross section set exists in the library
hid_t xs_grp;
if (object_exists(file_id, name)) {
xs_grp = open_group(file_id, name);
} else {
fatal_error("Data for " + std::string(name) + " does not exist in "
+ "provided MGXS Library");
}
Mgxs mg;
mg.from_hdf5(xs_grp, energy_groups, delayed_groups,
temperature, method, tolerance, max_order, legendre_to_tabular,
legendre_to_tabular_points);
nuclides_MG.push_back(mg);
}
bool query_fissionable_c(const int n_nuclides, const int i_nuclides[])
{
bool result = false;
for (int n = 0; n < n_nuclides; n++) {
if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
}
return result;
}
void create_macro_xs_c(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)
{
Mgxs macro;
if (n_temps > 0) {
// // Convert temps to a vector
double_1dvec temperature;
temperature.assign(temps, temps + n_temps);
// Convert atom_densities to a vector
double_1dvec atom_densities_vec;
atom_densities_vec.assign(atom_densities, atom_densities + n_nuclides);
// Build array of pointers to nuclides_MG's Mgxs objects needed for this
// material
std::vector<Mgxs*> mgxs_ptr(n_nuclides);
for (int n = 0; n < n_nuclides; n++) {
mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
}
macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
method, tolerance);
}
macro_xs.push_back(macro);
}
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
void calculate_xs_c(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 sample_scatter_c(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);
// adjust return value for fortran indexing
gout = gout_c + 1;
// Rotate the angle
rotate_angle_c(uvw, mu, nullptr);
}
void sample_fission_energy_c(const int i_mat, const int gin, int& dg, int& gout)
{
int dg_c = 0;
int gout_c = 0;
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
// adjust return values for fortran indexing
dg = dg_c + 1;
gout = gout_c + 1;
}
void get_name_c(const int index, int name_len, char* name)
{
// First blank out our input string
std::string str(name_len, ' ');
std::strcpy(name, str.c_str());
// Now get the data and copy to the C-string
str = nuclides_MG[index - 1].name;
std::strcpy(name, str.c_str());
// Finally, remove the null terminator
name[std::strlen(name)] = ' ';
}
double get_awr_c(const int index)
{
return nuclides_MG[index - 1].awr;
}
double get_nuclide_xs_c(const int index, const int xstype, const int gin,
int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
double get_macro_xs_c(const int index, const int xstype, const int gin,
int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
void set_nuclide_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3])
{
// Store the old
last_pol = nuclides_MG[index - 1].index_pol;
last_azi = nuclides_MG[index - 1].index_azi;
last_uvw[0] = nuclides_MG[index - 1].last_uvw[0];
last_uvw[1] = nuclides_MG[index - 1].last_uvw[1];
last_uvw[2] = nuclides_MG[index - 1].last_uvw[2];
// Update the values
nuclides_MG[index - 1].set_angle_index(uvw);
}
void reset_nuclide_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3])
{
nuclides_MG[index - 1].index_pol = last_pol;
nuclides_MG[index - 1].index_azi = last_azi;
nuclides_MG[index - 1].last_uvw[0] = last_uvw[0];
nuclides_MG[index - 1].last_uvw[1] = last_uvw[1];
nuclides_MG[index - 1].last_uvw[2] = last_uvw[2];
}
void set_macro_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3])
{
// Store the old
last_pol = macro_xs[index - 1].index_pol;
last_azi = macro_xs[index - 1].index_azi;
last_uvw[0] = macro_xs[index - 1].last_uvw[0];
last_uvw[1] = macro_xs[index - 1].last_uvw[1];
last_uvw[2] = macro_xs[index - 1].last_uvw[2];
// Update the values
macro_xs[index - 1].set_angle_index(uvw);
}
void reset_macro_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3])
{
macro_xs[index - 1].index_pol = last_pol;
macro_xs[index - 1].index_azi = last_azi;
macro_xs[index - 1].last_uvw[0] = last_uvw[0];
macro_xs[index - 1].last_uvw[1] = last_uvw[1];
macro_xs[index - 1].last_uvw[2] = last_uvw[2];
}
int set_nuclide_temperature_index_c(const int index, const double sqrtkT)
{
int old = nuclides_MG[index - 1].index_temp;
nuclides_MG[index - 1].set_temperature_index(sqrtkT);
return old;
}
void reset_nuclide_temperature_index_c(const int index, const int last_temp)
{
nuclides_MG[index - 1].index_temp = last_temp;
}
} // namespace openmc

66
src/mgxs_interface.h Normal file
View file

@ -0,0 +1,66 @@
//! \file mgxs_interface.h
//! A collection of C interfaces to the C++ Mgxs class
#ifndef MGXS_INTERFACE_H
#define MGXS_INTERFACE_H
#include "mgxs.h"
namespace openmc {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points);
extern "C" bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]);
extern "C" void create_macro_xs_c(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_c(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 sample_scatter_c(const int i_mat, const int gin, int& gout,
double& mu, double& wgt, double uvw[3]);
extern "C" void sample_fission_energy_c(const int i_mat, const int gin,
int& dg, int& gout);
extern "C" void get_name_c(const int index, int name_len, char* name);
extern "C" double get_awr_c(const int index);
extern "C" double get_nuclide_xs_c(const int index, const int xstype,
const int gin, int* gout, double* mu, int* dg);
extern "C" double get_macro_xs_c(const int index, const int xstype,
const int gin, int* gout, double* mu, int* dg);
extern "C" void set_nuclide_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3]);
extern "C" void reset_nuclide_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3]);
extern "C" void set_macro_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3]);
extern "C" void reset_macro_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3]);
extern "C" int set_nuclide_temperature_index_c(const int index,
const double sqrtkT);
extern "C" void reset_nuclide_temperature_index_c(const int index,
const int last_temp);
} // namespace openmc
#endif // MGXS_INTERFACE_H

6
src/output.F90
View file

@ -12,7 +12,7 @@ module output
use math, only: t_percentile
use mesh_header, only: RegularMesh, meshes
use message_passing, only: master, n_procs
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, Particle
use plot_header
@ -680,6 +680,7 @@ contains
character(36) :: score_name ! names of scoring function
! to be applied at write-time
type(TallyFilterMatch), allocatable :: matches(:)
character(MAX_WORD_LEN) :: temp_name
! Skip if there are no tallies
if (n_tallies == 0) return
@ -843,8 +844,9 @@ contains
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides(i_nuclide) % name)
else
call get_name_c(i_nuclide, len(temp_name), temp_name)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides_MG(i_nuclide) % obj % name)
trim(temp_name)
end if
end if

2
src/particle_restart.F90
View file

@ -6,7 +6,7 @@ module particle_restart
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_header, only: energy_bin_avg
use mgxs_interface, only: energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: print_particle
use particle_header, only: Particle

33
src/physics_mg.F90
View file

@ -8,13 +8,12 @@ module physics_mg
use material_header, only: Material, materials
use math, only: rotate_angle
use mesh_header, only: meshes
use mgxs_header
use mgxs_interface
use message_passing
use nuclide_header, only: material_xs
use particle_header, only: Particle
use physics_common
use random_lcg, only: prn
use scattdata_header
use settings
use simulation_header
use string, only: to_str
@ -22,20 +21,6 @@ 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
!===============================================================================
@ -157,14 +142,8 @@ 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)
! 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)
call sample_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)
@ -194,10 +173,6 @@ contains
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
class(Mgxs), pointer :: xs
! Get Pointers
xs => macro_xs(p % material) % obj
! TODO: Heat generation from fission
@ -277,7 +252,7 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
call xs % sample_fission_energy(p % g, bank_array(i) % uvw, dg, gout)
call sample_fission_energy_c(p % material, p % g, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg

46
src/scattdata.cpp
View file

@ -56,35 +56,38 @@ void ScattData::sample_energy(int gin, int& gout, int& i_gout)
}
double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu)
double ScattData::get_xs(const int xstype, int gin, int* gout, double* mu)
{
// Set the outgoing group offset index as needed
int i_gout = 0;
if (gout != nullptr) {
// short circuit the function if gout is from a zero portion of the
// scattering matrix
if ((*gout < gmin[gin]) || (*gout >= gmax[gin])) { // > gmax?
if ((*gout < gmin[gin]) || (*gout > gmax[gin])) { // > gmax?
return 0.;
}
i_gout = *gout - gmin[gin];
}
double val = 0.;
if (std::strcmp(xstype, "scatter")) {
switch(xstype) {
case MG_GET_XS_SCATTER:
if (gout != nullptr) {
val = scattxs[gin] * energy[gin][i_gout];
} else {
val = scattxs[gin];
}
} else if (std::strcmp(xstype, "scatter/mult")) {
break;
case MG_GET_XS_SCATTER_MULT:
if (gout != nullptr) {
val = scattxs[gin] * energy[gin][i_gout] / mult[gin][i_gout];
} else {
val = scattxs[gin] / std::inner_product(mult[gin].begin(),
mult[gin].end(),
energy[gin].begin(), 0.0);
val = scattxs[gin] /
std::inner_product(mult[gin].begin(), mult[gin].end(),
energy[gin].begin(), 0.0);
}
} else if (std::strcmp(xstype, "scatter*f_mu/mult")) {
break;
case MG_GET_XS_SCATTER_FMU_MULT:
if ((gout != nullptr) && (mu != nullptr)) {
val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu);
} else {
@ -92,7 +95,8 @@ double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu)
// group or mu is not useful
fatal_error("Invalid call to get_xs");
}
} else if (std::strcmp(xstype, "scatter*f_mu")) {
break;
case MG_GET_XS_SCATTER_FMU:
if ((gout != nullptr) && (mu != nullptr)) {
val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu) /
mult[gin][i_gout];
@ -101,6 +105,7 @@ double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu)
// group or mu is not useful
fatal_error("Invalid call to get_xs");
}
break;
}
return val;
}
@ -205,13 +210,11 @@ void ScattDataLegendre::update_max_val()
double ScattDataLegendre::calc_f(int gin, int gout, double mu)
{
// TODO: gout >= or gout >?
double f;
if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
// TODO: size() -1 or just size?
int i_gout = gout - gmin[gin]; //TODO: + 1?
int i_gout = gout - gmin[gin];
f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
dist[gin][i_gout].data(), mu);
}
@ -479,19 +482,18 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double ScattDataHistogram::calc_f(int gin, int gout, double mu)
{
// TODO: gout >= or gout >?
double f;
if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
// Find mu bin
int i_gout = gout - gmin[gin]; //TODO: + 1?
int i_gout = gout - gmin[gin];
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.);
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
f = fmu[gin][i_gout][imu];
@ -776,19 +778,18 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double ScattDataTabular::calc_f(int gin, int gout, double mu)
{
// TODO: gout >= or gout >?
double f;
if ((gout < gmin[gin]) || (gout >= gmax[gin])) {
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
f = 0.;
} else {
// Find mu bin
int i_gout = gout - gmin[gin]; //TODO: + 1?
int i_gout = gout - gmin[gin];
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.);
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
double r = (mu - this->mu[imu]) / (this->mu[imu + 1] - this->mu[imu]);
@ -1006,7 +1007,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg,
tab.dmu = 2. / (n_mu - 1);
tab.mu[0] = -1.;
for (int imu = 1; imu < n_mu - 1; imu++) {
tab.mu[imu] = -1. + (imu - 1) * tab.dmu;
tab.mu[imu] = -1. + imu * tab.dmu;
}
tab.mu[n_mu - 1] = 1.;
@ -1032,6 +1033,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg,
// Now re-normalize for numerical integration issues and to take care of
// the above negative fix-up. Also accrue the CDF
double norm = 0.;
tab.dist[gin][i_gout][0] = 0.;
for (int imu = 1; imu < n_mu; imu++) {
norm += 0.5 * tab.dmu * (tab.fmu[gin][i_gout][imu - 1] +
tab.fmu[gin][i_gout][imu]);

2
src/scattdata.h
View file

@ -38,7 +38,7 @@ class ScattData {
virtual void init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_mult, double_3dvec& coeffs) = 0;
void sample_energy(int gin, int& gout, int& i_gout);
double get_xs(const char* xstype, int gin, int* gout, double* mu);
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,

2
src/simulation.F90
View file

@ -20,7 +20,7 @@ module simulation
use geometry_header, only: n_cells
use material_header, only: n_materials, materials
use message_passing
use mgxs_header, only: energy_bins, energy_bin_avg
use mgxs_interface, only: energy_bins, energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: header, print_columns, &
print_batch_keff, print_generation, print_runtime, &

2
src/source.F90
View file

@ -14,7 +14,7 @@ module source
use hdf5_interface
use math
use message_passing, only: rank
use mgxs_header, only: rev_energy_bins, num_energy_groups
use mgxs_interface, only: rev_energy_bins, num_energy_groups
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream

11
src/state_point.F90
View file

@ -22,7 +22,7 @@ module state_point
use hdf5_interface
use mesh_header, only: RegularMesh, meshes, n_meshes
use message_passing
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header, only: nuclides
use output, only: time_stamp
use random_lcg, only: openmc_get_seed, openmc_set_seed
@ -73,6 +73,7 @@ contains
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
err = 0
if (present(filename)) then
@ -307,11 +308,13 @@ contains
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
else
i_xs = index(nuclides_MG(tally % nuclide_bins(j)) % obj % name, '.')
call get_name_c(tally % nuclide_bins(j), len(temp_name), &
temp_name)
i_xs = index(temp_name, '.')
if (i_xs > 0) then
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name(1 : i_xs-1)
str_array(j) = trim(temp_name(1 : i_xs-1))
else
str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name
str_array(j) = trim(temp_name)
end if
end if
else

30
src/string_functions.cpp Normal file
View file

@ -0,0 +1,30 @@
#include "string_functions.h"
namespace openmc {
std::string& strtrim(std::string& s)
{
const char* t = " \t\n\r\f\v";
s.erase(s.find_last_not_of(t) + 1);
s.erase(0, s.find_first_not_of(t));
return s;
}
char* strtrim(char* c_str)
{
std::string std_str;
std_str.assign(c_str);
strtrim(std_str);
int length = std_str.copy(c_str, std_str.size());
c_str[length] = '\0';
return c_str;
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
} // namespace openmc

31
src/string_functions.h
View file

@ -4,38 +4,15 @@
#ifndef STRING_FUNCTIONS_H
#define STRING_FUNCTIONS_H
// for string functions
#include <functional>
#include <cctype>
#include <locale>
#include <string>
namespace openmc {
std::string& strtrim(std::string& s)
{
const char* t = " \t\n\r\f\v";
s.erase(s.find_last_not_of(t) + 1);
s.erase(0, s.find_first_not_of(t));
return s;
}
std::string& strtrim(std::string& s);
char* strtrim(char* c_str);
char* strtrim(char* c_str)
{
std::string std_str;
std_str.assign(c_str);
strtrim(std_str);
int length = std_str.copy(c_str, std_str.size());
c_str[length] = '\0';
return c_str;
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
void to_lower(std::string& str);
} // namespace openmc
#endif // STRING_FUNCTIONS_H
#endif // STRING_FUNCTIONS_H

24
src/summary.F90
View file

@ -8,7 +8,7 @@ module summary
use material_header, only: Material, n_materials
use mesh_header, only: RegularMesh
use message_passing
use mgxs_header, only: nuclides_MG
use mgxs_interface
use nuclide_header
use output, only: time_stamp
use settings, only: run_CE
@ -94,7 +94,7 @@ contains
num_nuclides = 0
num_macros = 0
do i = 1, n_nuclides
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
if (get_awr_c(i) /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
@ -119,12 +119,14 @@ contains
nuc_names(i) = nuclides(i) % name
awrs(i) = nuclides(i) % awr
else
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(j) = nuclides_MG(i) % obj % name
awrs(j) = nuclides_MG(i) % obj % awr
if (get_awr_c(i) /= MACROSCOPIC_AWR) then
call get_name_c(i, len(nuc_names(j)), nuc_names(j))
nuc_names(j) = trim(nuc_names(j))
awrs(j) = get_awr_c(i)
j = j + 1
else
macro_names(k) = nuclides_MG(i) % obj % name
call get_name_c(i, len(macro_names(k)), macro_names(k))
macro_names(k) = trim(macro_names(k))
k = k + 1
end if
end if
@ -458,7 +460,7 @@ contains
num_nuclides = 0
num_macros = 0
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
@ -484,12 +486,14 @@ contains
k = 1
n = 1
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(k) = nuclides_MG(m % nuclide(j)) % obj % name
if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then
call get_name_c(m % nuclide(j), len(nuc_names(k)), nuc_names(k))
nuc_names(k) = trim(nuc_names(k))
nuc_densities(k) = m % atom_density(j)
k = k + 1
else
macro_names(n) = nuclides_MG(m % nuclide(j)) % obj % name
call get_name_c(m % nuclide(j), len(macro_names(n)), macro_names(n))
macro_names(n) = trim(macro_names(n))
n = n + 1
end if
end do

259
src/tallies/tally.F90
View file

@ -10,7 +10,7 @@ module tally
use math, only: t_percentile
use mesh_header, only: RegularMesh, meshes
use message_passing
use mgxs_header
use mgxs_interface
use nuclide_header
use output, only: header
use particle_header, only: LocalCoord, Particle
@ -1229,8 +1229,14 @@ contains
real(8) :: p_uvw(3) ! Particle's current uvw
integer :: p_g ! Particle group to use for getting info
! to tally with.
class(Mgxs), pointer :: matxs
class(Mgxs), pointer :: nucxs
! Storage of the indices the Mgxs object arrived with for resetting later
integer(C_INT) :: last_nuc_azi
integer(C_INT) :: last_nuc_pol
integer(C_INT) :: last_mat_azi
integer(C_INT) :: last_mat_pol
integer(C_INT) :: last_nuc_temp
real(C_DOUBLE) :: last_mat_uvw(3)
real(C_DOUBLE) :: last_nuc_uvw(3)
! Set the direction and group to use with get_xs
if (t % estimator == ESTIMATOR_ANALOG .or. &
@ -1268,13 +1274,15 @@ contains
! To significantly reduce de-referencing, point matxs to the
! macroscopic Mgxs for the material of interest
matxs => macro_xs(p % material) % obj
call set_macro_angle_index_c(p % material, p_uvw, last_mat_pol, &
last_mat_azi, last_mat_uvw)
! Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide > 0) then
nucxs => nuclides_MG(i_nuclide) % obj
! And since we haven't calculated this temperature index yet, do so now
call nucxs % find_temperature(p % sqrtkT)
last_nuc_temp = set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, p_uvw, last_nuc_pol, &
last_nuc_azi, last_nuc_uvw)
end if
i = 0
@ -1329,13 +1337,13 @@ contains
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('total', p_g, UVW=p_uvw) / &
matxs % get_xs('total', p_g, UVW=p_uvw) * flux
get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g) * flux
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('total', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
atom_density * flux
else
score = material_xs % total * flux
@ -1358,19 +1366,23 @@ contains
end if
if (i_nuclide > 0) then
score = score * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
/ matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
else
score = score * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) &
/ matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
end if
else
if (i_nuclide > 0) then
score = flux * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
score = flux * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
else
score = flux * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw)
score = flux * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
end if
end if
@ -1392,21 +1404,23 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu) / &
matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
else
! Get the scattering x/s and take away
! the multiplication baked in to sigS
score = flux * &
matxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
end if
end if
@ -1428,19 +1442,20 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('scatter*f_mu', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu) / &
matxs % get_xs('scatter*f_mu', p % last_g, p % g, &
UVW=p_uvw, MU=p % mu)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * &
atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g)
else
! Get the scattering x/s, which includes multiplication
score = matxs % get_xs('scatter', p_g, UVW=p_uvw) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
end if
end if
@ -1460,13 +1475,13 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('absorption', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('absorption', p_g, UVW=p_uvw) * &
atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
else
score = material_xs % absorption * flux
end if
@ -1491,19 +1506,19 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
end if
end if
@ -1529,12 +1544,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1547,17 +1562,17 @@ contains
score = keff * p % wgt_bank * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('nu-fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
end if
end if
@ -1583,12 +1598,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1602,17 +1617,17 @@ contains
/ real(p % n_bank, 8)) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
end if
end if
@ -1638,7 +1653,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1653,13 +1668,13 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed-nu-fission', &
p_g, UVW=p_uvw, dg=d) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('delayed-nu-fission', &
p_g, UVW=p_uvw, dg=d) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1669,11 +1684,13 @@ contains
else
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end if
end if
@ -1703,8 +1720,8 @@ contains
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1715,8 +1732,8 @@ contains
score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
end if
@ -1735,11 +1752,11 @@ contains
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) * atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1748,11 +1765,12 @@ contains
end select
else
if (i_nuclide > 0) then
score = nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
* atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
else
score = matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) &
* flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
end if
end if
end if
@ -1767,7 +1785,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (matxs % get_xs('absorption', p_g, UVW=p_uvw) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1782,17 +1800,15 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * nucxs % get_xs('decay rate', p_g, &
UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
p_g, UVW=p_uvw)
score = score * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * matxs % get_xs('decay rate', p_g, &
UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, &
UVW=p_uvw, dg=d) / matxs % get_xs('absorption', &
p_g, UVW=p_uvw)
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1809,15 +1825,15 @@ contains
! for all delayed groups.
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + p % absorb_wgt * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score + p % absorb_wgt * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score + p % absorb_wgt * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux
score = score + p % absorb_wgt * flux * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end do
end if
@ -1846,13 +1862,13 @@ contains
if (i_nuclide > 0) then
score = score + keff * atom_density * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * &
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
matxs % get_xs('fission', p_g, UVW=p_uvw) * flux
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * flux
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
end if
! if the delayed group filter is present, tally to corresponding
@ -1904,13 +1920,13 @@ contains
d = filt % groups(d_bin)
if (i_nuclide > 0) then
score = nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) * atom_density * flux
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, &
dg=d) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1927,12 +1943,12 @@ contains
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + atom_density * flux * &
nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = score + flux * &
matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * &
matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, dg=d)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
end do
end if
@ -1957,19 +1973,20 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / &
matxs % get_xs('absorption', p_g, UVW=p_uvw)
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
atom_density * flux
else
score = matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * flux
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
end if
end if
@ -1989,7 +2006,15 @@ contains
end do SCORE_LOOP
nullify(matxs, nucxs)
! Reset temporary Mgxs indices
call reset_macro_angle_index_c(p % material, last_mat_pol, last_mat_azi, &
last_mat_uvw);
if (i_nuclide > 0) then
call reset_nuclide_temperature_index_c(i_nuclide, last_nuc_temp)
call reset_nuclide_angle_index_c(i_nuclide, last_nuc_pol, last_nuc_azi, &
last_nuc_uvw)
end if
end subroutine score_general_mg
!===============================================================================

2
src/tallies/tally_filter_energy.F90
View file

@ -6,7 +6,7 @@ module tally_filter_energy
use constants
use error
use hdf5_interface
use mgxs_header, only: num_energy_groups, rev_energy_bins
use mgxs_interface, only: num_energy_groups, rev_energy_bins
use particle_header, only: Particle
use settings, only: run_CE
use string, only: to_str

21
src/tracking.F90
View file

@ -9,7 +9,7 @@ module tracking
check_cell_overlap
use material_header, only: materials, Material
use message_passing
use mgxs_header
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, Particle
use physics, only: collision
@ -29,21 +29,6 @@ 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
!===============================================================================
@ -129,10 +114,6 @@ 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)

15
src/xsdata.cpp
View file

@ -24,6 +24,8 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
if (fissionable) {
fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
nu_fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
prompt_nu_fission = double_3dvec(n_pol, double_2dvec(n_azi,
double_1dvec(energy_groups, 0.)));
kappa_fission = double_3dvec(n_pol, double_2dvec(n_azi,
@ -514,6 +516,17 @@ void XsData::_fissionable_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 p = 0; p < n_pol; p++) {
for (int a = 0; a < n_azi; a++) {
for (int gin = 0; gin < energy_groups; gin++) {
nu_fission[p][a][gin] =
std::accumulate(delayed_nu_fission[p][a][gin].begin(),
delayed_nu_fission[p][a][gin].end(),
prompt_nu_fission[p][a][gin]);
}
}
}
}
@ -668,6 +681,8 @@ void XsData::combine(std::vector<XsData*> those_xs, double_1dvec& scalars)
inverse_velocity[p][a][gin] +=
scalar * that->inverse_velocity[p][a][gin];
if (that->prompt_nu_fission.size() > 0) {
nu_fission[p][a][gin] +=
scalar * that->nu_fission[p][a][gin];
prompt_nu_fission[p][a][gin] +=
scalar * that->prompt_nu_fission[p][a][gin];
kappa_fission[p][a][gin] +=

1
src/xsdata.h
View file

@ -38,6 +38,7 @@ class XsData {
// [phi][theta][incoming group]
double_3dvec total;
double_3dvec absorption;
double_3dvec nu_fission;
double_3dvec prompt_nu_fission;
double_3dvec kappa_fission;
double_3dvec fission;