crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fix errors in conversion of Nuclide methods

Browse source
This commit is contained in:
Paul Romano 2018-12-21 08:22:50 -06:00
parent bc13a65711
commit 015360d477
10 changed files with 78 additions and 79 deletions
Download patch file Download diff file Expand all files Collapse all files

15
include/openmc/nuclide.h
View file

@ -98,10 +98,7 @@ public:
Nuclide(hid_t group, const double* temperature, int n, int i_nuclide);
//! Initialize logarithmic grid for energy searches
//! \param E_min Minimum energy in [eV]
//! \param E_max Maximum energy in [eV]
//! \param M Number of equally log-spaced bins
void init_grid(double E_min, double E_max, int M);
void init_grid();
void calculate_xs(int i_sab, double E, int i_log_union,
double sqrtkT, double sab_frac);
@ -160,11 +157,11 @@ public:
private:
void create_derived();
static constexpr int XS_TOTAL {0};
static constexpr int XS_ABSORPTION {1};
static constexpr int XS_FISSION {2};
static constexpr int XS_NU_FISSION {3};
static constexpr int XS_PHOTON_PROD {4};
static int XS_TOTAL;
static int XS_ABSORPTION;
static int XS_FISSION;
static int XS_NU_FISSION;
static int XS_PHOTON_PROD;
};
//==============================================================================

12
include/openmc/thermal.h
View file

@ -137,18 +137,6 @@ public:
std::vector<ThermalData> data_;
};
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature,
int n, int method, double tolerance, const double* minmax);
void sab_free(ThermalScattering* data);
bool sab_has_nuclide(ThermalScattering* data, const char* name);
double sab_threshold(ThermalScattering* data);
}
} // namespace openmc
#endif // OPENMC_THERMAL_H

3
src/input_xml.F90
View file

@ -2157,6 +2157,9 @@ contains
end if
! Set up logarithmic grid for nuclides
do i = 1, size(nuclides)
call nuclides(i) % init_grid()
end do
log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins
do i = 1, size(materials)

2
src/material_header.F90
View file

@ -443,7 +443,7 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
.or. i_sab /= micro_xs(i_nuclide) % index_sab &
.or. i_sab /= micro_xs(i_nuclide) % index_sab + 1 &
.or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
call nuclides(i_nuclide) % calculate_xs(i_sab, p % E, i_grid, &
p % sqrtkT, sab_frac)

35
src/nuclide.cpp
View file

