Updating SA-Diffusion term Discretisation by Bot-Enigma-0 · Pull Request #2508 · su2code/SU2

Bot-Enigma-0 · 2025-05-16T18:26:47Z

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 · 2025-05-20T09:23:01Z

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 · 2025-05-20T14:33:05Z

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 · 2025-05-20T16:32:50Z

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 · 2025-05-26T18:40:42Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+    SolverSpecificNumerics(iPoint, jPoint);
+
+    /*--- Compute residual, and Jacobians ---*/
+    numerics->ComputeResidual(Residual_i, Residual_j, Jacobian_ii, Jacobian_ij, Jacobian_ji, Jacobian_jj, const_cast<CConfig*>(config));


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));

SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp

pcarruscag · 2025-05-27T15:31:08Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+    auto jPoint = geometry->edges->GetNode(iEdge, 1);
+
+  /*--- Helper function to compute the flux ---*/
+  auto ComputeFlux = [&](unsigned long point_i, unsigned long point_j, const su2double* normal) {


Nice use of lambdas 👍

pcarruscag · 2025-05-27T15:32:32Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+  auto residual_ij = ComputeFlux(iPoint, jPoint, normal);
+  auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal);


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 · 2025-05-27T15:33:28Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+    EdgeFluxes.SubtractBlock(iEdge, residual_ij);  
+    EdgeFluxes.AddBlock(iEdge, residual_ji);  
+    if (implicit) {
+      Jacobian.UpdateBlocksSub(iEdge, residual_ij.jacobian_i, residual_ij.jacobian_j);
+      Jacobian.UpdateBlocks(iEdge, residual_ji.jacobian_i, residual_ji.jacobian_j);
+    }


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

pcarruscag · 2025-05-27T15:34:47Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+    LinSysRes.SubtractBlock(iPoint, residual_ij);  
+    LinSysRes.AddBlock(jPoint,  residual_ji);  
+    if (implicit) {
+      Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual_ij.jacobian_i, residual_ij.jacobian_j);
+      Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual_ji.jacobian_i, residual_ji.jacobian_j);


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

SU2_CFD/src/solvers/CTurbSASolver.cpp

pcarruscag · 2025-05-28T19:46:51Z

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

+    if (implicit) { 
+      Jacobian_i[0][0] = (0.5*Proj_Mean_GradScalarVar[0]-nu_ij*proj_vector_ij)/sigma;
+      Jacobian_j[0][0] = (0.5*Proj_Mean_GradScalarVar[0]+nu_ij*proj_vector_ij)/sigma;


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

pcarruscag · 2025-05-28T19:46:54Z

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

+    su2double term_1 = (nu_ij + (1 + cb2)*nu_tilde_ij*fn_i) * Proj_Mean_GradScalarVar[0]/sigma;
+    su2double term_2 = cb2*nu_tilde_i*fn_i* Proj_Mean_GradScalarVar[0]/sigma;


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

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.

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.

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

YairMO · 2025-05-29T16:09:06Z

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 · 2025-05-29T16:14:24Z

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 · 2025-05-29T16:16:52Z

Do you have a written draft?

Bot-Enigma-0 · 2025-05-29T16:29:08Z

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

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

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 · 2025-06-13T15:40:07Z

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

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


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 · 2025-06-13T15:40:23Z

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

+      if (config->GetKind_Solver() == MAIN_SOLVER::RANS) {
+        for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
+
+          LinSysRes.SetBlock_Zero(iPoint);
+
+          for (auto iEdge : geometry->nodes->GetEdges(iPoint)) {
+            if (iPoint == geometry->edges->GetNode(iEdge,0))
+              LinSysRes.AddBlock(iPoint, EdgeFluxes_ij.GetBlock(iEdge));
+            else
+              LinSysRes.SubtractBlock(iPoint, EdgeFluxes_ji.GetBlock(iEdge));
+          }
+        }        
+      }
+      else SumEdgeFluxes(geometry);


This can also be reverted

pcarruscag · 2025-06-13T15:40:35Z

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

+  if (ReducerStrategy) {
+    EdgeFluxes.Initialize(geometry.GetnEdge(), geometry.GetnEdge(), nVar, nullptr);
+    EdgeFluxes_ij.Initialize(geometry.GetnEdge(), geometry.GetnEdge(), nVar, nullptr); //initializing only when solver=rans, has no effect
+    EdgeFluxes_ji.Initialize(geometry.GetnEdge(), geometry.GetnEdge(), nVar, nullptr);
+    }


pcarruscag · 2025-06-13T15:42:13Z

SU2_CFD/include/solvers/CScalarSolver.hpp

  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). */


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 · 2025-06-13T15:45:31Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+    if (implicit) {
+      Jacobian.UpdateBlocksSub(iEdge, residual_ij.jacobian_i, residual_ij.jacobian_j);
+    }


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 · 2025-06-13T15:47:39Z

SU2_CFD/src/solvers/CTurbSASolver.cpp

+  }
+
+  /*--- Compute fluxes and jacobians j->i ---*/
+  su2double flipped_normal[3];


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 · 2025-06-20T16:03:42Z

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

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


Suggested change

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

