Implementation of EDDES length scale in su2 is different from the published paper

[ShiheJia]

Describe the bug

Hi, everyone. I'm trying to using hybrid rans/les method in my work. However, it seems that the implementation of SA_EDDES length scale in su2 is different from the paper, so I wonder whether it is a bug or just for the sake of experience.

The relevant code locate in the void SetDES_LengthScale() function in SU2_CFD/src/solvers/CTurbSASolver.cpp， in line 1480 to 1495 the vorticity-adaptive SGS is calculated. But different for the published function, the absolute value of edge vector between vertex i and j are used and such a form seems unphysical.

1480 ln_max = 0.0; 1481 deltaDDES = 0.0; 1482 for (iNeigh = 0;iNeigh < nNeigh; iNeigh++){ 1483 jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); 1484 coord_j = geometry->nodes->GetCoord(jPoint); 1485 for (iDim = 0; iDim < nDim; iDim++){ 1486 delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]); //why absolute value 1487 } 1488 deltaDDES = geometry->nodes->GetMaxLength(iPoint); 1489 ln[0] = delta[1]*ratioOmega[2] - delta[2]*ratioOmega[1]; 1490 ln[1] = delta[2]*ratioOmega[0] - delta[0]*ratioOmega[2]; 1491 ln[2] = delta[0]*ratioOmega[1] - delta[1]*ratioOmega[0]; 1492 aux_ln = sqrt(ln[0]*ln[0] + ln[1]*ln[1] + ln[2]*ln[2]); 1493 ln_max = max(ln_max,aux_ln); 1494 vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint); 1495 }

Any help would be appreciated.

[ShiheJia]

[bigfooted]

Hi, could you place a link to the paper that you are referring to?

[TobiKattmann]

To everyone searching for the code in latest develop (there also @bigfooted's paper is given maybe, also it might be another from Eduardo Molina?) I posted it below. It seems you are using an older version of the code so please compare it with the latest which is posted below. (I am too lazy to check myself, and I am also not knowledgeable in this stuff)

(@ShiheJia when posting code please dont include every line number, keep the formatting and post the code wrapped in '''c++ <code> ''' with ' being the forward pointing tick)

case SA_EDDES: {
        /*--- An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows.
         Shur et al.
         Flow Turbulence Combust - 2015
         ---*/

        su2double vortexTiltingMeasure = nodes->GetVortex_Tilting(iPoint);

        const su2double omega = GeometryToolbox::Norm(3, vorticity);

        su2double ratioOmega[MAXNDIM] = {};

        for (auto iDim = 0; iDim < 3; iDim++){
          ratioOmega[iDim] = vorticity[iDim]/omega;
        }

        const su2double deltaDDES = geometry->nodes->GetMaxLength(iPoint);

        su2double ln_max = 0.0;
        for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) {
          const auto coord_j = geometry->nodes->GetCoord(jPoint);
          su2double delta[MAXNDIM] = {};
          for (auto iDim = 0u; iDim < nDim; iDim++){
            delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);
          }
          su2double ln[3];
          ln[0] = delta[1]*ratioOmega[2] - delta[2]*ratioOmega[1];
          ln[1] = delta[2]*ratioOmega[0] - delta[0]*ratioOmega[2];
          ln[2] = delta[0]*ratioOmega[1] - delta[1]*ratioOmega[0];
          const su2double aux_ln = sqrt(ln[0]*ln[0] + ln[1]*ln[1] + ln[2]*ln[2]);
          ln_max = max(ln_max, aux_ln);
          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
        }
        vortexTiltingMeasure /= (nNeigh + 1);

        const su2double f_kh = max(f_min,
                                   min(f_max,
                                       f_min + ((f_max - f_min)/(a2 - a1)) * (vortexTiltingMeasure - a1)));

        const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0));
        const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0));

        su2double maxDelta = (ln_max/sqrt(3.0)) * f_kh;
        if (f_d < 0.999){
          maxDelta = deltaDDES;
        }

        const su2double distDES = constDES * maxDelta;
        lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES));

        break;
      }
    }

    nodes->SetDES_LengthScale(iPoint, lengthScale);

In case @EduardoMolina is still around, this might be interesting for him.

[ShiheJia]

Sorry for the late reply! Thanks for all the help, and i am afraid i am using a 7.2.0 version and the newly released version 7.2.1 remains the same code. The problem is just lying on "delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);" which is also in @TobiKattmann 's code post(thanks for your kind guidance ). In deed, this part of code should give the credit to @EduardoMolina, and the function is a part of his doctoral thesis(2018) which i just couldnot found a link or doi of. But in this paper(https://www.researchgate.net/publication/318143234 ), he gives the official definition as below without the implementation in above picture.

[ShiheJia]

[TobiKattmann]

Unfortunately I am not much of a help when it comes to deciding whether the paper or code version is the one to use. Is it possible to just implement the paper version and compare results? Not that is solves the discrepancy, but maybe there is a reason why the code version is different, and that might become more apparent when comparing results.

[EduardoMolina]

Hello guys!

@ShiheJia Could you please elaborate why do you think it is a bug? Did you tested the code?

[ShiheJia]

Thanks for the suggestion and I am preparing for a test. As far as I understand it, I just find that to get the absolute value of r_ij ( showed in the highest equation) in this part of code is unnecessary. The cross-product operation is to find the grid vector mostly parallel to the vorticity vector and the absolute value may cause a nonphysical recognition. @EduardoMolina, I don't know if I got it wrong and wish more guidance,.

for (auto iDim = 0u; iDim < nDim; iDim++){
      delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);
}

[HexFluid]

