Expand spatial distribution capabilities, notably independent x,y,z

distributions.
This commit is contained in:
Paul Romano 2015-12-17 16:41:12 -06:00
parent 8207028d8c
commit efc7d2ecbf
8 changed files with 372 additions and 308 deletions

View file

@ -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

View file

@ -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

View file

@ -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)

View file

@ -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.")

View file

@ -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 {

View file

@ -266,175 +266,165 @@
</optional>
<optional>
<element name="source">
<interleave>
<optional>
<element name="file">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
<optional>
<element name="space">
<interleave>
<grammar>
<start>
<interleave>
<optional>
<element name="file">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
<optional>
<element name="space">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="parameters">
<oneOrMore>
<data type="double"/>
</oneOrMore>
</element>
<attribute name="parameters">
<oneOrMore>
<data type="double"/>
</oneOrMore>
</attribute>
</choice>
</optional>
<optional>
<element name="x">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="y">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="z">
<ref name="distribution"/>
</element>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="angle">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="reference_uvw">
<data type="double"/>
<data type="double"/>
<data type="double"/>
</element>
<attribute name="reference_uvw">
<group>
<data type="double"/>
<data type="double"/>
<data type="double"/>
</group>
</attribute>
</choice>
</optional>
<optional>
<element name="mu">
<ref name="distribution"/>
</element>
</optional>
<optional>
<element name="phi">
<ref name="distribution"/>
</element>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="energy">
<ref name="distribution"/>
</element>
</optional>
<optional>
<choice>
<element name="type">
<element name="write_initial">
<data type="boolean"/>
</element>
<attribute name="write_initial">
<data type="boolean"/>
</attribute>
</choice>
</optional>
</interleave>
</start>
<define name="distribution">
<interleave>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">16</param>
<param name="maxLength">10</param>
</data>
</element>
<attribute name="type">
<attribute name="interpolation">
<data type="string">
<param name="maxLength">16</param>
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="angle">
<interleave>
</optional>
<optional>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="energy">
<interleave>
<choice>
<element name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</element>
<attribute name="type">
<data type="string">
<param name="maxLength">16</param>
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="parameters">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<choice>
<element name="write_initial">
<data type="boolean"/>
</element>
<attribute name="write_initial">
<data type="boolean"/>
</attribute>
</choice>
</optional>
</interleave>
</optional>
</interleave>
</define>
</grammar>
</element>
</optional>
<optional>

View file

@ -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()

View file

@ -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