[WIP] Wall Model LES #1120

EduardoMolina · 2020-12-03T02:24:24Z

Proposed Changes

This PR adds Wall Model LES (and also Wall Resolved LES) capabilities to SU2. Several wall models were included, including the Equilibrium Algebraic and Adverse Pressure Gradient Wall Models. The Vreman and WALE with constant model coefficients are used to model the subgrid-scales but others SGS models can be easily added.
It is worth noting that the exchange location interpolation for the wall model inputs is not included here (see feature_WallModelLES branch). Thus, the wall model will use the information of the 1st grid point off the wall.

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.

[X ] I am submitting my contribution to the develop branch.
[X ] 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).
[X ] My contribution is commented and consistent with SU2 style.
[X ] 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.

vdweide · 2020-12-03T06:46:43Z

@EduardoMolina, you mentioned that the wall model will use the information of the 1st grid point off the wall and it is not possible to specify the exchange location manually anymore in this PR. I realize the latter requires a more complicated search for the donor point, but wouldn't it be worthwhile to have this option, especially if you do not have a structured grid direction normal to the wall?

pcarruscag

Thank you for bringing this to develop Eduardo!
Even if this is still WIP I have a few comments, judging from parts of the code you were probably planning to address some of them, in any case here they are:

pcarruscag · 2020-12-03T10:32:26Z

SU2_CFD/include/sgs_model.hpp

+ * \class CSIGMAModel
+ * \brief Derived class for defining the SIGMA SGS model.
+ * \author: E. van der Weide, T. Economon, P. Urbanczyk
+ * \version 7.0.3 "Blackbird"


Update the version, or remove this line, these "local" version numbers are always out of date.

pcarruscag · 2020-12-03T10:33:15Z

SU2_CFD/include/sgs_model.hpp

+                                   const su2double drhodx,
+                                   const su2double drhody,
+                                   const su2double dudx,
+                                   const su2double dudy,
+                                   const su2double dvdx,
+                                   const su2double dvdy,
+                                   const su2double d2udx2,


These derivatives could be passed in a more elegant way.

Can you give an idea what you have in mind? Using a struct maybe or some other structure that bundles this together?

the math thing not the movie :)

pcarruscag · 2020-12-03T10:40:55Z

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

+  if (TauWallFlag_i && TauWallFlag_j){
+     UseWallFunction = false;
+     Mean_TauWall = 0.0;
+   }
+   else if (TauWallFlag_i && !TauWallFlag_j){
+     UseWallFunction = true;
+     Mean_TauWall = TauWall_i;
+   }
+   else if (!TauWallFlag_i && TauWallFlag_j){
+     UseWallFunction = true;
+     Mean_TauWall = TauWall_j;
+   }
+   else{
+     UseWallFunction = false;
+     Mean_TauWall = 0.0;
+   }


UseWallFunction = TauWallFlag_i ^ TauWallFlag_j; if (UseWallFunction) Mean_TauWall = TauWallFlag_i? TauWall_i : TauWall_j;

Nice :)
Maybe add a Mean_TauWall = 0.0 above or as an else... not soure what the value is otherwise

pcarruscag · 2020-12-03T10:46:39Z

SU2_CFD/src/solvers/CNSSolver.cpp

+  su2double *Normal     = new su2double[nDim];
+  su2double *UnitNormal = new su2double[nDim];


Use static arrays of MAXNDIM size

pcarruscag · 2020-12-03T10:48:08Z

SU2_CFD/src/solvers/CNSSolver.cpp

+      Area = 0.0;
+      for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim];
+      Area = sqrt (Area);


Use the Geometry toolbox please

pcarruscag · 2020-12-03T10:51:38Z

SU2_CFD/src/solvers/CNSSolver.cpp

+          su2double *Coord_i = geometry->nodes->GetCoord(iPoint);
+          su2double *Coord_j = geometry->nodes->GetCoord(Point_Normal);
+          su2double dist_ij = 0;
+          for (iDim = 0; iDim < nDim; iDim++)
+            dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]);
+          dist_ij = sqrt(dist_ij);


