diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 0ab11b7c54..b738dc43d8 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -12,6 +12,10 @@ import openmc.checkvalue as cv from openmc.stats.univariate import Univariate, Uniform +__all__ = ['UnitSphere', 'PolarAzimuthal', 'Isotropic', 'Monodirectional', + 'Spatial', 'CartesianIndependent', 'Box', 'Point'] + + @add_metaclass(ABCMeta) class UnitSphere(object): """Distribution of points on the unit sphere. diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index c0476ce45c..185b222047 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -11,6 +11,10 @@ import openmc.checkvalue as cv from openmc.mixin import EqualityMixin +__all__ = ['Univariate', 'Discrete', 'Uniform', 'Maxwell', 'Watt', 'Tabular', + 'Legendre', 'Mixture'] + + _INTERPOLATION_SCHEMES = ['histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'] diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 83cb4ea257..0bdadd87ee 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -3,9 +3,8 @@ module cross_section use algorithm, only: binary_search use constants use error, only: fatal_error - use global use list_header, only: ListElemInt - use material_header, only: Material + use material_header, only: Material, materials use math, only: faddeeva, w_derivative, broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, RM_RF, & MLBW_RT, MLBW_RX, MLBW_RA, MLBW_RF, FIT_T, FIT_A,& @@ -13,7 +12,9 @@ module cross_section use nuclide_header use particle_header, only: Particle use random_lcg, only: prn, future_prn, prn_set_stream - use sab_header, only: SAlphaBeta + use sab_header, only: SAlphaBeta, sab_tables + use settings + use simulation_header implicit none diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 588cfea387..5f99d1580a 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -5,13 +5,16 @@ module eigenvalue use algorithm, only: binary_search use constants, only: ZERO use error, only: fatal_error, warning - use global use math, only: t_percentile use mesh, only: count_bank_sites - use mesh_header, only: RegularMesh + use mesh_header, only: RegularMesh, meshes use message_passing use random_lcg, only: prn, set_particle_seed, advance_prn_seed + use settings + use simulation_header use string, only: to_str + use tally_header + use timer_header implicit none diff --git a/src/geometry.F90 b/src/geometry.F90 index e27ceeedb1..3a325dfd92 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -2,16 +2,17 @@ module geometry use constants use error, only: fatal_error, warning - use geometry_header, only: Cell, Universe, Lattice, & - &RectLattice, HexLattice - use global + use geometry_header use output, only: write_message use particle_header, only: LocalCoord, Particle use particle_restart_write, only: write_particle_restart + use simulation_header + use settings use surface_header use stl_vector, only: VectorInt use string, only: to_str use tally, only: score_surface_current + use tally_header implicit none diff --git a/src/global.F90 b/src/global.F90 index 09652f294a..999a5a6763 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -13,6 +13,7 @@ module global use sab_header use settings use simulation_header + use source_header use surface_header use tally_filter_header use tally_header @@ -131,13 +132,4 @@ contains end subroutine free_memory -!=============================================================================== -! OVERALL_GENERATION determines the overall generation number -!=============================================================================== - - pure function overall_generation() result(gen) - integer :: gen - gen = gen_per_batch*(current_batch - 1) + current_gen - end function overall_generation - end module global diff --git a/src/main.F90 b/src/main.F90 index fadf9f3bbd..d8e52643c1 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,11 +1,11 @@ program main use constants - use global use message_passing use openmc_api, only: openmc_init, openmc_finalize, openmc_run, & openmc_plot_geometry, openmc_calculate_volumes use particle_restart, only: run_particle_restart + use settings, only: run_mode implicit none diff --git a/src/mesh.F90 b/src/mesh.F90 index 785d7761c3..e997890a4a 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,8 +1,8 @@ module mesh use algorithm, only: binary_search + use bank_header, only: bank use constants - use global use mesh_header use message_passing diff --git a/src/multipole.F90 b/src/multipole.F90 index c59f947272..0282adb7d9 100644 --- a/src/multipole.F90 +++ b/src/multipole.F90 @@ -3,10 +3,11 @@ module multipole use hdf5 use constants - use global + use error, only: fatal_error use hdf5_interface use multipole_header, only: MultipoleArray, FIT_T, FIT_A, FIT_F, & MP_FISS, FORM_MLBW, FORM_RM + use nuclide_header, only: nuclides implicit none diff --git a/src/output.F90 b/src/output.F90 index e3d2b8ac14..34134168be 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -3,25 +3,30 @@ module output use, intrinsic :: ISO_C_BINDING use, intrinsic :: ISO_FORTRAN_ENV + use cmfd_header use constants use eigenvalue, only: openmc_get_keff use endf, only: reaction_name use error, only: fatal_error, warning - use geometry_header, only: Cell, Universe, Lattice, RectLattice, & - HexLattice - use global + use geometry_header use math, only: t_percentile - use mesh_header, only: RegularMesh + use mesh_header, only: RegularMesh, meshes use message_passing, only: master, n_procs + use mgxs_header, only: nuclides_MG use nuclide_header use particle_header, only: LocalCoord, Particle use plot_header use sab_header, only: SAlphaBeta + use settings + use simulation_header + use surface_header, only: surfaces use string, only: to_upper, to_str - use tally_header, only: TallyObject + use tally_header + use tally_derivative_header use tally_filter use tally_filter_mesh, only: MeshFilter use tally_filter_header, only: TallyFilterMatch + use timer_header implicit none diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 35b16224d1..53332c794f 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -4,11 +4,15 @@ module particle_restart use bank_header, only: Bank use constants - use global use hdf5_interface, only: file_open, file_close, read_dataset + use mgxs_header, only: energy_bin_avg + use nuclide_header, only: micro_xs, n_nuclides use output, only: write_message, print_particle use particle_header, only: Particle use random_lcg, only: set_particle_seed + use settings + use simulation_header + use tally_header, only: n_tallies use tracking, only: transport use hdf5, only: HID_T diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 index 08e80ea1b0..5addf6056f 100644 --- a/src/particle_restart_write.F90 +++ b/src/particle_restart_write.F90 @@ -1,9 +1,10 @@ module particle_restart_write - use bank_header, only: Bank - use global + use bank_header, only: Bank, source_bank use hdf5_interface use particle_header, only: Particle + use settings + use simulation_header use string, only: to_str use hdf5 diff --git a/src/physics_common.F90 b/src/physics_common.F90 index 98d240c077..63bf62ce62 100644 --- a/src/physics_common.F90 +++ b/src/physics_common.F90 @@ -1,9 +1,9 @@ module physics_common use constants - use global, only: weight_cutoff, weight_survive use particle_header, only: Particle use random_lcg, only: prn + use settings, only: weight_cutoff, weight_survive implicit none diff --git a/src/settings.F90 b/src/settings.F90 index add7f5b903..10b7ea958f 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -4,7 +4,6 @@ module settings use constants use set_header, only: SetInt - use source_header implicit none diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index e47c2dd045..3918307a14 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -2,6 +2,7 @@ module simulation_header use bank_header use constants + use settings, only: gen_per_batch implicit none @@ -67,4 +68,15 @@ module simulation_header !$omp threadprivate(trace, thread_id, current_work) +contains + +!=============================================================================== +! OVERALL_GENERATION determines the overall generation number +!=============================================================================== + + pure function overall_generation() result(gen) + integer :: gen + gen = gen_per_batch*(current_batch - 1) + current_gen + end function overall_generation + end module simulation_header diff --git a/src/source.F90 b/src/source.F90 index 517715b563..fe7a451cbb 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -103,14 +103,7 @@ contains integer :: i ! dummy loop index integer :: n_source ! number of source distributions real(8) :: c ! cumulative frequency - real(8) :: r(3) ! sampled coordinates - logical :: found ! Does the source particle exist within geometry? - type(Particle) :: p ! Temporary particle for using find_cell - integer, save :: n_accept = 0 ! Number of samples accepted - integer, save :: n_reject = 0 ! Number of samples rejected - - ! Set weight to one by default - site % wgt = ONE + real(8) :: xi ! Set the random number generator to the source stream. call prn_set_stream(STREAM_SOURCE) @@ -118,84 +111,18 @@ contains ! Sample from among multiple source distributions n_source = size(external_source) if (n_source > 1) then - r(1) = prn()*sum(external_source(:) % strength) + xi = prn()*sum(external_source(:) % strength) c = ZERO do i = 1, n_source c = c + external_source(i) % strength - if (r(1) < c) exit + if (xi < c) exit end do else i = 1 end if - ! Repeat sampling source location until a good site has been found - found = .false. - do while (.not. found) - ! Set particle defaults - call p % initialize() - - ! Sample spatial distribution - site % xyz(:) = external_source(i) % space % sample() - - ! Fill p with needed data - p % coord(1) % xyz(:) = site % xyz - p % coord(1) % uvw(:) = [ ONE, ZERO, ZERO ] - - ! Now search to see if location exists in geometry - call find_cell(p, found) - - ! Check if spatial site is in fissionable material - if (found) then - select type (space => external_source(i) % space) - type is (SpatialBox) - if (space % only_fissionable) then - if (p % material == MATERIAL_VOID) then - found = .false. - elseif (.not. materials(p % material) % fissionable) then - found = .false. - end if - end if - end select - end if - - ! Check for rejection - if (.not. found) then - n_reject = n_reject + 1 - if (n_reject >= EXTSRC_REJECT_THRESHOLD .and. & - real(n_accept, 8)/n_reject <= EXTSRC_REJECT_FRACTION) then - call fatal_error("More than 95% of external source sites sampled & - &were rejected. Please check your external source definition.") - end if - end if - end do - - ! Increment number of accepted samples - n_accept = n_accept + 1 - - call p % clear() - - ! Sample angle - site % uvw(:) = external_source(i) % angle % sample() - - ! Check for monoenergetic source above maximum neutron energy - select type (energy => external_source(i) % energy) - type is (Discrete) - if (any(energy % x >= energy_max_neutron)) then - call fatal_error("Source energy above range of energies of at least & - &one cross section table") - end if - end select - - do - ! Sample energy spectrum - site % E = external_source(i) % energy % sample() - - ! resample if energy is greater than maximum neutron energy - if (site % E < energy_max_neutron) exit - end do - - ! Set delayed group - site % delayed_group = 0 + ! Sample source site from i-th source distribution + site = external_source(i) % sample() ! If running in MG, convert site % E to group if (.not. run_CE) then diff --git a/src/source_header.F90 b/src/source_header.F90 index 3ae3aca960..76723212ca 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -1,15 +1,25 @@ module source_header + use, intrinsic :: ISO_C_BINDING + + use bank_header, only: Bank use constants use distribution_univariate use distribution_multivariate - use error, only: fatal_error + use error + use geometry, only: find_cell + use material_header, only: materials + use nuclide_header, only: energy_max_neutron + use particle_header, only: Particle use string, only: to_lower use xml_interface implicit none private + integer :: n_accept = 0 ! Number of samples accepted + integer :: n_reject = 0 ! Number of samples rejected + !=============================================================================== ! SOURCEDISTRIBUTION describes an external source of particles for a ! fixed-source problem or for the starting source in a k eigenvalue problem @@ -21,15 +31,19 @@ module source_header class(UnitSphereDistribution), allocatable :: angle ! angle distribution class(Distribution), allocatable :: energy ! energy distribution contains - procedure :: from_xml => source_from_xml + procedure :: from_xml + procedure :: sample end type SourceDistribution + ! Number of external source distributions + integer(C_INT32_T), public, bind(C) :: n_sources = 0 + ! External source distributions type(SourceDistribution), public, allocatable :: external_source(:) contains - subroutine source_from_xml(this, node, path_source) + subroutine from_xml(this, node, path_source) class(SourceDistribution), intent(inout) :: this type(XMLNode), intent(in) :: node character(MAX_FILE_LEN), intent(out) :: path_source @@ -196,6 +210,141 @@ contains end if end if - end subroutine source_from_xml + end subroutine from_xml + + function sample(this) result(site) + class(SourceDistribution), intent(in) :: this + type(Bank) :: site + + real(8) :: r(3) ! sampled coordinates + logical :: found ! Does the source particle exist within geometry? + type(Particle) :: p ! Temporary particle for using find_cell + + ! Set weight to one by default + site % wgt = ONE + + ! Repeat sampling source location until a good site has been found + found = .false. + do while (.not. found) + ! Set particle defaults + call p % initialize() + + ! Sample spatial distribution + site % xyz(:) = this % space % sample() + + ! Fill p with needed data + p % coord(1) % xyz(:) = site % xyz + p % coord(1) % uvw(:) = [ ONE, ZERO, ZERO ] + + ! Now search to see if location exists in geometry + call find_cell(p, found) + + ! Check if spatial site is in fissionable material + if (found) then + select type (space => this % space) + type is (SpatialBox) + if (space % only_fissionable) then + if (p % material == MATERIAL_VOID) then + found = .false. + elseif (.not. materials(p % material) % fissionable) then + found = .false. + end if + end if + end select + end if + + ! Check for rejection + if (.not. found) then + n_reject = n_reject + 1 + if (n_reject >= EXTSRC_REJECT_THRESHOLD .and. & + real(n_accept, 8)/n_reject <= EXTSRC_REJECT_FRACTION) then + call fatal_error("More than 95% of external source sites sampled & + &were rejected. Please check your external source definition.") + end if + end if + end do + + ! Increment number of accepted samples + n_accept = n_accept + 1 + + call p % clear() + + ! Sample angle + site % uvw(:) = this % angle % sample() + + ! Check for monoenergetic source above maximum neutron energy + select type (energy => this % energy) + type is (Discrete) + if (any(energy % x >= energy_max_neutron)) then + call fatal_error("Source energy above range of energies of at least & + &one cross section table") + end if + end select + + do + ! Sample energy spectrum + site % E = this % energy % sample() + + ! resample if energy is greater than maximum neutron energy + if (site % E < energy_max_neutron) exit + end do + + ! Set delayed group + site % delayed_group = 0 + + end function sample + +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_extend_sources(n, index_start, index_end) result(err) bind(C) + ! Extend the external_source array by n elements + integer(C_INT32_T), value, intent(in) :: n + integer(C_INT32_T), optional, intent(out) :: index_start + integer(C_INT32_T), optional, intent(out) :: index_end + integer(C_INT) :: err + + integer :: i + type(SourceDistribution), allocatable :: temp(:) ! temporary array + + if (n_sources == 0) then + ! Allocate external_source array + allocate(external_source(n)) + else + ! Allocate external_source array with increased size + allocate(temp(n_sources + n)) + + ! Copy original source array to temporary array + temp(1:n_sources) = external_source + + ! Move allocation from temporary array + call move_alloc(FROM=temp, TO=external_source) + end if + + ! Return indices in external_source array + if (present(index_start)) index_start = n_sources + 1 + if (present(index_end)) index_end = n_sources + n + n_sources = n_sources + n + + err = 0 + end function openmc_extend_sources + + function openmc_source_set_strength(index, strength) result(err) bind(C) + integer(C_INT32_T), value, intent(in) :: index + real(C_DOUBLE), value, intent(in) :: strength + integer(C_INT) :: err + + err = E_UNASSIGNED + if (index >= 1 .and. index <= n_sources) then + if (strength > ZERO) then + external_source(index) % strength = strength + err = 0 + end if + else + err = E_OUT_OF_BOUNDS + end if + end function openmc_source_set_strength + end module source_header diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 2a948a869a..bb228e892b 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -7,14 +7,19 @@ module tally use cross_section, only: multipole_deriv_eval use error, only: fatal_error use geometry_header - use global use math, only: t_percentile, calc_pn, calc_rn - use mesh_header, only: RegularMesh + use mesh_header, only: RegularMesh, meshes use message_passing + use mgxs_header + use nuclide_header use output, only: header use particle_header, only: LocalCoord, Particle + use settings + use simulation_header use string, only: to_str + use tally_derivative_header, only: tally_derivs use tally_filter + use tally_header implicit none