diff --git a/src/api.F90 b/src/api.F90 index 4a2dc8b681..901dcaeb6b 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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 !=============================================================================== diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 5a5cf2c26b..aaea8aa5be 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7bc345da4d..545c7e3510 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 & - ¬ 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 diff --git a/src/material_header.F90 b/src/material_header.F90 index e55272ccee..e716247fe0 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index efef69535c..a2f3d5022d 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 & + ¬ 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