Move many procedures out of api.F90 into modules where they belong. Also move

assign_0K_elastic_scattering and check_data_version to nuclide_header
This commit is contained in:
Paul Romano 2017-09-19 21:49:25 -05:00
parent 2e78305b02
commit 9674f1ce2e
5 changed files with 532 additions and 541 deletions

View file

@ -15,7 +15,6 @@ module openmc_api
use message_passing
use nuclide_header
use initialize, only: openmc_init
use input_xml, only: assign_0K_elastic_scattering, check_data_version
use particle_header, only: Particle
use plot, only: openmc_plot_geometry
use random_lcg, only: seed, initialize_prng
@ -78,116 +77,6 @@ module openmc_api
contains
!===============================================================================
! OPENMC_CELL_GET_ID returns the ID of a cell
!===============================================================================
function openmc_cell_get_id(index, id) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(cells)) then
id = cells(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_get_id
!===============================================================================
! OPENMC_CELL_SET_TEMPERATURE sets the temperature of a cell
!===============================================================================
function openmc_cell_set_temperature(index, T, instance) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index ! cell index in cells
real(C_DOUBLE), value, intent(in) :: T ! temperature
integer(C_INT32_T), optional, intent(in) :: instance ! cell instance
integer(C_INT) :: err ! error code
integer :: j ! looping variable
integer :: n ! number of cell instances
integer :: material_index ! material index in materials array
integer :: num_nuclides ! num nuclides in material
integer :: nuclide_index ! index of nuclide in nuclides array
real(8) :: min_temp ! min common-denominator avail temp
real(8) :: max_temp ! max common-denominator avail temp
real(8) :: temp ! actual temp we'll assign
logical :: outside_low ! lower than available data
logical :: outside_high ! higher than available data
outside_low = .false.
outside_high = .false.
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(cells)) then
! error if the cell is filled with another universe
if (cells(index) % fill /= NONE) then
err = E_CELL_NO_MATERIAL
else
! find which material is associated with this cell (material_index
! is the index into the materials array)
if (present(instance)) then
material_index = cells(index) % material(instance)
else
material_index = cells(index) % material(1)
end if
! number of nuclides associated with this material
num_nuclides = size(materials(material_index) % nuclide)
min_temp = ZERO
max_temp = INFINITY
do j = 1, num_nuclides
nuclide_index = materials(material_index) % nuclide(j)
min_temp = max(min_temp, minval(nuclides(nuclide_index) % kTs))
max_temp = min(max_temp, maxval(nuclides(nuclide_index) % kTs))
end do
! adjust the temperature to be within bounds if necessary
if (K_BOLTZMANN * T < min_temp) then
outside_low = .true.
temp = min_temp / K_BOLTZMANN
else if (K_BOLTZMANN * T > max_temp) then
outside_high = .true.
temp = max_temp / K_BOLTZMANN
else
temp = T
end if
associate (c => cells(index))
if (allocated(c % sqrtkT)) then
n = size(c % sqrtkT)
if (present(instance) .and. n > 1) then
if (instance >= 0 .and. instance < n) then
c % sqrtkT(instance + 1) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
else
c % sqrtkT(:) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
end if
end associate
! Assign error codes for outside of temperature bounds provided the
! temperature was changed correctly. This needs to be done after
! changing the temperature based on the logical structure above.
if (err == 0) then
if (outside_low) err = W_BELOW_MIN_BOUND
if (outside_high) err = W_ABOVE_MAX_BOUND
end if
end if
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_set_temperature
!===============================================================================
! OPENMC_FINALIZE frees up memory by deallocating arrays and resetting global
! variables
@ -309,76 +198,6 @@ contains
end function openmc_find
!===============================================================================
! OPENMC_GET_CELL returns the index in the cells array of a cell with a given ID
!===============================================================================
function openmc_get_cell(id, index) result(err) bind(C)
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has_key(id)) then
index = cell_dict % get_key(id)
err = 0
else
err = E_CELL_INVALID_ID
end if
else
err = E_CELL_NOT_ALLOCATED
end if
end function openmc_get_cell
!===============================================================================
! OPENMC_GET_MATERIAL returns the index in the materials array of a material
! with a given ID
!===============================================================================
function openmc_get_material(id, index) result(err) bind(C)
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(materials)) then
if (material_dict % has_key(id)) then
index = material_dict % get_key(id)
err = 0
else
err = E_MATERIAL_INVALID_ID
end if
else
err = E_MATERIAL_NOT_ALLOCATED
end if
end function openmc_get_material
!===============================================================================
! OPENMC_GET_NUCLIDE returns the index in the nuclides array of a nuclide
! with a given name
!===============================================================================
function openmc_get_nuclide(name, index) result(err) bind(C)
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT), intent(out) :: index
integer(C_INT) :: err
character(:), allocatable :: name_
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
if (allocated(nuclides)) then
if (nuclide_dict % has_key(to_lower(name_))) then
index = nuclide_dict % get_key(to_lower(name_))
err = 0
else
err = E_NUCLIDE_NOT_LOADED
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_get_nuclide
!===============================================================================
! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom
! generator state
@ -397,278 +216,6 @@ contains
call initialize_prng()
end subroutine openmc_hard_reset
!===============================================================================
! OPENMC_LOAD_NUCLIDE loads a nuclide from the cross section library
!===============================================================================
function openmc_load_nuclide(name) result(err) bind(C)
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT) :: err
integer :: i_library
integer :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
character(:), allocatable :: name_
real(8) :: minmax(2) = [ZERO, INFINITY]
type(VectorReal) :: temperature
type(Nuclide), allocatable :: new_nuclides(:)
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
err = 0
if (.not. nuclide_dict % has_key(to_lower(name_))) then
if (library_dict % has_key(to_lower(name_))) then
! allocate extra space in nuclides array
n = n_nuclides
allocate(new_nuclides(n + 1))
new_nuclides(1:n) = nuclides(:)
call move_alloc(FROM=new_nuclides, TO=nuclides)
n = n + 1
i_library = library_dict % get_key(to_lower(name_))
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
call check_data_version(file_id)
! Read nuclide data from HDF5
group_id = open_group(file_id, name_)
call nuclides(n) % from_hdf5(group_id, temperature, &
temperature_method, temperature_tolerance, minmax, &
master)
call close_group(group_id)
call file_close(file_id)
! Add entry to nuclide dictionary
call nuclide_dict % add_key(to_lower(name_), n)
n_nuclides = n
! Assign resonant scattering data
if (res_scat_on) call assign_0K_elastic_scattering(nuclides(n))
! Initialize nuclide grid
call nuclides(n) % init_grid(energy_min_neutron, &
energy_max_neutron, n_log_bins)
else
err = E_NUCLIDE_NOT_IN_LIBRARY
end if
end if
end function openmc_load_nuclide
!===============================================================================
! OPENMC_MATERIAL_ADD_NUCLIDE
!===============================================================================
function openmc_material_add_nuclide(index, name, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR) :: name(*)
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: j, k, n
real(8) :: awr
integer, allocatable :: new_nuclide(:)
real(8), allocatable :: new_density(:)
character(:), allocatable :: name_
name_ = to_f_string(name)
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! Check if nuclide is already in material
do j = 1, size(m % nuclide)
k = m % nuclide(j)
if (nuclides(k) % name == name_) then
awr = nuclides(k) % awr
m % density = m % density + density - m % atom_density(j)
m % density_gpcc = m % density_gpcc + (density - &
m % atom_density(j)) * awr * MASS_NEUTRON / N_AVOGADRO
m % atom_density(j) = density
err = 0
end if
end do
! If nuclide wasn't found, extend nuclide/density arrays
if (err /= 0) then
! If nuclide hasn't been loaded, load it now
err = openmc_load_nuclide(name)
if (err == 0) then
! Extend arrays
n = size(m % nuclide)
allocate(new_nuclide(n + 1))
new_nuclide(1:n) = m % nuclide
call move_alloc(FROM=new_nuclide, TO=m % nuclide)
allocate(new_density(n + 1))
new_density(1:n) = m % atom_density
call move_alloc(FROM=new_density, TO=m % atom_density)
! Append new nuclide/density
k = nuclide_dict % get_key(to_lower(name_))
m % nuclide(n + 1) = k
m % atom_density(n + 1) = density
m % density = m % density + density
m % density_gpcc = m % density_gpcc + &
density * nuclides(k) % awr * MASS_NEUTRON / N_AVOGADRO
m % n_nuclides = n + 1
end if
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_add_nuclide
!===============================================================================
! OPENMC_MATERIAL_GET_DENSITIES returns an array of nuclide densities in a
! material
!===============================================================================
function openmc_material_get_densities(index, nuclides, densities, n) &
result(err) bind(C)
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: nuclides
type(C_PTR), intent(out) :: densities
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
if (allocated(m % atom_density)) then
nuclides = C_LOC(m % nuclide(1))
densities = C_LOC(m % atom_density(1))
n = size(m % atom_density)
err = 0
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_densities
!===============================================================================
! OPENMC_MATERIAL_GET_ID returns the ID of a material
!===============================================================================
function openmc_material_get_id(index, id) result(err) bind(C)
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
id = materials(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_id
!===============================================================================
! OPENMC_MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm
!===============================================================================
function openmc_material_set_density(index, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
err = m % set_density(density)
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_density
!===============================================================================
! OPENMC_MATERIAL_SET_DENSITIES sets the densities for a list of nuclides in a
! material. If the nuclides don't already exist in the material, they will be
! added
!===============================================================================
function openmc_material_set_densities(index, n, name, density) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
type(C_PTR), intent(in) :: name(n)
real(C_DOUBLE), intent(in) :: density(n)
integer(C_INT) :: err
integer :: i
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: name_
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! If nuclide/density arrays are not correct size, reallocate
if (n /= size(m % nuclide)) then
deallocate(m % nuclide, m % atom_density, m % p0)
allocate(m % nuclide(n), m % atom_density(n), m % p0(n))
end if
do i = 1, n
! Convert C string to Fortran string
call c_f_pointer(name(i), string, [10])
name_ = to_lower(to_f_string(string))
if (.not. nuclide_dict % has_key(name_)) then
err = openmc_load_nuclide(string)
if (err < 0) return
end if
m % nuclide(i) = nuclide_dict % get_key(name_)
m % atom_density(i) = density(i)
end do
m % n_nuclides = n
! Set isotropic flags to flags
m % p0(:) = .false.
! Set total density to the sum of the vector
err = m % set_density(sum(density))
! Assign S(a,b) tables
call m % assign_sab_tables()
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_densities
!===============================================================================
! OPENMC_NUCLIDE_NAME returns the name of a nuclide with a given index
!===============================================================================
function openmc_nuclide_name(index, name) result(err) bind(C)
integer(C_INT), value, intent(in) :: index
type(c_ptr), intent(out) :: name
integer(C_INT) :: err
character(C_CHAR), pointer :: name_
err = E_UNASSIGNED
if (allocated(nuclides)) then
if (index >= 1 .and. index <= size(nuclides)) then
name_ => nuclides(index) % name(1:1)
name = C_LOC(name_)
err = 0
else
err = E_OUT_OF_BOUNDS
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_nuclide_name
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================

View file

@ -422,4 +422,132 @@ contains
end subroutine free_memory_geometry
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_get_cell(id, index) result(err) bind(C)
! Return the index in the cells array of a cell with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has_key(id)) then
index = cell_dict % get_key(id)
err = 0
else
err = E_CELL_INVALID_ID
end if
else
err = E_CELL_NOT_ALLOCATED
end if
end function openmc_get_cell
function openmc_cell_get_id(index, id) result(err) bind(C)
! Return the ID of a cell
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(cells)) then
id = cells(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_get_id
function openmc_cell_set_temperature(index, T, instance) result(err) bind(C)
! Set the temperature of a cell
integer(C_INT32_T), value, intent(in) :: index ! cell index in cells
real(C_DOUBLE), value, intent(in) :: T ! temperature
integer(C_INT32_T), optional, intent(in) :: instance ! cell instance
integer(C_INT) :: err ! error code
integer :: j ! looping variable
integer :: n ! number of cell instances
integer :: material_index ! material index in materials array
integer :: num_nuclides ! num nuclides in material
integer :: nuclide_index ! index of nuclide in nuclides array
real(8) :: min_temp ! min common-denominator avail temp
real(8) :: max_temp ! max common-denominator avail temp
real(8) :: temp ! actual temp we'll assign
logical :: outside_low ! lower than available data
logical :: outside_high ! higher than available data
outside_low = .false.
outside_high = .false.
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(cells)) then
! error if the cell is filled with another universe
if (cells(index) % fill /= NONE) then
err = E_CELL_NO_MATERIAL
else
! find which material is associated with this cell (material_index
! is the index into the materials array)
if (present(instance)) then
material_index = cells(index) % material(instance)
else
material_index = cells(index) % material(1)
end if
! number of nuclides associated with this material
num_nuclides = size(materials(material_index) % nuclide)
min_temp = ZERO
max_temp = INFINITY
do j = 1, num_nuclides
nuclide_index = materials(material_index) % nuclide(j)
min_temp = max(min_temp, minval(nuclides(nuclide_index) % kTs))
max_temp = min(max_temp, maxval(nuclides(nuclide_index) % kTs))
end do
! adjust the temperature to be within bounds if necessary
if (K_BOLTZMANN * T < min_temp) then
outside_low = .true.
temp = min_temp / K_BOLTZMANN
else if (K_BOLTZMANN * T > max_temp) then
outside_high = .true.
temp = max_temp / K_BOLTZMANN
else
temp = T
end if
associate (c => cells(index))
if (allocated(c % sqrtkT)) then
n = size(c % sqrtkT)
if (present(instance) .and. n > 1) then
if (instance >= 0 .and. instance < n) then
c % sqrtkT(instance + 1) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
else
c % sqrtkT(:) = sqrt(K_BOLTZMANN * temp)
err = 0
end if
end if
end associate
! Assign error codes for outside of temperature bounds provided the
! temperature was changed correctly. This needs to be done after
! changing the temperature based on the logical structure above.
if (err == 0) then
if (outside_low) err = W_BELOW_MIN_BOUND
if (outside_high) err = W_ABOVE_MAX_BOUND
end if
end if
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_cell_set_temperature
end module geometry_header