Bot-Enigma-0 added 4 commits May 16, 2025 18:43

Base model changes

ca7bf9e

freeze 'cross-production' from sources

fb0b376

something to search through code, delete later

457c590

update cell-centre nu_tilde

7d83805

bigfooted added the changelog:chore label May 20, 2025

Bot-Enigma-0 added 4 commits May 26, 2025 19:18

non-conservative class in scalar_diffusion

7e9d43a

evaluation of non-cons fluxes and jacobians for base model

8b61921

updated edge looping for SA Solver

bf50893

remove extra ;

0213d1f

pcarruscag reviewed May 26, 2025

View reviewed changes

github-advanced-security bot found potential problems May 26, 2025

View reviewed changes

SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp Fixed Show fixed Hide fixed

Bot-Enigma-0 added 3 commits May 27, 2025 16:26

revert turb_diff

0bd4013

compute fluxes twice in TurbSASolver

7049631

remove non-conservative class

8a5ca6c

pcarruscag reviewed May 27, 2025

View reviewed changes

Bot-Enigma-0 added 2 commits May 27, 2025 17:30

evaluate and accumulate each flux separately

10cc4dc

converging CTurbSASolver

57db642

github-advanced-security bot found potential problems May 27, 2025

View reviewed changes

SU2_CFD/src/solvers/CTurbSASolver.cpp Fixed Show fixed Hide fixed

SU2_CFD/src/solvers/CTurbSASolver.cpp Fixed Show fixed Hide fixed

SU2_CFD/src/solvers/CTurbSASolver.cpp Fixed Show fixed Hide fixed

SU2_CFD/src/solvers/CTurbSASolver.cpp Fixed Show fixed Hide fixed

negative model

823c335

pcarruscag reviewed May 28, 2025

View reviewed changes

exact jacobians for negative model

8c976d5

github-advanced-security bot found potential problems May 29, 2025

View reviewed changes

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp Fixed Show fixed Hide fixed

Bot-Enigma-0 added 2 commits June 3, 2025 14:53

clean up diffusion numerics

e590b37

clean up CTurbSASolver

42dec93

github-advanced-security bot found potential problems Jun 3, 2025

View reviewed changes

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp Fixed Show fixed Hide fixed

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp Fixed Show fixed Hide fixed

Bot-Enigma-0 added 7 commits June 6, 2025 18:33

remove exact jacobians from neg and add option to use exact for base

69d9885

remove cross-production calculations from sources

bcb0f94

small typo, does not affect results

1fcb8de

remove search script

53ee8ca

define vectors

75574d9

reducer strategy for SA

fe10797

update summation in FVM solver

947e807

pcarruscag reviewed Jun 13, 2025

View reviewed changes

Bot-Enigma-0 added 5 commits June 16, 2025 21:41

revert flow solver

5f52759

update scalar solver for additional EdgeFluxes_Diff

861574f

update SA solver

c0687bc

correct Reducer accumulation

5750252

update regression tests

77f797b

Bot-Enigma-0 mentioned this pull request Jun 18, 2025

Updating restart files for new SA diffusion implementation su2code/TestCases#172

Merged

pcarruscag reviewed Jun 20, 2025

View reviewed changes

change regression directory

bb7a0c4

This was referenced Jun 25, 2025

[WIP] updating cross-diffusion discretisation in sst #2521

Closed

Updating the cross-diffusion discretisation in SST #2522

Merged

Bot-Enigma-0 added 3 commits June 25, 2025 18:20

update regression tests

1148082

add option specific to turbulence for accurate jacobians

81a0620

update regression

344e939

pcarruscag changed the title ~~[WIP] Updating SA-Diffusion term Discretisation~~ Updating SA-Diffusion term Discretisation Jun 27, 2025

pcarruscag marked this pull request as ready for review June 27, 2025 20:02

pcarruscag added changelog:feature and removed changelog:chore labels Jun 27, 2025

pcarruscag changed the base branch from develop to feature_sa_diffusion June 28, 2025 19:05

pcarruscag merged commit 5b56b32 into su2code:feature_sa_diffusion Jun 28, 2025
217 of 233 checks passed

	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));

		auto residual_ij = ComputeFlux(iPoint, jPoint, normal);
		auto residual_ji = ComputeFlux(jPoint, iPoint, flipped_normal);

		su2double term_1 = (nu_ij + (1 + cb2)nu_tilde_ijfn_i) * Proj_Mean_GradScalarVar[0]/sigma;
		su2double term_2 = cb2nu_tilde_ifn_i* Proj_Mean_GradScalarVar[0]/sigma;

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

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

Conversation

Bot-Enigma-0 commented May 16, 2025

Proposed Changes

Related Work

PR Checklist

Uh oh!

YairMO commented May 20, 2025

Uh oh!

Bot-Enigma-0 commented May 20, 2025

Uh oh!

YairMO commented May 20, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

YairMO commented May 29, 2025

Uh oh!

Bot-Enigma-0 commented May 29, 2025

Uh oh!

YairMO commented May 29, 2025

Uh oh!

Bot-Enigma-0 commented May 29, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Bot-Enigma-0 commented May 29, 2025 •

edited

Loading