Skip to content

Commit 76508f8

Browse files
GuyStenGuySten
committed
added option to have cylindrical sources with general axis
1 parent b14eff4 commit 76508f8

7 files changed

Lines changed: 225 additions & 14 deletions

File tree

‎docs/source/io_formats/settings.rst‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -695,6 +695,7 @@ attributes/sub-elements:
695695

696696
For a "cylindrical" distribution, no parameters are specified. Instead,
697697
the ``r``, ``phi``, ``z``, and ``origin`` elements must be specified.
698+
Optionally, the ``r_dir`` and ``z_dir`` elements could be specified.
698699

699700
For a "spherical" distribution, no parameters are specified. Instead,
700701
the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified.
@@ -726,6 +727,10 @@ attributes/sub-elements:
726727
of a univariate probability distribution (see the description in
727728
:ref:`univariate`).
728729

730+
:r_dir:
731+
For "cylindrical" distributions, this element specifies the direction
732+
of the cylinder r-axis at phi=0, it defaults to (1.0, 0.0, 0.0).
733+
729734
:theta:
730735
For a "spherical" distribution, this element specifies the distribution
731736
of theta-coordinates. The necessary sub-elements/attributes are those of a
@@ -737,6 +742,10 @@ attributes/sub-elements:
737742
the distribution of phi-coordinates. The necessary
738743
sub-elements/attributes are those of a univariate probability
739744
distribution (see the description in :ref:`univariate`).
745+
746+
:phi_dir:
747+
For "cylindrical" distributions, this element specifies the direction
748+
of the cylinder phi-axis at phi=0, it defaults to (0.0, 1.0, 0.0).
740749

741750
:origin:
742751
For "cylindrical and "spherical" distributions, this element specifies

‎docs/source/pythonapi/stats.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,5 @@ Spatial Distributions
6767
:template: myfunction.rst
6868

6969
openmc.stats.spherical_uniform
70+
openmc.stats.cylindrical_uniform
71+
openmc.stats.ring_uniform

‎include/openmc/distribution_spatial.h‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -65,12 +65,18 @@ class CylindricalIndependent : public SpatialDistribution {
6565
Distribution* phi() const { return phi_.get(); }
6666
Distribution* z() const { return z_.get(); }
6767
Position origin() const { return origin_; }
68+
Direction r_dir() const { return r_dir_; }
69+
Direction phi_dir() const { return phi_dir_; }
70+
Direction z_dir() const { return z_dir_; }
6871

6972
private:
70-
UPtrDist r_; //!< Distribution of r coordinates
71-
UPtrDist phi_; //!< Distribution of phi coordinates
72-
UPtrDist z_; //!< Distribution of z coordinates
73-
Position origin_; //!< Cartesian coordinates of the cylinder center
73+
UPtrDist r_; //!< Distribution of r coordinates
74+
UPtrDist phi_; //!< Distribution of phi coordinates
75+
UPtrDist z_; //!< Distribution of z coordinates
76+
Position origin_; //!< Cartesian coordinates of the cylinder center
77+
Direction r_dir_; //!< Direction of r-axis at phi=0
78+
Direction phi_dir_; //!< Direction of phi-axis at phi=0
79+
Direction z_dir_; //!< Direction of z-axis
7480
};
7581

7682
//==============================================================================

‎openmc/stats/multivariate.py‎

Lines changed: 136 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212
import openmc.checkvalue as cv
1313
from .._xml import get_elem_list, get_text
1414
from ..mesh import MeshBase
15-
from .univariate import PowerLaw, Uniform, Univariate
15+
from .univariate import PowerLaw, Uniform, Univariate, delta_function
1616

1717

