-
Notifications
You must be signed in to change notification settings - Fork 37
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
Cleanup UniformCartesian #753
Cleanup UniformCartesian #753
Conversation
This comment was marked as outdated.
This comment was marked as outdated.
This comment was marked as outdated.
This comment was marked as outdated.
cb93802
to
72c6002
Compare
I ended up implementing a slightly different interface for Position of elements (previously
where *In particular with I've also dropped versions using Distance between elements (previously
Likewise, Integration elements (previously
These functions take all three indices to accommodate extension to spherical and cylindrical. They need
For each of these functions templated on
This convention avoids naming conflicts between simplified functions and runtime functions. We can implement these later as needed. We'll probably need to implement more coordinate functions as we need them for cylindrical and spherical coordinates. |
@par-hermes format |
I think these choices make sense. My one concern is regarding |
I encountered an error that I'm a little confused about as I tried out the new coordinates in AthenaPK. It occurred in this templated member function that was kind-of pointlessly templated on a coordinate class.
And when I call this function in the class
It causes this error with clang:
and this error with CUDA/11.7.0:
I don't think it's related to the Problematic AthenaPK code here: hydrostatic_equilibrium_sphere.cpp Work-around AthenaPK code here: |
I used this script to update AthenaPK to use the new coordinate names. It should work for any other derivative Parthenon code (or any Parthenon PR as well).
Edited once more for changes to API |
That's bizarre. I have no idea. |
I renamed these to |
I like being verbose in this case. 👍 |
👍 for the new naming scheme. I'm also in favor of the more verbose description and cleaned up interface! |
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.
LGTM but confused about template vs runtime call use in one spot, see comment
src/mesh/mesh_refinement.cpp
Outdated
const Real area00 = coords.faceArea(X1DIR, k, j, i); | ||
const Real area01 = coords.faceArea(X1DIR, k, j + 1, i); | ||
const Real area10 = coords.faceArea(X1DIR, k + 1, j, i); | ||
const Real area11 = coords.faceArea(X1DIR, k + 1, j + 1, 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.
Can this & below be FaceArea<XNDIR>()
? If not, is that something we can fix with constexpr?
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.
Similarly, is it worth standardizing all N->XNDIR using this PR, maybe with additional XNFACE variables?
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.
Can this & below be
FaceArea<XNDIR>()
? If not, is that something we can fix with constexpr?
These can be replaced, I just pushed that change along with a few other files.
Similarly, is it worth standardizing all N->XNDIR using this PR, maybe with additional XNFACE variables?
Maybe? It would mean a lot more function name bloat (i.e. x1f
->Xf<X1DIR,X1DIR>
) when it should be clear from context that the integers refer to directions and it wouldn't add any extra checks/code safety.
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.
Agreed it's verbose and unwieldy to use the macros. I just kept forgetting what Xf<1, 1> was referring to. Especially in terms of the argument order, since the convention switches by variable. Xf<X1FACE, X1DIR> or Dxf<X1DIR, X1FACE> is at least totally explicit. (I absolutely mixed up both those orders the first time).
It also wouldn't add any type checks, true. Which would make it hard to enforce, since using the macros would be purely a matter of convention and require much more typing.
...but that's given me a better (worse) idea: maybe we should type-check face locations vs edge locations vs directions? I mean, Xf<FaceLocation f, Direction d> seems like even more overkill, but at least then we get the compile-time check benefits and it's easy to enforce a readable standard.
Maybe the answer for now is just to use the Xf<int>
shortcut everywhere? That's unambiguous, resembles the old Athena convention, and it's all we need for UniformCartesian, right? And, it and leaves room for a future typed solution, since changing the full macro only involves changing the shortcut and not a bunch of calls scattered around the code.
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.
The Xf<int,int>
template is so that we can use Xf
for all references to locations, replacing the x*f
and x*s*
functions in Athena++.
That said, we're currently only using Xf<x,x>
in Parthenon (and Athena++) although that will change soonish when we add the edge-corrections for face centered fields on AMR blocks. We can add an overloaded template Xf<int>
when both the face and coordinate direction are the same.
For enforcing typing, we would have to change CoordinateDirection
from an enum
to enum class
I believe. Hopefully that wouldn't break anything. Anywhere we really need to covert from an int to this enum would be slightly more awkward. I would be fine with that change as long as people are onboard with the code bloat (in terms of characters added) that would create.
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 dunno if I want to move CoordinateDirection
to an enum class
... that's going to change everything everywhere.
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 definitely get why the <int, int> template is necessary. I was proposing something like:
template <int dirface, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const {
return Dxf<dirface, dirface>(args);
}
i.e., just overloading and calling through. I misread the diff and thought that a similar overload was already implemented. I think this would be simple enough to do in a general way, even if the specific code above doesn't work, but the difference is obviously minimal. If this isn't uniformly agreed to be more readable it's not worth it.
Re: making directions (and faces) an enum class, I am abstractly very in favor of this. Probably not in this PR though; it was just an idea for someday, and as Jonah says it will mean a lot of lines changed, even if there's always a clear way to translate the code.
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.
template <int dirface, class... Args> KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { return Dxf<dirface, dirface>(args); }
I see, that should be implemented, I've now added it in.
I dunno if I want to move
CoordinateDirection
to anenum class
... that's going to change everything everywhere.
That's mostly why I'm hesitant to make that change. But with just the enum, we don't have any type checking, so using the enum everywhere at the moment won't add much. In the few places where Xf<int,int>
will be used, it could be advantageous to use different enums to differentiate between coordinate direction and face direction.
Re: making directions (and faces) an enum class, I am abstractly very in favor of this. Probably not in this PR though; it was just an idea for someday, and as Jonah says it will mean a lot of lines changed, even if there's always a clear way to translate the code.
Yep, let's save it for another PR
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.
Overall I like the new naming scheme.
I was able to go through Phoebus and build and run with the new uniform-cartesian setup but it was not as trivial as I'd hoped. A number of locations in the code required manual intervention. I list a few questions/comments related to my manual interventions below.
IMO none of these are blocking, but I'd like to chat about them before we merge, so I'm requesting changes for now, and if you push back about the changes, I'll just click approve. :)
template <int dir, int face, class... Args> | ||
KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { | ||
assert(dir > 0 && dir < 4 && face > 0 && face < 4); | ||
return dx_[face - 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.
I'm not sure why... but this function in particular seemed to throw C++'s type deduction through a loop. A number of places, where I had
const auto &coords = pmb->coords;
I had to change to
const parthenon::Coordinates_t &coords = pmb->coords;
for template resolution of, e.g.,
coords.DXf<1,1>(i)
to resolve properly.
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.
Hmm, is this class... Args
causing the type deduction causing issues? Did it happen on multiple compilers?
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'm not sure... that's the likeliest culprit but I didn't see what was wrong. Only tried on gnu.
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 don't think this is urgent to fix, by the way.
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.
This issue came up when I switched from function argument functions to template arguments in the Poisson example.
/vast/home/forrestglines/code/parthenon-project/parthenon/example/poisson/poisson_package.cpp(233): error: expected an expression
/vast/home/forrestglines/code/parthenon-project/parthenon/example/poisson/poisson_package.cpp(162): error: expected an expression
2 errors detected in the compilation of "/vast/home/forrestglines/code/parthenon-project/parthenon/example/poisson/poisson_package.cpp".
Changing out auto
fixed the compilation issue
I also notice a number of tests seem to be failing... possibly related to the template type deduction issue I noticed above. |
What tests where are failing? In Parthenon or in Phoebus? The last merge introduced code using the old naming conventions so I need to fix those before it will compile. |
I just meant tests in parthenon. Phoebus tests passed after I fixed things. |
@jonahm-LANL This is ready for review again |
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.
Sorry for the delay---been under the weather. LGTM. Two minor nitpicks below. Neither blocking.
std::floor((x - x_min_) / coords_.Dx(CoordinateDirection::X1DIR))) + | ||
std::floor((x - x_min_) / coords_.CellWidthFA(CoordinateDirection::X1DIR))) + | ||
ib_s_; | ||
j = (ndim_ > 1) ? static_cast<int>(std::floor( | ||
(y - y_min_) / coords_.Dx(CoordinateDirection::X2DIR))) + | ||
jb_s_ | ||
: jb_s_; | ||
k = (ndim_ > 2) ? static_cast<int>(std::floor( | ||
(z - z_min_) / coords_.Dx(CoordinateDirection::X3DIR))) + | ||
kb_s_ | ||
: kb_s_; | ||
j = (ndim_ > 1) | ||
? static_cast<int>(std::floor( | ||
(y - y_min_) / coords_.CellWidthFA(CoordinateDirection::X2DIR))) + | ||
jb_s_ | ||
: jb_s_; | ||
k = (ndim_ > 2) | ||
? static_cast<int>(std::floor( | ||
(z - z_min_) / coords_.CellWidthFA(CoordinateDirection::X3DIR))) + | ||
kb_s_ | ||
: kb_s_; |
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've noticed lots of places where you've made this change: Dx -> CellWidthFA. I think the better replacement rule would be Dx -> Dxc<dir>
as that's probably the intent in all these places and it'll be the faster option down the road. For examples it also offers more concise and parthenonic code.
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.
If I recall correctly, Dx
was a run time determined function inherited from Athena++ and thus became CellWidthFA
from the Athena++ comments on its behavior vs. dx1v
. I think it might make a difference for some coordinate systems. Everywhere Dx
was used in Parthenon, however, would be better served by Dxc<dir>
so I'll make that change
@@ -146,10 +146,10 @@ void VTKOutput::WriteContainer(SimTime &tm, Mesh *pm, ParameterInput *pin, bool | |||
// write x1-coordinates as binary float in big endian order | |||
std::fprintf(pfile, "X_COORDINATES %d float\n", ncoord1); | |||
if (ncells1 == 1) { | |||
data[0] = static_cast<float>(pmb->coords.x1v(out_ib.s)); | |||
data[0] = static_cast<float>(pmb->coords.Xc<1>(out_ib.s)); |
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.
Does the vtk output even work? Not for this MR but maybe we should either revive it for real or kill it.
@forrestglines can you share the most recent version of your script to find/replace all the instances? |
The updated script in the comment above did the job for |
The script in the comment should work. I forgot that I had an old script in the description of the PR. I've now updated that to match the new interface |
Cool thanks, @forrestglines @pgrete |
PR Summary
Cleanup the current UniformCartesian coordinate class to have fewer functions, utilizing templates on coordinate direction to reduce the number of functions. Done in preparation to implement other coordinate system classes.
Here are the new functions and any differences highlighted
Position of elements (previously
x1v
,x1f
, andx1s2
*)where
1<=dir<=3
is which component of the direction being queried (i.e. x,y,z),face
andedge
specify which face/edge is being queried, andidx
is the index alongdir
. The addition offace
andedge
along with the functionsXe
andXn
are new. These allow for consistently getting the positions of faces, edges, and nodes with one function without mixing calls toXc
andXf
to get the position of off-cell-center locations.*In particular with
Xf
, whenface==dir
it returns the same asx1f
etc. did in the old version. Now whenface!=dir
, it returns the face-averaged position of the other directions, replacingx1s2
,x2s3
, etc.I've also dropped versions using
i,j,k
since the uniform cartesian, spherical, and cylindrical coordinates do not need them. If we want non-uniform coordinate systems later, we'll have to go back and change them.Distance between elements (previously
dx1v
,dx1f
etc. )Likewise,
face
,edge
,Dxe
, andDxn
are new to get distances between different elements along different dimensions.Integration elements (previously
Dx
,Area
,Edgelength
, andVolume
)These functions take all three indices to accommodate extension to spherical and cylindrical. They need
k
,j
, andi
for the sake of spherical and cylindrical coordinates but for now these arguments are unenforced in Parthenon.Area-averaged positions on face-centers (previouslyx1s2
,x3s1
etc. )These are used in the AMR machinery for face-centered fields, which is currently missing from Parthenon. However, I believe these are identical toXf<int dir, int face>
for the uniform coordinates proposed.Xs
might not be necessary anymore.*I removed
Xs
since the newXf
can give all the same information unambiguously.For each of these functions templated on
dir
,face
,edge
, etc, we can introduce undercase version with runtime parameters, i.e.This convention avoids naming conflicts between simplified functions and runtime functions. We can implement these later as needed.
We'll probably need to implement more coordinate functions as we need them for cylindrical and spherical coordinates.
Also, here is a script to do renaming of coordinates in downstream codes:
PR Checklist