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

[WIP] Implementation of eN transition model #1765

Draft
wants to merge 95 commits into
base: develop
Choose a base branch
from

Conversation

RichRoos
Copy link

Hello!

Proposed Changes

To continue current work being done on different transition models, I am working at the implementation of the eN model for low-turbulence external aerodynamics as proposed by James G. Goder and Mark D. Maughmer (10.2514/1.J052905). This model incorporates a new transport equation to describe the growth of the maximum Tollmien-Schlichting instability amplitude in the presence of a boundary layer.

Related Work

This work will be set-up in the same way as the LM model is in #1751 (which was also generally speaking how I had set up my own LM transition implementation) to make sure merging in a later stage will be easier.

Roadmap

I am hoping to achieve a working implementation of the eN transition method. This work will be done as finishing part of my master thesis. We will have to see if I have time to completely finish the V&V phase due to time restrictions. I have already done quite some work on implementing and verifying the LM method on a private branch and results show mixed results, especially at higher Re numbers (mentioned also by Menter et all.). Thus the step is made to see how these LM results will compare to eN results.

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 have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

@pr-triage pr-triage bot added the PR: draft label Sep 20, 2022
Comment on lines 6039 to 6042
switch (Kind_Trans_Model) {
case TURB_TRANS_MODEL::NONE: break;
case TURB_TRANS_MODEL::EN: cout << "Low-turbulence Transition model: eN 1 equation model (2014)" << endl; break;
}
Copy link
Contributor

Choose a reason for hiding this comment

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

Formatting.
Should we also include LM in this switch

Copy link
Author

@RichRoos RichRoos Sep 21, 2022

Choose a reason for hiding this comment

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

Thanks. Updated.

As for the LM option. There are a lot of locations that have this switch. For now I have left the 'case TURB_TRANS_MODEL::LM :' option out as it will be added to develop when the LM PR will be merged. I have done it in this way as in most locations the LM code is not there of course.

I should add that on these LM/EN switch locations the LM implementation has the same code set-up. The thought was that merging an extra case in the switch list would produce no problems?

trans_conv/diff and soures filled in. eN transition added to the ft2 term inside SA sources. I have added my T3A- mesh and .cfg that I use for testing to the Quickstart folder.

Next step is to debug eN itself. Regular run is not yet working.
Common/src/CConfig.cpp Outdated Show resolved Hide resolved
@lgtm-com
Copy link

lgtm-com bot commented Sep 22, 2022

This pull request introduces 1 alert when merging 56b93c1 into 2642f7f - view on LGTM.com

new alerts:

  • 1 for Comparison result is always the same

@lgtm-com
Copy link

lgtm-com bot commented Sep 22, 2022

This pull request introduces 1 alert when merging fd51593 into 2642f7f - view on LGTM.com

new alerts:

  • 1 for Comparison result is always the same

@pcarruscag
Copy link
Member

Don't add meshes to the code repo, we have TestCases for that.

@lgtm-com
Copy link

lgtm-com bot commented Sep 22, 2022

This pull request introduces 1 alert when merging 604aadb into 2642f7f - view on LGTM.com

new alerts:

  • 1 for Comparison result is always the same

