Merge pull request #2061 from cfichtlscherer/theta_distribution_spherical_independent

Uniform distribution of cos(theta) in SphericalIndependent
This commit is contained in:
Paul Romano 2022-05-27 06:46:00 -05:00 committed by GitHub
commit df50db22e1
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
24 changed files with 160 additions and 104 deletions

View file

@ -490,9 +490,10 @@ attributes/sub-elements:
independent distributions of r-, phi-, and z-coordinates where phi is the
azimuthal angle and the origin for the cylindrical coordinate system is
specified by origin. A "spherical" spatial distribution specifies
independent distributions of r-, theta-, and phi-coordinates where theta
is the angle with respect to the z-axis, phi is the azimuthal angle, and
the sphere is centered on the coordinate (x0,y0,z0).
independent distributions of r-, cos_theta-, and phi-coordinates where
cos_theta is the cosine of the angle with respect to the z-axis, phi is
the azimuthal angle, and the sphere is centered on the coordinate
(x0,y0,z0).
*Default*: None

View file

@ -70,7 +70,7 @@ private:
};
//==============================================================================
//! Distribution of points specified by spherical coordinates r,theta,phi
//! Distribution of points specified by spherical coordinates r,cos_theta,phi
//==============================================================================
class SphericalIndependent : public SpatialDistribution {
@ -83,15 +83,15 @@ public:
Position sample(uint64_t* seed) const;
Distribution* r() const { return r_.get(); }
Distribution* theta() const { return theta_.get(); }
Distribution* cos_theta() const { return cos_theta_.get(); }
Distribution* phi() const { return phi_.get(); }
Position origin() const { return origin_; }
private:
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist theta_; //!< Distribution of theta coordinates
UPtrDist phi_; //!< Distribution of phi coordinates
Position origin_; //!< Cartesian coordinates of the sphere center
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist cos_theta_; //!< Distribution of cos_theta coordinates
UPtrDist phi_; //!< Distribution of phi coordinates
Position origin_; //!< Cartesian coordinates of the sphere center
};
//==============================================================================

View file