View file

@ -22,6 +22,7 @@ module input_xml
use mgxs_data, only: create_macro_xs, read_mgxs
use mgxs_header
use multipole, only: multipole_read
use nuclide_header
use output, only: write_message, title, header, print_plot
use plot_header
use random_lcg, only: prn, seed, initialize_prng
@ -4455,7 +4456,7 @@ contains
call file_close(file_id)
! Assign resonant scattering data
if (res_scat_on) call assign_0K_elastic_scattering(nuclides(i_nuclide))
if (res_scat_on) call nuclides(i_nuclide) % assign_0K_elastic_scattering()
! Determine if minimum/maximum energy for this nuclide is greater/less
! than the previous
@ -4590,60 +4591,6 @@ contains
end do
end subroutine assign_temperatures
!===============================================================================
! ASSIGN_0K_ELASTIC_SCATTERING
!===============================================================================
subroutine assign_0K_elastic_scattering(nuc)
type(Nuclide), intent(inout) :: nuc
integer :: i
real(8) :: xs_cdf_sum
nuc % resonant = .false.
if (allocated(res_scat_nuclides)) then
! If resonant nuclides were specified, check the list explicitly
do i = 1, size(res_scat_nuclides)
if (nuc % name == res_scat_nuclides(i)) then
nuc % resonant = .true.
! Make sure nuclide has 0K data
if (.not. allocated(nuc % energy_0K)) then
call fatal_error("Cannot treat " // trim(nuc % name) // " as a &
&resonant scatterer because 0 K elastic scattering data is &
&not present.")
end if
exit
end if
end do
else
! Otherwise, assume that any that have 0 K elastic scattering data are
! resonant
nuc % resonant = allocated(nuc % energy_0K)
end if
if (nuc % resonant) then
! Build CDF for 0K elastic scattering
xs_cdf_sum = ZERO
allocate(nuc % xs_cdf(0:size(nuc % energy_0K)))
nuc % xs_cdf(0) = ZERO
associate (E => nuc % energy_0K, xs => nuc % elastic_0K)
do i = 1, size(E) - 1
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to zero.
if (xs(i) < ZERO) xs(i) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))&
/ TWO * (E(i+1) - E(i))
nuc % xs_cdf(i) = xs_cdf_sum
end do
end associate
end if
end subroutine assign_0K_elastic_scattering
!===============================================================================
! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the
! directory and loads it using multipole_read
@ -4693,31 +4640,6 @@ contains
end subroutine read_multipole_data
!===============================================================================
! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5
! files
!===============================================================================
subroutine check_data_version(file_id)
integer(HID_T), intent(in) :: file_id
integer, allocatable :: version(:)
if (attribute_exists(file_id, 'version')) then
call read_attribute(version, file_id, 'version')
if (version(1) /= HDF5_VERSION(1)) then
call fatal_error("HDF5 data format uses version " // trim(to_str(&
version(1))) // "." // trim(to_str(version(2))) // " whereas &
&your installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION(1))) // ".x data.")
end if
else
call fatal_error("HDF5 data does not indicate a version. Your &
&installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION(1))) // ".x data.")
end if
end subroutine check_data_version
!===============================================================================
! ADJUST_INDICES changes the values for 'surfaces' for each cell and the
! material index assigned to each to the indices in the surfaces and material

