diff --git a/src/material_header.F90 b/src/material_header.F90 index 03ecf69d9..9a1502493 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -172,7 +172,7 @@ contains found = .false. associate (sab => sab_tables(this % i_sab_tables(k))) FIND_NUCLIDE: do j = 1, size(this % nuclide) - if (any(sab % nuclides == nuclides(this % nuclide(j)) % name)) then + if (sab % has_nuclide(nuclides(this % nuclide(j)) % name)) then call i_sab_tables % push_back(this % i_sab_tables(k)) call i_sab_nuclides % push_back(j) call sab_fracs % push_back(this % sab_fracs(k)) @@ -326,7 +326,7 @@ contains ! If particle energy is greater than the highest energy for the ! S(a,b) table, then don't use the S(a,b) table - if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then + if (p % E > sab_tables(i_sab) % threshold()) then i_sab = 0 end if diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 6faaaec7f..0411ab102 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -198,7 +198,7 @@ module nuclide_header type(DictCharInt) :: nuclide_dict ! Cross section caches - type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide + type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide type(MaterialMacroXS) :: material_xs ! Cache for current material !$omp threadprivate(micro_xs, material_xs) @@ -1096,9 +1096,9 @@ contains real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - integer :: i_temp ! temperature index - real(8) :: inelastic ! S(a,b) inelastic cross section - real(8) :: elastic ! S(a,b) elastic cross section + integer(C_INT) :: i_temp ! temperature index + real(C_DOUBLE) :: inelastic ! S(a,b) inelastic cross section + real(C_DOUBLE) :: elastic ! S(a,b) elastic cross section ! Set flag that S(a,b) treatment should be used for scattering micro_xs % index_sab = i_sab diff --git a/src/physics.F90 b/src/physics.F90 index d5f6302a8..70c922b42 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -884,263 +884,16 @@ contains real(8), intent(inout) :: uvw(3) ! directional cosines real(8), intent(out) :: mu ! scattering cosine - integer :: i ! incoming energy bin - integer :: j ! outgoing energy bin - integer :: k ! outgoing cosine bin - integer :: i_temp ! temperature index - integer :: n_energy_out ! number of outgoing energy bins - real(8) :: f ! interpolation factor - real(8) :: r ! used for skewed sampling & continuous - real(8) :: E_ij ! outgoing energy j for E_in(i) - real(8) :: E_i1j ! outgoing energy j for E_in(i+1) - real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) - real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) - real(8) :: prob ! probability for sampling Bragg edge - ! Following are needed only for SAB_SECONDARY_CONT scattering - integer :: l ! sampled incoming E bin (is i or i + 1) - real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i - real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1 - real(8) :: E_1, E_J ! endpoints interpolated between i and i+1 - real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l - real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l - real(8) :: c_j, c_j1 ! cumulative probability - real(8) :: frac ! interpolation factor on outgoing energy - real(8) :: r1 ! RNG for outgoing energy - real(8) :: mu_left, mu_right ! adjacent mu values + real(C_DOUBLE) :: E_out + type(C_PTR) :: ptr - i_temp = micro_xs(i_nuclide) % index_temp_sab + ! Sample from C++ side + ptr = C_LOC(micro_xs(i_nuclide)) + call sab_tables(i_sab) % sample(ptr, E, E_out, mu) - ! Get pointer to S(a,b) table - associate (sab => sab_tables(i_sab) % data(i_temp)) - - ! Determine whether inelastic or elastic scattering will occur - if (prn() < micro_xs(i_nuclide) % thermal_elastic / & - micro_xs(i_nuclide) % thermal) then - ! elastic scattering - - ! Get index and interpolation factor for elastic grid - if (E < sab % elastic_e_in(1)) then - i = 1 - f = ZERO - else - i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E) - f = (E - sab%elastic_e_in(i)) / & - (sab%elastic_e_in(i+1) - sab%elastic_e_in(i)) - end if - - ! Select treatment based on elastic mode - if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then - ! With this treatment, we interpolate between two discrete cosines - ! corresponding to neighboring incoming energies. This is used for - ! data derived in the incoherent approximation - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_elastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % elastic_mu(k,i) - mu_i1jk = sab % elastic_mu(k,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - - elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then - ! This treatment is used for data derived in the coherent - ! approximation, i.e. for crystalline structures that have Bragg - ! edges. - - ! Sample a Bragg edge between 1 and i - prob = prn() * sab % elastic_P(i+1) - if (prob < sab % elastic_P(1)) then - k = 1 - else - k = binary_search(sab % elastic_P(1:i+1), i+1, prob) - end if - - ! Characteristic scattering cosine for this Bragg edge - mu = ONE - TWO*sab % elastic_e_in(k) / E - - end if - - ! Outgoing energy is same as incoming energy -- no need to do anything - - else - ! Perform inelastic calculations - - ! Get index and interpolation factor for inelastic grid - if (E < sab % inelastic_e_in(1)) then - i = 1 - f = ZERO - else - i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i)) / & - (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) - end if - - ! Now that we have an incoming energy bin, we need to determine the - ! outgoing energy bin. This will depend on the "secondary energy - ! mode". If the mode is 0, then the outgoing energy bin is chosen from a - ! set of equally-likely bins. If the mode is 1, then the first - ! two and last two bins are skewed to have lower probabilities than the - ! other bins (0.1 for the first and last bins and 0.4 for the second and - ! second to last bins, relative to a normal bin probability of 1). - ! Finally, if the mode is 2, then a continuous distribution (with - ! accompanying PDF and CDF is utilized) - - if ((sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) .or. & - (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED)) then - if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) then - ! All bins equally likely - - j = 1 + int(prn() * sab % n_inelastic_e_out) - elseif (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED) then - ! Distribution skewed away from edge points - - ! Determine number of outgoing energy and angle bins - n_energy_out = sab % n_inelastic_e_out - - r = prn() * (n_energy_out - 3) - if (r > ONE) then - ! equally likely N-4 middle bins - j = int(r) + 2 - elseif (r > 0.6_8) then - ! second to last bin has relative probability of 0.4 - j = n_energy_out - 1 - elseif (r > HALF) then - ! last bin has relative probability of 0.1 - j = n_energy_out - elseif (r > 0.1_8) then - ! second bin has relative probability of 0.4 - j = 2 - else - ! first bin has relative probability of 0.1 - j = 1 - end if - end if - - ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) - E_ij = sab % inelastic_e_out(j,i) - E_i1j = sab % inelastic_e_out(j,i+1) - - ! Outgoing energy - E = (1 - f)*E_ij + f*E_i1j - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % inelastic_mu(k,j,i) - mu_i1jk = sab % inelastic_mu(k,j,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - - else if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_CONT) then - ! Continuous secondary energy - this is to be similar to - ! Law 61 interpolation on outgoing energy - - ! Sample between ith and (i+1)th bin - r = prn() - if (f > r) then - l = i + 1 - else - l = i - end if - - ! Determine endpoints on grid i - n_energy_out = sab % inelastic_data(i) % n_e_out - E_i_1 = sab % inelastic_data(i) % e_out(1) - E_i_J = sab % inelastic_data(i) % e_out(n_energy_out) - - ! Determine endpoints on grid i + 1 - n_energy_out = sab % inelastic_data(i + 1) % n_e_out - E_i1_1 = sab % inelastic_data(i + 1) % e_out(1) - E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out) - - E_1 = E_i_1 + f * (E_i1_1 - E_i_1) - E_J = E_i_J + f * (E_i1_J - E_i_J) - - ! Determine outgoing energy bin - ! (First reset n_energy_out to the right value) - n_energy_out = sab % inelastic_data(l) % n_e_out - r1 = prn() - c_j = sab % inelastic_data(l) % e_out_cdf(1) - do j = 1, n_energy_out - 1 - c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1) - if (r1 < c_j1) exit - c_j = c_j1 - end do - - ! check to make sure k is <= n_energy_out - 1 - j = min(j, n_energy_out - 1) - - ! Get the data to interpolate between - E_l_j = sab % inelastic_data(l) % e_out(j) - p_l_j = sab % inelastic_data(l) % e_out_pdf(j) - - ! Next part assumes linear-linear interpolation in standard - E_l_j1 = sab % inelastic_data(l) % e_out(j + 1) - p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1) - - ! Find secondary energy (variable E) - frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j) - if (frac == ZERO) then - E = E_l_j + (r1 - c_j) / p_l_j - else - E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + & - TWO * frac * (r1 - c_j))) - p_l_j) / frac - end if - - ! Now interpolate between incident energy bins i and i + 1 - if (l == i) then - E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1) - else - E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) - end if - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Rather than use the sampled discrete mu directly, it is smeared over - ! a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the - ! discrete mu value itself. - associate (mu_l => sab % inelastic_data(l) % mu) - f = (r1 - c_j)/(c_j1 - c_j) - - ! Determine (k-1)th mu value - if (k == 1) then - mu_left = -ONE - else - mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j)) - end if - - ! Determine kth mu value - mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j)) - - ! Determine (k+1)th mu value - if (k == sab % n_inelastic_mu) then - mu_right = ONE - mu - else - mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu - end if - end associate - - ! Smear angle - mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF) - - end if ! (inelastic secondary energy treatment) - end if ! (elastic or inelastic) - end associate - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly - ! -1 or 1 - - if (abs(mu) > ONE) mu = sign(ONE,mu) - - ! change direction of particle + ! Set energy to outgoing, change direction of particle + E = E_out uvw = rotate_angle(uvw, mu) - end subroutine sab_scatter !=============================================================================== diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 615e7b281..ece262bd6 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -1,467 +1,169 @@ module sab_header use, intrinsic :: ISO_C_BINDING - use, intrinsic :: ISO_FORTRAN_ENV - use algorithm, only: find, sort, binary_search - use constants - use dict_header, only: DictIntInt, DictCharInt - use distribution_univariate, only: Tabular - use error, only: warning, fatal_error + use dict_header, only: DictCharInt use hdf5_interface - use random_lcg, only: prn - use secondary_correlated, only: CorrelatedAngleEnergy - use settings - use stl_vector, only: VectorInt, VectorReal - use string, only: to_str, str_to_int + use stl_vector, only: VectorReal + use string, only: to_c_string implicit none + private -!=============================================================================== -! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic -! thermal scattering collisions which utilize a continuous secondary energy -! representation. -!=============================================================================== - - type DistEnergySab - integer :: n_e_out - real(8), allocatable :: e_out(:) - real(8), allocatable :: e_out_pdf(:) - real(8), allocatable :: e_out_cdf(:) - real(8), allocatable :: mu(:,:) - end type DistEnergySab + public :: free_memory_sab !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc !=============================================================================== - type SabData - ! threshold for S(a,b) treatment (usually ~4 eV) - real(8) :: threshold_inelastic - real(8) :: threshold_elastic = ZERO - - ! Inelastic scattering data - integer :: n_inelastic_e_in ! # of incoming E for inelastic - integer :: n_inelastic_e_out ! # of outgoing E for inelastic - integer :: n_inelastic_mu ! # of outgoing angles for inelastic - real(8), allocatable :: inelastic_e_in(:) - real(8), allocatable :: inelastic_sigma(:) - ! The following are used only if secondary_mode is 0 or 1 - real(8), allocatable :: inelastic_e_out(:,:) - real(8), allocatable :: inelastic_mu(:,:,:) - ! The following is used only if secondary_mode is 3 - ! The different implementation is necessary because the continuous - ! representation has a variable number of outgoing energy points for each - ! incoming energy - type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein - - ! Elastic scattering data - integer :: elastic_mode ! elastic mode (discrete/exact) - integer :: n_elastic_e_in ! # of incoming E for elastic - integer :: n_elastic_mu ! # of outgoing angles for elastic - real(8), allocatable :: elastic_e_in(:) - real(8), allocatable :: elastic_P(:) - real(8), allocatable :: elastic_mu(:,:) - end type SabData - - type SAlphaBeta - character(150) :: name ! name of table, e.g. lwtr.10t - real(8) :: awr ! weight of nucleus in neutron masses - real(8), allocatable :: kTs(:) ! temperatures in eV (k*T) - character(10), allocatable :: nuclides(:) ! List of valid nuclides - integer :: secondary_mode ! secondary mode (equal/skewed/continuous) - - ! cross sections and distributions at each temperature - type(SabData), allocatable :: data(:) + type, public :: SAlphaBeta + type(C_PTR) :: ptr contains - procedure :: from_hdf5 => salphabeta_from_hdf5 - procedure :: calculate_xs => sab_calculate_xs + procedure :: from_hdf5 + procedure :: calculate_xs + procedure :: free + procedure :: has_nuclide + procedure :: sample + procedure :: threshold end type SAlphaBeta ! S(a,b) tables - type(SAlphaBeta), allocatable, target :: sab_tables(:) - integer(C_INT), bind(C) :: n_sab_tables - type(DictCharInt) :: sab_dict + type(SAlphaBeta), public, allocatable, target :: sab_tables(:) + integer(C_INT), public, bind(C) :: n_sab_tables + type(DictCharInt), public :: sab_dict + + interface + function sab_from_hdf5(group_id, temperature, n, method, & + tolerance, minmax) result(ptr) bind(C) + import HID_T, C_DOUBLE, C_INT, C_PTR + integer(HID_T), value :: group_id + real(C_DOUBLE), intent(in) :: temperature + integer(C_INT), value :: n + integer(C_INT), value :: method + real(C_DOUBLE), value :: tolerance + real(C_DOUBLE), intent(in) :: minmax(2) + type(C_PTR) :: ptr + end function + + subroutine sab_calculate_xs(ptr, E, sqrtkT, i_temp, elastic, & + inelastic) bind(C) + import C_PTR, C_DOUBLE, C_INT + type(C_PTR), value :: ptr + real(C_DOUBLE), value :: E + real(C_DOUBLE), value :: sqrtkT + integer(C_INT), intent(out) :: i_temp + real(C_DOUBLE), intent(out) :: elastic + real(C_DOUBLE), intent(out) :: inelastic + end subroutine + + subroutine sab_free(ptr) bind(C) + import C_PTR + type(C_PTR), value :: ptr + end subroutine + + function sab_has_nuclide(ptr, name) result(val) bind(C) + import C_PTR, C_CHAR, C_BOOL + type(C_PTR), value :: ptr + character(kind=C_CHAR), intent(in) :: name(*) + logical(C_BOOL) :: val + end function + + subroutine sab_sample(ptr, micro_xs, E_in, E_out, mu) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + type(C_PTR), value :: micro_xs + real(C_DOUBLE), value :: E_in + real(C_DOUBLE), intent(out) :: E_out + real(C_DOUBLE), intent(out) :: mu + end subroutine + + function sab_threshold(ptr) result(threshold) bind(C) + import C_PTR, C_double + type(C_PTR), value :: ptr + real(C_DOUBLE) :: threshold + end function + end interface contains - subroutine salphabeta_from_hdf5(this, group_id, temperature, method, & + subroutine from_hdf5(this, group_id, temperature, method, & tolerance, minmax) class(SAlphaBeta), intent(inout) :: this integer(HID_T), intent(in) :: group_id type(VectorReal), intent(in) :: temperature ! list of temperatures - integer, intent(in) :: method - real(8), intent(in) :: tolerance - real(8), intent(in) :: minmax(2) + integer(C_INT), intent(in) :: method + real(C_DOUBLE), intent(in) :: tolerance + real(C_DOUBLE), intent(in) :: minmax(2) - integer :: i, j - integer :: t - integer :: n_energy, n_energy_out, n_mu - integer :: i_closest - integer :: n_temperature - integer(HID_T) :: T_group - integer(HID_T) :: elastic_group - integer(HID_T) :: inelastic_group - integer(HID_T) :: dset_id - integer(HID_T) :: kT_group - integer(HSIZE_T) :: dims2(2) - integer(HSIZE_T) :: dims3(3) - real(8), allocatable :: temp(:,:) - character(20) :: type - type(CorrelatedAngleEnergy) :: correlated_dist + real(C_DOUBLE) :: dummy + integer(C_INT) :: n - character(MAX_WORD_LEN) :: temp_str - character(MAX_WORD_LEN), allocatable :: dset_names(:) - real(8), allocatable :: temps_available(:) ! temperatures available - real(8) :: temp_desired - real(8) :: temp_actual - type(VectorInt) :: temps_to_read - - ! Get name of table from group - this % name = get_name(group_id) - - ! Get rid of leading '/' - this % name = trim(this % name(2:)) - - call read_attribute(this % awr, group_id, 'atomic_weight_ratio') - call read_attribute(this % nuclides, group_id, 'nuclides') - call read_attribute(type, group_id, 'secondary_mode') - select case (type) - case ('equal') - this % secondary_mode = SAB_SECONDARY_EQUAL - case ('skewed') - this % secondary_mode = SAB_SECONDARY_SKEWED - case ('continuous') - this % secondary_mode = SAB_SECONDARY_CONT - end select - - ! Read temperatures - kT_group = open_group(group_id, 'kTs') - - ! Determine temperatures available - call get_datasets(kT_group, dset_names) - allocate(temps_available(size(dset_names))) - do i = 1, size(dset_names) - ! Read temperature value - call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) - temps_available(i) = temps_available(i) / K_BOLTZMANN - end do - call sort(temps_available) - - ! Determine actual temperatures to read -- start by checking whether a - ! temperature range was given, in which case all temperatures in the range - ! are loaded irrespective of what temperatures actually appear in the model - if (minmax(2) > ZERO) then - do i = 1, size(temps_available) - temp_actual = temps_available(i) - if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then - call temps_to_read % push_back(nint(temp_actual)) - end if - end do + n = temperature % size() + if (n > 0) then + this % ptr = sab_from_hdf5(group_id, temperature % data(1), n, method, tolerance, minmax) + else + ! In this case, temperatures % data(1) doesn't exist, so we just pass a + ! dummy value + this % ptr = sab_from_hdf5(group_id, dummy, n, method, tolerance, minmax) end if - - select case (method) - case (TEMPERATURE_NEAREST) - ! Determine actual temperatures to read - do i = 1, temperature % size() - temp_desired = temperature % data(i) - i_closest = minloc(abs(temps_available - temp_desired), dim=1) - temp_actual = temps_available(i_closest) - if (abs(temp_actual - temp_desired) < tolerance) then - if (find(temps_to_read, nint(temp_actual)) == -1) then - call temps_to_read % push_back(nint(temp_actual)) - end if - else - call fatal_error("Nuclear data library does not contain cross sections & - &for " // trim(this % name) // " at or near " // & - trim(to_str(nint(temp_desired))) // " K.") - end if - end do - - case (TEMPERATURE_INTERPOLATION) - ! If temperature interpolation or multipole is selected, get a list of - ! bounding temperatures for each actual temperature present in the model - TEMP_LOOP: do i = 1, temperature % size() - temp_desired = temperature % data(i) - - do j = 1, size(temps_available) - 1 - if (temps_available(j) <= temp_desired .and. & - temp_desired < temps_available(j + 1)) then - if (find(temps_to_read, nint(temps_available(j))) == -1) then - call temps_to_read % push_back(nint(temps_available(j))) - end if - if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then - call temps_to_read % push_back(nint(temps_available(j + 1))) - end if - cycle TEMP_LOOP - end if - end do - - call fatal_error("Nuclear data library does not contain cross sections & - &for " // trim(this % name) // " at temperatures that bound " // & - trim(to_str(nint(temp_desired))) // " K.") - end do TEMP_LOOP - - end select - - ! Sort temperatures to read - call sort(temps_to_read) - - n_temperature = temps_to_read % size() - allocate(this % kTs(n_temperature)) - allocate(this % data(n_temperature)) - - do t = 1, n_temperature - ! Get temperature as a string - temp_str = trim(to_str(temps_to_read % data(t))) // "K" - - ! Read exact temperature value - call read_dataset(this % kTs(t), kT_group, temp_str) - - ! Open group for temperature i - T_group = open_group(group_id, temp_str) - - ! Coherent elastic data - if (object_exists(T_group, 'elastic')) then - ! Read cross section data - elastic_group = open_group(T_group, 'elastic') - dset_id = open_dataset(elastic_group, 'xs') - call read_attribute(type, dset_id, 'type') - call get_shape(dset_id, dims2) - allocate(temp(dims2(1), dims2(2))) - call read_dataset(temp, dset_id) - call close_dataset(dset_id) - - ! Set cross section data and type - this % data(t) % n_elastic_e_in = int(dims2(1), 4) - allocate(this % data(t) % elastic_e_in(this % data(t) % n_elastic_e_in)) - allocate(this % data(t) % elastic_P(this % data(t) % n_elastic_e_in)) - this % data(t) % elastic_e_in(:) = temp(:, 1) - this % data(t) % elastic_P(:) = temp(:, 2) - select case (type) - case ('tab1') - this % data(t) % elastic_mode = SAB_ELASTIC_DISCRETE - case ('bragg') - this % data(t) % elastic_mode = SAB_ELASTIC_EXACT - end select - deallocate(temp) - - ! Set elastic threshold - this % data(t) % threshold_elastic = this % data(t) % elastic_e_in(& - this % data(t) % n_elastic_e_in) - - ! Read angle distribution - if (this % data(t) % elastic_mode /= SAB_ELASTIC_EXACT) then - dset_id = open_dataset(elastic_group, 'mu_out') - call get_shape(dset_id, dims2) - this % data(t) % n_elastic_mu = int(dims2(1), 4) - allocate(this % data(t) % elastic_mu(dims2(1), dims2(2))) - call read_dataset(this % data(t) % elastic_mu, dset_id) - call close_dataset(dset_id) - end if - - call close_group(elastic_group) - end if - - ! Inelastic data - if (object_exists(T_group, 'inelastic')) then - ! Read type of inelastic data - inelastic_group = open_group(T_group, 'inelastic') - - ! Read cross section data - dset_id = open_dataset(inelastic_group, 'xs') - call get_shape(dset_id, dims2) - allocate(temp(dims2(1), dims2(2))) - call read_dataset(temp, dset_id) - call close_dataset(dset_id) - - ! Set cross section data - this % data(t) % n_inelastic_e_in = int(dims2(1), 4) - allocate(this % data(t) % inelastic_e_in(this % data(t) % n_inelastic_e_in)) - allocate(this % data(t) % inelastic_sigma(this % data(t) % n_inelastic_e_in)) - this % data(t) % inelastic_e_in(:) = temp(:, 1) - this % data(t) % inelastic_sigma(:) = temp(:, 2) - deallocate(temp) - - ! Set inelastic threshold - this % data(t) % threshold_inelastic = this % data(t) % inelastic_e_in(& - this % data(t) % n_inelastic_e_in) - - if (this % secondary_mode /= SAB_SECONDARY_CONT) then - ! Read energy distribution - dset_id = open_dataset(inelastic_group, 'energy_out') - call get_shape(dset_id, dims2) - this % data(t) % n_inelastic_e_out = int(dims2(1), 4) - allocate(this % data(t) % inelastic_e_out(dims2(1), dims2(2))) - call read_dataset(this % data(t) % inelastic_e_out, dset_id) - call close_dataset(dset_id) - - ! Read angle distribution - dset_id = open_dataset(inelastic_group, 'mu_out') - call get_shape(dset_id, dims3) - this % data(t) % n_inelastic_mu = int(dims3(1), 4) - allocate(this % data(t) % inelastic_mu(dims3(1), dims3(2), dims3(3))) - call read_dataset(this % data(t) % inelastic_mu, dset_id) - call close_dataset(dset_id) - else - ! Read correlated angle-energy distribution - call correlated_dist % from_hdf5(inelastic_group) - - ! Convert to S(a,b) native format - n_energy = size(correlated_dist % energy) - allocate(this % data(t) % inelastic_data(n_energy)) - do i = 1, n_energy - associate (edist => correlated_dist % distribution(i)) - ! Get number of outgoing energies for incoming energy i - n_energy_out = size(edist % e_out) - this % data(t) % inelastic_data(i) % n_e_out = n_energy_out - allocate(this % data(t) % inelastic_data(i) % e_out(n_energy_out)) - allocate(this % data(t) % inelastic_data(i) % e_out_pdf(n_energy_out)) - allocate(this % data(t) % inelastic_data(i) % e_out_cdf(n_energy_out)) - - ! Copy outgoing energy distribution - this % data(t) % inelastic_data(i) % e_out(:) = edist % e_out - this % data(t) % inelastic_data(i) % e_out_pdf(:) = edist % p - this % data(t) % inelastic_data(i) % e_out_cdf(:) = edist % c - - do j = 1, n_energy_out - select type (adist => edist % angle(j) % obj) - type is (Tabular) - ! On first pass, allocate space for angles - if (j == 1) then - n_mu = size(adist % x) - this % data(t) % n_inelastic_mu = n_mu - allocate(this % data(t) % inelastic_data(i) % mu(& - n_mu, n_energy_out)) - end if - - ! Copy outgoing angles - this % data(t) % inelastic_data(i) % mu(:, j) = adist % x - end select - end do - end associate - end do - - ! Clear data on correlated angle-energy object - deallocate(correlated_dist % breakpoints) - deallocate(correlated_dist % interpolation) - deallocate(correlated_dist % energy) - deallocate(correlated_dist % distribution) - end if - - call close_group(inelastic_group) - end if - call close_group(T_group) - end do - - call close_group(kT_group) - end subroutine salphabeta_from_hdf5 + end subroutine from_hdf5 !=============================================================================== ! SAB_CALCULATE_XS determines the elastic and inelastic scattering ! cross-sections in the thermal energy range. !=============================================================================== - subroutine sab_calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic) + subroutine calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic) class(SAlphaBeta), intent(in) :: this ! S(a,b) object - real(8), intent(in) :: E ! energy - real(8), intent(in) :: sqrtkT ! temperature + real(C_DOUBLE), intent(in) :: E ! energy + real(C_DOUBLE), intent(in) :: sqrtkT ! temperature integer, intent(out) :: i_temp ! index in the S(a,b)'s temperature - real(8), intent(out) :: elastic ! thermal elastic cross section - real(8), intent(out) :: inelastic ! thermal inelastic cross section + real(C_DOUBLE), intent(out) :: elastic ! thermal elastic cross section + real(C_DOUBLE), intent(out) :: inelastic ! thermal inelastic cross section - integer :: i_grid ! index on S(a,b) energy grid - real(8) :: f ! interp factor on S(a,b) energy grid - real(8) :: kT + call sab_calculate_xs(this % ptr, E, sqrtkT, i_temp, elastic, inelastic) + end subroutine - ! Determine temperature for S(a,b) table - kT = sqrtkT**2 - if (temperature_method == TEMPERATURE_NEAREST) then - ! If using nearest temperature, do linear search on temperature - do i_temp = 1, size(this % kTs) - if (abs(this % kTs(i_temp) - kT) < & - K_BOLTZMANN*temperature_tolerance) exit - end do - else - ! Find temperatures that bound the actual temperature - do i_temp = 1, size(this % kTs) - 1 - if (this % kTs(i_temp) <= kT .and. & - kT < this % kTs(i_temp + 1)) exit - end do + subroutine free(this) + class(SAlphaBeta), intent(inout) :: this + call sab_free(this % ptr) + end subroutine - ! Randomly sample between temperature i and i+1 - f = (kT - this % kTs(i_temp)) / & - (this % kTs(i_temp + 1) - this % kTs(i_temp)) - if (f > prn()) i_temp = i_temp + 1 - end if + function has_nuclide(this, name) result(val) + class(SAlphaBeta), intent(in) :: this + character(len=*), intent(in) :: name + logical(C_BOOL) :: val + val = sab_has_nuclide(this % ptr, to_c_string(name)) + end function - ! Get pointer to S(a,b) table - associate (sab => this % data(i_temp)) + subroutine sample(this, micro_xs, E_in, E_out, mu) + class(SAlphaBeta), intent(in) :: this + type(C_PTR), value :: micro_xs + real(C_DOUBLE), value :: E_in + real(C_DOUBLE), intent(out) :: E_out + real(C_DOUBLE), intent(out) :: mu - ! Get index and interpolation factor for inelastic grid - if (E < sab % inelastic_e_in(1)) then - i_grid = 1 - f = ZERO - else - i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i_grid)) / & - (sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid)) - end if - - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + & - f * sab % inelastic_sigma(i_grid + 1) - - ! Check for elastic data - if (E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P - - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO - else - i_grid = binary_search(sab % elastic_e_in, & - sab % n_elastic_e_in, E) - elastic = sab % elastic_P(i_grid) / E - end if - else - ! Determine index on elastic energy grid - if (E < sab % elastic_e_in(1)) then - i_grid = 1 - else - i_grid = binary_search(sab % elastic_e_in, & - sab % n_elastic_e_in, E) - end if - - ! Get interpolation factor for elastic grid - f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - & - sab%elastic_e_in(i_grid)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE - f) * sab % elastic_P(i_grid) + & - f * sab % elastic_P(i_grid + 1) - end if - else - ! No elastic data - elastic = ZERO - end if - end associate - - end subroutine sab_calculate_xs + call sab_sample(this % ptr, micro_xs, E_in, E_out, mu) + end subroutine + function threshold(this) + class(SAlphaBeta), intent(in) :: this + real(C_DOUBLE) :: threshold + threshold = sab_threshold(this % ptr) + end function !=============================================================================== ! FREE_MEMORY_SAB deallocates global arrays defined in this module !=============================================================================== subroutine free_memory_sab() + integer :: i n_sab_tables = 0 + do i = 1, size(sab_tables) + call sab_tables(i) % free() + end do if (allocated(sab_tables)) deallocate(sab_tables) call sab_dict % clear() end subroutine free_memory_sab diff --git a/src/thermal.cpp b/src/thermal.cpp index 176980592..0738497b5 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -577,4 +577,41 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E, } +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +ThermalScattering* +sab_from_hdf5(hid_t group, const double* temperature, int n, + int method, double tolerance, const double* minmax) +{ + // Convert temperatures to a vector + std::vector T {temperature, temperature + n}; + + // Create new object and return it + return new ThermalScattering{group, T, method, tolerance, minmax}; +} + +void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT, + int* i_temp, double* elastic, double* inelastic) +{ + data->calculate_xs(E, sqrtkT, i_temp, elastic, inelastic); +} + +void sab_free(ThermalScattering* data) { delete data; } + +bool sab_has_nuclide(ThermalScattering* data, const char* name) +{ + return data->has_nuclide(name); +} + +void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs, + double E_in, double* E_out, double* mu) +{ + int i_temp = micro_xs->index_temp_sab; + data->data_[i_temp - 1].sample(micro_xs, E_in, E_out, mu); +} + +double sab_threshold(ThermalScattering* data) { return data->threshold(); } + } // namespace openmc diff --git a/src/thermal.h b/src/thermal.h index d33465bf8..e3571a0d5 100644 --- a/src/thermal.h +++ b/src/thermal.h @@ -111,6 +111,12 @@ public: //! \return Whether table applies to the nuclide bool has_nuclide(const char* name) const; + // Sample an outgoing energy and angle + void sample(const NuclideMicroXS* micro_xs, double E_in, + double* E_out, double* mu); + + double threshold() const { return data_[0].threshold_inelastic_; } + std::string name_; //!< name of table, e.g. "c_H_in_H2O" double awr_; //!< weight of nucleus in neutron masses std::vector kTs_; //!< temperatures in [eV] (k*T) @@ -120,6 +126,22 @@ public: std::vector 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_calculate_xs(ThermalScattering* data, double E, double sqrtkT, + int* i_temp, double* elastic, double* inelastic); + void sab_free(ThermalScattering* data); + bool sab_has_nuclide(ThermalScattering* data, const char* name); + void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs, + double E_in, double* E_out, double* mu); + double sab_threshold(ThermalScattering* data); +} + } // namespace openmc #endif // OPENMC_THERMAL_SCATTERING_H