Move external source reading to source_header and distribution modules

This commit is contained in:
Paul Romano 2017-08-24 10:56:12 -05:00
parent 1740d9fb8d
commit df60ae0dfd
3 changed files with 287 additions and 243 deletions

View file

@ -1,9 +1,11 @@
module distribution_multivariate
use constants, only: ONE, TWO, PI
use distribution_univariate, only: Distribution
use distribution_univariate
use error, only: fatal_error
use random_lcg, only: prn
use math, only: rotate_angle
use xml_interface
implicit none
@ -58,10 +60,17 @@ module distribution_multivariate
type, abstract :: SpatialDistribution
contains
procedure(spatial_distribution_from_xml_), deferred :: from_xml
procedure(spatial_distribution_sample_), deferred :: sample
end type SpatialDistribution
abstract interface
subroutine spatial_distribution_from_xml_(this, node)
import SpatialDistribution, XMLNode
class(SpatialDistribution), intent(inout) :: this
type(XMLNode), intent(in) :: node
end subroutine spatial_distribution_from_xml_
function spatial_distribution_sample_(this) result(xyz)
import SpatialDistribution
class(SpatialDistribution), intent(in) :: this
@ -74,6 +83,7 @@ module distribution_multivariate
class(Distribution), allocatable :: y
class(Distribution), allocatable :: z
contains
procedure :: from_xml => cartesian_independent_from_xml
procedure :: sample => cartesian_independent_sample
end type CartesianIndependent
@ -82,12 +92,14 @@ module distribution_multivariate
real(8) :: upper_right(3)
logical :: only_fissionable = .false.
contains
procedure :: from_xml => spatial_box_from_xml
procedure :: sample => spatial_box_sample
end type SpatialBox
type, extends(SpatialDistribution) :: SpatialPoint
real(8) :: xyz(3)
contains
procedure :: from_xml => spatial_point_from_xml
procedure :: sample => spatial_point_sample
end type SpatialPoint
@ -132,6 +144,55 @@ contains
uvw(:) = this % reference_uvw
end function monodirectional_sample
subroutine cartesian_independent_from_xml(this, node)
class(CartesianIndependent), intent(inout) :: this
type(XMLNode), intent(in) :: node
type(XMLNode) :: node_dist
! Read distribution for x coordinate
if (check_for_node(node, "x")) then
node_dist = node % child("x")
call distribution_from_xml(this % x, node_dist)
else
allocate(Discrete :: this % x)
select type (dist => this % 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, "y")) then
node_dist = node % child("y")
call distribution_from_xml(this % y, node_dist)
else
allocate(Discrete :: this % y)
select type (dist => this % 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, "z")) then
node_dist = node % child("z")
call distribution_from_xml(this % z, node_dist)
else
allocate(Discrete :: this % z)
select type (dist => this % z)
type is (Discrete)
allocate(dist % x(1), dist % p(1))
dist % x(1) = ZERO
dist % p(1) = ONE
end select
end if
end subroutine cartesian_independent_from_xml
function cartesian_independent_sample(this) result(xyz)
class(CartesianIndependent), intent(in) :: this
real(8) :: xyz(3)
@ -141,6 +202,26 @@ contains
xyz(3) = this % z % sample()
end function cartesian_independent_sample
subroutine spatial_box_from_xml(this, node)
class(SpatialBox), intent(inout) :: this
type(XMLNode), intent(in) :: node
real(8), allocatable :: temp_real(:)
! Make sure correct number of parameters are given
if (node_word_count(node, "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, "parameters", temp_real)
this % lower_left(:) = temp_real(1:3)
this % upper_right(:) = temp_real(4:6)
deallocate(temp_real)
end subroutine spatial_box_from_xml
function spatial_box_sample(this) result(xyz)
class(SpatialBox), intent(in) :: this
real(8) :: xyz(3)
@ -152,6 +233,20 @@ contains
xyz(:) = this % lower_left + r*(this % upper_right - this % lower_left)
end function spatial_box_sample
subroutine spatial_point_from_xml(this, node)
class(SpatialPoint), intent(inout) :: this
type(XMLNode), intent(in) :: node
! Make sure correct number of parameters are given
if (node_word_count(node, "parameters") /= 3) then
call fatal_error('Point spatial source must have &
&three parameters specified.')
end if
! Read location of point source
call get_node_array(node, "parameters", this % xyz)
end subroutine spatial_point_from_xml
function spatial_point_sample(this) result(xyz)
class(SpatialPoint), intent(in) :: this
real(8) :: xyz(3)

View file

@ -147,7 +147,6 @@ contains
integer(C_INT32_T) :: i_start, i_end
integer(C_INT) :: err
integer, allocatable :: temp_int_array(:)
real(8), allocatable :: temp_real(:)
integer :: n_tracks
logical :: file_exists
character(MAX_WORD_LEN) :: type
@ -155,10 +154,6 @@ contains
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_mode
type(XMLNode) :: node_source
type(XMLNode) :: node_space
type(XMLNode) :: node_angle
type(XMLNode) :: node_dist
type(XMLNode) :: node_cutoff
type(XMLNode) :: node_entropy
type(XMLNode) :: node_ufs
@ -400,241 +395,14 @@ contains
allocate(external_source(n))
end if
! Check if we want to write out source
if (check_for_node(root, "write_initial_source")) then
call get_node_value(root, "write_initial_source", write_initial_source)
end if
! Read each source
do i = 1, n
! Get pointer to source
node_source = node_source_list(i)
! Check if we want to write out source
if (check_for_node(node_source, "write_initial")) then
call get_node_value(node_source, "write_initial", write_initial_source)
end if
! Check for source strength
if (check_for_node(node_source, "strength")) then
call get_node_value(node_source, "strength", external_source(i)%strength)
else
external_source(i)%strength = ONE
end if
! Check for external source file
if (check_for_node(node_source, "file")) then
! Copy path of source file
call get_node_value(node_source, "file", path_source)
! Check if source file exists
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Source file '" // trim(path_source) &
// "' does not exist!")
end if
else
! Spatial distribution for external source
if (check_for_node(node_source, "space")) then
! Get pointer to spatial distribution
node_space = node_source % child("space")
! Check for type of spatial distribution
type = ''
if (check_for_node(node_space, "type")) &
call get_node_value(node_space, "type", type)
select case (to_lower(type))
case ('cartesian')
allocate(CartesianIndependent :: external_source(i)%space)
case ('box')
allocate(SpatialBox :: external_source(i)%space)
case ('fission')
allocate(SpatialBox :: external_source(i)%space)
select type(space => external_source(i)%space)
type is (SpatialBox)
space%only_fissionable = .true.
end select
case ('point')
allocate(SpatialPoint :: external_source(i)%space)
case default
call fatal_error("Invalid spatial distribution for external source: "&
// trim(type))
end select
select type (space => external_source(i)%space)
type is (CartesianIndependent)
! Read distribution for x coordinate
if (check_for_node(node_space, "x")) then
node_dist = node_space % child("x")
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
node_dist = node_space % child("y")
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
node_dist = node_space % child("z")
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 (node_word_count(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 (node_word_count(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
else
! If no spatial distribution specified, make it a point source
allocate(SpatialPoint :: external_source(i) % space)
select type (space => external_source(i) % space)
type is (SpatialPoint)
space % xyz(:) = [ZERO, ZERO, ZERO]
end select
end if
! Determine external source angular distribution
if (check_for_node(node_source, "angle")) then
! Get pointer to angular distribution
node_angle = node_source % child("angle")
! Check for type of angular distribution
type = ''
if (check_for_node(node_angle, "type")) &
call get_node_value(node_angle, "type", type)
select case (to_lower(type))
case ('isotropic')
allocate(Isotropic :: external_source(i)%angle)
case ('monodirectional')
allocate(Monodirectional :: external_source(i)%angle)
case ('mu-phi')
allocate(PolarAzimuthal :: external_source(i)%angle)
case default
call fatal_error("Invalid angular distribution for external source: "&
// trim(type))
end select
! Read reference directional unit vector
if (check_for_node(node_angle, "reference_uvw")) then
n = node_word_count(node_angle, "reference_uvw")
if (n /= 3) then
call fatal_error('Angular distribution reference direction must have &
&three parameters specified.')
end if
call get_node_array(node_angle, "reference_uvw", &
external_source(i)%angle%reference_uvw)
else
! By default, set reference unit vector to be positive z-direction
external_source(i)%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Read parameters for angle distribution
select type (angle => external_source(i)%angle)
type is (Monodirectional)
call get_node_array(node_angle, "reference_uvw", &
external_source(i)%angle%reference_uvw)
type is (PolarAzimuthal)
if (check_for_node(node_angle, "mu")) then
node_dist = node_angle % child("mu")
call distribution_from_xml(angle%mu, node_dist)
else
allocate(Uniform :: angle%mu)
select type (mu => angle%mu)
type is (Uniform)
mu%a = -ONE
mu%b = ONE
end select
end if
if (check_for_node(node_angle, "phi")) then
node_dist = node_angle % child("phi")
call distribution_from_xml(angle%phi, node_dist)
else
allocate(Uniform :: angle%phi)
select type (phi => angle%phi)
type is (Uniform)
phi%a = ZERO
phi%b = TWO*PI
end select
end if
end select
else
! Set default angular distribution isotropic
allocate(Isotropic :: external_source(i)%angle)
external_source(i)%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Determine external source energy distribution
if (check_for_node(node_source, "energy")) then
node_dist = node_source % child("energy")
call distribution_from_xml(external_source(i)%energy, node_dist)
else
! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
allocate(Watt :: external_source(i)%energy)
select type(energy => external_source(i)%energy)
type is (Watt)
energy%a = 0.988e6_8
energy%b = 2.249e-6_8
end select
end if
end if
call external_source(i) % from_xml(node_source_list(i), path_source)
end do
! Survival biasing

View file

@ -1,20 +1,201 @@
module source_header
use distribution_univariate, only: Distribution
use distribution_multivariate, only: UnitSphereDistribution, SpatialDistribution
use constants
use distribution_univariate
use distribution_multivariate
use error, only: fatal_error
use string, only: to_lower
use xml_interface
implicit none
private
!===============================================================================
! SOURCEDISTRIBUTION describes an external source of particles for a
! fixed-source problem or for the starting source in a k eigenvalue problem
!===============================================================================
type SourceDistribution
real(8) :: strength ! source strength
type, public :: SourceDistribution
real(8) :: strength = ONE ! source strength
class(SpatialDistribution), allocatable :: space ! spatial distribution
class(UnitSphereDistribution), allocatable :: angle ! angle distribution
class(Distribution), allocatable :: energy ! energy distribution
contains
procedure :: from_xml => source_from_xml
end type SourceDistribution
! External source distributions
type(SourceDistribution), public, allocatable :: external_source(:)
contains
subroutine source_from_xml(this, node, path_source)
class(SourceDistribution), intent(inout) :: this
type(XMLNode), intent(in) :: node
character(MAX_FILE_LEN), intent(out) :: path_source
integer :: n
logical :: file_exists
character(MAX_WORD_LEN) :: type
type(XMLNode) :: node_space
type(XMLNode) :: node_angle
type(XMLNode) :: node_dist
! Check for source strength
if (check_for_node(node, "strength")) then
call get_node_value(node, "strength", this % strength)
end if
! Check for external source file
if (check_for_node(node, "file")) then
! Copy path of source file
call get_node_value(node, "file", path_source)
! Check if source file exists
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Source file '" // trim(path_source) &
// "' does not exist!")
end if
else
! Spatial distribution for external source
if (check_for_node(node, "space")) then
! Get pointer to spatial distribution
node_space = node % child("space")
! Check for type of spatial distribution
type = ''
if (check_for_node(node_space, "type")) &
call get_node_value(node_space, "type", type)
select case (to_lower(type))
case ('cartesian')
allocate(CartesianIndependent :: this % space)
case ('box')
allocate(SpatialBox :: this % space)
case ('fission')
allocate(SpatialBox :: this % space)
select type(space => this % space)
type is (SpatialBox)
space % only_fissionable = .true.
end select
case ('point')
allocate(SpatialPoint :: this % space)
case default
call fatal_error("Invalid spatial distribution for external source: "&
// trim(type))
end select
! Read spatial distribution from XML
call this % space % from_xml(node_space)
else
! If no spatial distribution specified, make it a point source
allocate(SpatialPoint :: this % space)
select type (space => this % space)
type is (SpatialPoint)
space % xyz(:) = [ZERO, ZERO, ZERO]
end select
end if
! Determine external source angular distribution
if (check_for_node(node, "angle")) then
! Get pointer to angular distribution
node_angle = node % child("angle")
! Check for type of angular distribution
type = ''
if (check_for_node(node_angle, "type")) &
call get_node_value(node_angle, "type", type)
select case (to_lower(type))
case ('isotropic')
allocate(Isotropic :: this % angle)
case ('monodirectional')
allocate(Monodirectional :: this % angle)
case ('mu-phi')
allocate(PolarAzimuthal :: this % angle)
case default
call fatal_error("Invalid angular distribution for external source: "&
// trim(type))
end select
! Read reference directional unit vector
if (check_for_node(node_angle, "reference_uvw")) then
n = node_word_count(node_angle, "reference_uvw")
if (n /= 3) then
call fatal_error('Angular distribution reference direction must have &
&three parameters specified.')
end if
call get_node_array(node_angle, "reference_uvw", &
this % angle % reference_uvw)
else
! By default, set reference unit vector to be positive z-direction
this % angle % reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Read parameters for angle distribution
select type (angle => this % angle)
type is (Monodirectional)
call get_node_array(node_angle, "reference_uvw", &
this % angle % reference_uvw)
type is (PolarAzimuthal)
if (check_for_node(node_angle, "mu")) then
node_dist = node_angle % child("mu")
call distribution_from_xml(angle % mu, node_dist)
else
allocate(Uniform :: angle%mu)
select type (mu => angle%mu)
type is (Uniform)
mu % a = -ONE
mu % b = ONE
end select
end if
if (check_for_node(node_angle, "phi")) then
node_dist = node_angle % child("phi")
call distribution_from_xml(angle % phi, node_dist)
else
allocate(Uniform :: angle%phi)
select type (phi => angle%phi)
type is (Uniform)
phi % a = ZERO
phi % b = TWO*PI
end select
end if
end select
else
! Set default angular distribution isotropic
allocate(Isotropic :: this % angle)
this % angle % reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Determine external source energy distribution
if (check_for_node(node, "energy")) then
node_dist = node % child("energy")
call distribution_from_xml(this % energy, node_dist)
else
! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
allocate(Watt :: this % energy)
select type(energy => this % energy)
type is (Watt)
energy % a = 0.988e6_8
energy % b = 2.249e-6_8
end select
end if
end if
end subroutine source_from_xml
end module source_header