-
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
Boundary bugfix(es) and exposing FluxDiv_
interface
#558
Conversation
src/interface/update.hpp
Outdated
@@ -33,6 +33,23 @@ namespace parthenon { | |||
|
|||
namespace Update { | |||
|
|||
// Calculate the flux divergence for a specific component l of a variable v | |||
KOKKOS_FORCEINLINE_FUNCTION | |||
Real FluxDiv_(const int l, const int k, const int j, const int i, const int ndim, |
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 should consider renaming FluxDiv_
to FluxDiv
or FluxDivHelper
or something if it's publicly exposed.
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.
Done.
…ntially repeats work, but is more robust and allows the user to set only the interior in the problem generator
@pgrete Changing k-bounds initialized the uninitialized variables in the ghost zones at the initial time. I just pushed to your branch. |
What do you think about updating the advection example so that the background has a value of "1" versus "0" currently and then I'd add a |
I think this is a good idea, so long as the advected field now has whatever form it currently has + 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 great and great catch!
Since we're making changes to the boundary conditions file anyway, I thought this would be a good opportunity to reduce some code duplication. I've made one generic boundary condition function covering all 12 cases (X1, X2, X3) x (inner, outer) x (outflow, reflection). Hopefully, this would make it easier to avoid such copy&paste errors in the future (if the boundary condition code needs to change). Please let me know what you think.
In principal, I'm in favor of reduced code, but I worry that I can't quite parse whether or not the generic version covers all 12 specialized versions. I'm not sure we have an automated test to catch this. Perhaps we should add one. |
We should definitely have a test. However, I think it should be relatively straightforward to check. There's only a few lines that were different between the expanded 12 cases. That being said, I'm happy to move the code duplication reduction to a new PR if it would take too long to review as part of this PR. |
That's up to @pgrete I suppose. My preference would be to add an automated test before merging that ensures these are working as intended. Advecting off of each face in 3d would be a. reasonable regression test, I think. We've clearly just been burned by our lack of testing in this area. |
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 like the reduced code (see also comment on readability).
I also think that extending the test coverage is a good idea.
I have an idea to use dc reconstruction with a cfl of 1, so that we an check inflow, outflow, reflect exactly as as there'd be no diffusion.
Co-authored-by: Philipp Grete <gretephi@msu.edu>
Waiting on test for boundary condition before merging |
What kind of test are you looking for? simple unit test exercising the templated boundary condition function? |
Yes exactly. Just to make sure it behaves the same as the older non-templated code. (And/or is correct. ;) ) |
023ddf1
to
f287238
Compare
This is now ready for new comments/review.
test_reflect.mp4Regarding the user boundary conditions: |
Ping @Yurlungur for final re-review of latest changes as discussed on Thursday. |
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 @pgrete some minor comments below. Really fantastic to have this test in place, along with the reduced code.
if (fill_derived) { | ||
field_name = "one_minus_advected"; | ||
m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, | ||
std::vector<int>({num_vars})); | ||
pkg->AddField(field_name, m); | ||
pkg->AllFields(); |
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.
What is this line supposed to be doing? Doesn't this return the const reference to the dictionary? Why the empty call?
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.
Good catch. Looks like a debugging leftover. I'm removing it.
// The following block/logic is also just added for regression testing. | ||
// More specifically, the "smooth_gaussian" profile is initially != 0 everywhere, but | ||
// initialializes IndexDomain::interior. | ||
// FillDerived (here, SquareIt) is called after the ghost cells are exchanged and over | ||
// IndexDomain::entire. | ||
// Thus, no 0 (from initializing Kokkos views) should be left if all faces/corners/edges | ||
// are correct, which is what we check in the loop below. |
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.
High value comment 👍
v.flux(X1DIR, n, k, j, i) = | ||
ql(idx_v, i) > 0.0 ? ql(n, i) * ql(idx_v, i) : 0.0; | ||
v.flux(X1DIR, n, k, j, i) += | ||
qr(idx_v, i) < 0.0 ? qr(n, i) * qr(idx_v, i) : 0.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 like the +=
syntax here. Feels nicer than the if-else block.
// bvals testing. | ||
} else { | ||
v.flux(X2DIR, n, k, j, i) = | ||
ql(idx_v + 1, i) > 0.0 ? ql(n, i) * ql(idx_v + 1, i) : 0.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.
rather than +1
or +2
, one could do idx_v + XNDIR - 1
which might be a little cleaner. Not a big deal though.
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.
Done.
@@ -53,7 +58,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) { | |||
|
|||
auto profile_str = pin->GetOrAddString("Advection", "profile", "wave"); | |||
if (!((profile_str == "wave") || (profile_str == "smooth_gaussian") || | |||
(profile_str == "hard_sphere"))) { | |||
(profile_str == "hard_sphere") || (profile_str == "block"))) { |
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 love the name block
for the test profile. But I don't have a better name either.
class TestCase(utils.test_case.TestCaseAbs): | ||
def Prepare(self, parameters, step): | ||
|
||
# Step 1 (default vals in parameter file) is reflecting BC |
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 know it's default to reflecting, but I still might be explicit here, in case the parameter file changes.
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.
Good point. Will do.
Addressed comment and triggering automerge now. |
PR Summary
Also adjust
dt
calc as the original version resulted in anan
when it was initially set tomax
PR Checklist