Skip to content

Consistent mass fluxes between bulk flow and turbulence solvers #726

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.

Sign up for GitHub

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

Closed
pcarruscag wants to merge 3 commits into from
Closed

Consistent mass fluxes between bulk flow and turbulence solvers #726

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
17e5df0
Euler/NS solver stores fluxes that CSolverTurb can access
pcarruscag Jul 1, 2019
9c33db0
forgot to comment a function
pcarruscag Jul 1, 2019
b434c61
need to mark mass flux for preaccumulation
pcarruscag Jul 1, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions SU2_CFD/include/numerics_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ class CNumerics {
su2double *Enthalpy_formation;
su2double Prandtl_Lam; /*!< \brief Laminar Prandtl's number. */
su2double Prandtl_Turb; /*!< \brief Turbulent Prandtl's number. */
su2double MassFlux; /*!< \brief Mass flux across edge. */

public:

Expand Down Expand Up @@ -1484,10 +1485,15 @@ class CNumerics {
* \param[in] d: array that will hold the ordered eigenvalues
* \param[in] e: supplemental data structure
* \param[in] n: order of matrix V

*/
static void tql2(su2double **V, su2double *d, su2double *e, unsigned short n);

/*!
* \brief SetMassFlux
* \param[in] val_MassFlux: Mass flux across the edge
*/
inline void SetMassFlux(const su2double val_MassFlux) {MassFlux = val_MassFlux;}

};

/*!
Expand Down Expand Up @@ -2475,7 +2481,7 @@ class CUpwSca_TurbSST : public CUpwScalar {
/*!
* \brief Adds any extra variables to AD
*/
void ExtraADPreaccIn();
inline void ExtraADPreaccIn() {}
Copy link
Contributor

@clarkpede clarkpede Jul 8, 2019

Choose a reason for hiding this comment

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

The "inline" here is redundant. Functions defined inside a class definition are automatically "inline."

Copy link
Member Author

@pcarruscag pcarruscag Jul 8, 2019

Choose a reason for hiding this comment

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

In common practice yes but apparently that is not part of the standard (which I have not read) and so some compilers might complain about multiple definition.

Copy link
Contributor

@clarkpede clarkpede Jul 8, 2019

Choose a reason for hiding this comment

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

Ah ok. So no harm done in leaving it in.


/*!
* \brief SST specific steps in the ComputeResidual method
Expand Down
9 changes: 9 additions & 0 deletions SU2_CFD/include/solver_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,8 @@ class CSolver {
su2double **Smatrix, /*!< \brief Auxiliary structure for computing gradients by least-squares */
**Cvector; /*!< \brief Auxiliary structure for computing gradients by least-squares */

su2double *EdgeMassFluxes; /*!< \brief Mass fluxes across each edge, for discretization of passive scalars. */

int *Restart_Vars; /*!< \brief Auxiliary structure for holding the number of variables and points in a restart. */
int Restart_ExtIter; /*!< \brief Auxiliary structure for holding the external iteration offset from a restart. */
passivedouble *Restart_Data; /*!< \brief Auxiliary structure for holding the data values from a restart. */
Expand Down Expand Up @@ -4318,6 +4320,13 @@ class CSolver {
*/
virtual void ComputeVerificationError(CGeometry *geometry, CConfig *config);

/*!
* \brief Get the mass flux across an edge (computed and stored during the discretization of convective fluxes).
* \param[in] iEdge - Index of the edge.
* \return The mass flux across the edge.
*/
inline su2double GetEdgeMassFlux(const unsigned long iEdge) const {return EdgeMassFluxes[iEdge];}

protected:
/*!
* \brief Allocate the memory for the verification solution, if necessary.
Expand Down
75 changes: 24 additions & 51 deletions SU2_CFD/src/numerics_direct_turbulent.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,39 +65,17 @@ void CUpwScalar::ComputeResidual(su2double *val_residual,
CConfig *config) {

AD::StartPreacc();
AD::SetPreaccIn(Normal, nDim);
AD::SetPreaccIn(TurbVar_i, nVar); AD::SetPreaccIn(TurbVar_j, nVar);
if (grid_movement) {
AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim);
}
Copy link
Member Author

@pcarruscag pcarruscag Jul 1, 2019

Choose a reason for hiding this comment

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

By not computing the mass flux again we greatly reduce the number of input variables for preaccumulation.

AD::SetPreaccIn(TurbVar_i, nVar);
AD::SetPreaccIn(TurbVar_j, nVar);
AD::SetPreaccIn(MassFlux);

ExtraADPreaccIn();

Density_i = V_i[nDim+2];
Density_j = V_j[nDim+2];

q_ij = 0.0;
if (grid_movement) {
for (iDim = 0; iDim < nDim; iDim++) {
Velocity_i[iDim] = V_i[iDim+1] - GridVel_i[iDim];
Velocity_j[iDim] = V_j[iDim+1] - GridVel_j[iDim];
q_ij += 0.5*(Velocity_i[iDim]+Velocity_j[iDim])*Normal[iDim];
}
}
else {
for (iDim = 0; iDim < nDim; iDim++) {
Velocity_i[iDim] = V_i[iDim+1];
Velocity_j[iDim] = V_j[iDim+1];
q_ij += 0.5*(Velocity_i[iDim]+Velocity_j[iDim])*Normal[iDim];
}
}

a0 = 0.5*(q_ij+fabs(q_ij));
a1 = 0.5*(q_ij-fabs(q_ij));
a0 = max(0.0, MassFlux);
a1 = min(0.0, MassFlux);
Copy link
Member Author

@pcarruscag pcarruscag Jul 1, 2019

Choose a reason for hiding this comment

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

Here it is possible to set a0 = -a1 on the assumption that at convergence the summation of fluxes is 0, with that diagonal dominance is guaranteed (at least for the convective part) which improves the condition number of the assembled Jacobian matrix.


FinishResidualCalc(val_residual, val_Jacobian_i, val_Jacobian_j, config);


AD::SetPreaccOut(val_residual, nVar);
AD::EndPreacc();

Expand All @@ -113,16 +91,18 @@ CUpwSca_TurbSA::~CUpwSca_TurbSA(void) {
}

void CUpwSca_TurbSA::ExtraADPreaccIn() {
AD::SetPreaccIn(V_i, nDim+1); AD::SetPreaccIn(V_j, nDim+1);
AD::SetPreaccIn(V_i[nDim+2]); AD::SetPreaccIn(V_j[nDim+2]);
}

void CUpwSca_TurbSA::FinishResidualCalc(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) {

val_residual[0] = a0*TurbVar_i[0]+a1*TurbVar_j[0];
su2double OneOnAveRho = 2.0/(V_i[nDim+2]+V_j[nDim+2]);

Copy link
Member Author

@pcarruscag pcarruscag Jul 1, 2019

Choose a reason for hiding this comment

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

SA is funny and so density is still needed here...

val_residual[0] = (a0*TurbVar_i[0]+a1*TurbVar_j[0])*OneOnAveRho;

if (implicit) {
val_Jacobian_i[0][0] = a0;
val_Jacobian_j[0][0] = a1;
val_Jacobian_i[0][0] = a0*OneOnAveRho;
val_Jacobian_j[0][0] = a1*OneOnAveRho;
}
}

Expand Down Expand Up @@ -1087,27 +1067,20 @@ CUpwSca_TurbSST::CUpwSca_TurbSST(unsigned short val_nDim,
CUpwSca_TurbSST::~CUpwSca_TurbSST(void) {
}

void CUpwSca_TurbSST::ExtraADPreaccIn() {

AD::SetPreaccIn(V_i, nDim+3);
AD::SetPreaccIn(V_j, nDim+3);

}

void CUpwSca_TurbSST::FinishResidualCalc(su2double *val_residual,
su2double **val_Jacobian_i,
su2double **val_Jacobian_j,
CConfig *config) {

val_residual[0] = a0*Density_i*TurbVar_i[0]+a1*Density_j*TurbVar_j[0];
val_residual[1] = a0*Density_i*TurbVar_i[1]+a1*Density_j*TurbVar_j[1];
Copy link
Contributor

@clarkpede clarkpede Jul 8, 2019 •
edited
Loading

Choose a reason for hiding this comment

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

What's your rationale behind making these density changes?

I understand the problems with reconstruction and why the value of density is questionable. But it seems to me like removing density is incorrect. There's two reasons I think that:

  1. The SST variables are updated by the CVariable::AddConservativeSolution call. In that functions, the delta in the variables is assumed to be weighted by density, e.g. rho*k. Based on that you can conclude that the equations are equations for the density-weighted variables, and the residuals should also be density weighted.
  2. In the source term for the SST variables, (CSourcePieceWise_TurbSST) there are terms like (2/3*rho*k*diverg*Volume). These residuals are density weighted. To be consistent, all the residuals should be density weighted.

For a converged, steady state solution, I think these changes will only change the conditioning. Similarly, I expect no impact for a low-Mach test case. Have you tested these changes on a case like the RAE 2822?

Edit: I just saw the comments you made while I was typing out my review. The mass flux used here (in the convection terms) explains the first question. I still question the density changes in the source terms though.

Copy link
Member Author

@pcarruscag pcarruscag Jul 8, 2019

Choose a reason for hiding this comment

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

(I've been looking at the code a bit more to get an answer).
The rationale was: I did not want to divide the mass flux by an average density which would then break that "summation = 0" property we can use to improve conditioning a bit. Following that logic I linearized the other terms w.r.t. k and w instead of rho k and rho w, I did not change the source residual only its Jacobian. But I forgot to see how the solution was being updated so yeah in the end it is wrong.
Seeing how SA models need a volumetric flux anyway I will probably revert those changes to SST.

val_residual[0] = a0*TurbVar_i[0]+a1*TurbVar_j[0];
val_residual[1] = a0*TurbVar_i[1]+a1*TurbVar_j[1];

if (implicit) {
val_Jacobian_i[0][0] = a0; val_Jacobian_i[0][1] = 0.0;
val_Jacobian_i[1][0] = 0.0; val_Jacobian_i[1][1] = a0;
val_Jacobian_i[0][0] = a0; val_Jacobian_i[0][1] = 0.0;
val_Jacobian_i[1][0] = 0.0; val_Jacobian_i[1][1] = a0;

val_Jacobian_j[0][0] = a1; val_Jacobian_j[0][1] = 0.0;
val_Jacobian_j[1][0] = 0.0; val_Jacobian_j[1][1] = a1;
val_Jacobian_j[0][0] = a1; val_Jacobian_j[0][1] = 0.0;
val_Jacobian_j[1][0] = 0.0; val_Jacobian_j[1][1] = a1;
}
}

Expand Down Expand Up @@ -1156,11 +1129,11 @@ void CAvgGrad_TurbSST::FinishResidualCalc(su2double *val_residual, su2double **J

/*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
if (implicit) {
Jacobian_i[0][0] = -diff_kine*proj_vector_ij/Density_i; Jacobian_i[0][1] = 0.0;
Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_vector_ij/Density_i;
Jacobian_i[0][0] = -diff_kine*proj_vector_ij; Jacobian_i[0][1] = 0.0;
Jacobian_i[1][0] = 0.0; Jacobian_i[1][1] = -diff_omega*proj_vector_ij;

Jacobian_j[0][0] = diff_kine*proj_vector_ij/Density_j; Jacobian_j[0][1] = 0.0;
Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_vector_ij/Density_j;
Jacobian_j[0][0] = diff_kine*proj_vector_ij; Jacobian_j[0][1] = 0.0;
Jacobian_j[1][0] = 0.0; Jacobian_j[1][1] = diff_omega*proj_vector_ij;
}

}
Expand Down Expand Up @@ -1272,10 +1245,10 @@ void CSourcePieceWise_TurbSST::ComputeResidual(su2double *val_residual, su2doubl

/*--- Implicit part ---*/

val_Jacobian_i[0][0] = -beta_star*TurbVar_i[1]*Volume;
val_Jacobian_i[0][1] = -beta_star*TurbVar_i[0]*Volume;
val_Jacobian_i[0][0] = -beta_star*Density_i*TurbVar_i[1]*Volume;
val_Jacobian_i[0][1] = -beta_star*Density_i*TurbVar_i[0]*Volume;
val_Jacobian_i[1][0] = 0.0;
val_Jacobian_i[1][1] = -2.0*beta_blended*TurbVar_i[1]*Volume;
val_Jacobian_i[1][1] = -2.0*beta_blended*Density_i*TurbVar_i[1]*Volume;
}
Copy link
Contributor

@jayantmukho jayantmukho Jul 8, 2019

Choose a reason for hiding this comment

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

Damn, I've been staring at these lines of code for a while and I never realized that Density was missing from the Jacobian computation. Good catch!

Copy link
Member Author

@pcarruscag pcarruscag Jul 8, 2019

Choose a reason for hiding this comment

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

It was not missing, before the Jacobians were w.r.t. rhok and rhow, that is why in the diffusion part they had to be divided by density, I changed it to k and w (only for the SST model) since the mass flux already contains rho. I am not 100% sure this is correct though... In fact CTurbSolver updates the solution via node[iPoint]->AddConservativeSolution. So a better catch by you actually!


AD::SetPreaccOut(val_residual, nVar);
Expand Down
24 changes: 19 additions & 5 deletions SU2_CFD/src/solver_direct_mean.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,6 @@ CEulerSolver::CEulerSolver(void) : CSolver() {
CkInflow = NULL;
CkOutflow1 = NULL;
CkOutflow2 = NULL;

}

CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CSolver() {
Expand Down Expand Up @@ -887,7 +886,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, unsigned short

InitiateComms(geometry, config, SOLUTION);
CompleteComms(geometry, config, SOLUTION);

}

CEulerSolver::~CEulerSolver(void) {
Expand Down Expand Up @@ -3616,6 +3614,7 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta
bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
bool jst_scheme = ((config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0));
bool grid_movement = config->GetGrid_Movement();
bool rans = ((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS));

for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) {

Expand Down Expand Up @@ -3646,9 +3645,10 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta
numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[jPoint]->GetGridVel());
}

