Add material_add_nuclide function to C API (required a little restructuring)

This commit is contained in:
Paul Romano 2017-06-29 09:58:15 -05:00
parent 75af32e019
commit 76f99c2a00
10 changed files with 252 additions and 139 deletions

View file

@ -317,7 +317,6 @@ set(LIBOPENMC_FORTRAN_SRC
src/endf.F90
src/endf_header.F90
src/energy_distribution.F90
src/energy_grid.F90
src/error.F90
src/finalize.F90
src/geometry.F90

View file

@ -28,6 +28,9 @@ class _OpenMCLibrary(object):
self._dll.openmc_init.restype = None
self._dll.openmc_load_nuclide.argtypes = [c_char_p]
self._dll.openmc_load_nuclide.restype = c_int
self._dll.openmc_material_add_nuclide.argtypes = [
c_int, c_char_p, c_double]
self._dll.openmc_material_add_nuclide.restype = c_int
self._dll.openmc_material_get_densities.argtypes = [
c_int, _double_array]
self._dll.openmc_material_get_densities.restype = c_int
@ -126,6 +129,27 @@ class _OpenMCLibrary(object):
"""
return self._dll.openmc_load_nuclide(name.encode())
def material_add_nuclide(self, mat_id, name, density):
"""Add a nuclide to a material.
Parameters
----------
mat_id : int
ID of the material
name : str
Name of nuclide, e.g. 'U235'
density : float
Density in atom/b-cm
Returns
-------
int
Return status (negative if an error occurs).
"""
return self._dll.openmc_material_add_nuclide(
mat_id, name.encode(), density)
def material_get_densities(self, mat_id):
"""Get atom densities in a material.

View file

@ -24,6 +24,7 @@ module openmc_api
public :: openmc_finalize
public :: openmc_find
public :: openmc_init
public :: openmc_material_add_nuclide
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_load_nuclide
@ -36,44 +37,6 @@ module openmc_api
contains
function openmc_material_set_density(id, density) result(err) bind(C)
integer(C_INT), value, intent(in) :: id
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: i
err = -1
if (allocated(materials)) then
if (material_dict % has_key(id)) then
i = material_dict % get_key(id)
associate (m => materials(i))
err = m % set_density(density, nuclides)
end associate
end if
end if
end function openmc_material_set_density
function openmc_material_get_densities(id, ptr) result(n) bind(C)
integer(C_INT), intent(in), value :: id
type(C_PTR), intent(out) :: ptr
integer(C_INT) :: n
ptr = C_NULL_PTR
n = 0
if (allocated(materials)) then
if (material_dict % has_key(id)) then
i = material_dict % get_key(id)
associate (m => materials(i))
if (allocated(m % atom_density)) then
ptr = C_LOC(m % atom_density(1))
n = size(m % atom_density)
end if
end associate
end if
end if
end function openmc_material_get_densities
!===============================================================================
! OPENMC_CELL_SET_TEMPERATURE sets the temperature of a cell
!===============================================================================
@ -138,20 +101,16 @@ contains
character(kind=C_CHAR) :: name(*)
integer(C_INT) :: err
integer :: i, n
integer(HID_T) :: file_id, group_id
character(10) :: name_
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_ = ''
i = 1
do while (name(i) /= C_NULL_CHAR)
name_(i:i) = name(i)
i = i + 1
end do
name_ = to_f_string(name)
err = -1
if (.not. nuclide_dict % has_key(to_lower(name_))) then
@ -184,12 +143,136 @@ contains
! 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)
err = 0
else
err = -2
end if
end if
end function openmc_load_nuclide
!===============================================================================
! OPENMC_MATERIAL_ADD_NUCLIDE
!===============================================================================
function openmc_material_add_nuclide(id, name, density) result(err) bind(C)
integer(C_INT), value, intent(in) :: id
character(kind=C_CHAR) :: name(*)
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: i, j, k, n
integer :: err2
real(8) :: awr
integer, allocatable :: new_nuclide(:)
real(8), allocatable :: new_density(:)
character(:), allocatable :: name_
name_ = to_f_string(name)
err = -1
if (allocated(materials)) then
if (material_dict % has_key(id)) then
i = material_dict % get_key(id)
associate (m => materials(i))
! 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
err2 = openmc_load_nuclide(name)
if (err2 /= -2) 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
err = 0
end if
end if
end associate
end if
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(id, ptr) result(n) bind(C)
integer(C_INT), intent(in), value :: id
type(C_PTR), intent(out) :: ptr
integer(C_INT) :: n
ptr = C_NULL_PTR
n = 0
if (allocated(materials)) then
if (material_dict % has_key(id)) then
i = material_dict % get_key(id)
associate (m => materials(i))
if (allocated(m % atom_density)) then
ptr = C_LOC(m % atom_density(1))
n = size(m % atom_density)
end if
end associate
end if
end if
end function openmc_material_get_densities
!===============================================================================
! OPENMC_MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm
!===============================================================================
function openmc_material_set_density(id, density) result(err) bind(C)
integer(C_INT), value, intent(in) :: id
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
integer :: i
err = -1
if (allocated(materials)) then
if (material_dict % has_key(id)) then
i = material_dict % get_key(id)
associate (m => materials(i))
err = m % set_density(density, nuclides)
end associate
end if
end if
end function openmc_material_set_density
!===============================================================================
! OPENMC_RESET resets all tallies
!===============================================================================
@ -306,4 +389,23 @@ contains
end if
end subroutine openmc_tally_results
function to_f_string(c_string) result(f_string)
character(kind=C_CHAR), intent(in) :: c_string(*)
character(:), allocatable :: f_string
integer :: i, n
! Determine length of original string
n = 0
do while (c_string(n + 1) /= C_NULL_CHAR)
n = n + 1
end do
! Copy C string character by character
allocate(character(len=n) :: f_string)
do i = 1, n
f_string(i:i) = c_string(i)
end do
end function to_f_string
end module openmc_api

