Move most of external source sampling to source_header -- requiring removing a

lot of 'use global's, hooray!

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.

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']
 

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
 

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
 

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
 

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

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
 

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
 

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
 

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
 

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

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

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
 

src/settings.F90

@@ -4,7 +4,6 @@ module settings
 
   use constants
   use set_header,    only: SetInt
-  use source_header
 
   implicit none
 

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

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

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

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