Skip to content
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

Fix the neighbor-finding in CInterpolator::ReconstructBoundary #1346

Merged
merged 35 commits into from
Aug 16, 2021

Conversation

maxaehle
Copy link
Contributor

@maxaehle maxaehle commented Aug 4, 2021

Proposed Changes

CInterpolator::ReconstructBoundary builds arrays Buffer_Receive_StartLinkedNodes, _LinkedNodes,... that store the connectivity of the marker, i.e. the list of neighbor vertices for each vertex.

In the previous version of the code, two vertices were considered "neighbors" if they were connected through an edge of a volume element ("volume neighbors"). We propose to consider vertices as "neighbors" only if they are connected through an edge of a surface element ("surface neighbors").

These two notions are equivalent for smooth surfaces and clearly surface neighbors are always volume neighbors as well. However, at surfaces with sharp edges, two vertices can be volume neighbors without being surface neighbors. This happens when the edge of the volume element connecting the vertices lies inside the domain.

CInterpolator::ReconstructBoundary is called from CSlidingMesh::SetTransferCoeff, which then hands the connectivity structure to CSlidingMesh::Build_3D_surface_element. The latter function forms a circular ordering of the neighbors of each vertex s: e.g. if there are four neighbors a,b,c,d of s with connectivity a-b,c-d,b-c,a-d then it performs an iteration over e.g. a-b-c-d. If we understand "neighbors" as "surface neighbors", the neighbors of any vertex s can be brought into a "closed circle" or "open chain" ordering where each of them is neighbor only to its direct predecessor and successor. This is not true for the volume neighbor relation however, and in this case the iteration might enter an infinite loop e-a-b-c-d-a-b-... (until the access to element[ iElementNode ] segfaults). That's approximately how we (@maxaehle,@mcmorelli,@BeckettZhou,@EduardoMolina) noticed the bug.

Related Work

#1344

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags, or simply --warnlevel=2 when using meson).
  • My contribution is commented and consistent with SU2 style.
  • I have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.

Copy link
Member

@pcarruscag pcarruscag left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for looking into this, these interpolation functions can get quite complicated.
I think there is some room to simplify the new algorithm a bit, and to improve the old code if you have some time.
Just before merging this PR please run clang-format on CInterpolator.cpp.

Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
Co-authored-by: Pedro Gomes <[email protected]>
Copy link
Member

@pcarruscag pcarruscag left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, run clang-format on CInterpolator.cpp before merging please, it's a bit messy with unaligned comments, etc.

@TobiKattmann
Copy link
Contributor

That's really strange: The hybrid AD testcase discadj_incomp_turb_NACA0012_sst was ok after commit c68cf13, then I changed the stored residuals of unrelated testcases in 9d8f194 and the testcase fails... I'll run clang-format now, let's see what happens ;)

I had a similar thing over at the flamelet branch and a re-run of the regression tests solved that :)

Copy link
Contributor

@TobiKattmann TobiKattmann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @maxaehle 💐 ... just a little nagging in a few places.

Common/src/interface_interpolation/CSlidingMesh.cpp Outdated Show resolved Hide resolved
Common/src/interface_interpolation/CSlidingMesh.cpp Outdated Show resolved Hide resolved
Comment on lines 294 to 297
donor_iMidEdge_point[iDim] = ( DonorPoint_Coord[ donor_forward_point ][ iDim] +
DonorPoint_Coord[ donor_iPoint ][ iDim] ) / 2;
donor_jMidEdge_point[iDim] = ( DonorPoint_Coord[ donor_backward_point ][ iDim] +
DonorPoint_Coord[ donor_iPoint ][ iDim] ) / 2;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
donor_iMidEdge_point[iDim] = ( DonorPoint_Coord[ donor_forward_point ][ iDim] +
DonorPoint_Coord[ donor_iPoint ][ iDim] ) / 2;
donor_jMidEdge_point[iDim] = ( DonorPoint_Coord[ donor_backward_point ][ iDim] +
DonorPoint_Coord[ donor_iPoint ][ iDim] ) / 2;
donor_iMidEdge_point[iDim] = ( DonorPoint_Coord[donor_forward_point][iDim] +
DonorPoint_Coord[donor_iPoint][iDim] ) / 2;
donor_jMidEdge_point[iDim] = ( DonorPoint_Coord[donor_backward_point][iDim] +
DonorPoint_Coord[donor_iPoint][iDim] ) / 2;

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any reason to not use zero-spacing (like suggested) for []-access? I do understand how thats helpful for nested accesses [[][]]. But there is no consistency with the additions here [ iDim] [ iDim ] element[0][iDim] = coord[centralNode ][ iDim]. plz consolidate

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The right thing there would be operator ( , ) because Max replaced pointers with containers.
If we just want to clean the spaces then run clang-format and call it a day.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I forgot to replace by operator(i,j) here, fixed in 475c697

