[WIP] PGOMEGA Implementation #2354

TurboLdd · 2024-09-10T14:26:04Z

Proposed Changes

Hi all,
I have implemented the PGOmega S-A turbulence model(TM) which is proposed by Xiao He(2022). This TM enables us to get a more accurate prediction of stall margin which is very important in turbomachinery. This TM has been used in the simulations of ROTOR67, TUDarmstadt single stage compressor, and the validation work is shown in GPPS 2024 conference.
I also added the option related in configuration file.
The paper related is "He, X., Zhao, F., and Vahdati, M. (September 19, 2022). "A Turbo-Oriented Data-Driven Modification to the Spalart–Allmaras Turbulence Model."

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 --warnlevel=3 when using meson).
My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
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.

joshkellyjak · 2024-09-10T15:52:00Z

Do you have a testcase that demonstrates this implementation? Can you reproduce and compare to the validation plots published by Xiao He? I think you have shared some similar results with us before.

joshkellyjak

If this is a turbo-oriented feature I think you should highlight this in the documentation when you update it

SU2_CFD/include/numerics/CNumerics.hpp

joshkellyjak · 2024-09-10T15:52:52Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp


  /*--- List of auxiliary functions ---*/
  su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2;
+


Suggested change

joshkellyjak · 2024-09-10T15:53:17Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+          su2double RefReynold = 1.0e6;
+          dpomg = GeometryToolbox::DotProduct(3, PrimVar_Grad_i[idx.Pressure()], Vorticity_Rel) * laminar_viscosity *
+                  RefReynold / (omg * pow(velocityMag, 3) * pow(density, 2));  // reynold should be modified
+          //cout<<"dpomg="<<dpomg<<endl;


Remove comments

joshkellyjak · 2024-09-10T15:54:28Z

SU2_CFD/include/variables/CTurbSAVariable.hpp

@@ -42,6 +42,7 @@ class CTurbSAVariable final : public CTurbVariable {
  VectorType DES_LengthScale;
  VectorType Vortex_Tilting;

+


Suggested change

joshkellyjak · 2024-09-10T15:54:38Z

SU2_CFD/include/variables/CTurbSAVariable.hpp

@@ -87,4 +88,5 @@ class CTurbSAVariable final : public CTurbVariable {
   */
  inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }

+


Suggested change

joshkellyjak · 2024-09-10T15:54:47Z

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1413,6 +1413,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) {

  if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) {
    AddVolumeOutput("EDDY_VISCOSITY", "Eddy_Viscosity", "PRIMITIVE", "Turbulent eddy viscosity");
+


Suggested change

joshkellyjak · 2024-09-10T15:54:54Z

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1544,6 +1545,7 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
    SetVolumeOutputValue("EDDY_VISCOSITY", iPoint, Node_Flow->GetEddyViscosity(iPoint));
    SetVolumeOutputValue("TURB_DELTA_TIME", iPoint, Node_Turb->GetDelta_Time(iPoint));
    SetVolumeOutputValue("TURB_CFL", iPoint, Node_Turb->GetLocalCFL(iPoint));
+


Suggested change

joshkellyjak · 2024-09-10T16:05:21Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+        velocityMag = max(GeometryToolbox::Norm(3, Velocity_Rel), 0.001);
+        omg = max(GeometryToolbox::Norm(3, Vorticity_Rel), 0.001);
+
+        for (int iDim = 0; iDim < nDim; iDim++) dpds += PrimVar_Grad_i[idx.Pressure()][iDim] * Velocity_Rel[iDim];


Check capitialisation for consistency, should this be dP_ds?

Same with line 162

joshkellyjak · 2024-09-10T16:09:31Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+          Velocity_Rel[iDim] = V_i[idx.Velocity() + iDim] - RotationVelocityCrossR[iDim];
+          Vorticity_Rel[iDim] = Vorticity_i[iDim] - RotationalVelocity[iDim];
+
+        }


If this model only works for 3D flows (I assume you only use it in these cases) you should not allow users to select it with a 2D calculation

joshkellyjak · 2024-09-10T16:13:04Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+        su2double Velocity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double Vorticity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double RotationalVelocity[3] = {0.0, 0.0, 0.0};
+        su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};


Suggested change

su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};

su2double RotationVelocityCrossR[nDim] = {0.0};

or, since you are already doing static allocation:

Suggested change

su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};

su2double RotationVelocityCrossR[MAXNDIM];

and you do not rely on the initial value being 0 as far as I see

TobiKattmann

Hey @TurboLdd , I see this is your first PR, so congrats to that :)

I am with Josh there that it would be great to have the validation plots created with this PR-code, and in addition the setup of that added as a small Test or unit-tested.

Concerning the implementation, from my pov it is fine to first get a working version and then make it work with the existing structure - and for that I think you make it easier on yourself if you have a small test readily set up to iterate on, once you start refactoring.

TobiKattmann · 2024-09-11T07:05:10Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+        su2double Velocity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double Vorticity_Rel[3] = {0.0, 0.0, 0.0};
+        su2double RotationalVelocity[3] = {0.0, 0.0, 0.0};
+        su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};


