-
Notifications
You must be signed in to change notification settings - Fork 13
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
Reflecting BC #79
base: main
Are you sure you want to change the base?
Reflecting BC #79
Conversation
@pgrete @forrestglines I'd also like to get this merged in, since it's needed for my simulations. Let me know what would need to be polished to make it suitable for merging to main. |
So this is working and ready for review I'm guessing? I'll remove the WIP if that's the case. |
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 think this PR needs a test to exercise reflecting boundary conditions. A shock reflecting off a boundary should be pretty easy to add with the other Sod shock tube tests. I can do that next week if need be.
We can also add quick documentation.
I've also commented on a few code snippets that I think could be placed in Parthenon but I don't think those should hold up this PR. We can add those snippets to Parthenon later.
inline auto GetPhysicalZones(parthenon::MeshBlock *pmb, parthenon::IndexShape &bounds) | ||
-> std::tuple<parthenon::IndexRange, parthenon::IndexRange, parthenon::IndexRange> { |
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.
Any reason not to use:
inline auto GetPhysicalZones(parthenon::MeshBlock *pmb, parthenon::IndexShape &bounds) | |
-> std::tuple<parthenon::IndexRange, parthenon::IndexRange, parthenon::IndexRange> { | |
inline std::tuple<parthenon::IndexRange, parthenon::IndexRange, parthenon::IndexRange> | |
GetPhysicalZones(parthenon::MeshBlock *pmb, parthenon::IndexShape &bounds) { |
I'm not too used to this syntax, I only see it when the return type needs to be inferred with decltype
in a template function
void ReflectingInnerX3(std::shared_ptr<parthenon::MeshBlockData<Real>> &mbd, bool coarse); | ||
void ReflectingOuterX3(std::shared_ptr<parthenon::MeshBlockData<Real>> &mbd, bool coarse); |
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.
Why are these only defined for X3, don't we need reflecting boundary conditions for all dimensions?
That would be an easy template to 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.
Where are these used? In a pgen in another branch?
constexpr parthenon::IndexDomain domain = | ||
INNER ? (X1 ? parthenon::IndexDomain::inner_x1 | ||
: (X2 ? parthenon::IndexDomain::inner_x2 | ||
: parthenon::IndexDomain::inner_x3)) | ||
: (X1 ? parthenon::IndexDomain::outer_x1 | ||
: (X2 ? parthenon::IndexDomain::outer_x2 | ||
: parthenon::IndexDomain::outer_x3)); |
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.
Likewise this could be in Parthenon. An is_inner
function of IndexDomain
or like
return std::tuple<parthenon::IndexRange, parthenon::IndexRange, parthenon::IndexRange>{ | ||
parthenon::IndexRange{IsDomainBound(pmb, parthenon::BoundaryFace::inner_x1) | ||
? bounds.is(parthenon::IndexDomain::interior) | ||
: bounds.is(parthenon::IndexDomain::entire), | ||
IsDomainBound(pmb, parthenon::BoundaryFace::outer_x1) | ||
? bounds.ie(parthenon::IndexDomain::interior) | ||
: bounds.ie(parthenon::IndexDomain::entire)}, | ||
parthenon::IndexRange{IsDomainBound(pmb, parthenon::BoundaryFace::inner_x2) | ||
? bounds.js(parthenon::IndexDomain::interior) | ||
: bounds.js(parthenon::IndexDomain::entire), | ||
IsDomainBound(pmb, parthenon::BoundaryFace::outer_x2) | ||
? bounds.je(parthenon::IndexDomain::interior) | ||
: bounds.je(parthenon::IndexDomain::entire)}, | ||
parthenon::IndexRange{IsDomainBound(pmb, parthenon::BoundaryFace::inner_x3) | ||
? bounds.ks(parthenon::IndexDomain::interior) | ||
: bounds.ks(parthenon::IndexDomain::entire), | ||
IsDomainBound(pmb, parthenon::BoundaryFace::outer_x3) | ||
? bounds.ke(parthenon::IndexDomain::interior) | ||
: bounds.ke(parthenon::IndexDomain::entire)}}; |
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 snippet of code seems like it could belong in Parthenon. I'm not sure about the name, GetZonesInDomain
or something else referencing the problem domain might be better.
q(l, X3 ? offset - k : k, X2 ? offset - j : j, X1 ? offset - i : i); | ||
} else { | ||
q(l, k, j, i) = q(l, X3 ? ref : k, X2 ? ref : j, X1 ? ref : 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.
I might wonder about executing these ternary statements inside the kernel and wonder if the CUDA compiler is smart enough to compile them out.
However, this kernel is certainly memory bandwidth limited so the extra integer computes are irrelevant
//! \file bc.hpp | ||
// \brief Custom boundary conditions for AthenaPK | ||
// | ||
// Computes reflecting boundary conditions using AthenaPK's cons variable pack. |
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 this is only for reflecting boundary conditions, should this file be renamed? Should we have a bc
folder?
Or maybe this file can keep this name as a place holder for other commonly used boundary functions in AthenaPK
template <parthenon::CoordinateDirection DIR, BCSide SIDE, BCType TYPE> | ||
void ApplyBC(parthenon::MeshBlock *pmb, parthenon::VariablePack<Real> &q, | ||
parthenon::IndexRange &nvar, const bool is_normal, const bool coarse) { |
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 we instead template on:
template <parthenon::CoordinateDirection DIR, BCSide SIDE, BCType TYPE> | |
void ApplyBC(parthenon::MeshBlock *pmb, parthenon::VariablePack<Real> &q, | |
parthenon::IndexRange &nvar, const bool is_normal, const bool coarse) { | |
template <parthenon::BoundaryFace FACE, BCType TYPE> | |
void ApplyBC(parthenon::MeshBlock *pmb, parthenon::VariablePack<Real> &q, | |
parthenon::IndexRange &nvar, const bool is_normal, const bool coarse) { |
and remove the BCSide
enum? Would this work with your other 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.
I'm not quite sure how this function is used/should be used. In this PR it looks like it's only used to apply reflecting boundary conditions. Don't all the other pgen's have outflow boundary conditions applied as the default? Unless I've misunderstood something.
@BenWibking Is there a branch/pgen using the reflecting boundary conditions? If it's the precipitator boxes, could you link that here? |
|
User boundary conditions recently got an (internal) overhaul in Parthenon, which might make interfacing cleaner (and reduce the amount of code required on our end). |
Sure, if you want to create a new PR with that, I can close this one. |
New version of #58. (Renamed the branch and didn't realize the PR would be closed automatically.)