crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Support for tabular angular and energy distributions.

Browse source 
Implemented new polymorphic Distribution type which is used by
AngleDistribution, both of which appear in distribution_header.
This commit is contained in:
Paul Romano 2015-11-09 16:41:44 -06:00
parent c737c1f692
commit 2abc9f1e29
14 changed files with 563 additions and 257 deletions
Download patch file Download diff file Expand all files Collapse all files

60
docs/source/usersguide/input.rst
View file

@ -440,7 +440,6 @@ attributes/sub-elements:
has the following attributes:
:type:
The type of spatial distribution. Valid options are "box", "fission", and
"point". A "box" spatial distribution has coordinates sampled uniformly in
a parallelepiped. A "fission" spatial distribution samples locations from
@ -468,14 +467,23 @@ attributes/sub-elements:
has the following attributes:
:type:
The type of angular distribution. Valid options are "isotropic" and
"monodirectional". The angle of the particle emitted from a source site is
isotropic if the "isotropic" option is given. The angle of the particle
emitted from a source site is the direction specified in the <parameters>
attribute if "monodirectional" option is given.
The type of angular distribution. Valid options are "isotropic",
"monodirectional", and "tabular". The angle of the particle emitted from a
source site is isotropic if the "isotropic" option is given. The angle of
the particle emitted from a source site is the direction specified in the
<parameters> attribute if "monodirectional" option is given. The "tabular"
option produces directions with polar angles sampled from a tabulated
distribution.
*Default*: isotropic
:interpolation:
For a "tabular" angular distribution, ``interpolation`` can be set to
"histogram" or "linear-linear" thereby specifying how tabular points are
to be interpolated.
*Default*: histogram
:parameters:
For an "isotropic" angular distribution, ``parameters`` should not be
specified.
@ -484,6 +492,22 @@ attributes/sub-elements:
given as three real numbers which specify the angular cosines with respect
to each axis.
For a "tabular" angular distribution, ``parameters`` provides the
:math:`(\mu,p)` pairs defining the tabular distribution. All :math:`\mu`
points are given first followed by corresponding :math:`p` points. The
following example gives a histogram distribution with even probability of
selecting a polar angle in the range [-1,-0.5] and [0.5,1] (Note that the
last :math:`p` point is inconsequential):
.. code-block:: xml
<angle type="tabular" interpolation="histogram">
<parameters>
-1.0 -0.5 0.5 1.0
1.0 0.0 1.0 0.0
</parameters>
</angle>
*Default*: None
:energy:
@ -491,15 +515,23 @@ attributes/sub-elements:
has the following attributes:
:type:
The type of energy distribution. Valid options are "monoenergetic",
"watt", and "maxwell". The "monoenergetic" option produces source sites at
a single energy. The "watt" option produces source sites whose energy is
sampled from a Watt fission spectrum. The "maxwell" option produce source
sites whose energy is sampled from a Maxwell fission spectrum.
"watt", "maxwell", and "tabular". The "monoenergetic" option produces
source sites at a single energy. The "watt" option produces source sites
whose energy is sampled from a Watt fission spectrum. The "maxwell" option
produce source sites whose energy is sampled from a Maxwell fission
spectrum. The "tabular" option produces source sites whose energy is
sampled from a tabulated distribution.
*Default*: watt
:interpolation:
For a "tabular" angular distribution, ``interpolation`` can be set to
"histogram" or "linear-linear" thereby specifying how tabular points are
to be interpolated.
*Default*: histogram
:parameters:
For a "monoenergetic" energy distribution, ``parameters`` should be
given as the energy in MeV of the source sites.
@ -512,6 +544,10 @@ attributes/sub-elements:
real number :math:`a` that parameterizes the distribution :math:`p(E) dE =
c E e^{-E/a} dE`.
For a "tabular" energy distribution, ``parameters`` provides the
:math:`(E,p)` pairs defining the tabular distribution. All :math:`E`
points are given first followed by corresponding :math:`p` points.
*Default*: 0.988 2.249
:write_initial:
@ -519,7 +555,7 @@ attributes/sub-elements:
the beginning of the first batch. The output file is named
"initial_source.binary(h5)"
*Default*: false
*Default*: false
``<state_point>`` Element
-------------------------

