Updating SA-Diffusion term Discretisation

[Bot-Enigma-0]

Proposed Changes

Give a brief overview of your contribution here in a few sentences.
The current discretisation of the diffusion terms in the SA model are not consistent with the recommendation of S-A from their original paper [1]. The current implementation treats the non-linear term ($c_{b2}(\tilde{\nu})^2$) as a cross-production source term. By applying the product rule and simplifying, the diffusion term can be written as:

Then discretising gives:

Note that S-A assume that $\nabla\nu=0$, which does not hold in compressible cases.

The current changes are only for the Base model, I am extending it to the others. Would appreciate any feedback!

[1]: Spalart, Philippe, and Steven Allmaras. "A one-equation turbulence model for aerodynamic flows." 30th aerospace sciences meeting and exhibit. 1992.

Related Work

Resolve any issues (bug fix or feature request), note any related PRs, or mention interactions with the work of others, if any.

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.

[YairMO]

Hi,

Thank you for addressing this issue. The source diffusion (please note that this term is known as source diffusion rather than cross production source term) term in the SU2 code is currently treated explicitly (which may limit the CFL numbers of the SA model). One simple way is to adopt the original proposal by Spalart & Allmaras, as you pointed out. One should verify that the total sum of the coefficients is positive, resulting in a diffusive flux. Fortunately, using simple arithmetic averaging will result in a diffusive flux (as opposed to an anti-diffusive flux). Please bear in mind that for the negative version, some care should be taken: please follow the recent paper by Boris Diskin et al:

Revised USM3D-ME Implementation of Negative Version of Spalart–Allmaras Turbulence Model

AIAA Journal in press.

[Bot-Enigma-0]

Thank you for the paper. When extending this to the negative model I applied a very similar approach to the standard version. But after reading the positivity and boundedness discussion, I don't think the implementation I had would not be very stable. Rather the formulation the authors provide increases the stability, particularly due to the treatment of the $f_n$ term

[YairMO]

To make sure that the formulation is stable, it should be verified that the sum of the coefficients is always positive. Diskin addressed this issue in great detail.

[pcarruscag]

Residual_* and Jacobian_* are solver members that are not allocated for most solvers anymore.
The easiest way to test this is to compute the residual twice with the previous version you had, but swaping the order of ij (and flipping the normal direction) when setting things in the numerics class.

Suggested change

	numerics->ComputeResidual(Residual_i, Residual_j, Jacobian_ii, Jacobian_ij, Jacobian_ji, Jacobian_jj, const_cast<CConfig*>(config));
	numerics->ComputeResidual(Residual_i, Residual_j, Jacobian_ii, Jacobian_ij, Jacobian_ji, Jacobian_jj, const_cast<CConfig*>(config));

[pcarruscag]

Nice use of lambdas 👍

[pcarruscag]

The type returned by the numerics is just a view into their internal state. It's not a copy of the residuals and Jacobians.
So if you make these call back to back like this residual_ij and residual_ji will be pointing to the same values because they come from the same numerics class.
You need to compute one, use it to update the residual, then compute the other.

[pcarruscag]

Don't worry about this reducer stuff for now, it's not used unless you use OpenMP.

[pcarruscag]

Use only one of the residual results to update the Jacobians for now, it doesn't have to be super exact.

[pcarruscag]

This seems too approximate. nu_ij can be much smaller that the term multiplying the projected gradient.

[pcarruscag]

You can try limiting the overall diffusion coefficient (the value that multiplies Proj_Mean_GradScalarVar) to be positive.

[Bot-Enigma-0]

I exported the coefficient values to a buffer file, many of them appear to be negative. using the exact Jacobians reduces the negativeness, but still causes the solver to diverge.

[Bot-Enigma-0]

its the projected gradient (Proj_Mean_GradScalarVar[0]) which is causing the negativity. The diffusion coefficient always is positive but the product is not, which is causing the solver to diverge.

[YairMO]

Sorry, but I didn't follow through with all your work. Did you also modify the Jacobian due to the new additional diffusion flux (actually diffusion and anti-diffusion fluxes)?

[Bot-Enigma-0]

Sorry, but I didn't follow through with all your work. Did you also modify the Jacobian due to the new additional diffusion flux (actually diffusion and anti-diffusion fluxes)?

yes, i derived the exact Jacobians for the new discretisation

[YairMO]

Do you have a written draft?

[Bot-Enigma-0]

for the jacobians? i'll type them up and post them in a moment

@YairMO heres the derivation. I've written it in the way it appears in the code.

jacobain_derivation_of_neg_model.pdf

[pcarruscag]

You also need to modify CScalarSolver<VariableType>::SumEdgeFluxes
The convective fluxes still use EdgeFluxes, so one option is to store the difference between ij and ji in the second EdgeFlux vector.
Something like:

  for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
      if (iPoint == geometry->edges->GetNode(iEdge, 0)) {
        // Flux i-j.
        LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
      } else {
        // Flux j-i.
        LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
        if (EdgeFluxesDiff.size() > 0) {
          LinSysRes.SubtractBlock(iPoint, EdgeFluxesDiff.GetBlock(iEdge));
        }
      }
    }

[pcarruscag]

You don't need this change on the flow solver, just the scalar and turbulence.

Suggested change

	CSysVector<su2double> EdgeFluxes_ij; /*!< \brief Flux across each edge ij (non-conservative). */
	CSysVector<su2double> EdgeFluxes_ji; /*!< \brief Flux across each edge ji (non-conservative). */

[pcarruscag]

This can also be reverted

[pcarruscag]

same

[pcarruscag]

I think it's better to keep using EdgeFluxes as before and just add another one where needed (instead of 2).

Suggested change

	CSysVector<su2double> EdgeFluxes; /*!< \brief Flux across each edge. */
	CSysVector<su2double> EdgeFluxes_ij; /*!< \brief Flux across each edge ij (non-conservative). */
	CSysVector<su2double> EdgeFluxes_ji; /*!< \brief Flux across each edge ji (non-conservative). */
	CSysVector<su2double> EdgeFluxes; /*!< \brief Flux across each edge. */
	CSysVector<su2double> EdgeFluxes_ij; /*!< \brief Flux across each edge ij (non-conservative). */

[pcarruscag]

For the reducer strategy we need to compromise and use an "average" jacobian, i.e. assuming a conservative flux.
With what you have now the Jacobian is getting doubled

[pcarruscag]

Suggested change

	su2double flipped_normal[3];
	su2double flipped_normal[MAXNDIM];

[pcarruscag]

You can change the regression.yaml file to use your testcases branch with the updated restarts.

[pcarruscag]

Suggested change

#include "../scalar/scalar_diffusion.hpp"