mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Use radians in CylindricalMesh and SphericalMesh
This commit is contained in:
Paul Romano
parent
9df6820da1
commit
5dfc880d91
4 changed files with 39 additions and 46 deletions
|
|
@ -1,5 +1,6 @@
|
|||
from abc import ABC
|
||||
from collections.abc import Iterable
|
||||
from math import pi
|
||||
from numbers import Real, Integral
|
||||
import warnings
|
||||
from xml.etree import ElementTree as ET
|
||||
|
|
@ -705,7 +706,7 @@ class CylindricalMesh(MeshBase):
|
|||
Requirement is r >= 0.
|
||||
phi_grid : Iterable of float
|
||||
Mesh boundary points along the phi-axis.
|
||||
The default value is [0, 360], i.e. the full phi range.
|
||||
The default value is [0, 2π], i.e. the full phi range.
|
||||
z_grid : Iterable of float
|
||||
Mesh boundary points along the z-axis.
|
||||
indices : Iterable of tuple
|
||||
|
|
@ -718,7 +719,7 @@ class CylindricalMesh(MeshBase):
|
|||
super().__init__(mesh_id, name)
|
||||
|
||||
self._r_grid = None
|
||||
self._phi_grid = [0, 360]
|
||||
self._phi_grid = [0.0, 2*pi]
|
||||
self._z_grid = None
|
||||
|
||||
@property
|
||||
|
|
@ -761,12 +762,12 @@ class CylindricalMesh(MeshBase):
|
|||
@phi_grid.setter
|
||||
def phi_grid(self, grid):
|
||||
cv.check_type('mesh phi_grid', grid, Iterable, Real)
|
||||
self._phi_grid = np.array(grid)
|
||||
self._phi_grid = np.asarray(grid)
|
||||
|
||||
@z_grid.setter
|
||||
def z_grid(self, grid):
|
||||
cv.check_type('mesh z_grid', grid, Iterable, Real)
|
||||
self._z_grid = np.array(grid)
|
||||
self._z_grid = np.asarray(grid)
|
||||
|
||||
def __repr__(self):
|
||||
fmt = '{0: <16}{1}{2}\n'
|
||||
|
|
@ -796,7 +797,7 @@ class CylindricalMesh(MeshBase):
|
|||
# Read and assign mesh properties
|
||||
mesh = cls(mesh_id)
|
||||
mesh.r_grid = group['r_grid'][()]
|
||||
mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
|
||||
mesh.phi_grid = group['phi_grid'][()]
|
||||
mesh.z_grid = group['z_grid'][()]
|
||||
|
||||
return mesh
|
||||
|
|
@ -859,9 +860,9 @@ class CylindricalMesh(MeshBase):
|
|||
|
||||
"""
|
||||
|
||||
V_r = np.diff(np.array(self.r_grid)**2 / 2)
|
||||
V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0)
|
||||
V_z = np.diff(np.array(self.z_grid))
|
||||
V_r = np.diff(np.asarray(self.r_grid)**2 / 2)
|
||||
V_p = np.diff(self.phi_grid)
|
||||
V_z = np.diff(self.z_grid)
|
||||
|
||||
return np.multiply.outer(np.outer(V_r, V_p), V_z)
|
||||
|
||||
|
|
@ -891,10 +892,10 @@ class SphericalMesh(MeshBase):
|
|||
Requirement is r >= 0.
|
||||
theta_grid : Iterable of float
|
||||
Mesh boundary points along the theta-axis in degrees.
|
||||
The default value is [0, 180], i.e. the full theta range.
|
||||
The default value is [0, π], i.e. the full theta range.
|
||||
phi_grid : Iterable of float
|
||||
Mesh boundary points along the phi-axis in degrees.
|
||||
The default value is [0, 360], i.e. the full phi range.
|
||||
The default value is [0, 2π], i.e. the full phi range.
|
||||
indices : Iterable of tuple
|
||||
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
|
||||
(2, 1, 1), ...]
|
||||
|
|
@ -905,8 +906,8 @@ class SphericalMesh(MeshBase):
|
|||
super().__init__(mesh_id, name)
|
||||
|
||||
self._r_grid = None
|
||||
self._theta_grid = [0, 180]
|
||||
self._phi_grid = [0, 360]
|
||||
self._theta_grid = [0, pi]
|
||||
self._phi_grid = [0, 2*pi]
|
||||
|
||||
@property
|
||||
def dimension(self):
|
||||
|
|
@ -948,12 +949,12 @@ class SphericalMesh(MeshBase):
|
|||
@theta_grid.setter
|
||||
def theta_grid(self, grid):
|
||||
cv.check_type('mesh theta_grid', grid, Iterable, Real)
|
||||
self._theta_grid = np.array(grid)
|
||||
self._theta_grid = np.asarray(grid)
|
||||
|
||||
@phi_grid.setter
|
||||
def phi_grid(self, grid):
|
||||
cv.check_type('mesh phi_grid', grid, Iterable, Real)
|
||||
self._phi_grid = np.array(grid)
|
||||
self._phi_grid = np.asarray(grid)
|
||||
|
||||
def __repr__(self):
|
||||
fmt = '{0: <16}{1}{2}\n'
|
||||
|
|
@ -983,8 +984,8 @@ class SphericalMesh(MeshBase):
|
|||
# Read and assign mesh properties
|
||||
mesh = cls(mesh_id)
|
||||
mesh.r_grid = group['r_grid'][()]
|
||||
mesh.theta_grid = 180 / np.pi * group['theta_grid'][()]
|
||||
mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
|
||||
mesh.theta_grid = group['theta_grid'][()]
|
||||
mesh.phi_grid = group['phi_grid'][()]
|
||||
|
||||
return mesh
|
||||
|
||||
|
|
@ -1046,9 +1047,9 @@ class SphericalMesh(MeshBase):
|
|||
|
||||
"""
|
||||
|
||||
V_r = np.diff(np.array(self.r_grid)**3 / 3)
|
||||
V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0))
|
||||
V_p = np.diff(np.array(self.phi_grid) * np.pi / 180)
|
||||
V_r = np.diff(np.asarray(self.r_grid)**3 / 3)
|
||||
V_t = np.diff(-np.cos(self.theta_grid))
|
||||
V_p = np.diff(self.phi_grid)
|
||||
|
||||
return np.multiply.outer(np.outer(V_r, V_t), V_p)
|
||||
|
||||
|
|
|
|||
33
src/mesh.cpp
33
src/mesh.cpp
|
|
@ -1001,7 +1001,7 @@ double CylindricalMesh::find_r_crossing(
|
|||
double CylindricalMesh::find_phi_crossing(
|
||||
const Position& r, const Direction& u, double l, int shell) const
|
||||
{
|
||||
// Phi grid is [0, 360], thus there is no real surface to cross
|
||||
// Phi grid is [0, 2pi], thus there is no real surface to cross
|
||||
if (full_phi_ && (shape_[1] == 1))
|
||||
return INFTY;
|
||||
|
||||
|
|
@ -1104,18 +1104,14 @@ int CylindricalMesh::set_grid()
|
|||
"cylindrical meshes must start at phi >= 0.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
if (grid_[1].back() > 360) {
|
||||
if (grid_[1].back() > 2.0 * PI) {
|
||||
set_errmsg("phi-grids for "
|
||||
"cylindrical meshes must end with theta <= 360 degree.");
|
||||
"cylindrical meshes must end with theta <= 2*pi.");
|
||||
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 360.0);
|
||||
|
||||
// Transform phi-grid from degrees to radians
|
||||
std::transform(grid_[1].begin(), grid_[1].end(), grid_[1].begin(),
|
||||
[](double d) { return M_PI * d / 180.0; });
|
||||
full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 2.0 * PI);
|
||||
|
||||
lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
|
||||
upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
|
||||
|
|
@ -1225,7 +1221,7 @@ double SphericalMesh::find_r_crossing(
|
|||
double SphericalMesh::find_theta_crossing(
|
||||
const Position& r, const Direction& u, double l, int shell) const
|
||||
{
|
||||
// Theta grid is [0, 180], thus there is no real surface to cross
|
||||
// Theta grid is [0, pi], thus there is no real surface to cross
|
||||
if (full_theta_ && (shape_[1] == 1))
|
||||
return INFTY;
|
||||
|
||||
|
|
@ -1287,7 +1283,7 @@ double SphericalMesh::find_theta_crossing(
|
|||
double SphericalMesh::find_phi_crossing(
|
||||
const Position& r, const Direction& u, double l, int shell) const
|
||||
{
|
||||
// Phi grid is [0, 360], thus there is no real surface to cross
|
||||
// Phi grid is [0, 2pi], thus there is no real surface to cross
|
||||
if (full_phi_ && (shape_[2] == 1))
|
||||
return INFTY;
|
||||
|
||||
|
|
@ -1368,25 +1364,20 @@ int SphericalMesh::set_grid()
|
|||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
}
|
||||
if (grid_[1].back() > 180) {
|
||||
if (grid_[1].back() > PI) {
|
||||
set_errmsg("theta-grids for "
|
||||
"spherical meshes must end with theta <= 180 degree.");
|
||||
"spherical meshes must end with theta <= pi.");
|
||||
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
if (grid_[2].back() > 360) {
|
||||
if (grid_[2].back() > 2 * PI) {
|
||||
set_errmsg("phi-grids for "
|
||||
"spherical meshes must end with phi <= 180 degree.");
|
||||
"spherical meshes must end with phi <= 2*pi.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 180.0);
|
||||
full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 360.0);
|
||||
|
||||
// Transform theta- and phi-grid from degrees to radians
|
||||
for (int i = 1; i < 3; i++)
|
||||
std::transform(grid_[i].begin(), grid_[i].end(), grid_[i].begin(),
|
||||
[](double d) { return M_PI * d / 180.0; });
|
||||
full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == PI);
|
||||
full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 2 * PI);
|
||||
|
||||
lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
|
||||
upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
|
||||
|
|
|
|||
|
|
@ -55,13 +55,13 @@
|
|||
</mesh>
|
||||
<mesh id="5" type="cylindrical">
|
||||
<r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
|
||||
<phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
|
||||
<phi_grid>0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586</phi_grid>
|
||||
<z_grid>-7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5</z_grid>
|
||||
</mesh>
|
||||
<mesh id="6" type="spherical">
|
||||
<r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
|
||||
<theta_grid>0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0</theta_grid>
|
||||
<phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
|
||||
<theta_grid>0.0 0.39269908169872414 0.7853981633974483 1.1780972450961724 1.5707963267948966 1.9634954084936207 2.356194490192345 2.748893571891069 3.141592653589793</theta_grid>
|
||||
<phi_grid>0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586</phi_grid>
|
||||
</mesh>
|
||||
<filter id="1" type="mesh">
|
||||
<bins>1</bins>
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
import numpy as np
|
||||
from math import pi
|
||||
|
||||
import openmc
|
||||
import pytest
|
||||
|
|
@ -53,13 +54,13 @@ def model():
|
|||
|
||||
cyl_mesh = openmc.CylindricalMesh()
|
||||
cyl_mesh.r_grid = np.linspace(0, 7.5, 18)
|
||||
cyl_mesh.phi_grid = np.linspace(0, 360, 19)
|
||||
cyl_mesh.phi_grid = np.linspace(0, 2*pi, 19)
|
||||
cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17)
|
||||
|
||||
sph_mesh = openmc.SphericalMesh()
|
||||
sph_mesh.r_grid = np.linspace(0, 7.5, 18)
|
||||
sph_mesh.theta_grid = np.linspace(0, 180, 9)
|
||||
sph_mesh.phi_grid = np.linspace(0, 360, 19)
|
||||
sph_mesh.theta_grid = np.linspace(0, pi, 9)
|
||||
sph_mesh.phi_grid = np.linspace(0, 2*pi, 19)
|
||||
|
||||
# Create filters
|
||||
reg_filters = [
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue