From 76f99c2a006fc8f7a674875ed69097c041610929 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 29 Jun 2017 09:58:15 -0500 Subject: [PATCH] Add material_add_nuclide function to C API (required a little restructuring) --- CMakeLists.txt | 1 - openmc/capi.py | 24 +++++ src/api.F90 | 196 +++++++++++++++++++++++++++++++---------- src/cross_section.F90 | 1 - src/energy_grid.F90 | 67 -------------- src/global.F90 | 3 + src/initialize.F90 | 7 -- src/input_xml.F90 | 11 ++- src/nuclide_header.F90 | 39 ++++++++ src/simulation.F90 | 42 ++++++--- 10 files changed, 252 insertions(+), 139 deletions(-) delete mode 100644 src/energy_grid.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 3f8ce98e4..6dad4851c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/openmc/capi.py b/openmc/capi.py index 790fa1e20..a4ac38c32 100644 --- a/openmc/capi.py +++ b/openmc/capi.py @@ -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. diff --git a/src/api.F90 b/src/api.F90 index cba92da82..a7a7e3ab4 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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 diff --git a/src/cross_section.F90 b/src/cross_section.F90 index ce7d97133..28014f773 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 deleted file mode 100644 index 47408943e..000000000 --- a/src/energy_grid.F90 +++ /dev/null @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 72628f6d0..2690ed151 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 329455e81..98f02863e 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 09a93d118..597b4f4c0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index a55e9c841..c03ebee1a 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 diff --git a/src/simulation.F90 b/src/simulation.F90 index 1283752f6..87a17452d 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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