@ -1,6 +1,6 @@
from abc import ABC, abstractmethod
from collections.abc import Iterable
from math import pi
from math import pi, cos
from numbers import Real
from xml.etree import ElementTree as ET
@ -8,7 +8,7 @@ import numpy as np
import openmc.checkvalue as cv
from .._xml import get_text
from .univariate import Univariate, Uniform
from .univariate import Univariate, Uniform, PowerLaw
class UnitSphere(ABC):
@ -208,7 +208,6 @@ class Monodirectional(UnitSphere):
"""
def __init__(self, reference_uvw=[1., 0., 0.]):
super().__init__(reference_uvw)
@ -375,8 +374,8 @@ class SphericalIndependent(Spatial):
r"""Spatial distribution represented in spherical coordinates.
This distribution allows one to specify coordinates whose :math:`r`,
:math:`\theta`, and :math:`\phi` components are sampled independently from
one another and centered on the coordinates (x0, y0, z0).
:math:`\theta`, and :math:`\phi` components are sampled independently
from one another and centered on the coordinates (x0, y0, z0).
.. versionadded: 0.12
@ -385,9 +384,9 @@ class SphericalIndependent(Spatial):
r : openmc.stats.Univariate
Distribution of r-coordinates in a reference frame specified by
the origin parameter
theta : openmc.stats.Univariate
Distribution of theta-coordinates (angle relative to the z-axis) in a
reference frame specified by the origin parameter
cos_theta : openmc.stats.Univariate
Distribution of the cosine of the theta-coordinates (angle relative to
the z-axis) in a reference frame specified by the origin parameter
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle) in a reference frame
specified by the origin parameter
@ -399,9 +398,9 @@ class SphericalIndependent(Spatial):
----------
r : openmc.stats.Univariate
Distribution of r-coordinates in the local reference frame
theta : openmc.stats.Univariate
Distribution of theta-coordinates (angle relative to the z-axis) in the
local reference frame
cos_theta : openmc.stats.Univariate
Distribution of the cosine of the theta-coordinates (angle relative to
the z-axis) in the local reference frame
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle) in the local
reference frame
@ -411,9 +410,9 @@ class SphericalIndependent(Spatial):
"""
def __init__(self, r, theta, phi, origin=(0.0, 0.0, 0.0)):
def __init__(self, r, cos_theta, phi, origin=(0.0, 0.0, 0.0)):
self.r = r
self.theta = theta
self.cos_theta = cos_theta
self.phi = phi
self.origin = origin
@ -422,8 +421,8 @@ class SphericalIndependent(Spatial):
return self._r
@property
def theta(self):
return self._theta
def cos_theta(self):
return self._cos_theta
@property
def phi(self):
@ -438,10 +437,10 @@ class SphericalIndependent(Spatial):
cv.check_type('r coordinate', r, Univariate)
self._r = r
@theta.setter
def theta(self, theta):
cv.check_type('theta coordinate', theta, Univariate)
self._theta = theta
@cos_theta.setter
def cos_theta(self, cos_theta):
cv.check_type('cos_theta coordinate', cos_theta, Univariate)
self._cos_theta = cos_theta
@phi.setter
def phi(self, phi):
@ -466,7 +465,7 @@ class SphericalIndependent(Spatial):
element = ET.Element('space')
element.set('type', 'spherical')
element.append(self.r.to_xml_element('r'))
element.append(self.theta.to_xml_element('theta'))
element.append(self.cos_theta.to_xml_element('cos_theta'))
element.append(self.phi.to_xml_element('phi'))
element.set("origin", ' '.join(map(str, self.origin)))
return element
@ -487,10 +486,10 @@ class SphericalIndependent(Spatial):
"""
r = Univariate.from_xml_element(elem.find('r'))
theta = Univariate.from_xml_element(elem.find('theta'))
cos_theta = Univariate.from_xml_element(elem.find('cos_theta'))
phi = Univariate.from_xml_element(elem.find('phi'))
origin = [float(x) for x in elem.get('origin').split()]
return cls(r, theta, phi, origin=origin)
return cls(r, cos_theta, phi, origin=origin)
class CylindricalIndependent(Spatial):
@ -640,7 +639,6 @@ class Box(Spatial):
"""
def __init__(self, lower_left, upper_right, only_fissionable=False):
self.lower_left = lower_left
self.upper_right = upper_right
@ -779,3 +777,42 @@ class Point(Spatial):
"""
xyz = [float(x) for x in get_text(elem, 'parameters').split()]
return cls(xyz)
def spherical_uniform(r_outer, r_inner=0.0, thetas=(0., pi), phis=(0., 2*pi),
origin=(0., 0., 0.)):
"""Return a uniform spatial distribution over a spherical shell.
This function provides a uniform spatial distribution over a spherical
shell between `r_inner` and `r_outer`. Optionally, the range of angles
can be restricted by the `thetas` and `phis` arguments.
.. versionadded: 0.13.1
Parameters
----------
r_outer : float
Outer radius of the spherical shell in [cm]
r_inner : float, optional
Inner radius of the spherical shell in [cm]
thetas : iterable of float, optional
Starting and ending theta coordinates (angle relative to
the z-axis) in radius in a reference frame centered at `origin`
phis : iterable of float, optional
Starting and ending phi coordinates (azimuthal angle) in
radians in a reference frame centered at `origin`
origin: iterable of float, optional
Coordinates (x0, y0, z0) of the center of the spherical
reference frame for the distribution.
Returns
-------
openmc.stats.SphericalIndependent
Uniform distribution over the spherical shell
"""
r_dist = PowerLaw(r_inner, r_outer, 2)
cos_thetas_dist = Uniform(cos(thetas[1]), cos(thetas[0]))
phis_dist = Uniform(phis[0], phis[1])
return SphericalIndependent(r_dist, cos_thetas_dist, phis_dist, origin)

View file