View file

@ -12,11 +12,20 @@ module material_header
implicit none
private
public :: free_memory_material
public :: openmc_get_material
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_material_set_densities
!===============================================================================
! MATERIAL describes a material by its constituent nuclides
!===============================================================================
type Material
type, public :: Material
integer :: id ! unique identifier
character(len=104) :: name = "" ! User-defined name
integer :: n_nuclides ! number of nuclides
@ -54,12 +63,12 @@ module material_header
procedure :: assign_sab_tables => material_assign_sab_tables
end type Material
integer(C_INT32_T), bind(C) :: n_materials ! # of materials
integer(C_INT32_T), public, bind(C) :: n_materials ! # of materials
type(Material), allocatable, target :: materials(:)
type(Material), public, allocatable, target :: materials(:)
! Dictionary that maps user IDs to indices in 'materials'
type(DictIntInt) :: material_dict
type(DictIntInt), public :: material_dict
contains
@ -222,4 +231,200 @@ contains
call material_dict % clear()
end subroutine free_memory_material
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_get_material(id, index) result(err) bind(C)
! Returns the index in the materials array of a material with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(materials)) then
if (material_dict % has_key(id)) then
index = material_dict % get_key(id)
err = 0
else
err = E_MATERIAL_INVALID_ID
end if
else
err = E_MATERIAL_NOT_ALLOCATED
end if
end function openmc_get_material
function openmc_material_add_nuclide(index, name, density) result(err) bind(C)
! Add a nuclide at a specified density in atom/b-cm to a material
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR) :: name(*)
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: j, k, n
real(8) :: awr
integer, allocatable :: new_nuclide(:)
real(8), allocatable :: new_density(:)
character(:), allocatable :: name_
name_ = to_f_string(name)
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! Check if nuclide is already in material
do j = 1, size(m % nuclide)
k = m % nuclide(j)
if (nuclides(k) % name == name_) then
awr = nuclides(k) % awr
m % density = m % density + density - m % atom_density(j)
m % density_gpcc = m % density_gpcc + (density - &
m % atom_density(j)) * awr * MASS_NEUTRON / N_AVOGADRO
m % atom_density(j) = density
err = 0
end if
end do
! If nuclide wasn't found, extend nuclide/density arrays
if (err /= 0) then
! If nuclide hasn't been loaded, load it now
err = openmc_load_nuclide(name)
if (err == 0) then
! Extend arrays
n = size(m % nuclide)
allocate(new_nuclide(n + 1))
new_nuclide(1:n) = m % nuclide
call move_alloc(FROM=new_nuclide, TO=m % nuclide)
allocate(new_density(n + 1))
new_density(1:n) = m % atom_density
call move_alloc(FROM=new_density, TO=m % atom_density)
! Append new nuclide/density
k = nuclide_dict % get_key(to_lower(name_))
m % nuclide(n + 1) = k
m % atom_density(n + 1) = density
m % density = m % density + density
m % density_gpcc = m % density_gpcc + &
density * nuclides(k) % awr * MASS_NEUTRON / N_AVOGADRO
m % n_nuclides = n + 1
end if
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_add_nuclide
function openmc_material_get_densities(index, nuclides, densities, n) &
result(err) bind(C)
! returns an array of nuclide densities in a material
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: nuclides
type(C_PTR), intent(out) :: densities
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
if (allocated(m % atom_density)) then
nuclides = C_LOC(m % nuclide(1))
densities = C_LOC(m % atom_density(1))
n = size(m % atom_density)
err = 0
end if
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_densities
function openmc_material_get_id(index, id) result(err) bind(C)
! returns the ID of a material
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
id = materials(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_get_id
function openmc_material_set_density(index, density) result(err) bind(C)
! Set the total density of a material in atom/b-cm
integer(C_INT32_T), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
err = m % set_density(density)
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_density
function openmc_material_set_densities(index, n, name, density) result(err) bind(C)
! Sets the densities for a list of nuclides in a material. If the nuclides
! don't already exist in the material, they will be added
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
type(C_PTR), intent(in) :: name(n)
real(C_DOUBLE), intent(in) :: density(n)
integer(C_INT) :: err
integer :: i
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: name_
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
! If nuclide/density arrays are not correct size, reallocate
if (n /= size(m % nuclide)) then
deallocate(m % nuclide, m % atom_density, m % p0)
allocate(m % nuclide(n), m % atom_density(n), m % p0(n))
end if
do i = 1, n
! Convert C string to Fortran string
call c_f_pointer(name(i), string, [10])
name_ = to_lower(to_f_string(string))
if (.not. nuclide_dict % has_key(name_)) then
err = openmc_load_nuclide(string)
if (err < 0) return
end if
m % nuclide(i) = nuclide_dict % get_key(name_)
m % atom_density(i) = density(i)
end do
m % n_nuclides = n
! Set isotropic flags to flags
m % p0(:) = .false.
! Set total density to the sum of the vector
err = m % set_density(sum(density))
! Assign S(a,b) tables
call m % assign_sab_tables()
end associate
else
err = E_OUT_OF_BOUNDS
end if
end function openmc_material_set_densities
end module material_header

View file

@ -10,16 +10,16 @@ module nuclide_header
use dict_header, only: DictIntInt, DictCharInt
use endf, only: reaction_name, is_fission, is_disappearance
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error, only: fatal_error, warning
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape, get_datasets, &
object_exists, get_name, get_groups
use error
use hdf5_interface
use list_header, only: ListInt
use math, only: evaluate_legendre
use message_passing
use multipole_header, only: MultipoleArray
use product_header, only: AngleEnergyContainer
use reaction_header, only: Reaction
use secondary_uncorrelated, only: UncorrelatedAngleEnergy
use settings
use stl_vector, only: VectorInt, VectorReal
use string
use urr_header, only: UrrData
@ -95,6 +95,7 @@ module nuclide_header
class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas
contains
procedure :: assign_0K_elastic_scattering
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
@ -177,6 +178,60 @@ module nuclide_header
contains
!===============================================================================
! ASSIGN_0K_ELASTIC_SCATTERING
!===============================================================================
subroutine assign_0K_elastic_scattering(this)
class(Nuclide), intent(inout) :: this
integer :: i
real(8) :: xs_cdf_sum
this % resonant = .false.
if (allocated(res_scat_nuclides)) then
! If resonant nuclides were specified, check the list explicitly
do i = 1, size(res_scat_nuclides)
if (this % name == res_scat_nuclides(i)) then
this % resonant = .true.
! Make sure nuclide has 0K data
if (.not. allocated(this % energy_0K)) then
call fatal_error("Cannot treat " // trim(this % name) // " as a &
&resonant scatterer because 0 K elastic scattering data is &
&not present.")
end if
exit
end if
end do
else
! Otherwise, assume that any that have 0 K elastic scattering data are
! resonant
this % resonant = allocated(this % energy_0K)
end if
if (this % resonant) then
! Build CDF for 0K elastic scattering
xs_cdf_sum = ZERO
allocate(this % xs_cdf(0:size(this % energy_0K)))
this % xs_cdf(0) = ZERO
associate (E => this % energy_0K, xs => this % elastic_0K)
do i = 1, size(E) - 1
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to zero.
if (xs(i) < ZERO) xs(i) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))&
/ TWO * (E(i+1) - E(i))
this % xs_cdf(i) = xs_cdf_sum
end do
end associate
end if
end subroutine assign_0K_elastic_scattering
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide
!===============================================================================
@ -742,6 +797,31 @@ contains
end subroutine nuclide_init_grid
!===============================================================================
! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5
! files
!===============================================================================
subroutine check_data_version(file_id)
integer(HID_T), intent(in) :: file_id
integer, allocatable :: version(:)
if (attribute_exists(file_id, 'version')) then
call read_attribute(version, file_id, 'version')
if (version(1) /= HDF5_VERSION(1)) then
call fatal_error("HDF5 data format uses version " // trim(to_str(&
version(1))) // "." // trim(to_str(version(2))) // " whereas &
&your installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION(1))) // ".x data.")
end if
else
call fatal_error("HDF5 data does not indicate a version. Your &
&installation of OpenMC expects version " // trim(to_str(&
HDF5_VERSION(1))) // ".x data.")
end if
end subroutine check_data_version
!===============================================================================
! FREE_MEMORY_NUCLIDE deallocates global arrays defined in this module
!===============================================================================
@ -766,4 +846,113 @@ contains
end subroutine free_memory_nuclide
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_get_nuclide(name, index) result(err) bind(C)
! Return the index in the nuclides array of a nuclide with a given name
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT), intent(out) :: index
integer(C_INT) :: err
character(:), allocatable :: name_
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
if (allocated(nuclides)) then
if (nuclide_dict % has_key(to_lower(name_))) then
index = nuclide_dict % get_key(to_lower(name_))
err = 0
else
err = E_NUCLIDE_NOT_LOADED
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_get_nuclide
function openmc_load_nuclide(name) result(err) bind(C)
! Load a nuclide from the cross section library
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT) :: err
integer :: i_library
integer :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
character(:), allocatable :: name_
real(8) :: minmax(2) = [ZERO, INFINITY]
type(VectorReal) :: temperature
type(Nuclide), allocatable :: new_nuclides(:)
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
err = 0
if (.not. nuclide_dict % has_key(to_lower(name_))) then
if (library_dict % has_key(to_lower(name_))) then
! allocate extra space in nuclides array
n = n_nuclides
allocate(new_nuclides(n + 1))
new_nuclides(1:n) = nuclides(:)
call move_alloc(FROM=new_nuclides, TO=nuclides)
n = n + 1
i_library = library_dict % get_key(to_lower(name_))
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
call check_data_version(file_id)
! Read nuclide data from HDF5
group_id = open_group(file_id, name_)
call nuclides(n) % from_hdf5(group_id, temperature, &
temperature_method, temperature_tolerance, minmax, &
master)
call close_group(group_id)
call file_close(file_id)
! Add entry to nuclide dictionary
call nuclide_dict % add_key(to_lower(name_), n)
n_nuclides = n
! Assign resonant scattering data
if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering()
! Initialize nuclide grid
call nuclides(n) % init_grid(energy_min_neutron, &
energy_max_neutron, n_log_bins)
else
err = E_NUCLIDE_NOT_IN_LIBRARY
end if
end if
end function openmc_load_nuclide
function openmc_nuclide_name(index, name) result(err) bind(C)
! Return the name of a nuclide with a given index
integer(C_INT), value, intent(in) :: index
type(c_ptr), intent(out) :: name
integer(C_INT) :: err
character(C_CHAR), pointer :: name_
err = E_UNASSIGNED
if (allocated(nuclides)) then
if (index >= 1 .and. index <= size(nuclides)) then
name_ => nuclides(index) % name(1:1)
name = C_LOC(name_)
err = 0
else
err = E_OUT_OF_BOUNDS
end if
else
err = E_NUCLIDE_NOT_ALLOCATED
end if
end function openmc_nuclide_name
end module nuclide_header