@ -39,7 +39,14 @@ MaterialMacroXS material_xs;
// Nuclide implementation
//==============================================================================
int Nuclide::XS_TOTAL {0};
int Nuclide::XS_ABSORPTION {1};
int Nuclide::XS_FISSION {2};
int Nuclide::XS_NU_FISSION {3};
int Nuclide::XS_PHOTON_PROD {4};
Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
: i_nuclide_{i_nuclide}
{
// Get name of nuclide from group, removing leading '/'
name_ = object_name(group).substr(1);
@ -48,7 +55,6 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
read_attribute(group, "A", A_);
read_attribute(group, "metastable", metastable_);
read_attribute(group, "atomic_weight_ratio", awr_);
i_nuclide_ = i_nuclide;
// Determine temperatures available
hid_t kT_group = open_group(group, "kTs");
@ -382,8 +388,13 @@ void Nuclide::create_derived()
}
}
void Nuclide::init_grid(double E_min, double E_max, int M)
void Nuclide::init_grid()
{
int neutron = static_cast<int>(ParticleType::neutron) - 1;
double E_min = data::energy_min[neutron];
double E_max = data::energy_max[neutron];
int M = settings::n_log_bins;
// Determine equal-logarithmic energy spacing
double spacing = std::log(E_max/E_min)/M;
@ -398,10 +409,10 @@ void Nuclide::init_grid(double E_min, double E_max, int M)
// equal-logarithmic grid
int j = 0;
for (int k = 0; k <= M; ++k) {
while (std::log(grid.energy[j]/E_min) <= umesh(k)) {
while (std::log(grid.energy[j + 1]/E_min) <= umesh(k)) {
// Ensure that for isotopes where maxval(grid.energy) << E_max that
// there are no out-of-bounds issues.
if (j == grid.energy.size()) break;
if (j + 1 == grid.energy.size()) break;
++j;
}
grid.grid_index[k] = j;
@ -453,7 +464,7 @@ void Nuclide::calculate_elastic_xs() const
{
// Get temperature index, grid index, and interpolation factor
auto& micro = simulation::micro_xs[i_nuclide_];
int i_temp = micro.index_temp - 1;
int i_temp = micro.index_temp;
int i_grid = micro.index_grid - 1;
double f = micro.interp_factor;
@ -589,7 +600,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
int i_high = grid.grid_index[i_log_union + 1] + 1;
// Perform binary search over reduced range
i_grid = lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E) - 1;
i_grid = i_low + lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E);
}
// check for rare case where two energy points are the same
@ -600,7 +611,8 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
(grid.energy[i_grid + 1]- grid.energy[i_grid]);
micro_xs.index_temp = i_temp;
micro_xs.index_grid = i_grid;
// TODO: off-by-one
micro_xs.index_grid = i_grid + 1;
micro_xs.interp_factor = f;
// Calculate microscopic nuclide total cross section
@ -871,6 +883,8 @@ extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature,
return data::nuclides.back().get();
}
extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); }
extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx)
{
return nuc->reactions_[i_rx-1].get();
@ -882,10 +896,15 @@ extern "C" void nuclide_calculate_xs_c(Nuclide* nuc, int i_sab, double E,
nuc->calculate_xs(i_sab, E, i_log_union, sqrtkT, sab_frac);
}
extern "C" void nuclide_calculate_elastic_xs_c(Nuclide* nuc)
{
nuc->calculate_elastic_xs();
}
extern "C" void nuclides_clear() { data::nuclides.clear(); }
extern "C" NuclideMicroXS* micro_xs_ptr();
void set_micro_xs()
{
#pragma omp parallel

45
src/nuclide_header.F90
View file

@ -94,6 +94,7 @@ module nuclide_header
contains
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
procedure :: calculate_xs => nuclide_calculate_xs
@ -673,6 +674,17 @@ contains
end function nuclide_nu
subroutine nuclide_init_grid(this)
class(nuclide), intent(inout) :: this
interface
subroutine nuclide_init_grid_c(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call nuclide_init_grid_c(this % ptr)
end subroutine
!===============================================================================
! NUCLIDE_CALCULATE_XS determines microscopic cross sections for the nuclide
! at the energy of the given particle
@ -689,8 +701,9 @@ contains
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
interface
subroutine nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C)
import C_INT, C_DOUBLE
subroutine nuclide_calculate_xs_c(ptr, i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: i_sab
real(C_DOUBLE), value :: E
integer(C_INT), value :: i_log_union
@ -699,7 +712,7 @@ contains
end subroutine
end interface
call nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac)
call nuclide_calculate_xs_c(this % ptr, i_sab - 1, E, i_log_union, sqrtkT, sab_frac)
end subroutine nuclide_calculate_xs
!===============================================================================
@ -709,25 +722,15 @@ contains
! to calculate it early to adjust the total cross section correctly.
!===============================================================================
subroutine nuclide_calculate_elastic_xs(this, micro_xs)
subroutine nuclide_calculate_elastic_xs(this)
class(Nuclide), intent(in) :: this
type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
integer :: i_temp
integer :: i_grid
real(8) :: f
! Get temperature index, grid index, and interpolation factor
i_temp = micro_xs % index_temp
i_grid = micro_xs % index_grid
f = micro_xs % interp_factor
if (i_temp > 0) then
associate (rx => this % reactions(1))
micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + &
f * rx % xs(i_temp, i_grid + 1)
end associate
end if
interface
subroutine nuclide_calculate_elastic_xs_c(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call nuclide_calculate_elastic_xs_c(this % ptr)
end subroutine nuclide_calculate_elastic_xs
!===============================================================================

7
src/physics.cpp
View file

@ -518,7 +518,7 @@ Reaction* sample_fission(int i_nuclide, double E)
}
// Get grid index and interpolatoin factor and sample fission cdf
int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1;
int i_temp = simulation::micro_xs[i_nuclide].index_temp;
int i_grid = simulation::micro_xs[i_nuclide].index_grid;
double f = simulation::micro_xs[i_nuclide].interp_factor;
double cutoff = prn() * simulation::micro_xs[i_nuclide].fission;
@ -542,8 +542,7 @@ Reaction* sample_fission(int i_nuclide, double E)
void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product)
{
// Get grid index and interpolation factor and sample photon production cdf
// TODO: off-by-one
int i_temp = simulation::micro_xs[i_nuclide].index_temp - 1;
int i_temp = simulation::micro_xs[i_nuclide].index_temp;
int i_grid = simulation::micro_xs[i_nuclide].index_grid;
double f = simulation::micro_xs[i_nuclide].interp_factor;
double cutoff = prn() * simulation::micro_xs[i_nuclide].photon_prod;
@ -616,7 +615,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
// Get pointer to nuclide and grid index/interpolation factor
const auto& nuc {data::nuclides[i_nuclide]};
const auto& micro {simulation::micro_xs[i_nuclide]};
int i_temp = micro.index_temp - 1;
int i_temp = micro.index_temp;
int i_grid = micro.index_grid - 1;
double f = micro.interp_factor;

19
src/sab_header.F90
View file

@ -21,7 +21,6 @@ module sab_header
type(C_PTR) :: ptr
contains
procedure :: from_hdf5
procedure :: free
procedure :: has_nuclide
procedure :: threshold
end type SAlphaBeta
@ -44,9 +43,7 @@ module sab_header
type(C_PTR) :: ptr
end function
subroutine sab_free(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
subroutine sab_clear() bind(C)
end subroutine
function sab_has_nuclide(ptr, name) result(val) bind(C)
@ -87,11 +84,6 @@ contains
end if
end subroutine from_hdf5
subroutine free(this)
class(SAlphaBeta), intent(inout) :: this
call sab_free(this % ptr)
end subroutine
function has_nuclide(this, name) result(val)
class(SAlphaBeta), intent(in) :: this
character(len=*), intent(in) :: name
@ -111,14 +103,9 @@ contains
!===============================================================================
subroutine free_memory_sab()
integer :: i
n_sab_tables = 0
if (allocated(sab_tables)) then
do i = 1, size(sab_tables)
call sab_tables(i) % free()
end do
deallocate(sab_tables)
end if
call sab_clear()
if (allocated(sab_tables)) deallocate(sab_tables)
call sab_dict % clear()
end subroutine free_memory_sab

8
src/tallies/tally.F90
View file

@ -958,7 +958,7 @@ contains
else
if (i_nuclide > 0) then
if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
call nuclides(i_nuclide) % calculate_elastic_xs(micro_xs(i_nuclide))
call nuclides(i_nuclide) % calculate_elastic_xs()
end if
score = micro_xs(i_nuclide) % elastic * atom_density * flux
else
@ -971,7 +971,7 @@ contains
! Get index in nuclides array
i_nuc = materials(p % material) % nuclide(l)
if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then
call nuclides(i_nuc) % calculate_elastic_xs(micro_xs(i_nuc))
call nuclides(i_nuc) % calculate_elastic_xs()
end if
score = score + micro_xs(i_nuc) % elastic * atom_density_ * flux
@ -1126,7 +1126,7 @@ contains
if (m /= 0) then
! Retrieve temperature and energy grid index and interpolation
! factor
i_temp = micro_xs(i_nuclide) % index_temp
i_temp = micro_xs(i_nuclide) % index_temp + 1
if (i_temp > 0) then
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
@ -1159,7 +1159,7 @@ contains
if (m /= 0) then
! Retrieve temperature and energy grid index and
! interpolation factor
i_temp = micro_xs(i_nuc) % index_temp
i_temp = micro_xs(i_nuc) % index_temp + 1
if (i_temp > 0) then
i_energy = micro_xs(i_nuc) % index_grid
f = micro_xs(i_nuc) % interp_factor

11
src/thermal.cpp
View file

@ -590,7 +590,7 @@ ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
// Fortran compatibility functions
//==============================================================================
ThermalScattering*
extern "C" ThermalScattering*
sab_from_hdf5(hid_t group, const double* temperature, int n,
int method, double tolerance, const double* minmax)
{
@ -604,13 +604,16 @@ sab_from_hdf5(hid_t group, const double* temperature, int n,
return data::thermal_scatt.back().get();
}
void sab_free(ThermalScattering* data) { delete data; }
extern "C" void sab_clear() { data::thermal_scatt.clear(); }
bool sab_has_nuclide(ThermalScattering* data, const char* name)
extern "C" bool sab_has_nuclide(ThermalScattering* data, const char* name)
{
return data->has_nuclide(name);
}
double sab_threshold(ThermalScattering* data) { return data->threshold(); }
extern "C" double sab_threshold(ThermalScattering* data)
{
return data->threshold();
}
} // namespace openmc