1818
class UnitSphere(ABC):
@@ -532,6 +532,10 @@ class CylindricalIndependent(Spatial):
532532
origin: Iterable of float, optional
533533
coordinates (x0, y0, z0) of the center of the cylindrical reference
534534
frame. Defaults to (0.0, 0.0, 0.0)
535+
r_dir : Iterable of float, optional
536+
Unit vector of the cylinder r axis at phi=0. Defaults to (1.0, 0.0, 0.0).
537+
z_dir : Iterable of float, optional
538+
Unit vector of the cylinder z axis direction. Defaults to (0.0, 0.0, 1.0).
535539
536540
Attributes
537541
----------
@@ -545,14 +549,20 @@ class CylindricalIndependent(Spatial):
545549
origin: Iterable of float, optional
546550
coordinates (x0, y0, z0) of the center of the cylindrical reference
547551
frame. Defaults to (0.0, 0.0, 0.0)
552+
r_dir : Iterable of float, optional
553+
Unit vector of the cylinder r axis at phi=0. Defaults to (1.0, 0.0, 0.0).
554+
z_dir : Iterable of float, optional
555+
Unit vector of the cylinder z axis direction. Defaults to (0.0, 0.0, 1.0).
548556
549557
"""
550558

551-
def __init__(self, r, phi, z, origin=(0.0, 0.0, 0.0)):
559+
def __init__(self, r, phi, z, origin=(0.0, 0.0, 0.0), r_dir=(1.0, 0.0, 0.0), z_dir=(0.0, 0.0, 1.0)):
552560
self.r = r
553561
self.phi = phi
554562
self.z = z
555563
self.origin = origin
564+
self.z_dir = z_dir
565+
self.r_dir = r_dir
556566

557567
@property
558568
def r(self):
@@ -590,6 +600,33 @@ def origin(self, origin):
590600
cv.check_type('origin coordinates', origin, Iterable, Real)
591601
origin = np.asarray(origin)
592602
self._origin = origin
603+
604+
@property
605+
def z_dir(self):
606+
return self._z_dir
607+
608+
@z_dir.setter
609+
def z_dir(self, z_dir):
610+
cv.check_type('z-axis direction', z_dir, Iterable, Real)
611+
z_dir = np.asarray(z_dir)
612+
norm2 = np.dot(z_dir,z_dir)
613+
cv.check_greater_than('z-axis direction magnitude', norm2, 0.0)
614+
z_dir /= np.sqrt(norm2)
615+
self._z_dir = z_dir
616+
617+
@property
618+
def r_dir(self):
619+
return self._r_dir
620+
621+
@r_dir.setter
622+
def r_dir(self, r_dir):
623+
cv.check_type('r-axis direction', r_dir, Iterable, Real)
624+
r_dir = np.asarray(r_dir)
625+
r_dir -= np.dot(r_dir,self.z_dir)*self.z_dir
626+
norm2 = np.dot(r_dir,r_dir)
627+
cv.check_greater_than('r-axis direction magnitude', norm2,0.0)
628+
r_dir /= np.sqrt(norm2)
629+
self._r_dir = r_dir
593630

594631
def to_xml_element(self):
595632
"""Return XML representation of the spatial distribution
@@ -606,6 +643,8 @@ def to_xml_element(self):
606643
element.append(self.phi.to_xml_element('phi'))
607644
element.append(self.z.to_xml_element('z'))
608645
element.set("origin", ' '.join(map(str, self.origin)))
646+
element.set("r_dir", ' '.join(map(str, self.r_dir)))
647+
element.set("z_dir", ' '.join(map(str, self.z_dir)))
609648
return element
610649

611650
@classmethod
@@ -627,7 +666,9 @@ def from_xml_element(cls, elem: ET.Element):
627666
phi = Univariate.from_xml_element(elem.find('phi'))
628667
z = Univariate.from_xml_element(elem.find('z'))
629668
origin = get_elem_list(elem, "origin", float)
630-
return cls(r, phi, z, origin=origin)
669+
r_dir = get_elem_list(elem, "r_dir", float) or [1.0, 0.0, 0.0]
670+
z_dir = get_elem_list(elem, "z_dir", float) or [0.0, 0.0, 1.0]
671+
return cls(r, phi, z, origin=origin, r_dir=r_dir, z_dir=z_dir)
631672

632673