View file

@ -2,7 +2,6 @@ module cross_section
use algorithm, only: binary_search
use constants
use energy_grid, only: log_spacing
use error, only: fatal_error
use global
use list_header, only: ListElemInt

View file

@ -1,67 +0,0 @@
module energy_grid
use global
use list_header, only: ListReal
use output, only: write_message
implicit none
integer :: grid_method ! how to treat the energy grid
integer :: n_log_bins ! number of bins for logarithmic grid
real(8) :: log_spacing ! spacing on logarithmic grid
contains
!===============================================================================
! LOGARITHMIC_GRID determines a logarithmic mapping for energies to bounding
! indices on a nuclide energy grid
!===============================================================================
subroutine logarithmic_grid()
integer :: i, j, k ! Loop indices
integer :: t ! temperature index
integer :: M ! Number of equally log-spaced bins
real(8) :: E_max ! Maximum energy in MeV
real(8) :: E_min ! Minimum energy in MeV
real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid
! Set minimum/maximum energies
E_max = energy_max_neutron
E_min = energy_min_neutron
! Determine equal-logarithmic energy spacing
M = n_log_bins
log_spacing = log(E_max/E_min)/M
! Create equally log-spaced energy grid
allocate(umesh(0:M))
umesh(:) = [(i*log_spacing, i=0, M)]
do i = 1, n_nuclides_total
associate (nuc => nuclides(i))
do t = 1, size(nuc % grid)
! Allocate logarithmic mapping for nuclide
allocate(nuc % grid(t) % grid_index(0:M))
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(nuc % grid(t) % energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(nuc % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == size(nuc % grid(t) % energy)) exit
j = j + 1
end do
nuc % grid(t) % grid_index(k) = j
end do
end do
end associate
end do
deallocate(umesh)
end subroutine logarithmic_grid
end module energy_grid

View file

@ -112,6 +112,9 @@ module global
real(8) :: temperature_default = 293.6_8
real(8) :: temperature_range(2) = [ZERO, ZERO]
integer :: n_log_bins ! number of bins for logarithmic grid
real(8) :: log_spacing ! spacing on logarithmic grid
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES

View file

@ -11,7 +11,6 @@ module initialize
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use energy_grid, only: logarithmic_grid, grid_method
use error, only: fatal_error, warning
use geometry, only: neighbor_lists, count_instance, calc_offsets, &
maximum_levels
@ -109,12 +108,6 @@ contains
end if
if (run_mode /= MODE_PLOTTING) then
! Construct information needed for nuclear data
if (run_CE) then
! Construct log energy grid for cross-sections
call logarithmic_grid()
end if
! Allocate and setup tally stride, filter_matches, and tally maps
call configure_tallies()

View file

@ -7,7 +7,6 @@ module input_xml
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures, root_universe
@ -5240,9 +5239,6 @@ contains
allocate(nuclides(n_nuclides_total))
allocate(sab_tables(n_sab_tables))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! Read cross sections
do i = 1, size(materials)
@ -5294,6 +5290,13 @@ contains
end do
end do
! Set up logarithmic grid for nuclides
do i = 1, size(nuclides)
call nuclides(i) % init_grid(energy_min_neutron, &
energy_max_neutron, n_log_bins)
end do
log_spacing = log(energy_max_neutron/energy_min_neutron) / n_log_bins
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
if (.not. allocated(materials(i) % sab_names)) cycle

View file

@ -97,6 +97,7 @@ module nuclide_header
contains
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
end type Nuclide
@ -683,4 +684,42 @@ module nuclide_header
end function nuclide_nu
subroutine nuclide_init_grid(this, E_min, E_max, M)
class(Nuclide), intent(inout) :: this
real(8), intent(in) :: E_min ! Minimum energy in MeV
real(8), intent(in) :: E_max ! Maximum energy in MeV
integer, intent(in) :: M ! Number of equally log-spaced bins
integer :: i, j, k ! Loop indices
integer :: t ! temperature index
real(8) :: spacing
real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid
! Determine equal-logarithmic energy spacing
spacing = log(E_max/E_min)/M
! Create equally log-spaced energy grid
allocate(umesh(0:M))
umesh(:) = [(i*spacing, i=0, M)]
do t = 1, size(this % grid)
! Allocate logarithmic mapping for nuclide
allocate(this % grid(t) % grid_index(0:M))
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M
do while (log(this % grid(t) % energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(this % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == size(this % grid(t) % energy)) exit
j = j + 1
end do
this % grid(t) % grid_index(k) = j
end do
end do
end subroutine nuclide_init_grid
end module nuclide_header

View file

@ -42,17 +42,7 @@ contains
type(Particle) :: p
integer(8) :: i_work
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", 3)
if (verbosity >= 7) call print_columns()
end if
end if
call initialize_simulation()
! Turn on inactive timer
call time_inactive % start()
@ -380,12 +370,40 @@ contains
end subroutine replay_batch_history
!===============================================================================
! INITIALIZE_SIMULATION
!===============================================================================
subroutine initialize_simulation()
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", 3)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", 3)
if (verbosity >= 7) call print_columns()
end if
end if
end subroutine initialize_simulation
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
subroutine finalize_simulation
subroutine finalize_simulation()
!$omp parallel
deallocate(micro_xs)
!$omp end parallel
! Increment total number of generations
total_gen = total_gen + n_batches*gen_per_batch