Comment on lines 201 to 203
iPoint = geom->vertex[val_marker][iVertex]->GetNode();
if (geom->nodes->GetDomain(iPoint)) {
unsigned long iLocalVertex = iVertex_to_iLocalVertex[iVertex];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
iPoint = geom->vertex[val_marker][iVertex]->GetNode();
if (geom->nodes->GetDomain(iPoint)) {
unsigned long iLocalVertex = iVertex_to_iLocalVertex[iVertex];
const auto iPoint = geom->vertex[val_marker][iVertex]->GetNode();
if (geom->nodes->GetDomain(iPoint)) {
const auto iLocalVertex = iVertex_to_iLocalVertex[iVertex];

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not sure if the const for the iPoint is possible but that would be prefered.

This approach would be appreciated in the rest the blocks as well not just here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I replaced unsigned long by const auto at some locations


unsigned long iVertex, jVertex, kVertex;
unsigned long iVertex, kVertex;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

😞 not very gut

I always prefer to keep those loop variables local to their scope, independent of the number of occurrences in the function. This locality makes it easier imo to modify things because block are somewhat more self-contained and one can be sure that no-one does some shenanigans with break where the value has to leave its loop scope. I also like to do it with stuff like iPoint (see below) for the same reasons. I did not get hit by @pcarruscag with the performance-club yet so I expect that the compiler does efficient things with memory (de)allocation

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here i mean of course not the change you did, but keeping i,kVertex

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Those are either on the stack or probably don't even leave the registers.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I completely agree with @TobiKattmann that these big declaration blocks in the beginning of a function are not very meaningful to the reader. As according to @pcarruscag there is no performance concern, I improved this in CInterpolator::ReconstructBoundary and CSlidingMesh::Build_3D_surface_element.

Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
/*--- Get the number of domain vertices on the marker, and a mapping
* (iVertex) -> (iLocalVertex, non-domain points being ignored). ---*/
su2vector<unsigned long> iVertex_to_iLocalVertex(nVertex);
nLocalVertex = 0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

declare it here if it is needed here for the first time

Common/src/interface_interpolation/CInterpolator.cpp Outdated Show resolved Hide resolved
@pcarruscag
Copy link
Member

Hybrid AD tests are a bit flaky ATM 🤷 you can merge things without those tests passing, but if more than 1 fails let me know.

Update Common/src/interface_interpolation/CSlidingMesh.cpp

Co-authored-by: TobiKattmann <[email protected]>
maxaehle and others added 7 commits August 16, 2021 07:59
Update Common/src/interface_interpolation/CSlidingMesh.cpp

Co-authored-by: TobiKattmann <[email protected]>
Update Common/src/interface_interpolation/CInterpolator.cpp

Co-authored-by: TobiKattmann <[email protected]>
nNodes++;
/*--- Define the neighbors map. ---*/
for (unsigned long iElem = 0; iElem < nElems; iElem++) {
CPrimalGrid* elem = geom->bound[val_marker][iElem];
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This cannot be made const CPrimalGrid* elem at the moment because CPrimalGrid::GetNode etc are not marked as const.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess/assume you could make the various GetNode functions const :) ... and then const auto elem ... but this goes a little beyond PR scope so up to you

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

target_jMidEdge_point[iDim] = ( TargetPoint_Coord[ target_segment[1] ][ iDim ] +
target_geometry->nodes->GetCoord( target_iPoint , iDim) ) / 2;
target_iMidEdge_point[iDim] = ( TargetPoint_Coord(target_segment[0], iDim ) +
target_geometry->nodes->GetCoord( target_iPoint , iDim) ) / 2.;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am curious: Did you get a compiler notification or so on this 2. or was it just personal preference?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Personal preference. Once I spent a lot of time finding out that 1/2==0.

@maxaehle
Copy link
Contributor Author

That's really strange: The hybrid AD testcase discadj_incomp_turb_NACA0012_sst was ok after commit c68cf13, then I changed the stored residuals of unrelated testcases in 9d8f194 and the testcase fails... I'll run clang-format now, let's see what happens ;)

I had a similar thing over at the flamelet branch and a re-run of the regression tests solved that :)

Seems to be like that here as well:
After e95cdfa, only transonic_stator in hybrid_regression_AD.py failed.
After 475c697, only discadj_incomp_turb_NACA0012_sa in hybrid_regression_AD.py failed.
After b704f9e no test in hybrid_regression_AD.py failed.

Update Common/include/interface_interpolation/CInterpolator.hpp

Co-authored-by: Pedro Gomes <[email protected]>
@maxaehle
Copy link
Contributor Author

After 24d81f5, only ea_naca64206 in parallel_regression_AD.py failed

@pcarruscag pcarruscag merged commit 1a84301 into develop Aug 16, 2021
@pcarruscag pcarruscag deleted the fix_slidingmesh_reconstructboundary_pr2 branch August 16, 2021 13:45
@maxaehle maxaehle mentioned this pull request Aug 16, 2021
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants