From d4c8c3d34e4f377ca25eb14bdbc062a309f497d5 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Thu, 6 Dec 2018 19:28:26 -0500 Subject: [PATCH] completed conversion of Urr to C++ --- CMakeLists.txt | 1 - include/openmc/nuclide.h | 3 +- include/openmc/urr.h | 1 + src/nuclide.cpp | 12 +++++-- src/nuclide_header.F90 | 73 ++++++---------------------------------- src/urr.cpp | 3 +- src/urr_header.F90 | 72 --------------------------------------- 7 files changed, 25 insertions(+), 140 deletions(-) delete mode 100644 src/urr_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index da55400c4..eb2224867 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -361,7 +361,6 @@ add_library(libopenmc SHARED src/timer_header.F90 src/tracking.F90 src/track_output.F90 - src/urr_header.F90 src/vector_header.F90 src/volume_calc.F90 src/volume_header.F90 diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 662a6fa7c..67444a106 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -175,10 +175,9 @@ extern "C" void set_micro_xs(); extern "C" bool nuclide_wmp_present(int i_nuclide); extern "C" double nuclide_wmp_emin(int i_nuclide); extern "C" double nuclide_wmp_emax(int i_nuclide); -extern "C" void nuclide_calculate_urr_xs(const int i_nuclide, +extern "C" void nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide, const int i_temp, const double E); - } // namespace openmc #endif // OPENMC_NUCLIDE_H diff --git a/include/openmc/urr.h b/include/openmc/urr.h index acfe4d059..36ea07332 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -20,6 +20,7 @@ public: int inelastic_flag_; // inelastic competition flag int absorption_flag_; // other absorption flag bool multiply_smooth_; // multiply by smooth cross section? + int n_energy_; // number of energy points xt::xtensor energy_; // incident energies xt::xtensor prob_; // Actual probability tables diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 64c983fd3..727968a4d 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -576,9 +576,17 @@ void set_micro_xs() } extern "C" void -nuclide_calculate_urr_xs(const int i_nuclide, const int i_temp, const double E) +nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide, + const int i_temp, const double E) { - data::nuclides[i_nuclide - 1]->calculate_urr_xs(i_temp - 1, E); + Nuclide* nuc = data::nuclides[i_nuclide - 1].get(); + if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) { + if ((E > nuc->urr_data_[i_temp - 1].energy_(0)) && + (E < nuc->urr_data_[i_temp - 1].energy_( + nuc->urr_data_[i_temp - 1].n_energy_ - 1))) { + nuc->calculate_urr_xs(i_temp - 1, E); + } + } } } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index e05a70313..b556df665 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -23,7 +23,6 @@ module nuclide_header use stl_vector, only: VectorInt, VectorReal use string use simulation_header, only: need_depletion_rx - use urr_header, only: UrrData implicit none @@ -76,11 +75,6 @@ module nuclide_header integer, allocatable :: index_fission(:) ! indices in reactions class(Function1D), allocatable :: total_nu - ! Unresolved resonance data - logical :: urr_present = .false. - integer :: urr_inelastic - type(UrrData), allocatable :: urr_data(:) - ! Multipole data logical :: mp_present = .false. type(MultipoleArray), pointer :: multipole => null() @@ -199,12 +193,13 @@ module nuclide_header type(C_PTR) :: path end function - subroutine nuclide_calculate_urr_xs_c(i_nuclide, i_temp, E) & + subroutine nuclide_calculate_urr_xs_c(use_mp, i_nuclide, i_temp, E) & bind(C, name='nuclide_calculate_urr_xs') - import C_INT, C_DOUBLE - integer(C_INT), value, intent(in) :: i_nuclide - integer(C_INT), value, intent(in) :: i_temp - real(C_DOUBLE), value, intent(in) :: E + import C_BOOL, C_INT, C_DOUBLE + logical(C_BOOL), value, intent(in) :: use_mp + integer(C_INT), value, intent(in) :: i_nuclide + integer(C_INT), value, intent(in) :: i_temp + real(C_DOUBLE), value, intent(in) :: E end subroutine nuclide_calculate_urr_xs_c end interface @@ -261,7 +256,7 @@ contains integer :: i integer :: i_closest integer :: n_temperature - integer(HID_T) :: urr_group, nu_group + integer(HID_T) :: nu_group integer(HID_T) :: energy_group, energy_dset integer(HID_T) :: kT_group integer(HID_T) :: rxs_group @@ -446,47 +441,6 @@ contains end do call close_group(rxs_group) - ! Read unresolved resonance probability tables if present - if (object_exists(group_id, 'urr')) then - this % urr_present = .true. - allocate(this % urr_data(n_temperature)) - - do i = 1, n_temperature - ! Get temperature as a string - temp_str = trim(to_str(temps_to_read % data(i))) // "K" - - ! Read probability tables for i-th temperature - urr_group = open_group(group_id, 'urr/' // trim(temp_str)) - call this % urr_data(i) % from_hdf5(urr_group) - call close_group(urr_group) - - ! Check for negative values - if (any(this % urr_data(i) % prob < ZERO) .and. master) then - call warning("Negative value(s) found on probability table & - &for nuclide " // this % name // " at " // trim(temp_str)) - end if - end do - - ! if the inelastic competition flag indicates that the inelastic cross - ! section should be determined from a normal reaction cross section, we - ! need to get the index of the reaction - if (n_temperature > 0) then - if (this % urr_data(1) % inelastic_flag > 0) then - do i = 1, size(this % reactions) - if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then - this % urr_inelastic = i - end if - end do - - ! Abort if no corresponding inelastic reaction was found - if (this % urr_inelastic == NONE) then - call fatal_error("Could not find inelastic reaction specified on & - &unresolved resonance probability table.") - end if - end if - end if - end if - ! Check for nu-total if (object_exists(group_id, 'total_nu')) then nu_group = open_group(group_id, 'total_nu') @@ -784,7 +738,7 @@ contains real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - logical :: use_mp ! true if XS can be calculated with windowed multipole + logical(C_BOOL) :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index @@ -839,7 +793,7 @@ contains ! 1. physics.F90 - scatter - For inelastic scatter. ! 2. physics.F90 - sample_fission - For partial fissions. ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. - ! 4. nuclide_header.F90 - calculate_urr_xs - For unresolved purposes. + ! 4. nuclide.h - calculate_urr_xs - For unresolved purposes. ! It is worth noting that none of these occur in the resolved ! resonance range, so the value here does not matter. index_temp is ! set to -1 to force a segfault in case a developer messes up and tries @@ -982,15 +936,10 @@ contains call calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs) end if + ! If the particle is in the unresolved resonance range and there are ! probability tables, we need to determine cross sections from the table - - if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then - if (E > this % urr_data(i_temp) % energy(1) .and. E < this % & - urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then - call nuclide_calculate_urr_xs_c(this % i_nuclide, i_temp, E) - end if - end if + call nuclide_calculate_urr_xs_c(use_mp, this % i_nuclide, i_temp, E) micro_xs % last_E = E micro_xs % last_sqrtkT = sqrtkT diff --git a/src/urr.cpp b/src/urr.cpp index 1ecc8ca8c..1ed7bca39 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -27,6 +27,7 @@ UrrData::from_hdf5(hid_t group_id) interp_ = Interpolation::log_log; } + // read the metadata read_attribute(group_id, "inelastic", inelastic_flag_); read_attribute(group_id, "absorption", absorption_flag_); int temp_multiply_smooth; @@ -38,7 +39,7 @@ UrrData::from_hdf5(hid_t group_id) hsize_t dims[1]; get_shape(dset, dims); close_dataset(dset); - // n_energy_ = static_cast(dims[0]); + n_energy_ = static_cast(dims[0]); energy_ = xt::xtensor({dims[0]}, 0.); read_dataset_as_shape(group_id, "energy", energy_); diff --git a/src/urr_header.F90 b/src/urr_header.F90 deleted file mode 100644 index 0ec9e89de..000000000 --- a/src/urr_header.F90 +++ /dev/null @@ -1,72 +0,0 @@ -module urr_header - - use hdf5_interface, only: read_attribute, open_dataset, read_dataset, & - close_dataset, get_shape, HID_T, HSIZE_T - - implicit none - -!=============================================================================== -! URRDATA contains probability tables for the unresolved resonance range. -!=============================================================================== - - type UrrData - integer :: n_energy ! # of incident neutron energies - integer :: n_prob ! # of probabilities - integer :: interp ! inteprolation (2=lin-lin, 5=log-log) - integer :: inelastic_flag ! inelastic competition flag - integer :: absorption_flag ! other absorption flag - logical :: multiply_smooth ! multiply by smooth cross section? - real(8), allocatable :: energy(:) ! incident energies - real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - contains - procedure :: from_hdf5 => urr_from_hdf5 - end type UrrData - -contains - - subroutine urr_from_hdf5(this, group_id) - class(UrrData), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j, k - integer(HID_T) :: energy - integer(HID_T) :: table - integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: dims3(3) - real(8), allocatable :: temp(:,:,:) - - ! Read interpolation and other flags - call read_attribute(this % interp, group_id, 'interpolation') - call read_attribute(this % inelastic_flag, group_id, 'inelastic') - call read_attribute(this % absorption_flag, group_id, 'absorption') - call read_attribute(i, group_id, 'multiply_smooth') - this % multiply_smooth = (i == 1) - - ! Read energies at which tables exist - energy = open_dataset(group_id, 'energy') - call get_shape(energy, dims) - this % n_energy = int(dims(1), 4) - allocate(this % energy(this % n_energy)) - call read_dataset(this % energy, energy) - call close_dataset(energy) - - ! Read URR tables - table = open_dataset(group_id, 'table') - call get_shape(table, dims3) - this % n_prob = int(dims3(1), 4) - allocate(temp(this % n_prob, 6, this % n_energy)) - call read_dataset(temp, table) - call close_dataset(table) - - ! Swap first and last indices - allocate(this % prob(this % n_energy, 6, this % n_prob)) - do i = 1, this % n_energy - do j = 1, 6 - do k = 1, this % n_prob - this % prob(i, j, k) = temp(k, j, i) - end do - end do - end do - end subroutine urr_from_hdf5 - -end module urr_header