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

Add basic windowed multipole library support

Browse source 
This commit adds the ability to load in windowed multipole
temperature-dependent cross section data.  Given a windowed multipole
file named ZAID.h5 in the same directory as the *.xml files of a
simulation, OpenMC will now load this HDF5 file.  This will replace the
resolved resonance region with windowed multipole data, and all other
pointwise data by a small pointwise library to minimize memory
footprint.

Temperatures are tenatively set by adding a <temperature> tag inside of
a material in materials.xml.  Units are in Kelvin for now.

This commit also adds the Faddeeva.cc library from
http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package which requires
OpenMC to compile both Fortran and C code.  Modifications have been done
to most of CMakeLists.txt, but not all of it.
This commit is contained in:
Colin Josey 2016-01-15 19:30:10 -05:00
parent 39d2b8945f
commit 9b3e1b86c6
14 changed files with 3589 additions and 93 deletions
Download patch file Download diff file Expand all files Collapse all files

40
CMakeLists.txt
View file

@ -1,5 +1,5 @@
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)
project(openmc Fortran)
project(openmc Fortran C)
# Setup output directories
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
@ -109,47 +109,61 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
# GNU Fortran compiler options
list(APPEND f90flags -cpp -std=f2008 -fbacktrace)
list(APPEND cflags -cpp -std=c99)
if(debug)
if(NOT (GCC_VERSION VERSION_LESS 4.7))
list(APPEND f90flags -Wall)
list(APPEND cflags -Wall)
endif()
list(APPEND f90flags -g -pedantic -fbounds-check
-ffpe-trap=invalid,overflow,underflow)
list(APPEND cflags -g -pedantic -fbounds-check
-ffpe-trap=invalid,overflow,underflow)
list(APPEND ldflags -g)
endif()
if(profile)
list(APPEND f90flags -pg)
list(APPEND cflags -pg)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(APPEND f90flags -O3)
list(APPEND cflags -O3)
endif()
if(openmp)
list(APPEND f90flags -fopenmp)
list(APPEND cflags -fopenmp)
list(APPEND ldflags -fopenmp)
endif()
if(coverage)
list(APPEND f90flags -coverage)
list(APPEND cflags -coverage)
list(APPEND ldflags -coverage)
endif()
elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel)
# Intel Fortran compiler options
list(APPEND f90flags -fpp -std08 -assume byterecl -traceback)
list(APPEND cflags -std=c99)
if(debug)
list(APPEND f90flags -g -warn -ftrapuv -fp-stack-check
"-check all" -fpe0)
list(APPEND cflags -g -warn -ftrapuv -fp-stack-check
"-check all" -fpe0)
list(APPEND ldflags -g)
endif()
if(profile)
list(APPEND f90flags -pg)
list(APPEND cflags -pg)
list(APPEND ldflags -pg)
endif()
if(optimize)
list(APPEND f90flags -O3)
list(APPEND cflags -O3)
endif()
if(openmp)
list(APPEND f90flags -openmp)
list(APPEND cflags -openmp)
list(APPEND ldflags -openmp)
endif()
@ -248,7 +262,7 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
#===============================================================================
set(program "openmc")
file(GLOB source src/*.F90 src/xml/openmc_fox.F90)
file(GLOB source src/Faddeeva.c src/*.F90 src/xml/openmc_fox.F90)
add_executable(${program} ${source})
# target_include_directories was added in CMake 2.8.11 and is the recommended
@ -260,15 +274,19 @@ else()
endif()
# target_compile_options was added in CMake 2.8.12 and is the recommended way to
# set compile flags. Note that this sets the COMPILE_OPTIONS property (also
# available only in 2.8.12+) rather than the COMPILE_FLAGS property, which is
# deprecated. The former can handle lists whereas the latter cannot.
if(CMAKE_VERSION VERSION_LESS 4.8.12)
string(REPLACE ";" " " f90flags "${f90flags}")
set_property(TARGET ${program} PROPERTY COMPILE_FLAGS "${f90flags}")
else()
target_compile_options(${program} PUBLIC ${f90flags})
endif()
# set compile flag, but it doesn't seem to work with both Fortran and C. In
# theory the commented code below might work but I couldn't get it running.
# Maybe it requires a newer version of CMake?
#target_compile_options(${program} PUBLIC
# $<$<COMPILE_LANGUAGE:Fortran>:${f90flags}>
# $<$<COMPILE_LANGUAGE:C>:${cflags}>)
# This is not the recommended method, but it is a functional method. Convert
# the semicolon-delineated lists into space-delineated ones then set them as the
# compiler flags.
string(REPLACE ";" " " f90flags "${f90flags}")
string(REPLACE ";" " " cflags "${cflags}")
set(CMAKE_C_FLAGS "${cflags}")
set(CMAKE_Fortran_FLAGS "${f90flags}")
# Add HDF5 library directories to link line with -L
foreach(LIBDIR ${HDF5_LIBRARY_DIRS})

3
src/Faddeeva.c Normal file
View file

@ -0,0 +1,3 @@
/* The Faddeeva.cc file contains macros to let it compile as C code
(assuming C99 complex-number support), so just #include it. */
#include "Faddeeva.cc"

2516
src/Faddeeva.cc Normal file
View file

File diff suppressed because it is too large Load diff

68
src/Faddeeva.h Normal file
View file

@ -0,0 +1,68 @@
/* Copyright (c) 2012 Massachusetts Institute of Technology
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
/* Available at: http://ab-initio.mit.edu/Faddeeva
Header file for Faddeeva.c; see Faddeeva.cc for more information. */
#ifndef FADDEEVA_H
#define FADDEEVA_H 1
// Require C99 complex-number support
#include <complex.h>
#ifdef __cplusplus
extern "C"
{
#endif /* __cplusplus */
// compute w(z) = exp(-z^2) erfc(-iz) [ Faddeeva / scaled complex error func ]
extern double complex Faddeeva_w(double complex z,double relerr);
extern double Faddeeva_w_im(double x); // special-case code for Im[w(x)] of real x
// Various functions that we can compute with the help of w(z)
// compute erfcx(z) = exp(z^2) erfc(z)
extern double complex Faddeeva_erfcx(double complex z, double relerr);
extern double Faddeeva_erfcx_re(double x); // special case for real x
// compute erf(z), the error function of complex arguments
extern double complex Faddeeva_erf(double complex z, double relerr);
extern double Faddeeva_erf_re(double x); // special case for real x
// compute erfi(z) = -i erf(iz), the imaginary error function
extern double complex Faddeeva_erfi(double complex z, double relerr);
extern double Faddeeva_erfi_re(double x); // special case for real x
// compute erfc(z) = 1 - erf(z), the complementary error function
extern double complex Faddeeva_erfc(double complex z, double relerr);
extern double Faddeeva_erfc_re(double x); // special case for real x
// compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)
extern double complex Faddeeva_Dawson(double complex z, double relerr);
extern double Faddeeva_Dawson_re(double x); // special case for real x
#ifdef __cplusplus
}
#endif /* __cplusplus */
#endif // FADDEEVA_H

68
src/ace.F90
View file

@ -9,6 +9,8 @@ module ace
use global
use list_header, only: ListInt
use material_header, only: Material
use multipole, only: multipole_read
use multipole_header, only: max_L, max_poles, max_poly
use output, only: write_message
use set_header, only: SetChar
use string, only: to_str, to_lower
@ -105,6 +107,9 @@ contains
end do
end if
! Read multipole file into the appropriate entry on the nuclides array
call read_multipole_data(i_nuclide)
! Add name and alias to dictionary
call already_read % add(name)
call already_read % add(alias)
@ -408,6 +413,69 @@ contains
end subroutine read_ace_table
!===============================================================================
! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the
! directory and loads it using multipole_read
!===============================================================================
subroutine read_multipole_data(i_table)
integer, intent(in) :: i_table ! index in nuclides/sab_tables
logical :: file_exists ! does multipole library exist?
character(7) :: readable ! is multipole library readable?
character(6) :: zaid_string ! String of the ZAID
character(9) :: filename ! path to multipole cross section library
type(Nuclide), pointer :: nuc => null()
! For the time being, and I know this is a bit hacky, we just assume
! that the file will be zaid.h5.
nuc => nuclides(i_table)
write(zaid_string,'(I6.6)') nuc % zaid
filename = zaid_string // ".h5"
! Check if Multipole library exists and is readable
inquire(FILE=filename, EXIST=file_exists, READ=readable)
if (.not. file_exists) then
nuc % mp_present = .FALSE.
return
elseif (readable(1:3) == 'NO') then
call fatal_error("Multipole library '" // trim(filename) // "' is not readable! &
&Change file permissions with chmod command.")
end if
! display message
call write_message("Loading Multipole XS table: " // filename, 6)
allocate(nuc % multipole)
! Call the read routine
call multipole_read(filename, nuc % multipole, i_table)
nuc % mp_present = .TRUE.
! Update the maximum number of poles, l indices, and polynomial order
if(nuc % multipole % max_w > max_poles) then
max_poles = nuc % multipole % max_w
end if
if(nuc % multipole % num_l > max_L) then
max_L = nuc % multipole % num_l
end if
if(nuc % multipole % fit_order + 1 > max_poly) then
max_poly = nuc % multipole % fit_order + 1
end if
! Recreate nu-fission tables
if(nuc % fissionable) then
call generate_nu_fission(nuc)
end if
if(associated(nuc)) nullify(nuc)
end subroutine read_multipole_data
!===============================================================================
! READ_ESZ - reads through the ESZ block. This block contains the energy grid,
! total xs, absorption xs, elastic scattering xs, and heating numbers.

21
src/ace_header.F90
View file

@ -1,9 +1,10 @@
module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use stl_vector, only: VectorInt
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use multipole_header, only: MultipoleArray
use stl_vector, only: VectorInt
implicit none
@ -144,6 +145,10 @@ module ace_header
integer :: urr_inelastic
type(UrrData), pointer :: urr_data => null()
! Multipole data
logical :: mp_present
type(MultipoleArray), pointer :: multipole => null()
! Reactions
integer :: n_reaction ! # of reactions
type(Reaction), allocatable :: reactions(:)
@ -267,6 +272,9 @@ module ace_header
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
real(8) :: last_prn
! Information for Doppler broadening
real(8) :: last_sqrtkT = ZERO ! last temperature in sqrt(Boltzmann constant * temperature (MeV))
end type NuclideMicroXS
!===============================================================================
@ -338,6 +346,11 @@ module ace_header
deallocate(this % urr_data)
end if
if (associated(this % multipole)) then
call this % multipole % clear()
deallocate(this % multipole)
end if
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()

2
src/constants.F90
View file

@ -63,6 +63,7 @@ module constants
real(8), parameter :: &
PI = 3.1415926535898_8, & ! pi
SQRT_PI = 1.7724538509055_8, & ! square root of pi
MASS_NEUTRON = 1.008664916_8, & ! mass of a neutron in amu
MASS_NEUTRON_MEV = 939.565379_8, & ! mass of a neutron in MeV/c^2
MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu
@ -77,6 +78,7 @@ module constants
TWO = 2.0_8, &
THREE = 3.0_8, &
FOUR = 4.0_8
complex(8), parameter :: ONEI = (ZERO, ONE)
! ============================================================================
! GEOMETRY-RELATED CONSTANTS

381
src/cross_section.F90
View file

@ -1,19 +1,31 @@
module cross_section
use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData
use ace_header, only: Nuclide, SAlphaBeta, Reaction, UrrData
use constants
use energy_grid, only: grid_method, log_spacing
use error, only: fatal_error
use fission, only: nu_total
use energy_grid, only: grid_method, log_spacing
use error, only: fatal_error
use fission, only: nu_total
use global
use list_header, only: ListElemInt
use material_header, only: Material
use particle_header, only: Particle
use random_lcg, only: prn
use search, only: binary_search
use list_header, only: ListElemInt
use material_header, only: Material
use math, only: w, broaden_n_polynomials
use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, &
MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A, FIT_F, &
MultipoleArray, max_poly, max_L, max_poles
use particle_header, only: Particle
use random_lcg, only: prn
use search, only: binary_search
implicit none
! Allocatable arrays for multipole that are allocated once for speed purposes
complex(8), allocatable :: sigT_factor(:)
real(8), allocatable :: twophi(:)
real(8), allocatable :: broadened_polynomials(:)
logical :: mp_already_alloc = .FALSE.
!$omp threadprivate(sigT_factor, twophi, broadened_polynomials, mp_already_alloc)
contains
!===============================================================================
@ -93,10 +105,10 @@ contains
i_nuclide = mat % nuclide(i)
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
if (p % E /= micro_xs(i_nuclide) % last_E .or. mat % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT)
end if
! ========================================================================
@ -133,7 +145,7 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union, sqrtkT)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
@ -141,11 +153,13 @@ contains
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
real(8) :: f ! interp factor on nuclide energy grid
real(8) :: sigT, sigA, sigF ! Intermediate multipole variables
type(Nuclide), pointer :: nuc
type(Material), pointer :: mat
@ -153,53 +167,111 @@ contains
nuc => nuclides(i_nuclide)
mat => materials(i_mat)
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_MAT_UNION)
! If MP, don't interpolate, it's all already baked in.
if( nuc % mp_present .AND. &
(E >= nuc % multipole % start_E/1.0D6 .AND.&
E <= nuc % multipole % end_E/1.0D6)) then
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
! Call multipole kernel
call multipole_eval(nuc % multipole, E, sqrtkT, sigT, sigA, sigF)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
micro_xs(i_nuclide) % total = sigT
micro_xs(i_nuclide) % absorption = sigA
micro_xs(i_nuclide) % elastic = sigT - sigA
if (E < nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
if (nuc % fissionable) then
micro_xs(i_nuclide) % fission = sigF
micro_xs(i_nuclide) % nu_fission = sigF * nu_total(nuc,E)
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
i_low = nuc % grid_index(i_log_union)
i_high = nuc % grid_index(i_log_union + 1) + 1
! Perform binary search over reduced range
i_grid = binary_search(nuc % energy(i_low:i_high), &
i_high - i_low + 1, E) + i_low - 1
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
end if
case (GRID_NUCLIDE)
! Perform binary search on the nuclide energy grid in order to determine
! which points to interpolate between
! Ensure these values are set
! Note, the only time either is used is in one of 4 places:
! 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. cross_section.F90 - 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.
micro_xs(i_nuclide) % index_grid = 0
micro_xs(i_nuclide) % interp_factor = ZERO
else
! Determine index on nuclide energy grid
select case (grid_method)
case (GRID_MAT_UNION)
if (E <= nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
! the energy range over which a binary search needs to be performed
if (E < nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
i_low = nuc % grid_index(i_log_union)
i_high = nuc % grid_index(i_log_union + 1) + 1
! Perform binary search over reduced range
i_grid = binary_search(nuc % energy(i_low:i_high), &
i_high - i_low + 1, E) + i_low - 1
end if
case (GRID_NUCLIDE)
! Perform binary search on the nuclide energy grid in order to determine
! which points to interpolate between
if (E <= nuc % energy(1)) then
i_grid = 1
elseif (E > nuc % energy(nuc % n_grid)) then
i_grid = nuc % n_grid - 1
else
i_grid = binary_search(nuc % energy, nuc % n_grid, E)
end if
end select
! check for rare case where two energy points are the same
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
! calculate interpolation factor
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
micro_xs(i_nuclide) % index_grid = i_grid
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
+ f * nuc % total(i_grid+1)
! Calculate microscopic nuclide elastic cross section
micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) &
+ f * nuc % elastic(i_grid+1)
! Calculate microscopic nuclide absorption cross section
micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( &
i_grid) + f * nuc % absorption(i_grid+1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) &
+ f * nuc % fission(i_grid+1)
! Calculate microscopic nuclide nu-fission cross section
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
i_grid) + f * nuc % nu_fission(i_grid+1)
end if
end select
! check for rare case where two energy points are the same
if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1
! calculate interpolation factor
f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid))
micro_xs(i_nuclide) % index_grid = i_grid
micro_xs(i_nuclide) % interp_factor = f
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
@ -208,32 +280,6 @@ contains
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
+ f * nuc % total(i_grid+1)
! Calculate microscopic nuclide elastic cross section
micro_xs(i_nuclide) % elastic = (ONE - f) * nuc % elastic(i_grid) &
+ f * nuc % elastic(i_grid+1)
! Calculate microscopic nuclide absorption cross section
micro_xs(i_nuclide) % absorption = (ONE - f) * nuc % absorption( &
i_grid) + f * nuc % absorption(i_grid+1)
if (nuc % fissionable) then
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % fission = (ONE - f) * nuc % fission(i_grid) &
+ f * nuc % fission(i_grid+1)
! Calculate microscopic nuclide nu-fission cross section
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
i_grid) + f * nuc % nu_fission(i_grid+1)
end if
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
@ -525,6 +571,185 @@ contains
end function find_energy_index
!===============================================================================
! MULTIPOLE_EVAL_ALLOCATE allocates fixed-length arrays that vary based on
! what nuclides are loaded into the problem
!===============================================================================
subroutine multipole_eval_allocate()
allocate(sigT_factor(max_L))
allocate(twophi(max_L))
allocate(broadened_polynomials(max_poly))
mp_already_alloc = .TRUE.
end subroutine
!===============================================================================
! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross
! sections in the resolved resonance regions
!===============================================================================
subroutine multipole_eval(multipole, Emev, sqrtkT, sigT, sigA, sigF)
type(MultipoleArray), intent(in) :: multipole !< The windowed multipole object to process.
real(8), intent(in) :: Emev !< The energy at which to evaluate the cross section in MeV
real(8), intent(in) :: sqrtkT !< The temperature in the form sqrt(kT (in eV)), at which to evaluate the cross section.
real(8), intent(out) :: sigT !< Total cross section
real(8), intent(out) :: sigA !< Absorption cross section
real(8), intent(out) :: sigF !< Fission cross section
complex(8) :: PSIIKI
complex(8) :: CDUM1
complex(8) :: w_val
complex(8) :: Z
real(8) :: sqrtE
real(8) :: invE
real(8) :: DOPP
real(8) :: DOPP_ECOEF
real(8) :: temp
real(8) :: E
integer :: iP
integer :: iC
integer :: iW
integer :: startw
integer :: startw_1
integer :: startw_endw
integer :: endw
! Convert to eV
E = Emev * 1.0D6
sqrtE = sqrt(E)
invE = E**(-1)
if(mp_already_alloc .eqv. .FALSE.) then
call multipole_eval_allocate()
end if
! Locate us
iW = floor((sqrtE - sqrt(multipole % start_E))/multipole % spacing + 1.0_8)
startw = multipole % w_start(iW)
startw_1 = startw - 1 ! This is an index shift parameter.
endw = multipole % w_end(iW)
startw_endw = endw - startw + 1
! Fill in factors
if (startw <= endw) then
call fill_factors(multipole, sqrtE, sigT_factor, twophi, multipole % num_l)
end if
! Generate some doppler broadening parameters
! DOPP_ECOEF is inverse of dopp, divided by E, multiplied by sqrt(pi).
DOPP = multipole % sqrtAWR/sqrtKT
DOPP_ECOEF = DOPP*invE*SQRT_PI
sigT = 0.0_8
sigA = 0.0_8
sigF = 0.0_8
! Evaluate linefit first
if(sqrtkT /= 0 .AND. multipole % broaden_poly(iW) == 1) then ! Broaden the curvefit.
call broaden_n_polynomials(E, DOPP, multipole % fit_order + 1, broadened_polynomials)
do iC = 1,multipole % fit_order+1
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*broadened_polynomials(iC)
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*broadened_polynomials(iC)
if (multipole % fissionable) then
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*broadened_polynomials(iC)
end if
end do
else ! Evaluate as if it were a polynomial
temp = invE
do iC = 1,multipole % fit_order+1
sigT = sigT + multipole % curvefit(FIT_T, iC, iW)*temp
sigA = sigA + multipole % curvefit(FIT_A, iC, iW)*temp
if (multipole % fissionable) then
sigF = sigF + multipole % curvefit(FIT_F, iC, iW)*temp
end if
temp = temp * sqrtE
end do
end if
! Then get the poles we want and broaden them.
if (sqrtkT == 0.0_8) then
! If at 0K, use asymptotic form.
do iP = startw, endw
PSIIKI = -ONEI/(multipole % data(MP_EA, iP) - sqrtE)
CDUM1 = PSIIKI/E
if (multipole % formalism == FORM_MLBW) then
sigT = sigT + real(multipole % data(MLBW_RT, iP)*CDUM1* &
sigT_factor(multipole % l_value(iP))) &
+ real(multipole % data(MLBW_RX, iP)*CDUM1)
sigA = sigA + real(multipole % data(MLBW_RA, iP)*CDUM1)
sigF = sigF + real(multipole % data(MLBW_RF, iP)*CDUM1)
else if (multipole % formalism == FORM_RM) then
sigT = sigT + real(multipole % data(RM_RT, iP)*CDUM1* &
sigT_factor(multipole % l_value(iP)))
sigA = sigA + real(multipole % data(RM_RA, iP)*CDUM1)
sigF = sigF + real(multipole % data(RM_RF, iP)*CDUM1)
end if
end do
else
! At temperature, use Faddeeva function-based form.
if(endw >= startw) then
do iP = startw, endw
Z = (sqrtE - multipole % data(MP_EA, iP))*DOPP
call w(Z, w_val)
w_val = w_val*DOPP_ECOEF
if (multipole % formalism == FORM_MLBW) then
sigT = sigT + real((multipole % data(MLBW_RT, iP)* &
sigT_factor(multipole%l_value(iP)) + &
multipole % data(MLBW_RX, iP))*w_val)
sigA = sigA + real(multipole % data(MLBW_RA, iP)*w_val)
sigF = sigF + real(multipole % data(MLBW_RF, iP)*w_val)
else if (multipole % formalism == FORM_RM) then
sigT = sigT + real(multipole % data(RM_RT, iP)*w_val* &
sigT_factor(multipole % l_value(iP)))
sigA = sigA + real(multipole % data(RM_RA, iP)*w_val)
sigF = sigF + real(multipole % data(RM_RF, iP)*w_val)
end if
end do
end if
end if
end subroutine
!===============================================================================
! FILL_FACTORS calculates the value of phi, the hardsphere phase shift factor,
! and sigT_factor, a factor inside of the sigT equation not present in the
! sigA and sigF equations.
!===============================================================================
subroutine fill_factors(multipole, sqrtE, sigT_factor, twophi, max_L)
type(MultipoleArray), intent(in) :: multipole
real(8), intent(in) :: sqrtE
integer, intent(in) :: max_L
complex(8), intent(out) :: sigT_factor(max_L)
real(8), intent(out) :: twophi(max_L)
integer :: iL
real(8) :: arg
do iL = 1, max_L
twophi(iL) = multipole%pseudo_k0RS(iL)*sqrtE
if (iL == 2) then
twophi(iL) = twophi(iL) - atan(twophi(iL))
else if (iL == 3) then
arg = 3.0_8*twophi(iL)/(3.0_8-twophi(iL)**2)
twophi(iL) = twophi(iL) - atan(arg)
else if (iL == 4) then
arg = twophi(iL)*(15.0_8-twophi(iL)**2)/(15.0_8-6.0_8*twophi(iL)**2)
twophi(iL) = twophi(iL) - atan(arg)
end if
end do
twophi = 2.0_8 * twophi
sigT_factor = cmplx(cos(twophi),-sin(twophi), KIND=8)
end subroutine
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering

