crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Respond to @smharper comments on #556

Browse source
This commit is contained in:
Paul Romano 2016-01-14 06:46:32 -06:00
parent 84cf96a4be
commit ea9fb637f6
8 changed files with 68 additions and 44 deletions
Download patch file Download diff file Expand all files Collapse all files

2
openmc/source.py
View file

@ -95,6 +95,8 @@ class Source(object):
@energy.setter
def energy(self, energy):
cv.check_type('energy distribution', energy, Univariate)
if energy.name is None:
energy.name = 'energy'
self._energy = energy
@strength.setter

21
openmc/stats/multivariate.py
View file

@ -80,7 +80,7 @@ class PolarAzimuthal(UnitSphere):
mu : openmc.stats.Univariate
Distribution of the cosine of the polar angle
phi : openmc.stats.Univariate
Distribution of the azimuthal angle
Distribution of the azimuthal angle in radians
name : str, optional
Name of the distribution. Defaults to 'angle'.
reference_uvw : Iterable of Real
@ -92,21 +92,22 @@ class PolarAzimuthal(UnitSphere):
mu : openmc.stats.Univariate
Distribution of the cosine of the polar angle
phi : openmc.stats.Univariate
Distribution of the azimuthal angle
Distribution of the azimuthal angle in radians
"""
def __init__(self, mu=None, phi=None, name='angle',reference_uvw=[0., 0., 1.]):
def __init__(self, mu=None, phi=None, name='angle',
reference_uvw=[0., 0., 1.]):
super(PolarAzimuthal, self).__init__(name, reference_uvw)
if mu is not None:
self.mu = mu
else:
self.mu = Uniform('mu', -1., 1.)
self.mu = Uniform(-1., 1.)
if phi is not None:
self.phi = phi
else:
self.phi = Uniform('phi', 0., 2*pi)
self.phi = Uniform(0., 2*pi)
@property
def mu(self):
@ -119,11 +120,15 @@ class PolarAzimuthal(UnitSphere):
@mu.setter
def mu(self, mu):
cv.check_type('cosine of polar angle', mu, Univariate)
if mu.name is None:
mu.name = 'mu'
self._mu = mu
@phi.setter
def phi(self, phi):
cv.check_type('azimuthal angle', phi, Univariate)
if phi.name is None:
phi.name = 'phi'
self._phi = phi
def to_xml(self):
@ -271,16 +276,22 @@ class SpatialIndependent(Spatial):
@x.setter
def x(self, x):
cv.check_type('x coordinate', x, Univariate)
if x.name is None:
x.name = 'x'
self._x = x
@y.setter
def y(self, y):
cv.check_type('y coordinate', y, Univariate)
if y.name is None:
y.name = 'y'
self._y = y
@z.setter
def z(self, z):
cv.check_type('z coordinate', z, Univariate)
if z.name is None:
z.name = 'z'
self._z = z
def to_xml(self):

33
openmc/stats/univariate.py
View file

@ -30,8 +30,10 @@ class Univariate(object):
__metaclass__ = ABCMeta
def __init__(self, name):
self.name = name
def __init__(self, name=None):
self._name = None
if name is not None:
self.name = name
@property
def name(self):
@ -70,7 +72,7 @@ class Discrete(Univariate):
"""
def __init__(self, name, x, p):
def __init__(self, x, p, name=None):
super(Discrete, self).__init__(name)
self.x = x
self.p = p
@ -85,18 +87,25 @@ class Discrete(Univariate):
@x.setter
def x(self, x):
if cv._isinstance(x, Real):
x = [x]
cv.check_type('discrete values', x, Iterable, Real)
self._x = x
@p.setter
def p(self, p):
if cv._isinstance(p, Real):
p = [p]
cv.check_type('discrete probabilities', p, Iterable, Real)
for pk in p:
cv.check_greater_than('discrete probability', pk, 0.0, True)
self._p = p
def to_xml(self):
element = ET.Element(self.name)
if self.name is not None:
element = ET.Element(self.name)
else:
element = ET.Element('distribution')
element.set("type", "discrete")
params = ET.SubElement(element, "parameters")
@ -124,7 +133,7 @@ class Uniform(Univariate):
"""
def __init__(self, name, a=0.0, b=1.0):
def __init__(self, a=0.0, b=1.0, name=None):
super(Uniform, self).__init__(name)
self.a = a
self.b = b
@ -148,7 +157,10 @@ class Uniform(Univariate):
self._b = b
def to_xml(self):
element = ET.Element(self.name)
if self.name is not None:
element = ET.Element(self.name)
else:
element = ET.Element('distribution')
element.set("type", "uniform")
element.set("parameters", '{} {}'.format(self.a, self.b))
return element
@ -219,7 +231,7 @@ class Watt(Univariate):
"""
def __init__(self, a, b, name='energy'):
def __init__(self, a=0.988, b=2.249, name='energy'):
super(Watt, self).__init__(name)
self.a = a
self.b = b
@ -282,7 +294,7 @@ class Tabular(Univariate):
"""
def __init__(self, name, x, p, interpolation='linear-linear'):
def __init__(self, x, p, interpolation='linear-linear', name=None):
super(Tabular, self).__init__(name)
self.x = x
self.p = p
@ -319,7 +331,10 @@ class Tabular(Univariate):
self._interpolation = interpolation
def to_xml(self):
element = ET.Element(self.name)
if self.name is not None:
element = ET.Element(self.name)
else:
element = ET.Element('distribution')
element.set("type", "tabular")
element.set("interpolation", self.interpolation)

24
src/distribution_multivariate.F90
View file

@ -5,10 +5,12 @@ module distribution_multivariate
use math, only: rotate_angle
use random_lcg, only: prn
implicit none
!===============================================================================
! UNITSPHEREDISTRIBUTION type defines a probability density function for points
! on the unit sphere. Extensions of this type are used to sample angular
! distributions for starting soures
! distributions for starting sources
!===============================================================================
type, abstract :: UnitSphereDistribution
@ -99,13 +101,13 @@ contains
real(8) :: phi ! azimuthal angle
! Sample cosine of polar angle
mu = this%mu%sample()
mu = this % mu % sample()
if (mu == ONE) then
uvw(:) = this%reference_uvw
uvw(:) = this % reference_uvw
else
! Sample azimuthal angle
phi = this%phi%sample()
uvw(:) = rotate_angle(this%reference_uvw, mu, phi)
phi = this % phi % sample()
uvw(:) = rotate_angle(this % reference_uvw, mu, phi)
end if
end function polar_azimuthal_sample
@ -127,16 +129,16 @@ contains
class(Monodirectional), intent(in) :: this
real(8) :: uvw(3)
uvw(:) = this%reference_uvw
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()
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)
@ -147,14 +149,14 @@ contains
real(8) :: r(3)
r = [ (prn(), i = 1,3) ]
xyz(:) = this%lower_left + r*(this%upper_right - this%lower_left)
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
xyz(:) = this % xyz
end function spatial_point_sample
end module distribution_multivariate

10
src/distribution_univariate.F90
View file

@ -8,6 +8,8 @@ module distribution_univariate
use string, only: to_lower
use xml_interface
implicit none
!===============================================================================
! DISTRIBUTION type defines a probability density function
!===============================================================================
@ -82,6 +84,7 @@ contains
class(Discrete), intent(in) :: this
real(8) :: x
integer :: i ! loop counter
integer :: n ! size of distribution
real(8) :: c ! cumulative frequency
real(8) :: xi ! sampled CDF value
@ -197,10 +200,10 @@ contains
real(8), intent(in) :: p(:)
integer, intent(in) :: interp
integer :: i
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.')
@ -241,6 +244,7 @@ contains
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: temp_str
integer :: n
integer :: temp_int
real(8), allocatable :: temp_real(:)
@ -258,14 +262,14 @@ contains
! Allocate extension of Distribution
select case (to_lower(type))
case ('uniform')
allocate (Uniform :: dist)
allocate(Uniform :: dist)
if (n /= 2) then
call fatal_error('Uniform distribution must have two &
&parameters specified.')
end if
case ('maxwell')
allocate (Maxwell :: dist)
allocate(Maxwell :: dist)
if (n /= 1) then
call fatal_error('Maxwell energy distribution must have one &
&parameter specified.')

8
src/input_xml.F90
View file

@ -56,7 +56,6 @@ contains
character(MAX_LINE_LEN) :: temp_str
integer :: i
integer :: n
integer :: coeffs_reqd
integer :: temp_int
integer :: temp_int_array3(3)
integer, allocatable :: temp_int_array(:)
@ -502,13 +501,6 @@ contains
! Get pointer to angular distribution
call get_node_ptr(node_source, "angle", node_angle)
! Determine number of parameters specified
if (check_for_node(node_angle, "parameters")) then
n = get_arraysize_double(node_angle, "parameters")
else
n = 0
end if
! Check for type of angular distribution
type = ''
if (check_for_node(node_angle, "type")) &

2
src/source.F90
View file

@ -104,8 +104,6 @@ contains
integer :: n_source ! number of source distributions
real(8) :: c ! cumulative frequency
real(8) :: r(3) ! sampled coordinates
real(8) :: p_min(3) ! minimum coordinates of source
real(8) :: p_max(3) ! maximum coordinates of source
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

12
tests/test_source/test_source.py
View file

@ -36,15 +36,15 @@ class SourceTestHarness(PyAPITestHarness):
geometry_xml.export_to_xml()
# Create an array of different sources
x_dist = openmc.stats.Uniform('x', -3., 3.)
y_dist = openmc.stats.Discrete('y', [-4., -1., 3.], [0.2, 0.3, 0.5])
z_dist = openmc.stats.Tabular('z', [-2., 0., 2.], [0.2, 0.3, 0.2])
x_dist = openmc.stats.Uniform(-3., 3.)
y_dist = openmc.stats.Discrete([-4., -1., 3.], [0.2, 0.3, 0.5])
z_dist = openmc.stats.Tabular([-2., 0., 2.], [0.2, 0.3, 0.2])
spatial1 = openmc.stats.SpatialIndependent(x_dist, y_dist, z_dist)
spatial2 = openmc.stats.SpatialBox([-4., -4., -4.], [4., 4., 4.])
spatial3 = openmc.stats.SpatialPoint([1.2, -2.3, 0.781])
mu_dist = openmc.stats.Discrete('mu', [-1., 0., 1.], [0.5, 0.25, 0.25])
phi_dist = openmc.stats.Uniform('phi', 0., 6.28318530718)
mu_dist = openmc.stats.Discrete([-1., 0., 1.], [0.5, 0.25, 0.25])
phi_dist = openmc.stats.Uniform(0., 6.28318530718)
angle1 = openmc.stats.PolarAzimuthal(mu_dist, phi_dist)
angle2 = openmc.stats.Monodirectional(reference_uvw=[0., 1., 0.])
angle3 = openmc.stats.Isotropic()
@ -54,7 +54,7 @@ class SourceTestHarness(PyAPITestHarness):
p /= sum(np.diff(E)*p[:-1])
energy1 = openmc.stats.Maxwell(1.2895)
energy2 = openmc.stats.Watt(0.988, 2.249)
energy3 = openmc.stats.Tabular('energy', E, p, interpolation='histogram')
energy3 = openmc.stats.Tabular(E, p, interpolation='histogram')
source1 = Source(spatial1, angle1, energy1, strength=0.5)
source2 = Source(spatial2, angle2, energy2, strength=0.3)