const su2double* when possible and geometry toolbox for the distance too.

pcarruscag · 2020-12-03T10:53:11Z

SU2_CFD/src/solvers/CNSSolver.cpp

+              VelNormal = 0.0;
+              for (auto iDim = 0u; iDim < nDim; iDim++)
+                VelNormal += Vel[iDim] * UnitNormal[iDim];
+              for (auto iDim = 0u; iDim < nDim; iDim++)
+                VelTang[iDim] = Vel[iDim] - VelNormal*UnitNormal[iDim];
+
+              VelTangMod = 0.0;
+              for (auto iDim = 0u; iDim < nDim; iDim++)
+                VelTangMod += VelTang[iDim]*VelTang[iDim];
+              VelTangMod = sqrt(VelTangMod);
+
+              /*--- Compute normal distance of the interior point from the wall ---*/
+
+              for (auto iDim = 0u; iDim < nDim; iDim++)
+                WallDist[iDim] = (Coord[iDim] - Coord_Normal[iDim]);
+
+              WallDistMod = 0.0;
+              for (auto iDim = 0u; iDim < nDim; iDim++)
+                WallDistMod += WallDist[iDim]*WallDist[iDim];
+              WallDistMod = sqrt(WallDistMod);


More cases where the toolbox simplifies things to one-liners

pcarruscag · 2020-12-03T10:54:41Z

TestCases/wmles/channel/wmles_channel.cfg

+% SU2 configuration file                                                       %
+% Case description: Plane Channel Flow at Re_tau = 5200                         %
+% Author: Eduardo Molina                                                       %


Thank you for adding an example, I would make it a regression test too.

pcarruscag · 2020-12-03T10:57:29Z

SU2_CFD/include/solvers/CSolver.hpp

+  inline virtual su2double GetTauWall_WMLES(unsigned short val_marker, unsigned long val_vertex) const { return 0;}
+
+  /*!
+   * \brief A virtual member
+   * \param[in] val_marker - Surface marker where the coefficient is computed.
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the coefficient is evaluated.
+   * \return Value of the heat flux from the wall model.
+   */
+  inline virtual su2double GetHeatFlux_WMLES(unsigned short val_marker, unsigned long val_vertex) const { return 0;}
+


Are these methods really needed?
I think a lot of these exist due to flow solver - continuous adjoint solver interaction, but if these variables are only used within the NS solver, all these virtual functions can be avoided.

pcarruscag · 2020-12-03T11:02:36Z

SU2_CFD/include/solvers/CSolver.hpp

