-
Notifications
You must be signed in to change notification settings - Fork 95
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
Shiftedmetric updates #1176
Shiftedmetric updates #1176
Conversation
90fd318
to
4790016
Compare
0add534
to
e616269
Compare
This PR addresses a few things that aren't really related:
I think it's ready to merge, but should I split it up? |
@@ -609,12 +609,12 @@ class Mesh { | |||
typedef BoutReal (*flux_func)(stencil&, stencil &); | |||
|
|||
/// Transform a field into field-aligned coordinates | |||
const Field3D toFieldAligned(const Field3D &f) { | |||
return getParallelTransform().toFieldAligned(f); | |||
const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) { |
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.
Is there a reason for this choice of default region?
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 the reason I did this is that:
- we (usually, probably) want to calculate the result in the y-guard cells, to avoid interpolation and needing to apply parallel boundary conditions before taking derivatives or interpolating
- we don't want to shift in the corner guard cells, because they contain invalid data and I think I was getting exceptions thrown when including them
- we didn't have a region that excludes just the corner guard cells, but includes x- and y-guard cells
- since radial derivatives should probably be done in the orthogonal coordinates rather than the field-aligned ones, the x-guard cells shouldn't normally be needed.
include/field3d.hxx
Outdated
@@ -233,6 +233,11 @@ class Field3D : public Field, public FieldData { | |||
* Ensure that yup and ydown refer to this field | |||
*/ | |||
void mergeYupYdown(); | |||
|
|||
/*! | |||
* Clear all yup/ydown fields and field_fa |
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 do you mean by clear (delete, zero out, something else)?
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 mean delete (if they exist). It would make more sense to change the name...
src/field/field3d.cxx
Outdated
@@ -375,6 +378,7 @@ Field3D & Field3D::operator=(const Field3D &rhs) { | |||
|
|||
setLocation(rhs.location); | |||
|
|||
clearYupYdown(); |
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.
Would it be better to copy the yup/ydown (as we do for the data) in this 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.
Would it be better to copy the yup/ydown (as we do for the data) in this case?
I don't know. I guess it only makes sense to copy yup/ydown if we also calculate yup/ydown in arithmetic operators. But if we calculate yup/ydown everywhere possible by default, we'd be doing a lot of computations that might not be needed (if we don't take a y-derivative, etc.). It would be nice to give an option for including yup/ydown fields, but I don't see a nice way to do it...
src/mesh/index_derivs.cxx
Outdated
Field2D result(this); | ||
result.allocate(); // Make sure data allocated | ||
result.setLocation(loc); |
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 seems to (potentially) be changing the output location -- it's now certain to be at the requested loc
which sounds sensible, but by default I think loc
is CELL_DEFAULT
and I don't think this makes sense in a call to setLocation
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.
src/mesh/index_derivs.cxx
Outdated
@@ -872,6 +866,7 @@ const Field2D Mesh::applyYdiff(const Field2D &var, Mesh::deriv_func func, CELL_L | |||
|
|||
Field2D result(this); | |||
result.allocate(); // Make sure data allocated | |||
result.setLocation(diffloc); |
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 is different to the X functions which now use the passed loc
rather than diffloc
-- is there a reason for the difference?
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.
It was different because staggering was not implemented for applyYdiff(Field2D). I've implemented it now in 2a35d4d.
src/mesh/index_derivs.cxx
Outdated
Field3D result(this); | ||
result.allocate(); // Make sure data allocated | ||
result.setLocation(loc); |
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 is diffloc
in earlier y derivative function.
src/mesh/index_derivs.cxx
Outdated
result = applyYdiff(f, func, diffloc, region); | ||
|
||
result.setLocation(diffloc); // Set the result location | ||
Field3D result = applyYdiff(f, func, diffloc, region); |
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.
Could just directly return applyYdiff(f, func, diffloc, region);
I think.
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.
Yep, updated in ce76aa2.
src/mesh/parallel/shiftedmetric.cxx
Outdated
for(int jy=0;jy<mesh.LocalNy;jy++) { | ||
shiftZ(f(jx,jy), mesh.LocalNz, zangle(jx,jy), result(jx,jy)); | ||
} | ||
invalidateGuards(result); // Won't set x-guard cells, so allow checking to throw exception if they are used. |
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 shouldn't be needed. Creating a field and allocating should give you a field with invalid guards (if check is at the right level). Let me know if this isn't the case and I can look at why/fixing it!
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.
That's right. Sorry, I was trying to be over-cautious. Removed in b4934a5.
|
||
// We only use methods in ShiftedMetric to get fields for parallel operations | ||
// like interp_to or DDY. | ||
// Therefore we don't need x-guard cells, so do not set them. |
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 we might want to set the x-guard cells in some future situations. Is it easy enough to support this even if the default is to not do it?
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 it's not trivial at the moment, but should be when Region<T>
goodness is available #1245, if we define a region with everything but the corner guard cells, which is maybe what RGN_ALL
should be?
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 wanted to avoid using zShift in the guard cells, because in future at least, zShift on staggered grids might be calculated by interpolating and communicating, so the corner guard cells might not be initialized. At least I think that was the issue I was having before that made valgrind complain, maybe the better solution would just be to make sure the corner guard cells are initialized to something even when checking is turned off?
b4934a5
to
3e7b005
Compare
3e7b005
to
ab43f72
Compare
Remove the intermediate variable result_fa in branch of interp_to that transforms to field-aligned variables. Just store the intermediate result in 'result' instead, to save memory.
ab43f72
to
545aae5
Compare
Add region arguments (where it is valid to give RGN_NOX or RGN_NOBNDRY) so functions like toFieldAligned/fromFieldAligned can be used both when y-guard cells have been set and when they have not. e.g. when using fromFieldAligned on the results of interpolation or derivatives, y-guard cells cannot be included, but when using toFieldAligned on the inputs to interpolation or derivatives, y-guard cells must be included.
545aae5
to
79b5157
Compare
Ensure that no fields in the rhs function/method are using yup/ydown before calling communicate.
Avoids cost of new/delete.
Tests of splitYupYdown() and mergeYupYdown() need to call field.setHasValidYUpDown(true);
Fixing the low-level shiftZ method means that toFieldAligned/fromFieldAligned also do the right thing. Check that input to shiftZ is at same location as the shift angle zShift (i.e. CELL_CENTRE). Set location of result returned from zShift to be the same as the location of the input.
8d9ef31
to
d8de292
Compare
I've pulled out several of the commits, that were independent of the rest, into separate PRs (#1427, #1432, #1433). What's left is an attempt to delete, or mark as invalid, |
A collection of small updates:
These updates are in preparation for adding a shifted-metric type class that shifts variables to field-aligned coordinates (which is helpful when also using staggered grids), and for an update to ShiftedMetric to include 2 points up and down the field (if myg>1).