90
src/hdf5_interface.F90
View file

@ -62,6 +62,7 @@ module hdf5_interface
module procedure read_string_1D
module procedure read_tally_result_1D
module procedure read_tally_result_2D
module procedure read_complex_2D
end interface read_dataset
public :: write_dataset
@ -1921,4 +1922,93 @@ contains
mpio = (driver == H5FD_MPIO_F)
end function using_mpio_device
!===============================================================================
! READ_COMPLEX_2D reads double precision complex 2-D array data as output by
! the h5py HDF5 python module.
!===============================================================================
subroutine read_complex_2D(group_id, name, buffer, indep)
integer(HID_T), intent(in) :: group_id
character(*), intent(in) :: name ! name of data
complex(8), intent(inout), target :: buffer(:,:) ! data to write
logical, intent(in), optional :: indep ! independent I/O
integer(HSIZE_T) :: dims(2)
dims(:) = shape(buffer)
if (present(indep)) then
call read_complex_2D_explicit(group_id, dims, name, buffer, indep)
else
call read_complex_2D_explicit(group_id, dims, name, buffer)
end if
end subroutine read_complex_2D
subroutine read_complex_2D_explicit(group_id, dims, name, buffer, indep)
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T), intent(in) :: dims(2)
character(*), intent(in) :: name ! name of data
complex(8), intent(inout), target :: buffer(dims(1),dims(2))
logical, intent(in), optional :: indep ! independent I/O
real(8), target :: buffer_r(dims(1), dims(2))
real(8), target :: buffer_i(dims(1), dims(2))
integer(HSIZE_T) :: i, j
integer :: hdf5_err
integer :: data_xfer_mode
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
integer(HID_T) :: dset ! data set handle
type(c_ptr) :: f_ptr_r, f_ptr_i
! Components needed for complex type support
integer(HID_T) :: dtype_real
integer(HID_T) :: dtype_imag
integer(SIZE_T) :: size_double
integer :: error
! Create the complex type
call h5tget_size_f(H5T_NATIVE_DOUBLE, size_double, error)
! Insert the 'r' and 'i' identifiers
call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_real, error)
call h5tcreate_f(H5T_COMPOUND_F, size_double, dtype_imag, error)
call h5tinsert_f(dtype_real, "r", 0_8, H5T_NATIVE_DOUBLE, error)
call h5tinsert_f(dtype_imag, "i", 0_8, H5T_NATIVE_DOUBLE, error)
! Set up collective vs. independent I/O
data_xfer_mode = H5FD_MPIO_COLLECTIVE_F
if (present(indep)) then
if (indep) data_xfer_mode = H5FD_MPIO_INDEPENDENT_F
end if
call h5dopen_f(group_id, trim(name), dset, hdf5_err)
f_ptr_r = c_loc(buffer_r)
f_ptr_i = c_loc(buffer_i)
if (using_mpio_device(group_id)) then
#ifdef PHDF5
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, data_xfer_mode, hdf5_err)
call h5dread_f(dset, dtype_real, f_ptr_r, hdf5_err, xfer_prp=plist)
call h5dread_f(dset, dtype_imag, f_ptr_i, hdf5_err, xfer_prp=plist)
call h5pclose_f(plist, hdf5_err)
#endif
else
call h5dread_f(dset, dtype_real, f_ptr_r, hdf5_err)
call h5dread_f(dset, dtype_imag, f_ptr_i, hdf5_err)
end if
! Reconstitute the complex numbers
do i = 1,dims(1)
do j = 1,dims(2)
buffer(i,j) = cmplx(buffer_r(i,j), buffer_i(i,j), kind=8)
end do
end do
call h5dclose_f(dset, hdf5_err)
end subroutine read_complex_2D_explicit
end module hdf5_interface