+   inline virtual void SetEddyViscFirstPoint(CGeometry *geometry,
+                              CSolver** solver_container,
+                              CConfig* config) { }
+
+   /*!
+    * \brief A virtual member.
+    * \param[in] geometry - Geometrical definition of the problem.
+    * \param[in] solver_container - Container vector with all the solutions.
+    * \param[in] config - Definition of the particular problem.
+    */
+   inline virtual void Setmut_LES(CGeometry *geometry,


Similarly for these methods, if they are part of the preprocessing of the NS solver, they should be private methods of that solver.
The abstraction (i.e. the virtual function) we want in CSolver is "Preprocessing", everything that is physics-specific should concern only the derived classes (Euler, NS, etc.) and fall in the abstract categories of Preprocessing, postprocessing, iteration, etc. which are the only hooks we should use to augment the functionality of an abstract solver.

EduardoMolina · 2020-12-03T22:28:01Z

Hi @vdweide,
Thanks for the comment.
I am using the "geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor()" function. Thus, there isn't any search/interpolation of the donor points. I think we can add all the complicated search/interpolation stuff in a single pull request later.

TobiKattmann

Thanks @EduardoMolina for this addition
Although I don't understand much of the implementation I have some general comments towards the code. I like that you thoroughly commented a lot of it 👍

TobiKattmann · 2020-12-04T13:07:33Z

Common/src/wall_model.cpp

+  if (Temperature_Prescribed){
+    /* The Kader's law will be used to approximate the variations of the temperature inside the boundary layer.
+     */
+    const su2double y_plus = u_tau*h_wm/nu_wall;
+    const su2double lhs = - ((tExchange - TWall) * rho_wall * Cp * u_tau);
+    const su2double Gamma = - (0.01 * (pow(Pr_lam * y_plus,4.0))/(1.0 + 5.0*y_plus*pow(Pr_lam,3.0)));
+    const su2double rhs_1 = Pr_lam * y_plus * exp(Gamma);
+    const su2double rhs_2 = (2.12*log(1.0+y_plus) + pow((3.85*pow(Pr_lam,(1.0/3.0)) - 1.3),2.0) + 2.12*log(Pr_lam)) * exp(1./Gamma);
+    qWall = lhs/(rhs_1 + rhs_2);


I notice some code duplication of this and a bit more in 3 more places below. Is it possible to easily extract the common and specific parts to avoid that?

TobiKattmann · 2020-12-04T13:12:19Z

SU2_CFD/include/numerics/CNumerics.hpp

@@ -1307,6 +1307,13 @@ class CNumerics {
   */
  inline virtual void SetTauWall(su2double val_tauwall_i, su2double val_tauwall_j) { }

+  /*!
+   * \brief Set the value of the bollean flag to use (or not) the wall shear stress from the wall function.


TobiKattmann · 2020-12-04T13:17:51Z

SU2_CFD/include/sgs_model.hpp

+                                   const su2double drhodx,
+                                   const su2double drhody,
+                                   const su2double dudx,
+                                   const su2double dudy,
+                                   const su2double dvdx,
+                                   const su2double dvdy,
+                                   const su2double d2udx2,


Can you give an idea what you have in mind? Using a struct maybe or some other structure that bundles this together?

TobiKattmann · 2020-12-04T13:24:19Z

SU2_CFD/include/sgs_model.inl

+
+  if ( alpha1 < 0.0 ) alpha1 = abs(alpha1); // avoid sqrt of negative value
+


Do you maybe want to throw a personalized error in case of alpha1 being negative? Or have some other form of error handling

For performance avoid if statements like the plague, the code above is equivalent without the "if".

TobiKattmann · 2020-12-04T13:26:21Z

SU2_CFD/include/sgs_model.inl

+  cout << "CSIGMAModel::ComputeGradEddyViscosity_2D: Not implemented yet" << endl;
+  exit(1);


maybe use SU2_MPI::Error here

Always use SU2_MPI::Error

TobiKattmann · 2020-12-04T14:00:06Z

SU2_CFD/src/numerics/flow/flow_diffusion.cpp

+  if (TauWallFlag_i && TauWallFlag_j){
+     UseWallFunction = false;
+     Mean_TauWall = 0.0;
+   }
+   else if (TauWallFlag_i && !TauWallFlag_j){
+     UseWallFunction = true;
+     Mean_TauWall = TauWall_i;
+   }
+   else if (!TauWallFlag_i && TauWallFlag_j){
+     UseWallFunction = true;
+     Mean_TauWall = TauWall_j;
+   }
+   else{
+     UseWallFunction = false;
+     Mean_TauWall = 0.0;
+   }


Nice :)
Maybe add a Mean_TauWall = 0.0 above or as an else... not soure what the value is otherwise

TobiKattmann · 2020-12-04T14:14:38Z

SU2_CFD/src/solvers/CNSSolver.cpp

+    Alloc3D(nMarker, nVertex, nDim, VelTimeFilter_WMLES);
+
+    /*--- Check if the Wall models or Wall functions are unique. ---*/
+    /*--- OBS: All the markers must have the same wall model/function ---*/


TobiKattmann · 2020-12-04T14:21:48Z

SU2_CFD/src/solvers/CNSSolver.cpp

+  /*--- Compute the recovery factor ---*/
+  // su2double-check: laminar or turbulent Pr for this?
+  su2double Recovery = pow(config->GetPrandtl_Lam(),(1.0/3.0));


... when double-checking is not enough ;)

TobiKattmann · 2020-12-04T14:27:37Z

SU2_CFD/src/solvers/CNSSolver.cpp

+                cout << "WARNING: Y+ is too high (>1e4). The muT computation at the 1st point off the wall is disable." << endl;
+                cout << rank << " " << iPoint;
+//                for (auto iDim = 0u; iDim < nDim; iDim++)
+//                  cout << " " << Coord[iDim];
+                cout << endl;


disable'd' in the sentence
Add "Rank: " and "iPoint: " otherwise it will be quite confusing, or integrate it in sentence above
Remove commented debug output.

TobiKattmann · 2020-12-04T14:32:06Z

Common/include/option_structure.hpp

-  LOGARITHMIC_WALL_MODEL    = 6    /*!< \brief Logarithmic law-of-the-wall model for LES. */
+  EQUILIBRIUM_WALL_MODEL    = 2,   /*!< \brief Equilibrium wall model for LES. */
+  LOGARITHMIC_WALL_MODEL    = 3,   /*!< \brief Reichardt's law-of-the-wall model for LES. */
+  ALGEBRAIC_WALL_MODEL      = 4,    /*!< \brief Algebraic wall model for LES. */
+  APGLL_WALL_MODEL          = 5,     /*!< \brief Adverse Pressure Gradient Wall Model for LES. */
+  TEMPLATE_WALL_MODEL       = 6     /*!< \brief Template Wall Model */


Can you please add/update the config_template.cfg with the changes introduced in this PR.

If the normal point lies on the surface, skip the WMLES calculation.

stale · 2021-02-13T07:20:59Z

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. If this is still a relevant issue please comment on it to restart the discussion. Thank you for your contributions.

pcarruscag · 2021-02-13T12:09:55Z

Hi @EduardoMolina,
I tried to fix the conflicts since all of them were caused by my recent cleanups and fixes of the flow solvers. But I very likely broke compilation as I did it trough github...

I think this is an extremely valuable feature to have in SU2, more so for putting in place what is needed to finally have RANS wall functions.
At the same time it is hard to keep the stale bot away from these PR's that don't see any action for a long time. And we cannot stop doing work in other areas that need attention.
So, what can I do to help get this into develop? Even if you are still working on some part of the implementation, or tuning some aspects of the methods, we can always merge it and warn users about something that might be experimental or hide for expert-use only.
But it really is a lot easier to maintain code that is already in develop.

EduardoMolina · 2021-02-18T18:00:07Z

Hi @pcarruscag ,

Yes, I am still tweaking parameters while we are participating on the AIAA-HLPW. But I agree with you that we need to merge asap even if we do an update after the workshop. Let me fix the compilation and add a regression test case. I probably will need your help with the AD part.

Thanks.
Eduardo

bigfooted · 2021-02-21T08:21:36Z

Good to hear about the progress. In the mean time, I have had a look at the wall model implementation for RANS in your implementation and the current develop. I have a working version validated for SA and SST in the branch 'fix_wallfunctions', which was branched from develop.

… vectorized implementation

stale · 2021-06-02T16:07:51Z

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. If this is still a relevant issue please comment on it to restart the discussion. Thank you for your contributions.

* fixed wall functions for SA and SST * fix unused variables update * fix wall model constant config option * computing skinfriction when y+ < 5 * fix OMP command * update to enum class * Apply suggestions from code review * a little more refactoring * apply fix from #1120 (do not consider tau wall for edges along wall), vectorized implementation * some fixes and better warnings Co-authored-by: TobiKattmann <31306376+TobiKattmann@users.noreply.github.com> Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Co-authored-by: Pedro Gomes <pcarruscag@gmail.com>

pcarruscag · 2022-04-13T14:19:06Z

Re-opened because @GomerOfDoom is interested in bringing this up to speed with develop.

pcarruscag · 2022-05-02T11:28:34Z

Closing again based on #1622. @GomerOfDoom please have a look at the comments Tobi and I left here, might help with the re-implementation.

EduardoMolina added 12 commits November 24, 2020 09:02

Xcode update

5a2a44a

Start merging WMLES (not working)

9a822dc

Finish merging WMLES (need testing)

49c5931

Update CNSSolver

8c708cc

Update missing declarations (ready for test)

31c012e

Update .travis.yml

942cdfc

Add WMLES output

de689da

Update CNSSolver

0eada7a

Fix small typo in SetTime_Step

0417f3d

Added wmles test case

b682814

Merge remote-tracking branch 'su2code/develop' into feature_WMLES-Dev

aa55317

Update .travis.yml back to develop

bdea14e

pcarruscag reviewed Dec 3, 2020

View reviewed changes

TobiKattmann reviewed Dec 4, 2020

View reviewed changes

EduardoMolina added 8 commits December 6, 2020 12:32

Trying to fix isnan compilation

6310a6a

Update BC_WallModel

85bcfc2

Minor updates for PR

f9f7124

Add check for normal point in WMLES

554d94d

If the normal point lies on the surface, skip the WMLES calculation.

Add WMLES to incompressible solver (need testing)

f9e60a9

Merge remote-tracking branch 'su2code/develop' into feature_WMLES-Dev

7ef5c9a

SU2_Type for std:isnan

07470c9

Add WMLES test cases

db80002

stale bot added the stale label Feb 13, 2021

Merge branch 'develop' into feature_WMLES-Dev

b3743e9

stale bot removed the stale label Feb 13, 2021

pcarruscag added the changelog:feature label Feb 25, 2021

pcarruscag added this to In progress in RPSVV Working Group via automation Feb 25, 2021

pcarruscag mentioned this pull request May 29, 2021

Fix SA and SST wall functions #1204

Merged

3 tasks

pcarruscag added a commit that referenced this pull request May 30, 2021

apply fix from #1120 (do not consider tau wall for edges along wall),…

fb3c0cc

… vectorized implementation

stale bot added the stale label Jun 2, 2021

stale bot closed this Jun 9, 2021

RPSVV Working Group automation moved this from In progress to Done Jun 9, 2021

pcarruscag mentioned this pull request Jan 23, 2022

MAJOR ISSUES / DEVELOPMENT AREAS #1487

Open

pcarruscag reopened this Apr 13, 2022

RPSVV Working Group automation moved this from Done to In progress Apr 13, 2022

stale bot removed the stale label Apr 13, 2022

pcarruscag closed this May 2, 2022

RPSVV Working Group automation moved this from In progress to Done May 2, 2022

		su2double *Normal = new su2double[nDim];
		su2double *UnitNormal = new su2double[nDim];


		if ( alpha1 < 0.0 ) alpha1 = abs(alpha1); // avoid sqrt of negative value

		cout << "CSIGMAModel::ComputeGradEddyViscosity_2D: Not implemented yet" << endl;
		exit(1);

[WIP] Wall Model LES #1120

[WIP] Wall Model LES #1120

Conversation

EduardoMolina commented Dec 3, 2020

Proposed Changes

PR Checklist

vdweide commented Dec 3, 2020

pcarruscag left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pcarruscag Dec 3, 2020 • edited Loading

Choose a reason for hiding this comment

EduardoMolina commented Dec 3, 2020

TobiKattmann left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stale bot commented Feb 13, 2021

pcarruscag commented Feb 13, 2021

EduardoMolina commented Feb 18, 2021

bigfooted commented Feb 21, 2021

stale bot commented Jun 2, 2021

pcarruscag commented Apr 13, 2022

pcarruscag commented May 2, 2022

pcarruscag Dec 3, 2020 •

edited

Loading