/*--- Compute residuals, and Jacobians ---*/
/*--- Compute residuals, and Jacobians. Store mass flux if required. ---*/

numerics->ComputeResidual(Res_Conv, Jacobian_i, Jacobian_j, config);
if (rans && (iMesh == MESH_0)) EdgeMassFluxes[iEdge] = Res_Conv[0];

/*--- Update convective and artificial dissipation residuals ---*/

Expand Down Expand Up @@ -3689,6 +3689,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE;
bool low_mach_corr = config->Low_Mach_Correction();
unsigned short kind_dissipation = config->GetKind_RoeLowDiss();
bool rans = ((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS));

/*--- Loop over all the edges ---*/

Expand Down Expand Up @@ -3870,9 +3871,10 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
}
}

/*--- Compute the residual ---*/
/*--- Compute the residual and store the mass flux if required ---*/

numerics->ComputeResidual(Res_Conv, Jacobian_i, Jacobian_j, config);
if (rans && (iMesh == MESH_0)) EdgeMassFluxes[iEdge] = Res_Conv[0];

/*--- Update residual value ---*/

Expand Down Expand Up @@ -14524,7 +14526,9 @@ CNSSolver::CNSSolver(void) : CEulerSolver() {
DonorPrimVar = NULL; DonorGlobalIndex = NULL;

HeatConjugateVar = NULL;


/*--- Edge mass fluxes ---*/
EdgeMassFluxes = NULL;
}