or, since you are already doing static allocation:

Suggested change

su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};

su2double RotationVelocityCrossR[MAXNDIM];

and you do not rely on the initial value being 0 as far as I see

TobiKattmann · 2024-09-11T07:15:31Z

Common/include/option_structure.hpp

@@ -1113,6 +1113,7 @@ enum class SA_OPTIONS {
  ROT,      /*!< \brief Rotation correction. */
  BC,       /*!< \brief Bas-Cakmakcioclu transition. */
  EXP,      /*!< \brief Allow experimental combinations of options (according to TMR). */
+  PGOMGA,   /*!< \brief Use PGOMGA term. */


Please add the DOI address of your reference somewhere (maybe here), i.e. some https://doi.org/<the DOI string>address

TobiKattmann · 2024-09-11T07:22:58Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -283,7 +323,7 @@ struct Omega {
 struct Bsl {
  template <class MatrixType>
  static void get(const su2double* vorticity, unsigned short, const MatrixType&, CSAVariables& var) {
-    var.Omega = GeometryToolbox::Norm(3, vorticity);
+    var.Omega = GeometryToolbox::Norm(3, vorticity)*(1.0 + var.beta_PGO);


I am not familiar with the PGOMGA model but my guess is that we can better fit it into the current modular structure than it is done here. E.g. under the struct Omega { you could create you own struct PGOMGA { and put down you omega contribution there.

Currently also in the ComputeResidual routine there is specific PGOMGA contributions which would be great to move into specialized PGOMGA contributions following the current structure

TobiKattmann · 2024-09-11T07:32:41Z

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+        for (int iDim = 0; iDim < nDim; iDim++) dpds += PrimVar_Grad_i[idx.Pressure()][iDim] * Velocity_Rel[iDim];
+
+        if (dpds >= 0.0) {
+          su2double RefReynold = 1.0e6;


Suggested change

su2double RefReynold = 1.0e6;

const su2double RefReynold = 1.0e6;

const all the things possible :)

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

+          var.beta_PGO = var.cpw1 * tanh(var.cpw2 * pow(dpomg, var.cpw3));
+        }
+
+        //var.Omega = omg * (1.0 + var.beta_PGO);


TurboLdd added 2 commits September 10, 2024 22:10

PGOMEGA Implementation

2f41cb1

Merge branch 'develop' into PGOMEGA

881e2d5

TurboLdd marked this pull request as draft September 10, 2024 14:35

TurboLdd marked this pull request as ready for review September 10, 2024 14:35

joshkellyjak reviewed Sep 10, 2024

View reviewed changes

TobiKattmann reviewed Sep 11, 2024

View reviewed changes

TobiKattmann added the changelog:feature label Sep 11, 2024

github-advanced-security bot found potential problems Sep 12, 2024

View reviewed changes

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

var.beta_PGO = var.cpw1 * tanh(var.cpw2 * pow(dpomg, var.cpw3));

}

//var.Omega = omg * (1.0 + var.beta_PGO);

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

pcarruscag changed the title ~~PGOMEGA Implementation~~ [WIP] PGOMEGA Implementation Sep 14, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] PGOMEGA Implementation #2354

[WIP] PGOMEGA Implementation #2354

TurboLdd commented Sep 10, 2024

joshkellyjak commented Sep 10, 2024

joshkellyjak left a comment

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

joshkellyjak Sep 10, 2024

TobiKattmann Sep 11, 2024

TobiKattmann left a comment

TobiKattmann Sep 11, 2024

TobiKattmann Sep 11, 2024

TobiKattmann Sep 11, 2024

TobiKattmann Sep 11, 2024


		/--- List of auxiliary functions ---/
		su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2;

		@@ -42,6 +42,7 @@ class CTurbSAVariable final : public CTurbVariable {
		VectorType DES_LengthScale;
		VectorType Vortex_Tilting;

		@@ -87,4 +88,5 @@ class CTurbSAVariable final : public CTurbVariable {
		*/
		inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); }

		@@ -1413,6 +1413,7 @@ void CFlowOutput::SetVolumeOutputFieldsScalarPrimitive(const CConfig* config) {

		if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) {
		AddVolumeOutput("EDDY_VISCOSITY", "Eddy_Viscosity", "PRIMITIVE", "Turbulent eddy viscosity");

	su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};
	su2double RotationVelocityCrossR[nDim] = {0.0};

	su2double RotationVelocityCrossR[3] = {0.0, 0.0, 0.0};
	su2double RotationVelocityCrossR[MAXNDIM];

	su2double RefReynold = 1.0e6;
	const su2double RefReynold = 1.0e6;

[WIP] PGOMEGA Implementation #2354

Are you sure you want to change the base?

[WIP] PGOMEGA Implementation #2354

Conversation

TurboLdd commented Sep 10, 2024

Proposed Changes

PR Checklist

joshkellyjak commented Sep 10, 2024

joshkellyjak 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

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