diff --git a/src/constants.F90 b/src/constants.F90 index 8eabdaf60..d51b039bb 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -357,28 +357,6 @@ module constants integer, parameter :: STREAM_TALLIES = 2 integer, parameter :: STREAM_SOURCE = 3 - ! ============================================================================ - ! EXTERNAL SOURCE PARAMETERS - - ! Source spatial distribution types - integer, parameter :: & - SRC_SPACE_BOX = 1, & ! Source in a rectangular prism - SRC_SPACE_POINT = 2, & ! Source at a single point - SRC_SPACE_FISSION = 3 ! Source in prism filtered by fissionable mats - - ! Source angular distribution types - integer, parameter :: & - SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular - SRC_ANGLE_MONO = 2, & ! Monodirectional source - SRC_ANGLE_TABULAR = 3 ! Tabular distribution - - ! Source energy distribution types - integer, parameter :: & - SRC_ENERGY_MONO = 1, & ! Monoenergetic source - SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum - SRC_ENERGY_WATT = 3, & ! Watt fission spectrum - SRC_ENERGY_TABULAR = 4 ! Tabular distribution - ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index d5a2d1f05..970311d7b 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -49,6 +49,46 @@ module distribution_multivariate procedure :: sample => monodirectional_sample end type Monodirectional +!=============================================================================== +! SPATIALDISTRIBUTION type defines a probability density function for arbitrary +! points in Euclidean space. +!=============================================================================== + + type, abstract :: SpatialDistribution + contains + procedure(iSampleSpatial), deferred :: sample + end type SpatialDistribution + + abstract interface + function iSampleSpatial(this) result(xyz) + import SpatialDistribution + class(SpatialDistribution), intent(in) :: this + real(8) :: xyz(3) + end function iSampleSpatial + end interface + + type, extends(SpatialDistribution) :: SpatialIndependent + class(Distribution), allocatable :: x + class(Distribution), allocatable :: y + class(Distribution), allocatable :: z + contains + procedure :: sample => spatial_independent_sample + end type SpatialIndependent + + type, extends(SpatialDistribution) :: SpatialBox + real(8) :: lower_left(3) + real(8) :: upper_right(3) + logical :: only_fissionable = .false. + contains + procedure :: sample => spatial_box_sample + end type SpatialBox + + type, extends(SpatialDistribution) :: SpatialPoint + real(8) :: xyz(3) + contains + procedure :: sample => spatial_point_sample + end type SpatialPoint + contains function polar_azimuthal_sample(this) result(uvw) @@ -90,4 +130,31 @@ contains uvw(:) = this%reference_uvw end function monodirectional_sample + function spatial_independent_sample(this) result(xyz) + class(SpatialIndependent), intent(in) :: this + real(8) :: xyz(3) + + xyz(1) = this%x%sample() + xyz(2) = this%y%sample() + xyz(3) = this%z%sample() + end function spatial_independent_sample + + function spatial_box_sample(this) result(xyz) + class(SpatialBox), intent(in) :: this + real(8) :: xyz(3) + + integer :: i + real(8) :: r(3) + + r = [ (prn(), i = 1,3) ] + xyz(:) = this%lower_left + r*(this%upper_right - this%lower_left) + end function spatial_box_sample + + function spatial_point_sample(this) result(xyz) + class(SpatialPoint), intent(in) :: this + real(8) :: xyz(3) + + xyz(:) = this%xyz + end function spatial_point_sample + end module distribution_multivariate diff --git a/src/global.F90 b/src/global.F90 index e8d20a72b..50974282c 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -450,12 +450,9 @@ contains if (allocated(micro_xs)) deallocate(micro_xs) ! Deallocate external source - if (allocated(external_source % params_space)) & - deallocate(external_source % params_space) - if (allocated(external_source % angle)) & - deallocate(external_source % angle) - if (allocated(external_source % energy)) & - deallocate(external_source % energy) + if (allocated(external_source % space)) deallocate(external_source % space) + if (allocated(external_source % angle)) deallocate(external_source % angle) + if (allocated(external_source % energy)) deallocate(external_source % energy) ! Deallocate k and entropy if (allocated(k_generation)) deallocate(k_generation) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 11374b9cf..f4211bfb3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -70,6 +70,7 @@ contains type(Node), pointer :: doc => null() type(Node), pointer :: node_mode => null() type(Node), pointer :: node_source => null() + type(Node), pointer :: node_space => null() type(Node), pointer :: node_angle => null() type(Node), pointer :: node_dist => null() type(Node), pointer :: node_cutoff => null() @@ -376,46 +377,106 @@ contains if (check_for_node(node_source, "space")) then ! Get pointer to spatial distribution - call get_node_ptr(node_source, "space", node_dist) + call get_node_ptr(node_source, "space", node_space) ! Check for type of spatial distribution type = '' - if (check_for_node(node_dist, "type")) & - call get_node_value(node_dist, "type", type) + if (check_for_node(node_space, "type")) & + call get_node_value(node_space, "type", type) select case (to_lower(type)) + case ('independent') + allocate(SpatialIndependent :: external_source%space) + case ('box') - external_source % type_space = SRC_SPACE_BOX - coeffs_reqd = 6 + allocate(SpatialBox :: external_source%space) + case ('fission') - external_source % type_space = SRC_SPACE_FISSION - coeffs_reqd = 6 + allocate(SpatialBox :: external_source%space) + select type(space => external_source%space) + type is (SpatialBox) + space%only_fissionable = .true. + end select + case ('point') - external_source % type_space = SRC_SPACE_POINT - coeffs_reqd = 3 + allocate(SpatialPoint :: external_source%space) + case default call fatal_error("Invalid spatial distribution for external source: "& - &// trim(type)) + // trim(type)) end select - ! Determine number of parameters specified - if (check_for_node(node_dist, "parameters")) then - n = get_arraysize_double(node_dist, "parameters") - else - n = 0 - end if + select type (space => external_source%space) + type is (SpatialIndependent) + ! Read distribution for x coordinate + if (check_for_node(node_space, "x")) then + call get_node_ptr(node_space, "x", node_dist) + call distribution_from_xml(space%x, node_dist) + else + allocate(Discrete :: space%x) + select type (dist => space%x) + type is (Discrete) + allocate(dist%x(1), dist%p(1)) + dist%x(1) = ZERO + dist%p(1) = ONE + end select + end if + + ! Read distribution for y coordinate + if (check_for_node(node_space, "y")) then + call get_node_ptr(node_space, "y", node_dist) + call distribution_from_xml(space%y, node_dist) + else + allocate(Discrete :: space%y) + select type (dist => space%y) + type is (Discrete) + allocate(dist%x(1), dist%p(1)) + dist%x(1) = ZERO + dist%p(1) = ONE + end select + end if + + if (check_for_node(node_space, "z")) then + call get_node_ptr(node_space, "z", node_dist) + call distribution_from_xml(space%z, node_dist) + else + allocate(Discrete :: space%z) + select type (dist => space%z) + type is (Discrete) + allocate(dist%x(1), dist%p(1)) + dist%x(1) = ZERO + dist%p(1) = ONE + end select + end if + + type is (SpatialBox) + ! Make sure correct number of parameters are given + if (get_arraysize_double(node_space, "parameters") /= 6) then + call fatal_error('Box/fission spatial source must have & + &six parameters specified.') + end if + + ! Read lower-right/upper-left coordinates + allocate(temp_real(6)) + call get_node_array(node_space, "parameters", temp_real) + space%lower_left(:) = temp_real(1:3) + space%upper_right(:) = temp_real(4:6) + deallocate(temp_real) + + type is (SpatialPoint) + ! Make sure correct number of parameters are given + if (get_arraysize_double(node_space, "parameters") /= 3) then + call fatal_error('Point spatial source must have & + &three parameters specified.') + end if + + ! Read location of point source + allocate(temp_real(3)) + call get_node_array(node_space, "parameters", temp_real) + space%xyz(:) = temp_real + deallocate(temp_real) + + end select - ! Read parameters for spatial distribution - if (n < coeffs_reqd) then - call fatal_error("Not enough parameters specified for spatial & - &distribution of external source.") - elseif (n > coeffs_reqd) then - call fatal_error("Too many parameters specified for spatial & - &distribution of external source.") - elseif (n > 0) then - allocate(external_source % params_space(n)) - call get_node_array(node_dist, "parameters", & - external_source % params_space) - end if else call fatal_error("No spatial distribution specified for external & &source.") diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index e29b4e22b..6bbd187b5 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -64,33 +64,37 @@ element settings { element seed { xsd:positiveInteger }? & + + element source { - element file { xsd:string { maxLength = "255" } }? & - element space { - (element type { xsd:string { maxLength = "16" } } | - attribute type { xsd:string { maxLength = "16" } }) & - (element interpolation { xsd:string { maxLength = "10" } } | - attribute interpolation { xsd:string { maxLength = "10" } })? & - (element parameters { list { xsd:double+ } } | - attribute parameters { list { xsd:double+ } })? - }? & - element angle { - (element type { xsd:string { maxLength = "16" } } | - attribute type { xsd:string { maxLength = "16" } }) & - (element interpolation { xsd:string { maxLength = "10" } } | - attribute interpolation { xsd:string { maxLength = "10" } })? & - (element parameters { list { xsd:double+ } } | - attribute parameters { list { xsd:double+ } })? - }? & - element energy { - (element type { xsd:string { maxLength = "16" } } | - attribute type { xsd:string { maxLength = "16" } }) & - (element interpolation { xsd:string { maxLength = "10" } } | - attribute interpolation { xsd:string { maxLength = "10" } })? & - (element parameters { list { xsd:double+ } } | - attribute parameters { list { xsd:double+ } })? - }? & - (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? + grammar { + start = + element file { xsd:string { maxLength = "255" } }? & + element space { + (element type { xsd:string } | attribute type { xsd:string }) & + (element parameters { xsd:double+ } | + attribute parameters { xsd:double+ })? & + element x { distribution }? & + element y { distribution }? & + element z { distribution }? + }? & + element angle { + (element type { xsd:string } | attribute type { xsd:string }) & + (element reference_uvw { xsd:double, xsd:double, xsd:double } | + attribute reference_uvw { xsd:double, xsd:double, xsd:double })? & + element mu { distribution }? & + element phi { distribution }? + }? & + element energy { distribution }? & + (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? + distribution = + (element type { xsd:string { maxLength = "16" } } | + attribute type { xsd:string { maxLength = "16" } }) & + (element interpolation { xsd:string { maxLength = "10" } } | + attribute interpolation { xsd:string { maxLength = "10" } })? & + (element parameters { list { xsd:double+ } } | + attribute parameters { list { xsd:double+ } })? + } }? & element state_point { diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 3536f3bd3..9fd2dcf2a 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -266,175 +266,165 @@ - - - - - 255 - - - - - - + + + + + + + 255 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + 16 + + + + + 16 + + + + + + - 16 + 10 - + - 16 + 10 - - - - - 10 - - - - - 10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - - - 16 - + + + + + + - - - 16 - + + + + + + - - - - - 10 - - - - - 10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 16 - - - - - 16 - - - - - - - - 10 - - - - - 10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + diff --git a/src/source.F90 b/src/source.F90 index ae1b7e108..2ac68a314 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -3,6 +3,7 @@ module source use bank_header, only: Bank use constants use distribution_univariate, only: Discrete + use distribution_multivariate, only: SpatialBox use error, only: fatal_error use geometry, only: find_cell use geometry_header, only: BASE_UNIVERSE @@ -113,76 +114,43 @@ contains ! Set the random number generator to the source stream. call prn_set_stream(STREAM_SOURCE) - ! Sample position - select case (external_source%type_space) - case (SRC_SPACE_BOX) - ! Set particle defaults - call p%initialize() - ! Repeat sampling source location until a good site has been found - found = .false. - do while (.not.found) - ! Coordinates sampled uniformly over a box - p_min = external_source%params_space(1:3) - p_max = external_source%params_space(4:6) - r = (/ (prn(), i = 1,3) /) - site%xyz = p_min + r*(p_max - p_min) + ! Set particle defaults + call p%initialize() - ! Fill p with needed data - p%coord(1)%xyz = site%xyz - p%coord(1)%uvw = [ ONE, ZERO, ZERO ] + ! Repeat sampling source location until a good site has been found + found = .false. + do while (.not.found) + ! Sample spatial distribution + site%xyz(:) = external_source%space%sample() - ! Now search to see if location exists in geometry - call find_cell(p, found) - if (.not. found) then - num_resamples = num_resamples + 1 - if (num_resamples == MAX_EXTSRC_RESAMPLES) then - call fatal_error("Maximum number of external source spatial & - &resamples reached!") + ! 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) + if (.not. found) then + num_resamples = num_resamples + 1 + if (num_resamples == MAX_EXTSRC_RESAMPLES) then + call fatal_error("Maximum number of external source spatial & + &resamples reached!") + end if + end if + + ! Check if spatial site is in fissionable material + select type (space => external_source%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 do - call p%clear() + end select + end do - case (SRC_SPACE_FISSION) - ! Repeat sampling source location until a good site has been found - found = .false. - do while (.not.found) - ! Set particle defaults - call p%initialize() - - ! Coordinates sampled uniformly over a box - p_min = external_source%params_space(1:3) - p_max = external_source%params_space(4:6) - r = (/ (prn(), i = 1,3) /) - site%xyz = p_min + r*(p_max - p_min) - - ! 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) - if (.not. found) then - num_resamples = num_resamples + 1 - if (num_resamples == MAX_EXTSRC_RESAMPLES) then - call fatal_error("Maximum number of external source spatial & - &resamples reached!") - end if - cycle - end if - if (p%material == MATERIAL_VOID) then - found = .false. - cycle - end if - if (.not. materials(p%material)%fissionable) found = .false. - end do - call p%clear() - - case (SRC_SPACE_POINT) - ! Point source - site%xyz = external_source%params_space - - end select + call p%clear() ! Sample angle site%uvw(:) = external_source%angle%sample() diff --git a/src/source_header.F90 b/src/source_header.F90 index f52afb919..d87c90589 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -1,7 +1,7 @@ module source_header use distribution_univariate, only: Distribution - use distribution_multivariate, only: UnitSphereDistribution + use distribution_multivariate, only: UnitSphereDistribution, SpatialDistribution implicit none @@ -11,10 +11,9 @@ module source_header !=============================================================================== type ExtSource - integer :: type_space ! spacial distribution, e.g. 'box' or 'point' - real(8), allocatable :: params_space(:) ! parameters for spatial distribution + class(SpatialDistribution), allocatable :: space ! spatial distribution class(UnitSphereDistribution), allocatable :: angle ! angle distribution - class(Distribution), allocatable :: energy ! energy distribution + class(Distribution), allocatable :: energy ! energy distribution end type ExtSource end module source_header