@ -132,15 +132,15 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
r_ = make_unique<Discrete>(x, p, 1);
}
// Read distribution for theta-coordinate
if (check_for_node(node, "theta")) {
pugi::xml_node node_dist = node.child("theta");
theta_ = distribution_from_xml(node_dist);
// Read distribution for cos_theta-coordinate
if (check_for_node(node, "cos_theta")) {
pugi::xml_node node_dist = node.child("cos_theta");
cos_theta_ = distribution_from_xml(node_dist);
} else {
// If no distribution was specified, default to a single point at theta=0
// If no distribution was specified, default to a single point at cos_theta=0
double x[] {0.0};
double p[] {1.0};
theta_ = make_unique<Discrete>(x, p, 1);
cos_theta_ = make_unique<Discrete>(x, p, 1);
}
// Read distribution for phi-coordinate
@ -171,11 +171,12 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
Position SphericalIndependent::sample(uint64_t* seed) const
{
double r = r_->sample(seed);
double theta = theta_->sample(seed);
double cos_theta = cos_theta_->sample(seed);
double phi = phi_->sample(seed);
double x = r * sin(theta) * cos(phi) + origin_.x;
double y = r * sin(theta) * sin(phi) + origin_.y;
double z = r * cos(theta) + origin_.z;
// sin(theta) by sin**2 + cos**2 = 1
double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;
double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
double z = r * cos_theta + origin_.z;
return {x, y, z};
}

View file

@ -53,9 +53,9 @@
<source strength="0.1">
<space origin="1.0 1.0 0.0" type="spherical">
<r parameters="2.0 3.0" type="uniform" />
<theta type="discrete">
<parameters>0.7853981633974483 1.5707963267948966 2.356194490192345 0.3 0.4 0.3</parameters>
</theta>
<cos_theta type="discrete">
<parameters>0.7071067811865476 0.0 -0.7071067811865475 0.3 0.4 0.3</parameters>
</cos_theta>
<phi parameters="0.0 6.283185307179586" type="uniform" />
</space>
<angle type="isotropic" />
@ -102,9 +102,9 @@
<source strength="0.1">
<space origin="1.0 1.0 0.0" type="spherical">
<r parameters="2.0 3.0 2.0" type="powerlaw" />
<theta type="discrete">
<parameters>0.7853981633974483 1.5707963267948966 2.356194490192345 0.3 0.4 0.3</parameters>
</theta>
<cos_theta type="discrete">
<parameters>0.7071067811865476 0.0 -0.7071067811865475 0.3 0.4 0.3</parameters>
</cos_theta>
<phi parameters="0.0 6.283185307179586" type="uniform" />
</space>
<angle type="isotropic" />

View file

@ -1,4 +1,4 @@
from math import pi
from math import pi, cos
import numpy as np
import openmc
@ -32,19 +32,19 @@ class SourceTestHarness(PyAPITestHarness):
r_dist = openmc.stats.Uniform(2., 3.)
r_dist1 = openmc.stats.PowerLaw(2., 3., 1.)
r_dist2 = openmc.stats.PowerLaw(2., 3., 2.)
theta_dist = openmc.stats.Discrete([pi/4, pi/2, 3*pi/4],
[0.3, 0.4, 0.3])
cos_theta_dist = openmc.stats.Discrete([cos(pi/4), 0.0, cos(3*pi/4)],
[0.3, 0.4, 0.3])
phi_dist = openmc.stats.Uniform(0.0, 2*pi)
spatial1 = openmc.stats.CartesianIndependent(x_dist, y_dist, z_dist)
spatial2 = openmc.stats.Box([-4., -4., -4.], [4., 4., 4.])
spatial3 = openmc.stats.Point([1.2, -2.3, 0.781])
spatial4 = openmc.stats.SphericalIndependent(r_dist, theta_dist,
spatial4 = openmc.stats.SphericalIndependent(r_dist, cos_theta_dist,
phi_dist,
origin=(1., 1., 0.))
spatial5 = openmc.stats.CylindricalIndependent(r_dist, phi_dist,
z_dist,
origin=(1., 1., 0.))
spatial6 = openmc.stats.SphericalIndependent(r_dist2, theta_dist,
spatial6 = openmc.stats.SphericalIndependent(r_dist2, cos_theta_dist,
phi_dist,
origin=(1., 1., 0.))
spatial7 = openmc.stats.CylindricalIndependent(r_dist1, phi_dist,

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -50,9 +50,9 @@
<source strength="1.0">
<space origin="0.0 0.0 0.0" type="spherical">
<r parameters="0.0 0.0" type="uniform" />
<theta type="discrete">
<parameters>0.0 1.0</parameters>
</theta>
<cos_theta type="discrete">
<parameters>1.0 1.0</parameters>
</cos_theta>
<phi type="discrete">
<parameters>0.0 1.0</parameters>
</phi>

View file

@ -228,10 +228,10 @@ def test_unstructured_mesh(test_opts):
# source setup
r = openmc.stats.Uniform(a=0.0, b=0.0)
theta = openmc.stats.Discrete(x=[0.0], p=[1.0])
cos_theta = openmc.stats.Discrete(x=[1.0], p=[1.0])
phi = openmc.stats.Discrete(x=[0.0], p=[1.0])
space = openmc.stats.SphericalIndependent(r, theta, phi)
space = openmc.stats.SphericalIndependent(r, cos_theta, phi)
energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
source = openmc.Source(space=space, energy=energy)
settings.source = source

View file

@ -1,6 +1,6 @@
import openmc
import openmc.stats
from math import pi
def test_source():
space = openmc.stats.Point()
@ -26,6 +26,22 @@ def test_source():
assert src.strength == 1.0
def test_spherical_uniform():
r_outer = 2.0
r_inner = 1.0
thetas = (0.0, pi/2)
phis = (0.0, pi)
origin = (0.0, 1.0, 2.0)
sph_indep_function = openmc.stats.spherical_uniform(r_outer,
r_inner,
thetas,
phis,
origin)
assert isinstance(sph_indep_function, openmc.stats.SphericalIndependent)
def test_source_file():
filename = 'source.h5'
src = openmc.Source(filename=filename)
@ -35,6 +51,7 @@ def test_source_file():
assert 'strength' in elem.attrib
assert 'file' in elem.attrib
def test_source_dlopen():
library = './libsource.so'
src = openmc.Source(library=library)