39
src/input_xml.F90
View file

@ -1854,6 +1854,7 @@ contains
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
logical :: temp_known ! Has the temperature yet been defined?
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
character(MAX_WORD_LEN) :: units ! units on density
@ -1866,6 +1867,7 @@ contains
type(Node), pointer :: doc => null()
type(Node), pointer :: node_mat => null()
type(Node), pointer :: node_dens => null()
type(Node), pointer :: node_temp => null()
type(Node), pointer :: node_nuc => null()
type(Node), pointer :: node_ele => null()
type(Node), pointer :: node_sab => null()
@ -1936,6 +1938,17 @@ contains
cycle
end if
! =======================================================================
! READ AND PARSE <temperature> TAG
if (check_for_node(node_mat, "temperature")) then
call get_node_ptr(node_mat, "temperature", node_temp)
call get_node_value(node_temp, "value", temp_dble)
mat % sqrtkT = sqrt(temp_dble * K_BOLTZMANN * 1.0D6)
temp_known = .TRUE.
else
temp_known = .FALSE.
end if
! =======================================================================
! READ AND PARSE <density> TAG
@ -2042,6 +2055,16 @@ contains
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
! If needed, look up temperature
if (temp_known .eqv. .FALSE.) then
! Find xs_listing and set the name/alias according to the listing
index_list = xs_listing_dict % get_key(to_lower(name))
if(xs_listings(index_list) % kT /= 0.0_8) then
mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6)
temp_known = .TRUE.
end if
end if
! save name and density to list
call list_names % append(name)
@ -2130,6 +2153,16 @@ contains
temp_str = "data"
end if
! If still needed, look up temperature
if (temp_known .eqv. .FALSE.) then
! Find xs_listing and set kT
index_list = xs_listing_dict % get_key(to_lower(list_names % tail % data))
if(xs_listings(index_list) % kT /= 0.0_8) then
mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6)
temp_known = .TRUE.
end if
end if
! Set ace or iso-in-lab scattering for each nuclide in element
do k = 1, n_nuc_ele
if (adjustl(to_lower(temp_str)) == "iso-in-lab") then
@ -2144,6 +2177,12 @@ contains
end do NATURAL_ELEMENTS
! If still undefined, set the temperature to zero
if (temp_known .eqv. .FALSE.) then
mat % sqrtkT = 0.0_8
temp_known = .TRUE.
end if
! ========================================================================
! COPY NUCLIDES TO ARRAYS IN MATERIAL

