Allow BC_Sym_Plane to be used as Slip wall for viscous flow by TobiKattmann · Pull Request #740 · su2code/SU2

TobiKattmann · 2019-07-23T09:47:10Z

Hi all,

Proposed Changes

BC_Sym_Plane was written with only straight lines or planes as symmetry boundaries in mind. Therefore a constant unit normal was used. If symmetry as a slip wall in viscous flows is used the boundary in MARKER_SYM can have some curvature which results in a non-constant unit normal. This is fixed in this PR. Additionally viscous terms are removed from Euler wall and MARKER_EULER can not be used in viscous flows.

Affected Regression tests:
a lot ... see a list in the Conversation below

Related Work

This PR is meant to resolve issue #735 opened by @EduardoMolina which was introduced by PR #657 . More details are in the issue

PR Checklist

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).
My contribution is commented and consistent with SU2 style.
I have added a test case that demonstrates my contribution, if necessary.

… flow. Removed viscous terms from Euler wall.

viscous flow.

viscous flow anymore.

…iny bit on symmetry marker, thus the change.

EduardoMolina · 2019-07-25T16:41:11Z

Hi @TobiKattmann
Is this PR ready for review?

Eduardo

TobiKattmann · 2019-07-25T21:41:12Z

Hey @EduardoMolina ,
Yes, I guess it would be good if you can test if this actually fixes your validation Testcase. Maybe also adding the case as regression/validation case could be an idea as there is currently none to my knowledge.
MfG, Tobi

talbring · 2019-07-26T13:33:29Z

Maybe instead of adding a new one, we can modify an existing one to check that feature. Or we add a new one and remove an existing one. We have to somehow reduce the number of tests ...

TobiKattmann · 2019-07-26T15:50:18Z

For clarification: you mean removing a Testcase from the active regression tests to get faster turnaround from travis, the Testcase-files can stay probably.
The person who put it in place most likely had some feature which (s)he wanted to test. The path can give the purpose of the case away, but not necessarily as it is purely the files. If you (or someone else) have a Testcase in mind which can go out or can be modified I would be thankful. Choosing myself is rather arbitrary.

Just as an additional idea: I removed README.md from the gitignore list of the Testcases in a branch and my plan is to shortly explain what the cases are intended to test specifically Link (Note: The Readme there is not done at all, but it shows the idea kinda)

pcarruscag

It should be possible to accomplish this while keeping the initial assumptions for true symmetry planes, and not requiring the users to use different names for the same boundary condition depending on flow type.

pcarruscag · 2019-07-26T16:28:51Z

+    /*---                The "Normal" in SU2 is an Area-Normal and is most likely not constant  ---*/
+    /*---                on the symmetry-plane.                                                 ---*/
+    /*--- Edit July 2019: In order to use this BC_Sym_Plane method for slip walls in viscous    ---*/
+    /*---                 flow the unit-normal & tangent computation is moved into the loop over---*/


I would just edit your initial comment, mention the issue number if you want, but the historical context should not be required to understand what any code currently does.

Yep, Going to change that, once the implementation itself is approved

pcarruscag · 2019-08-01T09:56:11Z

@TobiKattmann is the incompressible solver the only one that needs the special implementation?
Why don't you augment the BC_Sym function to take in a parameter telling it whether the normal is constant or not (default true) and then implement BC_EulerWall of the incompressible solver as a call to BC_Sym specifying varying normal?
I think this would save you the preprocessing of measuring flatness and keep the implementation of BC_Sym equal for the comp and incomp solvers, so that we may one day put it in common and remove duplication.

TobiKattmann · 2019-08-01T10:37:57Z

Hi @pcarruscag ,
compressible needs this as well. There is some NICF (Chi_b, Kappa_b,...) stuff in the Euler_Wall implementation which I need some help with to get it right (that's why it is not done yet ). That is also the only compressible-exclusive part that could be a reason to not have one single code for inc+comp. But of course one could do it anyway with an if(compressible) conditional to avoid code duplication :)

Towards input boolean: I thought about that as well, but with that preprocessing step it is impossible to get it wrong imo. I personally agree with BC_Sym always being flat and Euler allowing to be curvy, but as seen in issue #735 this opinion is not undisputed. By paying the price of checking straightness once, one eliminates the use of Marker_Sym (assumes const normal after your proposal) on curvy slip walls (happened in #735 ), or the possibly correct use of Marker_Euler on flat boundaries which would be a lot more expensive than checking once, right.

pcarruscag · 2019-08-01T11:08:23Z

I see... well every commercial code I have used to date complained about curvy symmetries, we would need to ask Euclid but I think symmetries are flat too.
At any rate, please no if(compressible) :). And on the issue of the compressible implementation couldn't CNSSolver::BC_Euler_Wall be implemented as:

{
    ...
    CEulerSolver::BC_Euler_Wall(...);
    "viscous pixie dust"
}

?

pcarruscag

Question: What happens in moving mesh cases if the surface starts out flat and then gets deformed?

pcarruscag · 2019-08-02T09:39:01Z

 	CVertex*** vertex;		/*!< \brief Boundary Vertex vector (dual grid information). */
  CTurboVertex**** turbovertex; /*!< \brief Boundary Vertex vector ordered for turbomachinery calculation(dual grid information). */
  unsigned long *nVertex;	/*!< \brief Number of vertex for each marker. */
+  bool *bound_is_straight;	/*!< \brief Bool if boundary-marker is straight(2D)/plane(3D) for each marker. */


Can I interest you in std::vector<bool>? It is optimum w.r.t. size and you no longer need to worry about memory management.

pcarruscag · 2019-08-02T09:42:18Z

+          if(RefUnitNormal_defined) {
+            for (iDim = 0; iDim < nDim; iDim++) {
+              if( abs(RefUnitNormal[iDim] - UnitNormal[iDim]) > epsilon )
+                bound_is_straight[iMarker] = false; break;


Are you missing curly brackets here?

Yes! thanks

pcarruscag · 2019-08-02T09:43:46Z

+      } //while iVertex
+
+    /*--- Print results on screen. ---*/
+    if(rank == MASTER_NODE && print_on_screen) {


Is this routine only run by the master node? If not, I think you need a reduction over all partitions that may have gotten a piece of the marker.

Yeah this one was a bit more involved because of the good'ol local vs global marker stuff.

In this previous version each process had its own bound_is_straight array and accessed that in the BC_sym routine -> therefore the computation was correct i think. It would have happened that on the Master-process the value is true and on another false for the same global marker (e.g. picture the 2D wedge cases split in the middle), but as only the master-process prints its value to screen ...the info on screen would be wrong.

I changed it now: there is an Allreduce in the end such that each process has the same globally correct info. I acknowledge that it is prob not the most elegant implementation, so if there is a hint on how to do it better...

MPI wise looks ok, idk if there is a better way in SU2 to get a iMarker to iMarkerGlobal correspondence than comparing strings though.
Maybe let each rank decide if its patch should be straight or not and communicate only the messages you want to print to master and don't do the reduction.
At least convert back from global to local in the preprocessing to avoid comparing strings inside the actual boundary condition.

Reverted implementation back to local vector<bool> bound_is_straight . But I kept a global version in the geometry->ComputeSurf_Straightness to print the result on screen.

pcarruscag · 2019-08-02T10:02:30Z

+          }
+          case 3: {
+            /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/
+            unsigned short Largest, Arbitrary, Zero;


I think the method in the top voted answer should allow you to do this with less logic.
https://math.stackexchange.com/questions/137362/how-to-find-perpendicular-vector-to-another-vector

Well, 1. It has to be a unit normal, therefore the division by sqrt(...) and 2. this Largest, Arbitratry, Zero stuff (which I understand bothers you) is there to avoid having the unit vector build upon numerical zeros which I found to give wrong results. Example: normal is n=( 1e-13, 1e-13, 1) then I could just take t=(n2, -n1, 0) (which is in fact the way v_1 from that top voted answer is build). Then t is no unit vector so I divide each entry by the norm of t which ends up with inaccurate tangentials. That problem occurred in my tests for #657 . I hope that makes it clearer why it is like that. Search the largest entry of the vector and do that switch (like here first and second value t=(n2, -n1, 0) with that largest value

Or did I misinterpret your suggestion

Yeah I got the purpose, but he does a linear combination of two trivially perpendicular vectors.
Which, if the coefficients of the combination are arbitrary, can still give you a small vector.
So one has to choose the coefficients that maximize the magnitude of the potentially zero element (that's one "if").
Worst case scenario this keeps the largest and smallest element in the tangent vector, but if one adds a third trivially perpendicular vector to the result that is avoided:

n = ai + bj +ck if |b| > |c| t = bi + (c-a)j - bk else t = ci - cj + (b-a)k

Which can actually be written without the conditional. But anyway maybe this is not even expensive at all, I just puzzles.

Ok, got it. And I agree that it is the more elegant solution so I took it.

pcarruscag · 2019-08-02T10:04:22Z


-void CTurbSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
+void CTurbSolver::BC_Sym_Plane(CGeometry      *geometry, 
+                               CSolver        **solver_container, 


@talbring should we start using the short container names everywhere now?

TobiKattmann · 2019-08-05T16:48:48Z

Question: What happens in moving mesh cases if the surface starts out flat and then gets deformed?

@pcarruscag I now explicitly stated that: config->GetKind_GridMovement() != RIGID_MOTION ->then recompute normals.

pcarruscag

Thanks @TobiKattmann , I won't nag anymore about implementation details.

pcarruscag · 2019-08-06T07:20:15Z

   */
-  void BC_Euler_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config,
-                     unsigned short val_marker);
+   void BC_Euler_Wall(CGeometry      *geometry, 


more indentation here and there

Went over (hopefully) everything and set correct indentation. Good'ol tabs vs spaces problems in IDE 😐

pcarruscag · 2019-08-06T08:15:21Z

+          }
+          case 3: {
+            /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/
+            unsigned short Largest, Arbitrary, Zero;


Yeah I got the purpose, but he does a linear combination of two trivially perpendicular vectors.
Which, if the coefficients of the combination are arbitrary, can still give you a small vector.
So one has to choose the coefficients that maximize the magnitude of the potentially zero element (that's one "if").
Worst case scenario this keeps the largest and smallest element in the tangent vector, but if one adds a third trivially perpendicular vector to the result that is avoided:

n = ai + bj +ck if |b| > |c| t = bi + (c-a)j - bk else t = ci - cj + (b-a)k

Which can actually be written without the conditional. But anyway maybe this is not even expensive at all, I just puzzles.

pcarruscag · 2019-08-06T08:25:10Z


+    if (iVertex == 0 || 
+        geometry->bound_is_straight[val_marker_Global] != true ||
+        config->GetKind_GridMovement() != RIGID_MOTION) {


Is the default grid movement RIGID_MOTION or NONE? I would bake this check into bound_is_straight during the straightness preprocessing, which lets you return early from that function without doing any computation / MPI and reduces the number of places you need to maintain should (when) the grid motions changes.

Good catch again. As suggested i moved the check into the preprocessing routine.

economon

Thanks, @TobiKattmann, looking good so far. Some comments below

economon · 2019-08-13T06:35:24Z

+      /*---                The "Normal" in SU2 is an Area-Normal and is most likely not constant  ---*/
+      /*---                on the symmetry-plane.                                                 ---*/
+      /*--- Edit July 2019: In order to use this BC_Sym_Plane method for slip walls in viscous    ---*/
+      /*---                 flow the unit-normal & tangent computation is moved into the loop over---*/


Can you please remove the old comment in favor of the new? No need to keep the history, since git takes care of this for us.. might confuse folks more than it helps

economon · 2019-08-13T06:37:29Z

+                                    CConfig        *config,
+                                    unsigned short val_marker) {
+
+  BC_Sym_Plane(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);


👍 nice and easy to maintain

economon · 2019-08-13T06:42:31Z

-  }
-  delete [] Jacobian_b;
-  delete [] DubDu;
+  BC_Sym_Plane(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker);


So, I understand that this is the only change needing some additional attention to make sure that the non-ideal cases are handled for inviscid flows or slip walls in viscous problems, correct? @AlbertoGuardone: if you or anyone in your group has a moment, can you please take a look?

Ciao @economon, ciao @TobiKattmann, we looked into it and it sounds good.
Implementation is correct: it uses a ghost node with an inverted velocity which does not changes the thermodynamic state. The Jacobian is OK as well. Green light from a NICFD perspective.
Thanks!

Thanks Alberto for checking

TobiKattmann · 2019-09-30T15:35:23Z

Alright everyone, this PR is now ready to get merged from my point of view. If there are further concerns regarding particular Testcases please let me know. @pcarruscag did you check airfoil_fsi_rbf and discadj_fsi_airfoil? If yes, your (dis-)approval for this PR would be appreciated (of course from everyone else as well).

I would like to get this integrated before #790 and I would recommend @economon not to start changing reg-test values before this is merged.

The code still fixes the initial Issue #735

economon · 2019-09-30T21:05:18Z

LGTM. Thanks for the fix and updating all regressions, @TobiKattmann

Final question: in the end, the Euler and symmetry BCs are identical implementations, so do we have a practical guideline for when to use one or the other (or some error check), or will we just carry both and allow them to be used interchangeably?

TobiKattmann · 2019-10-01T07:00:43Z

@economon Right now, Euler & Sym can be used interchangeably in both (incomp. & comp.) solvers. One could safely remove one but keeping both is reasonable I guess as both BC are expected by users and maybe the implementation deviates in the future.
Both will give the same results, independent of the 'curvyness' of the boundary (as that is checked in a preprocessing step).

Bonus info: Comp & Incomp implementations are identical, so if a higher-level FlowSolver is implemented, EulerWall and SymBC could move up there.

talbring · 2019-10-01T16:06:54Z

Thanks @TobiKattmann for that fix and going through the hassle of changing all the reg. test values.

pcarruscag · 2019-10-01T18:28:04Z

@TobiKattmann,
In the Airfoil_2D case, why are the surfaces detected as straight? Also the messages are a bit strange:
Boundary marker leading_edge is a single straight.
Boundary marker 0.0 is a single straight.
Boundary marker pressure_side is a single straight.
Boundary marker 0.0 is a single straight.
Boundary marker suction_side is a single straight.
Boundary marker 0.0 is a single straight.

TobiKattmann · 2019-10-01T20:00:40Z

Hi @pcarruscag ,
You probably mean airfoil_fsi_rbf from serial_regression (?)
I get the written on screen:

Boundary marker leading_edge is NOT a single straight.
Boundary marker 0.0 is a single straight.
Boundary marker pressure_side is NOT a single straight.
Boundary marker 0.0 is a single straight.
Boundary marker suction_side is NOT a single straight.
Boundary marker 0.0 is a single straight.

The diff of my local branch with (origin/)develop is empty, so can you double check for me if you are using the latest develop or fix_SlipWall branch?

For the 0.0 markers: In the configFlow.cfg you have
MARKER_EULER= ( leading_edge, 0.0, pressure_side, 0.0, suction_side, 0.0)
I dont understand why there are zeros as nothing has to be prescribed on the Euler Wall. So the 0.0 is interpreted as a marker tag. Removing the 0.0 does not change anything for me in Residuals as well. In the mesh there are also no boundaries defined which are called 0.0 :)

Does that help or am I on the wrong track?

TobiKattmann · 2019-10-01T20:06:44Z

Maybe also an explanation why I specifically state single straight:
Take a case with two symmetry planes on either side of a channel -> it could be reasonable to put both in the same Marker in the su2 mesh -> both planes are straight for themselves but as I simply loop over all nodes in a marker I then have 2 different unit normals for the same marker -> thats why the specific single is used

pcarruscag · 2019-10-01T21:40:45Z

I mean the case in disc_adj_fsi/airfoil_2d for which you updated a filediff regression.

You are probably right on the 0's, those were RANS cases...

TobiKattmann · 2019-10-02T07:46:08Z

For the fluid zone0 the boundaries are detected as NOT straight (captial letters because it is written on screen like that ;) ) and for the solid zone1 the fluid markers are still in the Global marker list but not in the local list... thats why they are written as straight on screen (explanation in next paragraph)

In parallel cases this case (process knows the global marker but does not own any nodes of it i.e. it is also not in the local marker list) the process has to assume that the boundary is straight for the mpi communication otherwise we would never have straight boundary-straightness predictions in high process-count simulations. Here the situation is the same... except for it is a serial case

I have to admit I am not really sure how to best tackle that. Can you split the cfg files like in the Airfoil_RBF case (not sure if this is an old vs new driver thing)? I could exclude that structral solver from the surf_straightness check in CDriver::Geometrical_Preprocessing but I am not sure if that possible and if, what boolean to use.

pcarruscag · 2019-10-02T09:29:11Z

Is the solver doing the right thing then? That is all I care about.

That case cannot be split for now, but since that driver will be retired in the near future I would say you don't need to worry (so long as the solver is doing the right thing).

TobiKattmann · 2019-10-02T14:33:40Z

As each zone has its own geometry_container in which the bound_is_straight data is stored, the false data for the solid zone should not interfere with the fluid zone. And as the fluid zone contains the correct straightness info everything should be fine (except for the screen output of course).

BC_Sym_Plane was written with only straight lines or planes as symmetry boundaries in mind. Therefore a constant unit normal was used. If symmetry as a slip wall in viscous flows is used the boundary in MARKER_SYM can have some curvature which results in a non-constant unit normal. This is fixed in this PR. + MARKER_SYM now allows non-straight walls with local symmetry enforced + For viscous flows, MARKER_EULER now accounts for viscosity + MARKER_EULER and MARKER_SYM (symmetry planes and Euler walls) now share their implementation Right now, Euler & Sym can be used interchangeably in both (incomp. & comp.) solvers. One could safely remove one but keeping both is reasonable as both BC are expected by users and maybe the implementation deviates in the future. Both will give the same results, independent of the 'curvyness' of the boundary (as that is checked in a preprocessing step). See the following issues/PRs for more details: su2code#735 su2code#740

TobiKattmann added 5 commits July 23, 2019 11:24

See Issue 735. Allow BC_Sym_Plane to be used as slip wall for viscous…

1816c3b

… flow. Removed viscous terms from Euler wall.

Added information in config_template to use MARKER_SYM for slip walls in

18905e7

viscous flow.

Change Regression test values. MARKER_EULER not used for slip wall in

1436847

viscous flow anymore.

Remove Unused Variables to not generate more warnings.

060797e

Change Regressiontest value of euler ONERAM6. Unit Normal changes a t…

b7ee8d5

…iny bit on symmetry marker, thus the change.

TobiKattmann marked this pull request as ready for review July 25, 2019 21:41

pcarruscag requested changes Jul 26, 2019

View reviewed changes

TobiKattmann added 3 commits July 31, 2019 13:41

Merged EulerWall and SymPlane into SlipWall.

2af4024

Merge remote-tracking branch 'origin/develop' into fix_SlipWall

779d353

Finalize Marker slip-wall for incompressible flow.

dc4cd06

Merge remote-tracking branch 'origin/develop' into fix_SlipWall

6bfe577

TobiKattmann and others added 2 commits August 1, 2019 15:44

BC_Sym_Plane used for Euler_wall in inc+comp flow.

49f3c1c

Revert unnecessary changes.

1ca346c

pcarruscag reviewed Aug 2, 2019

View reviewed changes

TobiKattmann added 2 commits August 5, 2019 17:29

Merge remote-tracking branch 'origin/develop' into fix_SlipWall

3e0283e

Added parallel comm for bound_is_straight boolean vector.

f24eb4b

pcarruscag approved these changes Aug 6, 2019

View reviewed changes

economon reviewed Aug 13, 2019

View reviewed changes

TobiKattmann added 4 commits August 13, 2019 13:49

Req Changes for PR 740

91be1f6

Merge remote-tracking branch 'origin/develop' into fix_SlipWall

07d1f21

Merge remote-tracking branch 'origin/develop' into fix_SlipWall

c1a2d8a

FIxes for PR 740.

74902e0

TobiKattmann added 5 commits September 30, 2019 10:07

PR740 adapting antother chunk of Reg.test values.

1f1fa29

PR740 Reg.test filediff changes.

c85e740

Fix little error in serial_reg_AD.py.

1de4655

PR740 Reg.test adaption. Wrong end-of-line of filediffs due to windows.

114e21b

Revert .travis.yml

db24032

TobiKattmann changed the title ~~WIP Allow BC_Sym_Plane to be used as Slip wall for viscous flow~~ Allow BC_Sym_Plane to be used as Slip wall for viscous flow Sep 30, 2019

pr-triage Bot added the PR: unreviewed label Sep 30, 2019

TobiKattmann mentioned this pull request Sep 30, 2019

Possible bug in develop compared to master branch. #735

Closed

WallyMaier added PR: reviewed-approved and removed PR: unreviewed labels Sep 30, 2019

economon approved these changes Sep 30, 2019

View reviewed changes

talbring approved these changes Oct 1, 2019

View reviewed changes

talbring merged commit 5241613 into develop Oct 1, 2019

talbring deleted the fix_SlipWall branch October 1, 2019 16:07

pr-triage Bot added PR: merged and removed PR: reviewed-approved labels Oct 1, 2019

talbring added the changelog:fix label Nov 6, 2019

TobiKattmann mentioned this pull request Dec 6, 2019

Different results on the same testcase of unsteady ONERA (v6.0.1 vs v6.2.0) #828

Closed

Uh oh!

Conversation

TobiKattmann commented Jul 23, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Related Work

PR Checklist

Uh oh!

EduardoMolina commented Jul 25, 2019

Uh oh!

TobiKattmann commented Jul 25, 2019

Uh oh!

talbring commented Jul 26, 2019

Uh oh!

TobiKattmann commented Jul 26, 2019

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pcarruscag commented Aug 1, 2019

Uh oh!

TobiKattmann commented Aug 1, 2019

Uh oh!

pcarruscag commented Aug 1, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

TobiKattmann commented Aug 5, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

TobiKattmann commented Jul 23, 2019 •

edited

Loading

pcarruscag commented Aug 1, 2019 •

edited

Loading

TobiKattmann commented Aug 5, 2019 •

edited

Loading

TobiKattmann Aug 21, 2019 •

edited

Loading

TobiKattmann commented Oct 1, 2019 •

edited

Loading

TobiKattmann commented Oct 2, 2019 •

edited

Loading