Comment on lines 714 to 716
virtual void SetAmplificationFactor(su2double val_amplification_factor_i, su2double val_amplification_factor_j) {
amplification_factor_i = val_amplification_factor_i;
amplification_factor_j = val_amplification_factor_j;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
virtual void SetAmplificationFactor(su2double val_amplification_factor_i, su2double val_amplification_factor_j) {
amplification_factor_i = val_amplification_factor_i;
amplification_factor_j = val_amplification_factor_j;
void SetAmplificationFactor(su2double val_amplification_factor_i, su2double val_amplification_factor_j) {
amplification_factor_i = val_amplification_factor_i;
amplification_factor_j = val_amplification_factor_j;

class CUpwSca_TransEN final : public CUpwScalar<FlowIndices> {
private:
using Base = CUpwScalar<FlowIndices>;
using Base::nDim;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
using Base::nDim;

/*!
* \brief Adds any extra variables to AD.
*/
void ExtraADPreaccIn() override {}
Copy link
Member

Choose a reason for hiding this comment

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

Densities are missing here

* \param[in] config - Definition of the particular problem.
*/
void FinishResidualCalc(const CConfig* config) override {
Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0];
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0];
Flux[0] = a0*V_i[idx.Density()]*ScalarVar_i[0] + a1*V_j[idx.Density()]*ScalarVar_j[0];

Comment on lines 77 to 82
* \param[in] val_nDim - Number of dimensions of the problem.
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwSca_TransEN(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
* \param[in] val_nDim - Number of dimensions of the problem.
* \param[in] val_nVar - Number of variables of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwSca_TransEN(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config)
: CUpwScalar<FlowIndices>(val_nDim, val_nVar, config) {}
* \param[in] val_nDim - Number of dimensions of the problem.
* \param[in] config - Definition of the particular problem.
*/
CUpwSca_TransEN(unsigned short val_nDim, const CConfig* config)
: CUpwScalar<FlowIndices>(val_nDim, 1, config) {}

Comment on lines 278 to 290
var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji;
static void get(const su2double& n_hat, CSAVariables& var) {
if (var.transEN == false){
const su2double xsi2 = pow(var.Ji, 2);
Copy link
Member

Choose a reason for hiding this comment

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

Please format the code... this is hard to read like this

@@ -103,6 +103,18 @@ class CSolverFactory {
*/
static CSolver* CreateTurbSolver(TURB_MODEL kindTurbModel, CSolver **solver, CGeometry *geometry, CConfig *config, int iMGLevel, int adjoint);

/*!
* \brief Create a turbulent solver
Copy link
Member

Choose a reason for hiding this comment

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

turbulent -> transition, check other places please

Comment on lines 51 to 54
/*!
* \brief Destructor of the class.
*/
~CTransENSolver() = default;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
/*!
* \brief Destructor of the class.
*/
~CTransENSolver() = default;

Comment on lines 1367 to 1370
if ( NONE) {
SU2_MPI::Error("Invalid convective scheme for the transition equations.", CURRENT_FUNCTION);
}
else if (EN) numerics[iMGlevel][TRANS_SOL][conv_term] = new CUpwSca_TransEN<Indices>(nDim, nVar_Trans, config);
Copy link
Member

Choose a reason for hiding this comment

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

I don't understand the if (NONE) parts, NONE is a constant and always 0

Comment on lines 43 to 48
void CTransENVariable::SetAmplificationFactor(unsigned long iPoint, su2double val_AmplificationFactor) {
AD::StartPreacc();
AmplificationFactor(iPoint) = val_AmplificationFactor;
AD::SetPreaccOut(AmplificationFactor(iPoint)); //Still need to be fixed?
AD::EndPreacc();
}
Copy link
Member

Choose a reason for hiding this comment

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

This should be done inline and preaccumulation does not need to be used for simple assignments like these

Suggested change
void CTransENVariable::SetAmplificationFactor(unsigned long iPoint, su2double val_AmplificationFactor) {
AD::StartPreacc();
AmplificationFactor(iPoint) = val_AmplificationFactor;
AD::SetPreaccOut(AmplificationFactor(iPoint)); //Still need to be fixed?
AD::EndPreacc();
}

@sun5k
Copy link
Contributor

sun5k commented Nov 25, 2022

Hi @RichRoos.
Hi. I was looking forward to implementing your e^N model. Relatively recently, I read #1820.
Can you share the problem you are facing?

AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim);

su2double rho = V_i[idx.Density()];
su2double p = V_i[idx.Pressure()];

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable p is not used.

const su2double vel_mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20);

su2double rhoInf = config->GetDensity_FreeStreamND();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable rhoInf is not used.
const su2double vel_mag = max(sqrt(vel_u * vel_u + vel_v * vel_v + vel_w * vel_w), 1e-20);

su2double rhoInf = config->GetDensity_FreeStreamND();
su2double pInf = config->GetPressure_FreeStreamND();

Check notice

Code scanning / CodeQL

Unused local variable Note

Variable pInf is not used.
// var.ft2 = var.ct3 * (1 - exp(2*(var.amplification - var.Ncrit)) ) * exp(-var.ct4 * xsi2);
var.ft2 = var.ct3 * (1 - exp(var.modifiedintermittency) );
var.d_ft2 = 0.0;
// var.d_ft2 = -2.0 * var.ct4 * var.Ji * var.ft2 * var.d_Ji;

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
const su2double xsi2 = pow(var.Ji, 2);

if (var.transEN == true){
// var.ft2 = var.ct3 * (1 - exp(2*(var.amplification - var.Ncrit)) ) * exp(-var.ct4 * xsi2);

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +168 to +169
// if (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) var.ct4 = 0.025;
// else var.ct4 = 0.05;

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +163 to +164
// var.Ncrit = -8.43 - 2.4*log(config->GetTurbulenceIntensity_FreeStream()/100);
// var.amplification = min(amplification_factor_i, var.Ncrit);

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
Comment on lines +171 to +188
// void FinishResidualCalc(const CConfig* config) override {
// const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
//
// /*--- Compute mean effective dynamic viscosity ---*/
// const su2double diff_i_amplification = (Laminar_Viscosity_i + Eddy_Viscosity_i)/sigma_n;
// const su2double diff_j_amplification = (Laminar_Viscosity_j + Eddy_Viscosity_j)/sigma_n;
//
// const su2double diff_amplification = 0.5*(diff_i_amplification + diff_j_amplification);
//
// Flux[0] = diff_amplification*Proj_Mean_GradScalarVar[0];
//
// /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
//
// if (implicit) {
// Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-diff_amplification*proj_vector_ij);
// Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+diff_amplification*proj_vector_ij);
// }
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
const auto index = counter * Restart_Vars[1] + skipVars;
for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]);
nodes ->SetAmplificationFactor(iPoint_Local, Restart_Data[index + 2]);
// nodes ->SetModifiedIntermittency(iPoint_Local, Restart_Data[index + 3]);

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
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.

None yet

5 participants