From ea9fb637f63f9374c7436456141afa850b84acf9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Jan 2016 06:46:32 -0600 Subject: [PATCH] Respond to @smharper comments on #556 --- openmc/source.py | 2 ++ openmc/stats/multivariate.py | 21 +++++++++++++++----- openmc/stats/univariate.py | 33 ++++++++++++++++++++++--------- src/distribution_multivariate.F90 | 24 +++++++++++----------- src/distribution_univariate.F90 | 10 +++++++--- src/input_xml.F90 | 8 -------- src/source.F90 | 2 -- tests/test_source/test_source.py | 12 +++++------ 8 files changed, 68 insertions(+), 44 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index cac7d1e58c..5e4d9e53df 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -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 diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 52dc6b96ab..ea02d7f3fd 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -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): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index a0395e3431..01419420ea 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -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) diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index 970311d7b8..71aee646ab 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -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 diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 index 6d85038212..b3d77e2f82 100644 --- a/src/distribution_univariate.F90 +++ b/src/distribution_univariate.F90 @@ -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 & ¶meters specified.') end if case ('maxwell') - allocate (Maxwell :: dist) + allocate(Maxwell :: dist) if (n /= 1) then call fatal_error('Maxwell energy distribution must have one & ¶meter specified.') diff --git a/src/input_xml.F90 b/src/input_xml.F90 index d8224742b0..fa0ca10345 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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")) & diff --git a/src/source.F90 b/src/source.F90 index 3700884868..406a4bf5ad 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -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 diff --git a/tests/test_source/test_source.py b/tests/test_source/test_source.py index 9c2bfe737f..e2451ba3fc 100644 --- a/tests/test_source/test_source.py +++ b/tests/test_source/test_source.py @@ -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)