diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index fb55f9e51..c14767e46 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -423,12 +423,16 @@ attributes/sub-elements: :type: The type of spatial distribution. Valid options are "box", "fission", - "point", and "cartesian". A "box" spatial distribution has coordinates - sampled uniformly in a parallelepiped. A "fission" spatial distribution - samples locations from a "box" distribution but only locations in - fissionable materials are accepted. A "point" spatial distribution has - coordinates specified by a triplet. An "cartesian" spatial distribution - specifies independent distributions of x-, y-, and z-coordinates. + "point", "cartesian", and "spherical". A "box" spatial distribution has + coordinates sampled uniformly in a parallelepiped. A "fission" spatial + distribution samples locations from a "box" distribution but only + locations in fissionable materials are accepted. A "point" spatial + distribution has coordinates specified by a triplet. An "cartesian" + spatial distribution specifies independent distributions of x-, y-, and + z-coordinates. 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). *Default*: None @@ -446,6 +450,9 @@ attributes/sub-elements: For an "cartesian" distribution, no parameters are specified. Instead, the ``x``, ``y``, and ``z`` elements must be specified. + For a "spherical" distribution, no parameters are specified. Instead, + the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified. + *Default*: None :x: @@ -466,6 +473,28 @@ attributes/sub-elements: univariate probability distribution (see the description in :ref:`univariate`). + :r: + For a "spherical" distribution, this element specifies the distribution + of r-coordinates. The necessary sub-elements/attributes are those of a + univariate probability distribution (see the description in + :ref:`univariate`). + + :theta: + For a "spherical" distribution, this element specifies the distribution + of theta-coordinates. The necessary sub-elements/attributes are those of a + univariate probability distribution (see the description in + :ref:`univariate`). + + :phi: + For a "spherical" distribution, this element specifies the distribution + of phi-coordinates. The necessary sub-elements/attributes are those of a + univariate probability distribution (see the description in + :ref:`univariate`). + + :origin: + For a "spherical" distribution, this element specifies the coordinates of + the center of the sphere. + :angle: An element specifying the angular distribution of source sites. This element has the following attributes: diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index 6ad0e6b03..cacbab31b 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -46,5 +46,6 @@ Spatial Distributions openmc.stats.Spatial openmc.stats.CartesianIndependent + openmc.stats.SphericalIndependent openmc.stats.Box openmc.stats.Point diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 248f03588..6cf6f52d6 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -37,6 +37,24 @@ private: UPtrDist z_; //!< Distribution of z coordinates }; +//============================================================================== +//! Distribution of points specified by spherical coordinates r,theta,phi +//============================================================================== + +class SphericalIndependent : public SpatialDistribution { +public: + explicit SphericalIndependent(pugi::xml_node node); + + //! Sample a position from the distribution + //! \return Sampled position + Position sample() const; +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 +}; + //============================================================================== //! Uniform distribution of points over a box //============================================================================== diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 35afc21dd..6a97e6192 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -272,6 +272,8 @@ class Spatial(metaclass=ABCMeta): distribution = get_text(elem, 'type') if distribution == 'cartesian': return CartesianIndependent.from_xml_element(elem) + elif distribution == 'spherical': + return SphericalIndependent.from_xml_element(elem) elif distribution == 'box' or distribution == 'fission': return Box.from_xml_element(elem) elif distribution == 'point': @@ -281,7 +283,7 @@ class Spatial(metaclass=ABCMeta): class CartesianIndependent(Spatial): """Spatial distribution with independent x, y, and z distributions. - This distribution allows one to specify a coordinates whose x-, y-, and z- + This distribution allows one to specify coordinates whose x-, y-, and z- components are sampled independently from one another. Parameters @@ -304,7 +306,6 @@ class CartesianIndependent(Spatial): """ - def __init__(self, x, y, z): super().__init__() self.x = x @@ -375,6 +376,122 @@ class CartesianIndependent(Spatial): return cls(x, y, z) +class SphericalIndependent(Spatial): + """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). + + Parameters + ---------- + r : openmc.stats.Univariate + Distribution of r-coordinates + theta : openmc.stats.Univariate + Distribution of theta-coordinates (angle relative to the z-axis) + phi : openmc.stats.Univariate + Distribution of phi-coordinates (azimuthal angle) + origin: Iterable of float, optional + coordinates (x0, y0, z0) of the center of the sphere. Defaults to + (0.0, 0.0, 0.0) + + Attributes + ---------- + r : openmc.stats.Univariate + Distribution of r-coordinates + theta : openmc.stats.Univariate + Distribution of theta-coordinates (angle relative to the z-axis) + phi : openmc.stats.Univariate + Distribution of phi-coordinates (azimuthal angle) + origin: Iterable of float, optional + coordinates (x0, y0, z0) of the center of the sphere. Defaults to + (0.0, 0.0, 0.0) + + """ + + def __init__(self, r, theta, phi, origin=(0.0, 0.0, 0.0)): + super().__init__() + self.r = r + self.theta = theta + self.phi = phi + self.origin = origin + + @property + def r(self): + return self._r + + @property + def theta(self): + return self._theta + + @property + def phi(self): + return self._phi + + @property + def origin(self): + return self._origin + + @r.setter + def r(self, r): + 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 + + @phi.setter + def phi(self, phi): + cv.check_type('phi coordinate', phi, Univariate) + self._phi = phi + + @origin.setter + def origin(self, origin): + cv.check_type('origin coordinates', origin, Iterable, Real) + origin = np.asarray(origin) + self._origin = origin + + def to_xml_element(self): + """Return XML representation of the spatial distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing spatial distribution data + + """ + 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.phi.to_xml_element('phi')) + element.set("origin", ' '.join(map(str, self.origin))) + return element + + @classmethod + def from_xml_element(cls, elem): + """Generate spatial distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.SphericalIndependent + Spatial distribution generated from XML element + + """ + r = Univariate.from_xml_element(elem.find('r')) + theta = Univariate.from_xml_element(elem.find('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) + + class Box(Spatial): """Uniform distribution of coordinates in a rectangular cuboid. diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 15098407b..e92601016 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -51,6 +51,73 @@ Position CartesianIndependent::sample() const return {x_->sample(), y_->sample(), z_->sample()}; } +//============================================================================== +// SphericalIndependent implementation +//============================================================================== + +SphericalIndependent::SphericalIndependent(pugi::xml_node node) +{ + // Read distribution for r-coordinate + if (check_for_node(node, "r")) { + pugi::xml_node node_dist = node.child("r"); + r_ = distribution_from_xml(node_dist); + } else { + // If no distribution was specified, default to a single point at r=0 + double x[] {0.0}; + double p[] {1.0}; + r_ = std::make_unique(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); + } else { + // If no distribution was specified, default to a single point at theta=0 + double x[] {0.0}; + double p[] {1.0}; + theta_ = std::make_unique(x, p, 1); + } + + // Read distribution for phi-coordinate + if (check_for_node(node, "phi")) { + pugi::xml_node node_dist = node.child("phi"); + phi_ = distribution_from_xml(node_dist); + } else { + // If no distribution was specified, default to a single point at phi=0 + double x[] {0.0}; + double p[] {1.0}; + phi_ = std::make_unique(x, p, 1); + } + + // Read sphere center coordinates + if (check_for_node(node, "origin")) { + auto origin = get_node_array(node, "origin"); + if (origin.size() == 3) { + origin_ = origin; + } else { + std::stringstream err_msg; + err_msg << "Origin for spherical source distribution must be length 3"; + fatal_error(err_msg); + } + } else { + // If no coordinates were specified, default to (0, 0, 0) + origin_ = {0.0, 0.0, 0.0}; + } + +} + +Position SphericalIndependent::sample() const +{ + double r = r_->sample(); + double theta = theta_->sample(); + double phi = phi_->sample(); + 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; + return {x, y, z}; +} + //============================================================================== // SpatialBox implementation //============================================================================== diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 731aa70a9..eb27d8591 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -75,7 +75,11 @@ element settings { attribute parameters { list { xsd:double+ } })? & element x { distribution }? & element y { distribution }? & - element z { distribution }? + element z { distribution }? & + element r { distribution }? & + element theta { distribution }? & + element phi { distribution }? & + element origin { list { xsd:double, xsd:double, xsd:double } }? }? & element angle { (element type { xsd:string } | attribute type { xsd:string }) & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 3ed15c3ed..0c825d0d8 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -348,6 +348,30 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/source.cpp b/src/source.cpp index 6b5cb8aba..b1156dc08 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -86,6 +86,8 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) type = get_node_value(node_space, "type", true, true); if (type == "cartesian") { space_ = UPtrSpace{new CartesianIndependent(node_space)}; + } else if (type == "spherical") { + space_ = UPtrSpace{new SphericalIndependent(node_space)}; } else if (type == "box") { space_ = UPtrSpace{new SpatialBox(node_space)}; } else if (type == "fission") { diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index 7eeefbc00..6f05eb3a7 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -41,7 +41,7 @@ - + 1.2 -2.3 0.781 @@ -50,4 +50,17 @@ 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + + 0.7853981633974483 1.5707963267948966 2.356194490192345 0.3 0.4 0.3 + + + + + + 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 6fe76c8cd..8f24644f0 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -2.800827E-01 7.360163E-03 +3.004769E-01 3.944044E-03 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index cd468c5b8..5fd06286b 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -29,9 +29,16 @@ class SourceTestHarness(PyAPITestHarness): 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]) + r_dist = openmc.stats.Uniform(2., 3.) + theta_dist = openmc.stats.Discrete([pi/4, pi/2, 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, + phi_dist, + origin=(1.0, 1.0, 0.0)) mu_dist = openmc.stats.Discrete([-1., 0., 1.], [0.5, 0.25, 0.25]) phi_dist = openmc.stats.Uniform(0., 6.28318530718) @@ -48,13 +55,14 @@ class SourceTestHarness(PyAPITestHarness): source1 = openmc.Source(spatial1, angle1, energy1, strength=0.5) source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3) - source3 = openmc.Source(spatial3, angle3, energy3, strength=0.2) + source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1) + source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1) settings = openmc.Settings() settings.batches = 10 settings.inactive = 5 settings.particles = 1000 - settings.source = [source1, source2, source3] + settings.source = [source1, source2, source3, source4] settings.export_to_xml()