1
src/material_header.F90
View file

@ -13,6 +13,7 @@ module material_header
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8) :: sqrtkT ! sqrt(kT), kT in eV
! Energy grid information
integer :: n_grid ! # of union material grid points

102
src/math.F90
View file

@ -2,9 +2,25 @@ module math
use constants
use random_lcg, only: prn
use ISO_C_BINDING
implicit none
!===============================================================================
! FADDEEVA_W evaluates the scaled complementary error function. This
! interfaces with the MIT C library
!===============================================================================
interface
COMPLEX (C_DOUBLE_COMPLEX) FUNCTION faddeeva_w &
(Z, RELERR) BIND(C, NAME='Faddeeva_w')
use ISO_C_BINDING
implicit none
COMPLEX(C_DOUBLE_COMPLEX), value :: Z
REAL(C_DOUBLE) , value :: RELERR
end function faddeeva_w
end interface
contains
!===============================================================================
@ -659,4 +675,90 @@ contains
end function watt_spectrum
!===============================================================================
! W acts as a front end to the MIT Faddeeva function, Faddeeva_w.
!===============================================================================
subroutine w(Z,wv)
use ISO_C_BINDING
complex(C_DOUBLE_COMPLEX), intent(inout) :: Z ! The point to evaluate Z at
complex(8), intent(out) :: wv ! The resulting w(Z) value
real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT Faddeeva
logical :: mangle ! Do we need to perform the special mangling step?
! There is some special mangling the true W function in WHOPPER does. I
! replicate it here.
if(aimag(Z) < 0) then
Z = conjg(Z)
mangle = .TRUE.
else
mangle = .FALSE.
end if
! Calculate Faddeeva
relerr = 0.0_8
wv = faddeeva_w(Z,relerr)
! Finish mangling the results
if(mangle .eqv. .TRUE.) then
wv = -conjg(wv)
end if
end subroutine
!===============================================================================
! BROADEN_N_POLYNOMIALS doppler broadens polynomials of the form
! a/En + b/sqrt(En) + c + d sqrt(En) ...
! exactly and quickly.
!===============================================================================
subroutine broaden_n_polynomials(En, DOPP, n, factors)
real(8), intent(in) :: En ! Energy to evaluate at
real(8), intent(in) :: DOPP ! sqrt(atomic weight ratio / kT), kT given in eV.
integer, intent(in) :: n ! number of components to polynomial
real(8), intent(out):: factors(n) ! output leading coefficient
integer :: i
real(8) :: sqrtE
real(8) :: beta
real(8) :: halfinvDOPP2
real(8) :: quarterinvDOPP4
real(8) :: erfbeta
real(8) :: exp_m_beta2
sqrtE = sqrt(En)
beta = sqrtE*DOPP
halfinvDOPP2 = 0.5_8/DOPP**2
quarterinvDOPP4 = 0.25_8/DOPP**4
if (beta > 6.0_8) then
! Save time, ERF(6) is 1 to machine precision.
! beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
erfBeta = 1.0_8
exp_m_beta2 = 0.0_8
else
erfBeta = erf(beta)
exp_m_beta2 = exp(-beta**2)
end if
! Assume that, for sure, we'll use a second order (1/E, 1/V, const)
! fit, and no less.
factors(1) = erfbeta/En
factors(2) = 1.0_8/sqrtE
factors(3) = factors(1)*(halfinvDOPP2 + En) + exp_m_beta2/(beta*SQRT_PI)
! Perform recursive broadening of high order components
do i = 1,n-3
if (i /= 1) then
factors(i+3) = -factors(i-1)*(i - 1.0_8)*i*quarterinvDOPP4 + &
factors(i+1)*(En + (1.0_8 + 2.0_8*i)*halfinvDOPP2)
else
! Although it's mathematically identical, factors(0) will contain
! nothing, and we don't want to have to worry about memory.
factors(i+3) = factors(i+1)*(En + (1.0_8 + 2.0_8*i)*halfinvDOPP2)
end if
end do
end subroutine broaden_n_polynomials
end module math