633674
class MeshSpatial(Spatial):
@@ -1086,3 +1127,95 @@ def spherical_uniform(
10861127
phis_dist = Uniform(phis[0], phis[1])
10871128

10881129
return SphericalIndependent(r_dist, cos_thetas_dist, phis_dist, origin)
1130+
1131+
1132+
def cylindrical_uniform(
1133+
r_outer: float,
1134+
height: float,
1135+
r_inner: float = 0.0,
1136+
phis: Sequence[float] = (0., 2*pi),
1137+
origin: Sequence[float] = (0., 0., 0.),
1138+
r_dir: Sequence[float] = (1., 0., 0.),
1139+
z_dir: Sequence[float] = (0., 0., 1.),
1140+
):
1141+
"""Return a uniform spatial distribution over a cylindrical shell.
1142+
1143+
This function provides a uniform spatial distribution over a cylindrical
1144+
shell between `r_inner` and `r_outer`. Optionally, the range of angles
1145+
can be restricted by the `phis` arguments.
1146+
1147+
1148+
Parameters
1149+
----------
1150+
r_outer : float
1151+
Outer radius of the cylindrical shell in [cm]
1152+
height : float
1153+
Height of the cylindrical shell in [cm]
1154+
r_inner : float
1155+
Inner radius of the cylindrical shell in [cm]
1156+
phis : iterable of float
1157+
Starting and ending phi coordinates (azimuthal angle) in
1158+
radians in a reference frame centered at `origin`
1159+
origin: iterable of float
1160+
Coordinates (x0, y0, z0) of the center of the cylindrical
1161+
reference frame for the distribution. Defaults to (0.0, 0.0, 0.0)
1162+
r_dir : iterable of float
1163+
Direction of the r-axis at phi=0. Defaults to (1.0, 0.0, 0.0)
1164+
z_dir :
1165+
Direction of the z-axis. Defaults to (0.0, 0.0, 1.0)
1166+
1167+
Returns
1168+
-------
1169+
openmc.stats.CylindricalIndependent
1170+
Uniform distribution over the cylindrical shell
1171+
"""
1172+
1173+
r_dist = PowerLaw(r_inner, r_outer, 1)
1174+
phis_dist = Uniform(phis[0], phis[1])
1175+
z_dist = Uniform(-height/2, height/2)
1176+
1177+
return CylindricalIndependent(r_dist, phis_dist, z_dist, origin, r_dir, z_dir)
1178+
1179+
def ring_uniform(
1180+
r_outer: float,
1181+
r_inner: float = 0.0,
1182+
phis: Sequence[float] = (0., 2*pi),
1183+
origin: Sequence[float] = (0., 0., 0.),
1184+
r_dir: Sequence[float] = (1., 0., 0.),
1185+
z_dir: Sequence[float] = (0., 0., 1.),
1186+
):
1187+
"""Return a uniform spatial distribution over a ring.
1188+
1189+
This function provides a uniform spatial distribution over a ring
1190+
shell between `r_inner` and `r_outer`. Optionally, the range of angles
1191+
can be restricted by the `phis` arguments.
1192+
1193+
1194+
Parameters
1195+
----------
1196+
r_outer : float
1197+
Outer radius of the ring in [cm]
1198+
r_inner : float
1199+
Inner radius of the ring in [cm]
1200+
phis : iterable of float
1201+
Starting and ending phi coordinates (azimuthal angle) in
1202+
radians in a reference frame centered at `origin`
1203+
origin: iterable of float
1204+
Coordinates (x0, y0, z0) of the center of the cylindrical
1205+
reference frame for the distribution. Defaults to (0.0, 0.0, 0.0)
1206+
r_dir : iterable of float
1207+
Direction of the r-axis at phi=0. Defaults to (1.0, 0.0, 0.0)
1208+
z_dir :
1209+
Direction of the z-axis. Defaults to (0.0, 0.0, 1.0)
1210+
1211+
Returns
1212+
-------
1213+
openmc.stats.CylindricalIndependent
1214+
Uniform distribution over the ring
1215+
"""
1216+
1217+
r_dist = PowerLaw(r_inner, r_outer, 1)
1218+
phis_dist = Uniform(phis[0], phis[1])
1219+
z_dist = delta_function(np.dot(origin,z_dir))
1220+
1221+
return CylindricalIndependent(r_dist, phis_dist, z_dist, origin, r_dir, z_dir)

‎src/distribution_spatial.cpp‎

Lines changed: 37 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -137,16 +137,49 @@ CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
137137
// If no coordinates were specified, default to (0, 0, 0)
138138
origin_ = {0.0, 0.0, 0.0};
139139
}
140+
141+
// Read cylinder z_dir
142+
if (check_for_node(node, "z_dir")) {
143+
auto z_dir = get_node_array<double>(node, "z_dir");
144+
if (z_dir.size() == 3) {
145+
z_dir_ = z_dir;
146+
z_dir_ /= z_dir_.norm();
147+
} else {
148+
fatal_error("z_dir for cylindrical source distribution must be length 3");
149+
}
150+
} else {
151+
// If no z_dir was specified, default to (0, 0, 1)
152+
z_dir_ = {0.0, 0.0, 1.0};
153+
}
154+
155+
// Read cylinder r_dir
156+
if (check_for_node(node, "r_dir")) {
157+
auto r_dir = get_node_array<double>(node, "r_dir");
158+
if (r_dir.size() == 3) {
159+
r_dir_ = r_dir;
160+
r_dir_ /= r_dir_.norm();
161+
} else {
162+
fatal_error("r_dir for cylindrical source distribution must be length 3");
163+
}
164+
} else {
165+
// If no r_dir was specified, default to (1, 0, 0)
166+
r_dir_ = {1.0, 0.0, 0.0};
167+
}
168+
169+
if (r_dir_.dot(z_dir_) > 1e-12)
170+
fatal_error("r_dir must be perpendicular to z_dir");
171+
172+
auto phi_dir = z_dir_.cross(r_dir_);
173+
phi_dir /= phi_dir.norm();
174+
phi_dir_ = phi_dir;
140175
}
141176

