crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Corrections to application of origin for sphere and cyl mesh

Browse source
This commit is contained in:
Patrick Shriwise 2023-01-26 14:13:39 -06:00
parent f0721ac535
commit 11eee756dd
2 changed files with 16 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

3
include/openmc/mesh.h
View file

@ -241,6 +241,7 @@ public:
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
Position origin_ {0.0, 0.0, 0.0};
protected:
};
@ -359,7 +360,6 @@ public:
void to_hdf5(hid_t group) const override;
array<vector<double>, 3> grid_;
Position origin_;
int set_grid();
@ -414,7 +414,6 @@ public:
void to_hdf5(hid_t group) const override;
array<vector<double>, 3> grid_;
Position origin_;
int set_grid();

27
src/mesh.cpp
View file

@ -412,7 +412,6 @@ template<class T>
void StructuredMesh::raytrace_mesh(
Position r0, Position r1, const Direction& u, T tally) const
{
// TODO: when c++-17 is available, use "if constexpr ()" to compile-time
// enable/disable tally calls for now, T template type needs to provide both
// surface and track methods, which might be empty. modern optimizing
@ -445,6 +444,11 @@ void StructuredMesh::raytrace_mesh(
return;
}
// translate start and end positions,
// this needs to come after the get_indices call because it does its own translation
r0 -= origin_;
r1 -= origin_;
// Calculate initial distances to next surfaces in all three dimensions
std::array<MeshDistance, 3> distances;
for (int k = 0; k < n; ++k) {
@ -972,13 +976,12 @@ std::string CylindricalMesh::get_mesh_type() const
StructuredMesh::MeshIndex CylindricalMesh::get_indices(
Position r, bool& in_mesh) const
{
Position mapped_r;
r -= origin_;
Position mapped_r;
mapped_r[0] = std::hypot(r.x, r.y);
mapped_r[2] = r[2];
mapped_r -= origin_;
if (mapped_r[0] < FP_PRECISION) {
mapped_r[1] = 0.0;
} else {
@ -1090,22 +1093,23 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
double l) const
{
Position r = r0 - origin_;
if (i == 0) {
return std::min(
MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
MeshDistance(ijk[i] + 1, true, find_r_crossing(r, u, l, ijk[i])),
MeshDistance(ijk[i] - 1, false, find_r_crossing(r, u, l, ijk[i] - 1)));
} else if (i == 1) {
return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
find_phi_crossing(r0, u, l, ijk[i])),
find_phi_crossing(r, u, l, ijk[i])),
MeshDistance(sanitize_phi(ijk[i] - 1), false,
find_phi_crossing(r0, u, l, ijk[i] - 1)));
find_phi_crossing(r, u, l, ijk[i] - 1)));
} else {
return find_z_crossing(r0, u, l, ijk[i]);
return find_z_crossing(r, u, l, ijk[i]);
}
}
@ -1209,12 +1213,11 @@ std::string SphericalMesh::get_mesh_type() const
StructuredMesh::MeshIndex SphericalMesh::get_indices(
Position r, bool& in_mesh) const
{
r -= origin_;
Position mapped_r;
mapped_r[0] = r.norm();
mapped_r -= origin_;
if (mapped_r[0] < FP_PRECISION) {
mapped_r[1] = 0.0;
mapped_r[2] = 0.0;