CNSSolver::CNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CEulerSolver() {
Expand Down Expand Up @@ -15305,6 +15309,15 @@ CNSSolver::CNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh)
InitiateComms(geometry, config, SOLUTION);
CompleteComms(geometry, config, SOLUTION);

/*--- Allocate and initialize edge mass flux array if RANS and for the finest grid only ---*/

if (rans && (iMesh == MESH_0)) {

EdgeMassFluxes = new su2double [geometry->GetnEdge()];

for(unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); ++iEdge)
EdgeMassFluxes[iEdge] = 0.0;
}
}

CNSSolver::~CNSSolver(void) {
Expand Down Expand Up @@ -15376,6 +15389,7 @@ CNSSolver::~CNSSolver(void) {
delete [] Buffet_Sensor;
}

if (EdgeMassFluxes != NULL) delete [] EdgeMassFluxes;
}

void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) {
Expand Down
21 changes: 19 additions & 2 deletions SU2_CFD/src/solver_direct_mean_inc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1839,6 +1839,7 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
bool jst_scheme = ((config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0));
bool grid_movement = config->GetGrid_Movement();
bool rans = ((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS));

for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) {

Expand Down Expand Up @@ -1869,9 +1870,10 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[jPoint]->GetGridVel());
}

/*--- Compute residuals, and Jacobians ---*/
/*--- Compute residuals, and Jacobians. Store mass flux if required. ---*/

