diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index dae3e74b25..b724d7f109 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -725,10 +725,12 @@ a material default temperature. ```` Element -------------------------------- -The ```` element has an accepted value of "nearest" or -"interpolation". A value of "nearest" indicates that for each cell, the nearest -temperature at which cross sections are given is to be applied, within a given -tolerance (see :ref:`temperature_tolerance`). A value of "multipole" indicates +The ```` element has an accepted value of "nearest", +"interpolation", or "multipole". A value of "nearest" indicates that for each +cell, the nearest temperature at which cross sections are given is to be +applied, within a given tolerance (see :ref:`temperature_tolerance`). A value of +"interpolation" indicates that cross sections are to be interpolated between +temperatures at which nuclear data are present. A value of "multipole" indicates that the windowed multipole method should be used to evaluate temperature-dependent cross sections in the resolved resonance range (a :ref:`windowed multipole library ` must also be available). diff --git a/src/cross_section.F90 b/src/cross_section.F90 index e2cc31d5ea..4bb323c517 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -162,15 +162,33 @@ contains ! temperature. Note that there is no tolerance here, so this ! temperature could be very far off! kT = sqrtkT**2 + i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) end if else - ! If not using multipole data, do a linear search on temperature kT = sqrtkT**2 - do i_temp = 1, size(nuclides(i_nuclide) % kTs) - if (abs(nuclides(i_nuclide) % kTs(i_temp) - kT) < & - K_BOLTZMANN*temperature_tolerance) exit - end do + + select case (temperature_method) + case (TEMPERATURE_NEAREST) + ! If using nearest temperature, do linear search on temperature + do i_temp = 1, size(nuc % kTs) + if (abs(nuc % kTs(i_temp) - kT) < K_BOLTZMANN * & + temperature_tolerance) exit + end do + case (TEMPERATURE_INTERPOLATION) + ! Find temperatures that bound the actual temperature + do i_temp = 1, size(nuc % kTs) - 1 + if (nuc % kTs(i_temp) <= kT .and. kT < nuc % kTs(i_temp + 1)) exit + end do + + ! Randomly sample between temperature i and i+1 + f = (kT - nuc % kTs(i_temp)) / & + (nuc % kTs(i_temp + 1) - nuc % kTs(i_temp)) + if (f > prn()) i_temp = i_temp + 1 + case (TEMPERATURE_MULTIPOLE) + i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1) + end select + end if ! Evaluate multipole or interpolate @@ -317,10 +335,25 @@ contains ! Determine temperature for S(a,b) table kT = sqrtkT**2 - do i_temp = 1, size(sab_tables(i_sab) % kTs) - if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < & - K_BOLTZMANN*temperature_tolerance) exit - end do + if (temperature_method == TEMPERATURE_NEAREST) then + ! If using nearest temperature, do linear search on temperature + do i_temp = 1, size(sab_tables(i_sab) % kTs) + if (abs(sab_tables(i_sab) % kTs(i_temp) - kT) < & + K_BOLTZMANN*temperature_tolerance) exit + end do + else + ! Find temperatures that bound the actual temperature + do i_temp = 1, size(sab_tables(i_sab) % kTs) - 1 + if (sab_tables(i_sab) % kTs(i_temp) <= kT .and. & + kT < sab_tables(i_sab) % kTs(i_temp + 1)) exit + end do + + ! Randomly sample between temperature i and i+1 + f = (kT - sab_tables(i_sab) % kTs(i_temp)) / & + (sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp)) + if (f > prn()) i_temp = i_temp + 1 + end if + ! Get pointer to S(a,b) table associate (sab => sab_tables(i_sab) % data(i_temp)) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index eeeae865b6..88eb0e39ca 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4948,7 +4948,7 @@ contains call names % push_back('Ga0') call densities % push_back(density) else - call names % push_back('Ha69') + call names % push_back('Ga69') call densities % push_back(density * 0.60108_8) call names % push_back('Ga71') call densities % push_back(density * 0.39892_8) @@ -5840,7 +5840,7 @@ contains file_id = file_open(libraries(i_library) % path, 'r') group_id = open_group(file_id, name) call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), & - temperature_tolerance) + temperature_method, temperature_tolerance) call close_group(group_id) call file_close(file_id) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 1fbd691ee8..6b3965917e 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -241,11 +241,13 @@ module nuclide_header 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 select case (method) case (TEMPERATURE_NEAREST) - ! Determine actual temperatures to read - TEMP_LOOP: do i = 1, temperature % size() + ! Find nearest temperatures + 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) @@ -265,11 +267,31 @@ module nuclide_header &for " // trim(this % name) // " at or near " // & trim(to_str(nint(temp_desired))) // " K.") end if - end do TEMP_LOOP + end do case (TEMPERATURE_INTERPOLATION) - ! TODO: Get bounding temperatures - call fatal_error("Temperature interpolation not yet implemented") + ! 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 case (TEMPERATURE_MULTIPOLE) ! Add first available temperature diff --git a/src/sab_header.F90 b/src/sab_header.F90 index 147643065a..21a2ab0322 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -80,10 +80,11 @@ module sab_header contains - subroutine salphabeta_from_hdf5(this, group_id, temperature, tolerance) + subroutine salphabeta_from_hdf5(this, group_id, temperature, method, tolerance) 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 integer :: i, j @@ -142,25 +143,55 @@ contains 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 - TEMP_LOOP: 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)) + 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 - 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 TEMP_LOOP + end do - ! TODO: If using interpolation, add a block to add bounding temperatures for - ! each + 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 + + case (TEMPERATURE_MULTIPOLE) + ! Add first available temperature + call temps_to_read % push_back(nint(temps_available(1))) + + end select ! Sort temperatures to read call sort(temps_to_read) @@ -301,6 +332,12 @@ contains 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)