Thanks for the suggestion and I am preparing for a test. As far as I understand it, I just find that to get the absolute value of r_ij ( showed in the highest equation) in this part of code is unnecessary. The cross-product operation is to find the grid vector mostly parallel to the vorticity vector and the absolute value may cause a nonphysical recognition. @EduardoMolina, I don't know if I got it wrong and wish more guidance,.

for (auto iDim = 0u; iDim < nDim; iDim++){
      delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);
}

Hi Shihe,

I checked the implementation and I think it is correctly done. You do need the absolute value (i.e., delta has a unit of [m] or equivalent) to keep the correct dimension of nu_t based on a Smagorinsky-type SGS model.

You may find the appendix of this paper useful for your understanding of delta_omg: https://doi.org/10.1007/s00162-011-0240-z

Also note that delta_omg does not always outperforms its peers - vorticity may not be aligned to the rotation axis of a local vortex (e.g., in rotating reference frame, in attached boundary layer, to name a few), in which case the physical meaning of delta_omg becomes vague. See also my work for a brief review of DES-type methods and some applications: https://doi.org/10.1115/1.4052019

[pcarruscag]

Disclaimer: I'm not a turbulence specialist.

With that said, that definition of delta looks dodgy, what is the physical meaning of "rectifying" a distance vector to have all components positive? What happens to my results if I compute them in a rotated domain? (suppose I like my channels at 45 degrees...) The magnitude of the cross product is going to change.

[ShiheJia]

Well, I cannot agree with @pcarruscag any more. I try to introduce some comments into this part of code as below.

//  loop all the neighboorhood point of point i
for (iNeigh = 0;iNeigh < nNeigh; iNeigh++){

          // get the point number and coordinates
          jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
          coord_j = geometry->nodes->GetCoord(jPoint);

          // loop dimesion and calculate the vector r_ij
         //  the core of the issue
          for (iDim = 0; iDim < nDim; iDim++){
            delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]);
          }

          // get delta_max used for DES or DDES that is out of concern here
          deltaDDES = geometry->nodes->GetMaxLength(iPoint);

          // do the cross product of r_ij and n_omega
          ln[0] = delta[1]*ratioOmega[2] - delta[2]*ratioOmega[1];
          ln[1] = delta[2]*ratioOmega[0] - delta[0]*ratioOmega[2];
          ln[2] = delta[0]*ratioOmega[1] - delta[1]*ratioOmega[0];

          // get the length of the cross product result
          aux_ln = sqrt(ln[0]*ln[0] + ln[1]*ln[1] + ln[2]*ln[2]);

          // compare to get a max length value
          ln_max = max(ln_max,aux_ln);

          // used for another function of EDDES  that is also out of concern here
          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
        }

And the function is shown here. Just as @pcarruscag said, if we previously make all components of r_ij positive, it may just be unsuitable for may cases. @HexFluid thanks for the recommand of the papers but perhaps I just cannot get you, could you please explain a little more?

[HexFluid]

Well, I cannot agree with @pcarruscag any more. I try to introduce some comments into this part of code as below.

//  loop all the neighboorhood point of point i
for (iNeigh = 0;iNeigh < nNeigh; iNeigh++){

          // get the point number and coordinates
          jPoint = geometry->nodes->GetPoint(iPoint, iNeigh);
          coord_j = geometry->nodes->GetCoord(jPoint);

          // loop dimesion and calculate the vector r_ij
         //  the core of the issue
         // !!! should not use absolute here as it changes the vector direction
          for (iDim = 0; iDim < nDim; iDim++){
            delta[iDim] = fabs(coord_j[iDim] - coord_i[iDim]); 
          }

          // get delta_max used for DES or DDES that is out of concern here 
          // !!! we do need this here. It will be used when fd<0.999 later.
          deltaDDES = geometry->nodes->GetMaxLength(iPoint);

          // do the cross product of r_ij and n_omega
          // !!! this follows the Eqns. designed for structured solver. We need to upgrade this to the unstructure version.
          ln[0] = delta[1]*ratioOmega[2] - delta[2]*ratioOmega[1];
          ln[1] = delta[2]*ratioOmega[0] - delta[0]*ratioOmega[2];
          ln[2] = delta[0]*ratioOmega[1] - delta[1]*ratioOmega[0];

          // get the length of the cross product result
          aux_ln = sqrt(ln[0]*ln[0] + ln[1]*ln[1] + ln[2]*ln[2]);

          // compare to get a max length value
          ln_max = max(ln_max,aux_ln);

          // used for another function of EDDES  that is also out of concern here
          vortexTiltingMeasure += nodes->GetVortex_Tilting(jPoint);
        }

And the function is shown here. Just as @pcarruscag said, if we previously make all components of r_ij positive, it may just be unsuitable for may cases. @HexFluid thanks for the recommand of the papers but perhaps I just cannot get you, could you please explain a little more?

!!! The Eq. below is originally designed for solvers using structured grid, and it is the current version of implementation in SU2. For unstructure solvers please refer to the appendix of https://doi.org/10.1007/s00162-011-0240-z

Ah now I see what you mean now. There seems to be more than the issues you identified in line 1486. See my comments (starting with !!!) above.

I am not sure who is working on the dev of this script but I am happy to provide some knowledge on improving it. I have this feature implemented in my own code and is easily transferable to SU2.

[pcarruscag]

WIP #1622

[stale[bot]]

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]

@GomerOfDoom is this something you can look into? Maybe there is an easy fix to make this consistent with literature.

[pcarruscag]

@rois1995 since you are looking into DDES stuff can you take a look at this issue please?

[rois1995]

Sure! The implementation looks suspicious. Now that I have the results for the NACA0021 at 60deg I can try play with it a little bit.