numerics->ComputeResidual(Res_Conv, Jacobian_i, Jacobian_j, config);
if (rans && (iMesh == MESH_0)) EdgeMassFluxes[iEdge] = Res_Conv[0];

/*--- Update convective and artificial dissipation residuals ---*/

Expand Down Expand Up @@ -1905,6 +1907,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter());
bool grid_movement = config->GetGrid_Movement();
bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE;
bool rans = ((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS));

/*--- Loop over all the edges ---*/

Expand Down Expand Up @@ -1978,9 +1981,10 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont

}

/*--- Compute the residual ---*/
/*--- Compute the residual and store the mass flux if required ---*/

numerics->ComputeResidual(Res_Conv, Jacobian_i, Jacobian_j, config);
if (rans && (iMesh == MESH_0)) EdgeMassFluxes[iEdge] = Res_Conv[0];

/*--- Update residual value ---*/

Expand Down Expand Up @@ -6802,6 +6806,8 @@ CIncNSSolver::CIncNSSolver(void) : CIncEulerSolver() {
SlidingState = NULL;
SlidingStateNodes = NULL;

/*--- Edge mass fluxes ---*/
EdgeMassFluxes = NULL;
}

CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CIncEulerSolver() {
Expand All @@ -6820,6 +6826,7 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short
string filename_ = config->GetSolution_FlowFileName();

unsigned short direct_diff = config->GetDirectDiff();
bool rans = ((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS));

/*--- Store the multigrid level. ---*/
MGLevel = iMesh;
Expand Down Expand Up @@ -7375,7 +7382,16 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short

InitiateComms(geometry, config, SOLUTION);
CompleteComms(geometry, config, SOLUTION);

/*--- Allocate and initialize edge mass flux array if RANS and for the finest grid only ---*/

if (rans && (iMesh == MESH_0)) {

EdgeMassFluxes = new su2double [geometry->GetnEdge()];

for(unsigned long iEdge = 0; iEdge < geometry->GetnEdge(); ++iEdge)
EdgeMassFluxes[iEdge] = 0.0;
}
}

CIncNSSolver::~CIncNSSolver(void) {
Expand Down Expand Up @@ -7440,6 +7456,7 @@ CIncNSSolver::~CIncNSSolver(void) {
delete [] HeatConjugateVar;
}

if (EdgeMassFluxes != NULL) delete [] EdgeMassFluxes;
}

void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) {
Expand Down
Loading