142177
Position CylindricalIndependent::sample(uint64_t* seed) const
143178
{
144179
double r = r_->sample(seed);
145180
double phi = phi_->sample(seed);
146-
double x = r * cos(phi) + origin_.x;
147-
double y = r * sin(phi) + origin_.y;
148-
double z = z_->sample(seed) + origin_.z;
149-
return {x, y, z};
181+
double z = z_->sample(seed);
182+
return r * (cos(phi) * r_dir_ + sin(phi) * phi_dir_) + z * z_dir_ + origin_;
150183
}
151184

152185
//==============================================================================

‎tests/regression_tests/source/inputs_true.dat‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@
6363
</energy>
6464
</source>
6565
<source type="independent" strength="0.1" particle="neutron">
66-
<space type="cylindrical" origin="1.0 1.0 0.0">
66+
<space type="cylindrical" origin="1.0 1.0 0.0" r_dir="1.0 0.0 0.0" z_dir="0.0 0.0 1.0">
6767
<r type="uniform" parameters="2.0 3.0"/>
6868
<phi type="uniform" parameters="0.0 6.283185307179586"/>
6969
<z type="tabular" interpolation="linear-linear">
@@ -76,7 +76,7 @@
7676
</energy>
7777
</source>
7878
<source type="independent" strength="0.1" particle="neutron">
79-
<space type="cylindrical" origin="1.0 1.0 0.0">
79+
<space type="cylindrical" origin="1.0 1.0 0.0" r_dir="1.0 0.0 0.0" z_dir="0.0 0.0 1.0">
8080
<r type="uniform" parameters="2.0 3.0"/>
8181
<phi type="uniform" parameters="0.0 6.283185307179586"/>
8282
<z type="tabular" interpolation="linear-linear">
@@ -123,7 +123,7 @@
123123
<time type="uniform" parameters="2 5"/>
124124
</source>
125125
<source type="independent" strength="0.1" particle="neutron">
126-
<space type="cylindrical" origin="1.0 1.0 0.0">
126+
<space type="cylindrical" origin="1.0 1.0 0.0" r_dir="1.0 0.0 0.0" z_dir="0.0 0.0 1.0">
127127
<r type="powerlaw" parameters="2.0 3.0 1.0"/>
128128
<phi type="uniform" parameters="0.0 6.283185307179586"/>
129129
<z type="tabular" interpolation="linear-linear">

‎tests/unit_tests/test_source.py‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,34 @@ def test_spherical_uniform():
4949
origin)
5050

5151
assert isinstance(sph_indep_function, openmc.stats.SphericalIndependent)
52+
53+
def test_cylindrical_uniform():
54+
r_outer = 2.0
55+
r_inner = 1.0
56+
height = 1.0
57+
phis = (0.0, pi)
58+
origin = (0.0, 1.0, 2.0)
59+
60+
cyl_indep_function = openmc.stats.cylindrical_uniform(r_outer,
61+
height,
62+
r_inner,
63+
phis,
64+
origin)
65+
66+
assert isinstance(cyl_indep_function, openmc.stats.CylindricalIndependent)
67+
68+
def test_ring_uniform():
69+
r_outer = 2.0
70+
r_inner = 1.0
71+
phis = (0.0, pi)
72+
origin = (0.0, 1.0, 2.0)
73+
74+
ring_indep_function = openmc.stats.ring_uniform(r_outer,
75+
r_inner,
76+
phis,
77+
origin)
78+
79+
assert isinstance(ring_indep_function, openmc.stats.CylindricalIndependent)
5280

5381
def test_point_cloud():
5482
positions = [(1, 0, 2), (0, 1, 0), (0, 0, 3), (4, 9, 2)]

0 commit comments

Comments
 (0)