78
openmc/settings.py
View file

@ -139,8 +139,10 @@ class SettingsFile(object):
self._source_space_type = None
self._source_space_params = None
self._source_angle_type = None
self._source_angle_interpolation = None
self._source_angle_params = None
self._source_energy_type = None
self._source_energy_interpolation = None
self._source_energy_params = None
self._confidence_intervals = None
@ -243,6 +245,10 @@ class SettingsFile(object):
def source_angle_type(self):
return self._source_angle_type
@property
def source_angle_interpolation(self):
return self._source_angle_interpolation
@property
def source_angle_params(self):
return self._source_angle_params
@ -251,6 +257,10 @@ class SettingsFile(object):
def source_energy_type(self):
return self._source_energy_type
@property
def source_energy_interpolation(self):
return self._source_energy_interpolation
@property
def source_energy_params(self):
return self._source_energy_params
@ -511,17 +521,19 @@ class SettingsFile(object):
self._source_space_type = stype
self._source_space_params = params
def set_source_angle(self, stype, params=[]):
def set_source_angle(self, stype, params=[], interp='histogram'):
"""Defined the angular distribution of the external/starting source.
Parameters
----------
stype : str
The type of angular distribution. Valid options are "isotropic" and
"monodirectional". The angle of the particle emitted from a source
site is isotropic if the "isotropic" option is given. The angle of
the particle emitted from a source site is the direction specified
in ``params`` if the "monodirectional" option is given.
The type of angular distribution. Valid options are "isotropic",
"monodirectional", and "tabular". The angle of the particle emitted
from a source site is isotropic if the "isotropic" option is
given. The angle of the particle emitted from a source site is the
direction specified in ``params`` if the "monodirectional" option is
given. The "tabular" option produces directions with polar angles
sampled from a tabulated distribution.
params : Iterable of float
For an "isotropic" angular distribution, ``params`` should not
be specified.
@ -530,11 +542,20 @@ class SettingsFile(object):
be given as three floats which specify the angular cosines
with respect to each axis.
For a "tabular" angular distribution, ``parameters`` provides the
:math:`(\mu,p)` pairs defining the tabular distribution. All
:math:`\mu` points are given first followed by corresponding
:math:`p` points.
interp : { 'histogram', 'linear-linear' }
For a "tabular" angular distribution, ``interpolation`` can be set
to "histogram" or "linear-linear" thereby specifying how tabular
points are to be interpolated.
"""
check_type('source angle type', stype, basestring)
check_value('source angle type', stype,
['isotropic', 'monodirectional'])
['isotropic', 'monodirectional', 'tabular'])
check_type('source angle parameters', params, Iterable, Real)
if stype == 'isotropic' and params is not None:
msg = 'Unable to set source angle parameters since they are not ' \
@ -543,22 +564,30 @@ class SettingsFile(object):
elif stype == 'monodirectional':
check_length('source angle parameters for a monodirectional '
'source', params, 3)
elif stype == 'tabular':
check_type('source angle interpolation', interp, basestring)
check_value('source angle interpolation', interp,
['histogram', 'linear-linear'])
self._source_angle_interpolation = interp
self._source_angle_type = stype
self._source_angle_params = params
def set_source_energy(self, stype, params=[]):
def set_source_energy(self, stype, params=[], interp='histogram'):
"""Defined the energy distribution of the external/starting source.
Parameters
----------
stype : str
The type of energy distribution. Valid options are "monoenergetic",
"watt", and "maxwell". The "monoenergetic" option produces source
sites at a single energy. The "watt" option produces source sites
whose energy is sampled from a Watt fission spectrum. The "maxwell"
option produce source sites whose energy is sampled from a Maxwell
fission spectrum.
"watt", "maxwell", and "tabular". The "monoenergetic" option
produces source sites at a single energy. The "watt" option produces
source sites whose energy is sampled from a Watt fission
spectrum. The "maxwell" option produce source sites whose energy is
sampled from a Maxwell fission spectrum. The "tabular" option
produces source sites whose energy is sampled from a tabulated
distribution.
params : Iterable of float
For a "monoenergetic" energy distribution, ``params`` should be
given as the energy in MeV of the source sites.
@ -571,11 +600,19 @@ class SettingsFile(object):
one real number :math:`a` that parameterizes the distribution
:math:`p(E) dE = c E e^{-E/a} dE`.
For a "tabular" energy distribution, ``parameters`` provides the
:math:`(E,p)` pairs defining the tabular distribution. All :math:`E`
points are given first followed by corresponding :math:`p` points.
interp : { 'histogram', 'linear-linear' }
For a "tabular" energy distribution, ``interpolation`` can be set
to "histogram" or "linear-linear" thereby specifying how tabular
points are to be interpolated.
"""
check_type('source energy type', stype, basestring)
check_value('source energy type', stype,
['monoenergetic', 'watt', 'maxwell'])
['monoenergetic', 'watt', 'maxwell', 'tabular'])
check_type('source energy parameters', params, Iterable, Real)
if stype in ['monoenergetic', 'maxwell']:
check_length('source energy parameters for a monoenergetic '
@ -583,6 +620,11 @@ class SettingsFile(object):
elif stype == 'watt':
check_length('source energy parameters for a Watt source',
params, 2)
elif stype == 'tabular':
check_type('source energy interpolation', interp, basestring)
check_value('source energy interpolation', interp,
['histogram', 'linear-linear'])
self._source_energy_interpolation = interp
self._source_energy_type = stype
self._source_energy_params = params
@ -949,6 +991,9 @@ class SettingsFile(object):
element = ET.SubElement(self._source_subelement, "angle")
element.set("type", self._source_angle_type)
if self.source_angle_interpolation is not None:
element.set("interpolation", self.source_angle_interpolation)
subelement = ET.SubElement(element, "parameters")
subelement.text = ' '.join(map(str, self._source_angle_params))
@ -961,7 +1006,10 @@ class SettingsFile(object):
element = ET.SubElement(self._source_subelement, "energy")
element.set("type", self._source_energy_type)
subelement = ET.SubElement(element, "parameters")
if self.source_energy_interpolation is not None:
element.set("interpolation", self.source_energy_interpolation)
subelement = ET.SubElement(element, "parameters")
subelement.text = ' '.join(map(str, self._source_energy_params))
def _create_output_subelement(self):

216
src/distribution_header.F90 Normal file
View file

@ -0,0 +1,216 @@
module distribution_header
use constants, only: ZERO, HALF, ONE, HISTOGRAM, LINEAR_LINEAR
use error, only: fatal_error
use math, only: rotate_angle, maxwell_spectrum, watt_spectrum
use random_lcg, only: prn
!===============================================================================
! DISTRIBUTION type defines a probability density function
!===============================================================================
type, abstract :: Distribution
contains
procedure(iSample), deferred :: sample
end type Distribution
abstract interface
function iSample(this) result(x)
import Distribution
class(Distribution), intent(in) :: this
real(8) :: x
end function iSample
end interface
!===============================================================================
! Derived classes of Distribution
!===============================================================================
! delta function at a single point
type, extends(Distribution) :: Delta
real(8) :: x0
contains
procedure :: sample => delta_sample
end type Delta
! Uniform distribution over the interval [a,b]
type, extends(Distribution) :: Uniform
real(8) :: a
real(8) :: b
contains
procedure :: sample => uniform_sample
end type Uniform
! Maxwellian distribution of form c*E*exp(-E/a)
type, extends(Distribution) :: Maxwell
real(8) :: theta
contains
procedure :: sample => maxwell_sample
end type Maxwell
! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E))
type, extends(Distribution) :: Watt
real(8) :: a
real(8) :: b
contains
procedure :: sample => watt_sample
end type Watt
! Histogram or linear-linear interpolated tabular distribution
type, extends(Distribution) :: Tabular
integer :: interpolation
real(8), allocatable :: x(:) ! tabulated independent variable
real(8), allocatable :: p(:) ! tabulated probability density
real(8), allocatable, private :: c(:) ! cumulative distribution at tabulated values
contains
procedure :: sample => tabular_sample
procedure :: initialize => tabular_initialize
end type Tabular
!===============================================================================
! AngleDistribution
!===============================================================================
type :: AngleDistribution
real(8) :: reference_uvw(3)
class(Distribution), allocatable :: mu
contains
procedure :: sample => angle_sample
end type AngleDistribution
contains
function delta_sample(this) result(x)
class(Delta), intent(in) :: this
real(8) :: x
x = this%x0
end function delta_sample
function uniform_sample(this) result(x)
class(Uniform), intent(in) :: this
real(8) :: x
x = this%a + prn()*(this%b - this%a)
end function uniform_sample
function maxwell_sample(this) result(x)
class(Maxwell), intent(in) :: this
real(8) :: x
x = maxwell_spectrum(this%theta)
end function maxwell_sample
function watt_sample(this) result(x)
class(Watt), intent(in) :: this
real(8) :: x
x = watt_spectrum(this%a, this%b)
end function watt_sample
function tabular_sample(this) result(x)
class(Tabular), intent(in) :: this
real(8) :: x
integer :: i
real(8) :: c ! sampled cumulative frequency
real(8) :: m ! slope of PDF
real(8) :: x_i, x_i1 ! i-th and (i+1)th x values
real(8) :: c_i, c_i1 ! i-th and (i+1)th cumulative distribution values
real(8) :: p_i, p_i1 ! i-th and (i+1)th probability density values
! Sample value of CDF
c = prn()
! Find first CDF bin which is above the sampled value
c_i = this%c(1)
do i = 1, size(this%c) - 1
c_i1 = this%c(i + 1)
if (c <= c_i1) exit
c_i = c_i1
end do
! Determine bounding PDF values
x_i = this%x(i)
p_i = this%p(i)
if (this%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_i > ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i
end if
else
! Linear-linear interpolation
x_i1 = this%x(i + 1)
p_i1 = this%p(i + 1)
m = (p_i1 - p_i)/(x_i1 - x_i)
if (m == ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i + (sqrt(max(ZERO, p_i*p_i + 2*m*(c - c_i))) - p_i)/m
end if
end if
end function tabular_sample
subroutine tabular_initialize(this, x, p, interp)
class(Tabular), intent(inout) :: this
real(8), intent(in) :: x(:)
real(8), intent(in) :: p(:)
integer, intent(in) :: interp
integer :: n
! Check interpolation parameter
if (interp /= HISTOGRAM .and. interp /= LINEAR_LINEAR) then
call fatal_error('Only histogram and linear-linear interpolation for tabular &
&distribution is supported.')
end if
! Check length of x, p arrays
if (size(x) /= size(p)) then
call fatal_error('Tabulated probabilities not of same length as &
&independent variable.')
end if
! Copy probability density function and interpolation parameter
n = size(x)
allocate(this%x(n), this%p(n), this%c(n))
this%interpolation = interp
this%x(:) = x(:)
this%p(:) = p(:)
! Calculate cumulative distribution function
this%c(1) = ZERO
do i = 2, n
if (this%interpolation == HISTOGRAM) then
this%c(i) = this%c(i-1) + this%p(i-1)*(this%x(i) - this%x(i-1))
elseif (this%interpolation == LINEAR_LINEAR) then
this%c(i) = this%c(i-1) + HALF*(this%p(i-1) + this%p(i)) * &
(this%x(i) - this%x(i-1))
end if
end do
! Normalize density and distribution functions
this%p(:) = this%p(:)/this%c(n)
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
function angle_sample(this) result(uvw)
class(AngleDistribution), intent(in) :: this
real(8) :: uvw(3)
real(8) :: mu
mu = this%mu%sample()
if (mu == ONE) then
uvw(:) = this%reference_uvw
else
uvw(:) = rotate_angle(this%reference_uvw, mu)
end if
end function angle_sample
end module distribution_header

15
src/error.F90
View file

@ -3,8 +3,6 @@ module error
use, intrinsic :: ISO_FORTRAN_ENV
use constants
use global
#ifdef MPI
use message_passing
#endif
@ -87,6 +85,9 @@ contains
integer :: line_wrap ! length of line
integer :: length ! length of message
integer :: indent ! length of indentation
#ifdef MPI
integer :: mpi_err
#endif
! set default error code
@ -136,16 +137,6 @@ contains
end if
end do
! Write information on current batch, generation, and particle
if (current_batch > 0) then
write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch
write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen
write(ERROR_UNIT,*)
end if
! Release memory from all allocatable arrays
call free_memory()
#ifdef MPI
! Abort MPI
call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err)

11
src/global.F90
View file

@ -289,9 +289,6 @@ module global
character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart
character(MAX_FILE_LEN) :: path_output = '' ! Path to output directory
! Random number seed
integer(8) :: seed = 1_8
! The verbosity controls how much information will be printed to the
! screen and in logs
integer :: verbosity = 7
@ -455,10 +452,10 @@ contains
! Deallocate external source
if (allocated(external_source % params_space)) &
deallocate(external_source % params_space)
if (allocated(external_source % params_angle)) &
deallocate(external_source % params_angle)
if (allocated(external_source % params_energy)) &
deallocate(external_source % params_energy)
if (allocated(external_source % angle % mu)) &
deallocate(external_source % angle % mu)
if (allocated(external_source % energy)) &
deallocate(external_source % energy)
! Deallocate k and entropy
if (allocated(k_generation)) deallocate(k_generation)

214
src/input_xml.F90
View file

@ -3,6 +3,7 @@ module input_xml
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
use distribution_header
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
@ -11,7 +12,7 @@ module input_xml
use mesh_header, only: RegularMesh
use output, only: write_message
use plot_header
use random_lcg, only: prn
use random_lcg, only: prn, seed
use surface_header
use stl_vector, only: VectorInt
use string, only: to_lower, to_str, str_to_int, str_to_real, &
@ -59,6 +60,7 @@ contains
integer :: temp_int_array3(3)
integer, allocatable :: temp_int_array(:)
integer(8) :: temp_long
real(8), allocatable :: temp_real(:)
integer :: n_tracks
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable
@ -423,24 +425,6 @@ contains
! Get pointer to angular distribution
call get_node_ptr(node_source, "angle", node_dist)
! Check for type of angular distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('isotropic')
external_source % type_angle = SRC_ANGLE_ISOTROPIC
coeffs_reqd = 0
case ('monodirectional')
external_source % type_angle = SRC_ANGLE_MONO
coeffs_reqd = 3
case ('tabular')
external_source % type_angle = SRC_ANGLE_TABULAR
case default
call fatal_error("Invalid angular distribution for external source: "&
&// trim(type))
end select
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
@ -448,21 +432,76 @@ contains
n = 0
end if
! Check for type of angular distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('isotropic')
allocate(Uniform :: external_source%angle%mu)
case ('monodirectional')
allocate(Delta :: external_source%angle%mu)
if (n /= 3) then
call fatal_error('Monodirectional angular distribution must have &
&three parameters specified.')
end if
case ('tabular')
allocate(Tabular :: external_source%angle%mu)
case default
call fatal_error("Invalid angular distribution for external source: "&
// trim(type))
end select
! Set reference unit vector to be positive z-direction
external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
! Read parameters for angle distribution
if (n < coeffs_reqd) then
call fatal_error("Not enough parameters specified for angle &
&distribution of external source.")
elseif (n > coeffs_reqd) then
call fatal_error("Too many parameters specified for angle &
&distribution of external source.")
elseif (n > 0) then
allocate(external_source % params_angle(n))
select type (mu => external_source%angle%mu)
type is (Uniform)
mu%a = -ONE
mu%b = ONE
type is (Delta)
mu%x0 = ONE
call get_node_array(node_dist, "parameters", &
external_source % params_angle)
end if
external_source%angle%reference_uvw)
type is (Tabular)
! Read interpolation
if (check_for_node(node_source, "interpolation")) then
call get_node_value(node_source, "interpolation", temp_str)
select case(to_lower(temp_str))
case ('histogram')
temp_int = HISTOGRAM
case ('linear-linear')
temp_int = LINEAR_LINEAR
case default
call fatal_error("Unknown interpolation type for source &
&angular distribution: " // trim(temp_str))
end select
else
temp_int = HISTOGRAM
end if
! Read and initialize tabular distribution
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call mu%initialize(temp_real(1:n), temp_real(n+1:2*n), temp_int)
deallocate(temp_real)
end select
else
! Set default angular distribution isotropic
external_source % type_angle = SRC_ANGLE_ISOTROPIC
allocate(Uniform :: external_source%angle%mu)
select type(mu => external_source%angle%mu)
type is (Uniform)
mu%a = -ONE
mu%b = ONE
end select
external_source%angle%reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Determine external source energy distribution
@ -471,27 +510,6 @@ contains
! Get pointer to energy distribution
call get_node_ptr(node_source, "energy", node_dist)
! Check for type of energy distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('monoenergetic')
external_source % type_energy = SRC_ENERGY_MONO
coeffs_reqd = 1
case ('maxwell')
external_source % type_energy = SRC_ENERGY_MAXWELL
coeffs_reqd = 1
case ('watt')
external_source % type_energy = SRC_ENERGY_WATT
coeffs_reqd = 2
case ('tabular')
external_source % type_energy = SRC_ENERGY_TABULAR
case default
call fatal_error("Invalid energy distribution for external source: " &
&// trim(type))
end select
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = get_arraysize_double(node_dist, "parameters")
@ -499,23 +517,87 @@ contains
n = 0
end if
! Check for type of energy distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('monoenergetic')
allocate(Delta :: external_source%energy)
if (n /= 1) then
call fatal_error('Monoenergetic energy distribution must have one &
&parameter specified.')
end if
case ('maxwell')
allocate(Maxwell :: external_source%energy)
if (n /= 1) then
call fatal_error('Maxwell energy distribution must have one &
&parameter specified.')
end if
case ('watt')
allocate(Watt :: external_source%energy)
if (n /= 2) then
call fatal_error('Watt energy distribution must have two &
&parameter specified.')
end if
case ('tabular')
allocate(Tabular :: external_source%energy)
case default
call fatal_error("Invalid energy distribution for external source: " &
// trim(type))
end select
! Read parameters for energy distribution
if (n < coeffs_reqd) then
call fatal_error("Not enough parameters specified for energy &
&distribution of external source.")
elseif (n > coeffs_reqd) then
call fatal_error("Too many parameters specified for energy &
&distribution of external source.")
elseif (n > 0) then
allocate(external_source % params_energy(n))
call get_node_array(node_dist, "parameters", &
external_source % params_energy)
end if
select type(energy => external_source%energy)
type is (Delta)
call get_node_value(node_dist, "parameters", energy%x0)
type is (Maxwell)
call get_node_value(node_dist, "parameters", energy%theta)
type is (Watt)
allocate(temp_real(2))
call get_node_array(node_dist, "parameters", temp_real)
energy%a = temp_real(1)
energy%b = temp_real(2)
deallocate(temp_real)
type is (Tabular)
! Read interpolation
if (check_for_node(node_dist, "interpolation")) then
call get_node_value(node_dist, "interpolation", temp_str)
select case(to_lower(temp_str))
case ('histogram')
temp_int = HISTOGRAM
case ('linear-linear')
temp_int = LINEAR_LINEAR
case default
call fatal_error("Unknown interpolation type for source &
&angular distribution: " // trim(temp_str))
end select
else
temp_int = HISTOGRAM
end if
! Read and initialize tabular distribution
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call energy%initialize(temp_real(1:n), temp_real(n+1:2*n), temp_int)
deallocate(temp_real)
end select
else
! Set default energy distribution to Watt fission spectrum
external_source % type_energy = SRC_ENERGY_WATT
allocate(external_source % params_energy(2))
external_source % params_energy = (/ 0.988_8, 2.249_8 /)
! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
allocate(Watt :: external_source%energy)
select type(energy => external_source%energy)
type is (Watt)
energy%a = 0.988_8
energy%b = 2.249_8
end select
end if
end if

49
src/math.F90
View file

@ -557,6 +557,55 @@ contains
end function calc_rn
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi)
phi = TWO * PI * prn()
! Precompute factors to save flops
sinphi = sin(phi)
cosphi = cos(phi)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
!===============================================================================
! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based
! on a direct sampling scheme. The probability distribution function for a

51
src/physics.F90
View file

@ -9,7 +9,7 @@ module physics
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use math, only: maxwell_spectrum, watt_spectrum
use math, only: rotate_angle, maxwell_spectrum, watt_spectrum
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: Particle
@ -1517,55 +1517,6 @@ contains
end function sample_angle
!===============================================================================
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
! mu and through an azimuthal angle sampled uniformly. Note that this is done
! with direct sampling rather than rejection as is done in MCNP and SERPENT.
!===============================================================================
function rotate_angle(uvw0, mu) result(uvw)
real(8), intent(in) :: uvw0(3) ! directional cosine
real(8), intent(in) :: mu ! cosine of angle in lab or CM
real(8) :: uvw(3) ! rotated directional cosine
real(8) :: phi ! azimuthal angle
real(8) :: sinphi ! sine of azimuthal angle
real(8) :: cosphi ! cosine of azimuthal angle
real(8) :: a ! sqrt(1 - mu^2)
real(8) :: b ! sqrt(1 - w^2)
real(8) :: u0 ! original cosine in x direction
real(8) :: v0 ! original cosine in y direction
real(8) :: w0 ! original cosine in z direction
! Copy original directional cosines
u0 = uvw0(1)
v0 = uvw0(2)
w0 = uvw0(3)
! Sample azimuthal angle in [0,2pi)
phi = TWO * PI * prn()
! Precompute factors to save flops
sinphi = sin(phi)
cosphi = cos(phi)
a = sqrt(max(ZERO, ONE - mu*mu))
b = sqrt(max(ZERO, ONE - w0*w0))
! Need to treat special case where sqrt(1 - w**2) is close to zero by
! expanding about the v component rather than the w component
if (b > 1e-10) then
uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
uvw(3) = mu*w0 - a*b*cosphi
else
b = sqrt(ONE - v0*v0)
uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
uvw(2) = mu*v0 - a*b*cosphi
uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
end if
end function rotate_angle
!===============================================================================
! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary
! neutron from a collision or for a prompt/delayed fission neutron

5
src/random_lcg.F90
View file

@ -7,6 +7,9 @@ module random_lcg
private
save
! Random number seed
integer(8), public :: seed = 1_8
integer(8) :: prn_seed0 ! original seed
integer(8) :: prn_seed(N_STREAMS) ! current seed
integer(8) :: prn_mult ! multiplication factor, g
@ -56,8 +59,6 @@ contains
subroutine initialize_prng()
use global, only: seed
integer :: i
prn_seed0 = seed

9
src/relaxng/settings.rnc
View file

@ -69,27 +69,24 @@ element settings {
element space {
(element type { xsd:string { maxLength = "16" } } |
attribute type { xsd:string { maxLength = "16" } }) &
(element length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { 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 length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { 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 length { xsd:int } | attribute length { xsd:int })? &
(element interpolation { xsd:string { maxLength = "10" } } |
attribute interplation { xsd:string { maxLength = "10" } })? &
attribute interpolation { xsd:string { maxLength = "10" } })? &
(element parameters { list { xsd:double+ } } |
attribute parameters { list { xsd:double+ } })?
}? &

36
src/relaxng/settings.rng
View file

@ -289,16 +289,6 @@
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
@ -306,7 +296,7 @@
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
@ -349,16 +339,6 @@
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
@ -366,7 +346,7 @@
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>
@ -409,16 +389,6 @@
</data>
</attribute>
</choice>
<optional>
<choice>
<element name="length">
<data type="int"/>
</element>
<attribute name="length">
<data type="int"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="interpolation">
@ -426,7 +396,7 @@
<param name="maxLength">10</param>
</data>
</element>
<attribute name="interplation">
<attribute name="interpolation">
<data type="string">
<param name="maxLength">10</param>
</data>

67
src/source.F90
View file

@ -2,6 +2,7 @@ module source
use bank_header, only: Bank
use constants
use distribution_header, only: Delta
use error, only: fatal_error
use geometry, only: find_cell
use geometry_header, only: BASE_UNIVERSE
@ -100,12 +101,8 @@ contains
integer :: i ! dummy loop index
real(8) :: r(3) ! sampled coordinates
real(8) :: phi ! azimuthal angle
real(8) :: mu ! cosine of polar angle
real(8) :: p_min(3) ! minimum coordinates of source
real(8) :: p_max(3) ! maximum coordinates of source
real(8) :: a ! Arbitrary parameter 'a'
real(8) :: b ! Arbitrary parameter 'b'
logical :: found ! Does the source particle exist within geometry?
type(Particle) :: p ! Temporary particle for using find_cell
integer, save :: num_resamples = 0 ! Number of resamples encountered
@ -188,58 +185,28 @@ contains
end select
! Sample angle
select case (external_source%type_angle)
case (SRC_ANGLE_ISOTROPIC)
! Sample isotropic distribution
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
site%uvw(1) = mu
site%uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site%uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
site%uvw(:) = external_source%angle%sample()
case (SRC_ANGLE_MONO)
! Monodirectional source
site%uvw = external_source%params_angle
case default
call fatal_error("No angle distribution specified for external source!")
end select
! Sample energy distribution
select case (external_source%type_energy)
case (SRC_ENERGY_MONO)
! Monoenergtic source
site%E = external_source%params_energy(1)
if (site%E >= energy_max_neutron) then
! Check for monoenergetic source above maximum neutron energy
select type (energy => external_source%energy)
type is (Delta)
if (energy%x0 >= energy_max_neutron) then
call fatal_error("Source energy above range of energies of at least &
&one cross section table")
end if
case (SRC_ENERGY_MAXWELL)
a = external_source%params_energy(1)
do
! Sample Maxwellian fission spectrum
site%E = maxwell_spectrum(a)
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
end do
case (SRC_ENERGY_WATT)
a = external_source%params_energy(1)
b = external_source%params_energy(2)
do
! Sample Watt fission spectrum
site%E = watt_spectrum(a, b)
! resample if energy is greater than maximum neutron energy
if (site%E < energy_max_neutron) exit
end do
case default
call fatal_error("No energy distribution specified for external source!")
end select
do
! Sample energy spectrum
site%E = external_source%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
! Set the random number generator back to the tracking stream.
call prn_set_stream(STREAM_TRACKING)

8
src/source_header.F90
View file

@ -1,5 +1,7 @@
module source_header
use distribution_header, only: Distribution, AngleDistribution
implicit none
!===============================================================================
@ -9,11 +11,9 @@ module source_header
type ExtSource
integer :: type_space ! spacial distribution, e.g. 'box' or 'point'
integer :: type_angle ! angle distribution, e.g. 'isotropic'
integer :: type_energy ! energy distribution, e.g. 'Watt'
real(8), allocatable :: params_space(:) ! parameters for spatial distribution
real(8), allocatable :: params_angle(:) ! parameters for angle distribution
real(8), allocatable :: params_energy(:) ! parameters for energy distribution
type(AngleDistribution) :: angle ! angle distribution
class(Distribution), allocatable :: energy ! energy distribution
end type ExtSource
end module source_header

1
src/state_point.F90
View file

@ -22,6 +22,7 @@ module state_point
use tally_header, only: TallyObject
use mesh_header, only: RegularMesh
use dict_header, only: ElemKeyValueII, ElemKeyValueCI
use random_lcg, only: seed
#ifdef MPI
use message_passing