Correct inconsistency regarding the density used in conservative scalar solver

[emaberman]

Proposed Changes

In SU2 the flow equations and the RANS turbulence model equations are solved in a decoupled manner. The flow equations are first solved, and their solution is advanced to the next time level, 'n+1'. The advanced flow variables, specifically the density and velocity (of time level 'n+1'), are used to solve the RANS turbulence model equations.

Since the SST model solves the conservative variables (density * TKE, density * Omega), it's residuals and linearization of the implicit operator are taken at the mixed time level of density^(n+1) * (TKE, Omega)^n.

In a similar manner for consistency , the updated conservative variables (of the SST model) of time level 'n+1' should be advanced as follows:
density^(n+1) * (TKE, Omega)^(n+1) = density^(n+1) * (TKE,OMEGA)^n + Delta(densityTKE,densityOMEGA)

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.

[pcarruscag]

In theory we are doing exactly what you describe, we need density old because of the underrelaxation and clipping.
Have you tried printing some density values from inside the clipped update function to check if they are consistent or not?

[emaberman]

I'm not sure I understand your answer. The density_old is not used for under-relaxation or clipping. Instead, it is used when advancing the time step, here is the function:

` /*!

• \brief Update the variables using a conservative format.
• \param[in] iPoint - Point index.
• \param[in] iVar - Index of the variable.
• \param[in] solution - Value of the solution change.
• \param[in] lowerlimit - Lower value for Solution clipping.
• \param[in] upperlimit - Upper value for Solution clipping.
• \param[in] Sol2Conservative - Factor multiplied to Solution to get transported variable.
• \param[in] Sol2Conservative_old - Factor multiplied to Solution to get transported variable, of the previous Iteration.
  */
  inline void AddClippedSolution(unsigned long iPoint, unsigned long iVar, su2double solution,
  su2double lowerlimit, su2double upperlimit,
  su2double Sol2Conservative = 1.0, su2double Sol2Conservative_old = 1.0) {

su2double val_new = (Solution_Old(iPoint,iVar)*Sol2Conservative_old + solution)/Sol2Conservative;
Solution(iPoint,iVar) = min(max(val_new, lowerlimit), upperlimit);

}
`

I have validated this inconsistency by printing out the density values. Moreover we have identified that the current implementation, which uses the old density values, leads in cases to negative values in the conserved SST variables (before clipping). When using the updated density values, these negative values do not appear.

[pcarruscag]

Ok, I suppose what you mean is that we should not use density old because when we compute the residuals to advance the SST equations we are already using the density at n+1?
Using density at n would only be correct if we used that density to compute the SST residuals?
Is the convergence improved in those cases with negative values? Can you share some comparison of the convergence?

[emaberman]

Yes that is exactly my point.
I did not check for improved convergence of this specific change alone.
We encountered this while implementing a method to ensure positivity of the turbulence models (apart from negative values introduced due to the linear solver). This full method (which includes this change among others) is currently in testing and initial results show that using this method better stability and robustness of the solver will be introduced (we will of course commit feature and share results when complete).

While the change proposed here is critical for our feature, it is in our understanding that regardless of our work, the current implementation (even if not catastrophic because of the small differences between the two density values) is incorrect, and should be changed regardless of better convergence.

[YairMO]

Hi All,

Let F denote the convective flux, and F^(n) is the flux at the current time level n and is given as follows:

F^n = r^(n+1) * v^(n+1) * phi^(n)

Where 'r' is the density, v is the velocity, and phi is the turbulence variables (TKE or Omega). It should be noted that 'r' and 'v' are taken at time level (n+1) because the flow equations were solved and updated before solving the turbulence model equations.

Similarly, F at time level (n+1) is given as follows:

F^(n+1) = r^(n+1) * v^(n+1) * phi^(n+1)

Next, a linearization is of F^(n+1) is conducted as follows:

F^(n+1) = F^n + DF/D(r*phi) * Delta(r * phi) = r^(n+1) * v^(n+1) * phi^(n) + v^(n+1) * Delta(r * phi)

Now to find the definition of Delta(r * phi) we simply substract F^n from F^(n+1) as follows:

F^(n+1) - F^n = v^(n+1) * Delta(r * phi)
divide the above equation by v^(n+1) we get the correct definition of Delta(r*phi)

r^(n+1) * phi^(n+1) - r^n+1 * phi^n = Delta(r*phi)

This may clarify the inconsistency present in the code.
Thank you

[YairMO]

I should have emphasized that the current situation in the code is different. The present coding is given as follows:

r^(n+1) * phi^(n+1) - r^n * phi^n = Delta(r*phi)

Note that the r of the second term on the left-hand side is given at time level n, which should be at time level n+1.

[pcarruscag]

Thank you for the details, those will be helpful for someone reading the PR later.

As you can see a lot of testcases need to be updated to be able to merge this PR https://github.com/su2code/SU2/actions/runs/13782145187/job/38543199068?pr=2458

If you want to do that you can check the values on GitHub and carry out the update of the test scripts (.py) inside the TestCases folder of the code.
If you think this is not critical for convergence and you have other changes to SST, you could bundle this fix with your other improvements and save the tedious update of testcase values.

[emaberman]

Thanks for your response and feedback.

Our work includes a non trivial modification to the implicit treatment of the problem and will therefore will be added as an option for the user. Unless adopted as the default treatment in SU2 our modifications should have no impact on existing SST model results. To ensure this I thought it was more correct to divide the PR into two different ones.

The easier option for me would indeed be to wrap this as part of our treatment, however it is my understanding that this fix should be available to all regardless of our treatment, and even if it does not donate to convergence in most/all cases. I see this as my minute contribution to the SU2 community.

Though it may take me a few days I will correct the regression tests..

[bigfooted]

Looking at some of the new results from the regression tests it looks like the impact on convergence is minor. Sometimes a little better, sometimes not.
Really great work!

[rois1995]

Looking at some of the new results from the regression tests it looks like the impact on convergence is minor. Sometimes a little better, sometimes not. Really great work!

Maybe it is due to the fact that the regression cases only run for very few iterations. Then, when a full simulation is considered, we may have visible gains in robustness.

[pcarruscag]

Thank you for updating the tests 👍