-
Notifications
You must be signed in to change notification settings - Fork 35
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
WIP: Another stab at more general coordinate support #1137
base: develop
Are you sure you want to change the base?
Conversation
template <int dir> | ||
KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { | ||
static_assert(dir > 0 && dir < 4); | ||
return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
looks like you're shadowing these in derived classes. Should they be marked overload
or something? Or is that not necessary due to CRTP?
|
||
template <int dir> | ||
KOKKOS_FORCEINLINE_FUNCTION | ||
Real hx(const Real r, const Real th, const Real phi) const { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
so this is the diagonal metric h
factor but what is the x
for in the name?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
And do we need any additional curvature terms we should track?
if (dir == X1DIR) return static_cast<const System *>(this)->template Dxc<X1DIR>(k, j, i); | ||
else if (dir == X2DIR) return static_cast<const System *>(this)->template Dxc<X2DIR>(k, j, i); | ||
else if (dir == X3DIR) return static_cast<const System *>(this)->template Dxc<X3DIR>(k, j, i); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are the if
s necessary here?
template <int dir> | ||
KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { | ||
static_assert(dir > 0 && dir < 4); | ||
if constexpr (dir == X1DIR) return static_cast<const System *>(this)->template Dxc<dir>(i); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
FWIW, could just have one line return static_cast<const System *>(this)->template Dxc<dir>(std::get<dir - 1>(std::make_tuple(i, j, k)));
// Xf: Positions on Faces | ||
//---------------------------------------- | ||
template <int dir> | ||
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this still be supported?
} | ||
|
||
template <int dir, TopologicalElement el> | ||
KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this still be supported?
using TE = TopologicalElement; | ||
if (cl == CellLevel::same) { | ||
if (el == TE::CC) { | ||
return cell_volume_; | ||
} else if (el == TE::F1) { | ||
return area_[X1DIR - 1]; | ||
} else if (el == TE::F2) { | ||
return area_[X2DIR - 1]; | ||
} else if (el == TE::F3) { | ||
return area_[X3DIR - 1]; | ||
} else if (el == TE::E1) { | ||
return dx_[X1DIR - 1]; | ||
} else if (el == TE::E2) { | ||
return dx_[X2DIR - 1]; | ||
} else if (el == TE::E3) { | ||
return dx_[X3DIR - 1]; | ||
} else if (el == TE::NN) { | ||
return 1.0; | ||
} | ||
} else if (cl == CellLevel::fine) { | ||
if (el == TE::CC) { | ||
return cell_volume_ / 8.0; | ||
} else if (el == TE::F1) { | ||
return area_[X1DIR - 1] / 4.0; | ||
} else if (el == TE::F2) { | ||
return area_[X2DIR - 1] / 4.0; | ||
} else if (el == TE::F3) { | ||
return area_[X3DIR - 1] / 4.0; | ||
} else if (el == TE::E1) { | ||
return dx_[X1DIR - 1] / 2.0; | ||
} else if (el == TE::E2) { | ||
return dx_[X2DIR - 1] / 2.0; | ||
} else if (el == TE::E3) { | ||
return dx_[X3DIR - 1] / 2.0; | ||
} else if (el == TE::NN) { | ||
return 1.0; | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess this still needs to be updated.
// Xc: Positions at cell centroids | ||
//---------------------------------------- | ||
template <int dir> | ||
KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
should this take k, j, i
?
else if (el == TE::E3) return EdgeLength<X3DIR>(k, j, i); | ||
else if (el == TE::NN) return 1.0; | ||
} else { | ||
PARTHENON_FAIL("Have not yet implemented fine fields for UniformSpherical coordinates."); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We could just have a coordinates object own another coordinates object that has twice the number of zones and then just call that for fine fields.
PR Summary
Add support for spherical and cylindrical coordinates, as well as supporting downstream defined coordinate systems.
PR Checklist