199
src/multipole.F90 Normal file
View file

@ -0,0 +1,199 @@
module multipole
use constants
use global
use hdf5
use hdf5_interface
use multipole_header, only: MultipoleArray, FIT_T, FIT_A, FIT_F, max_L, &
max_poles, max_poly, MP_FISS, FORM_MLBW, FORM_RM
implicit none
contains
!---------------------------------------------------------------------------
! MULTIPOLE_READ Reads in a multipole HDF5 file with the original API
! specification. Subject to change as the library format matures.
!---------------------------------------------------------------------------
subroutine multipole_read(filename, multipole, i_table)
character(len=*), intent(in) :: filename ! Filename of the multipole library to load
type(MultipoleArray), intent(out), target :: multipole ! The object to fill
integer, intent(in) :: i_table ! index in nuclides/sab_tables
type(Nuclide), pointer :: nuc => null()
integer(HID_T) :: file_id
integer(HID_T) :: group_id
! Intermediate loading components
integer :: NMT
integer :: i, j, k
integer, allocatable :: MT(:)
logical :: accumulated_fission
character(len=3) :: MT_string
character(len=24) :: MT_n ! Takes the form '/nuclide/reactions/MT???'
integer :: is_fissionable
nuc => nuclides(i_table)
! Open file for reading and move into the /isotope group
file_id = file_open(filename, 'r', parallel=.true.)
group_id = open_group(file_id, "/nuclide")
! Load in all the array size scalars
call read_dataset(group_id, "length", multipole % length)
call read_dataset(group_id, "windows", multipole % windows)
call read_dataset(group_id, "num_l", multipole % num_l)
call read_dataset(group_id, "fit_order", multipole % fit_order)
call read_dataset(group_id, "max_w", multipole % max_w)
call read_dataset(group_id, "fissionable", is_fissionable)
if (is_fissionable == MP_FISS) then
multipole % fissionable = .true.
else
multipole % fissionable = .false.
end if
call read_dataset(group_id, "formalism", multipole % formalism)
call read_dataset(group_id, "spacing", multipole % spacing)
call read_dataset(group_id, "sqrtAWR", multipole % sqrtAWR)
call read_dataset(group_id, "start_E", multipole % start_E)
call read_dataset(group_id, "end_E", multipole % end_E)
! Allocate the multipole array components
call multipole % allocate()
! Read in arrays
call read_dataset(group_id, "data", multipole % data)
call read_dataset(group_id, "pseudo_K0RS", multipole % pseudo_k0RS)
call read_dataset(group_id, "l_value", multipole % l_value)
call read_dataset(group_id, "w_start", multipole % w_start)
call read_dataset(group_id, "w_end", multipole % w_end)
call read_dataset(group_id, "broaden_poly", multipole % broaden_poly)
call read_dataset(group_id, "curvefit", multipole % curvefit)
! Delete ACE pointwise data
call read_dataset(group_id, "n_grid", nuc % n_grid)
deallocate(nuc % energy)
deallocate(nuc % total)
deallocate(nuc % elastic)
deallocate(nuc % fission)
deallocate(nuc % nu_fission)
deallocate(nuc % absorption)
allocate(nuc % energy(nuc % n_grid))
allocate(nuc % total(nuc % n_grid))
allocate(nuc % elastic(nuc % n_grid))
allocate(nuc % fission(nuc % n_grid))
allocate(nuc % nu_fission(nuc % n_grid))
allocate(nuc % absorption(nuc % n_grid))
nuc % total = ZERO
nuc % absorption = ZERO
nuc % fission = ZERO
! Read in new energy axis (converting eV to MeV)
call read_dataset(group_id, "energy_points", nuc % energy)
nuc % energy = nuc % energy / 1.0D6
! Get count and list of MT tables
call read_dataset(group_id, "MT_count", NMT)
allocate(MT(NMT))
call read_dataset(group_id, "MT_list", MT)
call close_group(group_id)
accumulated_fission = .false.
! Loop over each MT entry and load it into a reaction.
do i = 1, NMT
write(MT_string, '(I3.3)') MT(i)
MT_n = "/nuclide/reactions/MT" // MT_string
group_id = open_group(file_id, MT_n)
! Each MT needs to be treated slightly differently.
select case (MT(i))
case(ELASTIC)
call read_dataset(group_id, "MT_sigma", nuc % elastic)
nuc % total = nuc % total + nuc % elastic
case(N_FISSION)
call read_dataset(group_id, "MT_sigma", nuc % fission)
nuc % total = nuc % total + nuc % fission
nuc % absorption = nuc % absorption + nuc % fission
accumulated_fission = .true.
case default
! Search through all of our secondary reactions
do j = 1, nuc % n_reaction
if (nuc % reactions(j) % MT == MT(i)) then
! Match found
! Individual Fission components exist, so remove the combined
! fission cross section.
if ( (MT(i) == N_F .OR. MT(i) == N_NF .OR. MT(i) == N_2NF &
.OR. MT(i) == N_3NF) .AND. (accumulated_fission .eqv. .TRUE.)) then
nuc % total = nuc % total - nuc % fission
nuc % absorption = nuc % absorption - nuc % fission
nuc % fission = 0.0_8
accumulated_fission = .FALSE.
end if
deallocate(nuc % reactions(j) % sigma)
allocate(nuc % reactions(j) % sigma(nuc % n_grid))
call read_dataset(group_id, "MT_sigma", nuc % reactions(j) % sigma)
call read_dataset(group_id, "Q_value", nuc % reactions(j) % Q_value)
call read_dataset(group_id, "threshold", nuc % reactions(j) % threshold)
nuc % reactions(j) % threshold = 1 ! TODO: reconsider implications.
nuc % reactions(j) % Q_value = nuc % reactions(j) % Q_value / 1.0D6
! Update edist values TODO: does this do anything when
! the version of data is identical?
if (associated(nuc % reactions(j) % edist)) then
if (nuc % reactions(j) % edist % law == 3) then
! Search for first XS /= 0 for the threshold.
searchQ: do k = 1, nuc % n_grid
if (nuc % reactions(j) % sigma(k) /= ZERO) then
if (k /= 1) then
nuc % reactions(j) % edist % data(1) = nuc % energy(k-1)
end if
exit searchQ
end if
end do searchQ
end if
end if
! Accumulate total
if (MT(i) /= N_LEVEL .AND. MT(i) <= N_DA) then
nuc % total = nuc % total + nuc % reactions(j) % sigma
end if
! Accumulate absorption
if (MT(i) >= N_GAMMA .and. MT(i) <= N_DA) then
nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
end if
! Accumulate fission (if needed)
if ( (MT(i) == N_F .OR. MT(i) == N_NF .OR. MT(i) == N_2NF &
.OR. MT(i) == N_3NF) ) then
nuc % fission = nuc % fission + nuc % reactions(j) % sigma
nuc % absorption = nuc % absorption + nuc % reactions(j) % sigma
end if
end if
end do
end select
call close_group(group_id)
end do
! Close file
call file_close(file_id)
! Close nuc
if(associated(nuc)) nullify(nuc)
end subroutine multipole_read
end module multipole

152
src/multipole_header.F90 Normal file
View file

@ -0,0 +1,152 @@
module multipole_header
implicit none
!========================================================================
! Multipole related constants
! Formalisms
integer, parameter :: FORM_MLBW = 2, &
FORM_RM = 3, &
FORM_LRM = 7
! Constants that determine which value to access
integer, parameter :: MP_EA = 1 ! Pole
! Reich-Moore indices
integer, parameter :: RM_RT = 2, & ! Residue total
RM_RA = 3, & ! Residue absorption
RM_RF = 4 ! Residue fission
! Multi-level Breit Wigner indices
integer, parameter :: MLBW_RT = 2, & ! Residue total
MLBW_RX = 3, & ! Residue compettitive
MLBW_RA = 4, & ! Residue absorption
MLBW_RF = 5 ! Residue fission
! Polynomial fit indices
integer, parameter :: FIT_T = 1, & ! Total
FIT_A = 2, & ! Absorption
FIT_F = 3 ! Fission
! Value of 'true' when checking if nuclide is fissionable
integer, parameter :: MP_FISS = 1
! These variables store the maximum value from every nuclide in order
! to preallocate some arrays to improve performance.
integer :: max_poly ! Maximum number of polynomials we expect
integer :: max_poles ! Maximum number of poles in the problem for allocation
integer :: max_L ! Maximum L value for allocation
!===============================================================================
! MULTIPOLE contains all the components needed for the windowed multipole
! temperature dependent cross section libraries for the resolved resonance
! region.
!===============================================================================
type MultipoleArray
!=========================================================================
! Isotope Properties
logical :: fissionable = .false. ! Is this isotope fissionable?
integer :: length ! Number of poles
integer, allocatable :: l_value(:) ! The l index of the pole
real(8), allocatable :: pseudo_k0RS(:) ! pseudo_k0RS for each of l
complex(8), allocatable :: data(:,:) ! Contains all of the pole-residue data
real(8) :: sqrtAWR ! Square root of the atomic weight ratio
!=========================================================================
! Windows
integer :: windows ! Number of windows
integer :: fit_order ! Order of the fit. 1 linear, 2 quadratic, etc.
real(8) :: start_E ! Start energy for the windows
real(8) :: end_E ! End energy for the windows
real(8) :: spacing ! The actual spacing in sqrt(E) space.
! spacing = sqrt(multipole_w%endE - multipole_w%startE)/multipole_w%windows
integer, allocatable :: w_start(:) ! Contains the index of the pole at the start of the window
integer, allocatable :: w_end(:) ! Contains the index of the pole at the end of the window
real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function. (reaction type, coeff index, window index)
integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't.
!=========================================================================
! Storage Helpers
integer :: num_l
integer :: max_w
integer :: formalism
contains
procedure :: clear => multipole_clear ! Deallocates Multipole
procedure :: allocate => multipole_allocate ! Allocates Multipole
end type MultipoleArray
contains
!===============================================================================
! MULTIPOLE_CLEAR resets and deallocates data in Multipole.
!===============================================================================
subroutine multipole_clear(this)
class(MultipoleArray), intent(inout) :: this
if (allocated(this % data)) then
deallocate(this % data)
deallocate(this % w_start)
deallocate(this % w_end)
deallocate(this % curvefit)
deallocate(this % l_value)
deallocate(this % pseudo_k0RS)
end if
end subroutine
!===============================================================================
! MULTIPOLE_ALLOCATE allocates necessary data for Multipole.
!===============================================================================
subroutine multipole_allocate(multipole)
class(MultipoleArray), intent(inout) :: multipole ! Multipole object to allocate.
! This function assumes length, numL, fissionable, windows, fitorder,
! and formalism are known
! Allocate the pole-residue storage.
! MLBW has one more pole than Reich-Moore, and fissionable nuclides
! have further one more.
if (multipole % formalism == FORM_MLBW) then
if (multipole % fissionable) then
allocate(multipole % data(5, multipole % length))
else
allocate(multipole % data(4, multipole % length))
end if
else if (multipole % formalism == FORM_RM) then
if (multipole % fissionable) then
allocate(multipole % data(4, multipole % length))
else
allocate(multipole % data(3, multipole % length))
end if
end if
! Allocate the l value table for each pole-residue set.
allocate(multipole % l_value(multipole % length))
! Allocate the table of pseudo_k0RS values at each l.
allocate(multipole % pseudo_k0RS(multipole % num_l))
! Allocate window start, window end
allocate(multipole % w_start(multipole % windows))
allocate(multipole % w_end(multipole % windows))
! Allocate broaden_poly
allocate(multipole % broaden_poly(multipole % windows))
! Allocate curvefit
if(multipole % fissionable) then
allocate(multipole % curvefit(FIT_F, multipole % fit_order+1, multipole % windows))
else
allocate(multipole % curvefit(FIT_A, multipole % fit_order+1, multipole % windows))
end if
end subroutine
end module multipole_header