Skip to content

Commit 8e75b33

Browse files
WallyMaierWallyMaier
committed
addressing PR comments, LGTM alerts
1 parent 88e9e9f commit 8e75b33

8 files changed

Lines changed: 105 additions & 77 deletions

File tree

‎SU2_CFD/include/numerics/CNumerics.hpp‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1400,6 +1400,32 @@ class CNumerics {
14001400
* \param[out] val_Jacobian_i - Jacobian of the source terms
14011401
*/
14021402
inline virtual ResidualType<> ComputeChemistry(const CConfig* config) { return ResidualType<>(nullptr,nullptr,nullptr); }
1403+
1404+
/*!
1405+
* \brief Check if residual constains a NaN value
1406+
* \param[in] config - Definition of the particular problem.
1407+
* \param[in] val_residual - residual of the numeric function.
1408+
* \param[out] ERR - Presencse of NaN in vector
1409+
*/
1410+
static bool CheckResidualNaNs(bool implicit, int nVar, ResidualType<>(residual)){
1411+
1412+
bool ERR = false;
1413+
bool jac_j = false;
1414+
if (residual.jacobian_j != nullptr) jac_j =true;
1415+
1416+
for (auto iVar = 0; iVar<nVar; iVar++){
1417+
if (residual[iVar] != residual[iVar]) ERR = true;
1418+
1419+
if (implicit) {
1420+
for (auto jVar = 0; jVar < nVar; jVar++){
1421+
if (residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) ERR = true;
1422+
if ((jac_j) && (residual.jacobian_j[iVar][jVar] != residual.jacobian_j[iVar][jVar])) ERR = true;
1423+
}
1424+
}
1425+
}
1426+
return ERR;
1427+
}
1428+
14031429
/*!
14041430
* \brief Set intermittency for numerics (used in SA with LM transition model)
14051431
*/

‎SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -308,7 +308,7 @@ class CFVMFlowSolverBase : public CSolver {
308308
* \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over a nonlinear
309309
* iteration for stability.
310310
*/
311-
void ComputeUnderRelaxationFactor(const CConfig* config);
311+
virtual void ComputeUnderRelaxationFactor(const CConfig* config);
312312

313313
/*!
314314
* \brief General implementation to load a flow solution from a restart file.

‎SU2_CFD/include/solvers/CFVMFlowSolverBase.inl‎

Lines changed: 1 addition & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -540,41 +540,17 @@ void CFVMFlowSolverBase<V, R>::ComputeUnderRelaxationFactor(const CConfig* confi
540540
for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
541541
su2double localUnderRelaxation = 1.0;
542542

543-
su2double num = 0.0;
544-
su2double denom = 0.0;
545-
546543
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
547544
/* We impose a limit on the maximum percentage that the
548545
density and energy can change over a nonlinear iteration. */
549546

550-
unsigned short tmp = 1; //DELETE ME TODO
551-
if (config->GetnSpecies() != 0) tmp = 2;
552-
if ((tmp!=2)&&((iVar == 0) || (iVar >= (nVar - tmp)))) {
547+
if ((iVar == 0) || (iVar >= nVar - 1)) {
553548
const unsigned long index = iPoint * nVar + iVar;
554549
su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
555550
if (ratio > allowableRatio) {
556551
localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
557552
}
558553
}
559-
if (tmp==2){
560-
const unsigned long index = iPoint * nVar + iVar;
561-
if (iVar < config->GetnSpecies()) {
562-
num += fabs(LinSysSol[index]);
563-
denom += fabs(nodes->GetSolution(iPoint, iVar));
564-
if (iVar == (config ->GetnSpecies()-1)){
565-
su2double ratio = (num/(denom+EPS));
566-
if (ratio > allowableRatio) {
567-
localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
568-
}
569-
}
570-
}
571-
if (iVar == (nVar-tmp)){
572-
su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
573-
if (ratio > allowableRatio) {
574-
localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
575-
}
576-
}
577-
}
578554
}
579555

580556
/* Threshold the relaxation factor in the event that there is

‎SU2_CFD/include/solvers/CNEMOEulerSolver.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -199,6 +199,13 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
199199
virtual unsigned long SetPrimitive_Variables(CSolver **solver_container,
200200
CConfig *config, bool Output);
201201

202+
/*!
203+
* \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
204+
* a nonlinear iteration for stability.
205+
* \param[in] config - Definition of the particular problem.
206+
*/
207+
void ComputeUnderRelaxationFactor(const CConfig *config) override;
208+
202209
/*!
203210
* \brief Set the fluid solver nondimensionalization.
204211
* \param[in] geometry - Geometrical definition of the problem.

‎SU2_CFD/src/numerics/NEMO/convection/ausm.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -209,7 +209,7 @@ CNumerics::ResidualType<> CUpwAUSM_NEMO::ComputeResidual(const CConfig *config)
209209
if (mF >= 0) {
210210

211211
/*--- Jacobian contribution: dFc terms ---*/
212-
for (unsigned short iVar = 0; iVar < nSpecies+nDim; iVar++) {
212+
for (auto iVar = 0u; iVar < nSpecies+nDim; iVar++) {
213213
for (unsigned short jVar = 0; jVar < nVar; jVar++) {
214214
Jacobian_i[iVar][jVar] += mF * FcL[iVar]/a_i * daL[jVar];
215215
}
@@ -292,7 +292,7 @@ CNumerics::ResidualType<> CUpwAUSM_NEMO::ComputeResidual(const CConfig *config)
292292
if (mF < 0) {
293293

294294
/*--- Jacobian contribution: dFc terms ---*/
295-
for (unsigned short iVar = 0; iVar < nSpecies+nDim; iVar++) {
295+
for (auto iVar = 0u; iVar < nSpecies+nDim; iVar++) {
296296
for (unsigned short jVar = 0; jVar < nVar; jVar++) {
297297
Jacobian_j[iVar][jVar] += mF * FcR[iVar]/a_j * daR[jVar];
298298
}

‎SU2_CFD/src/numerics/NEMO/convection/ausmplusup2.cpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -69,6 +69,8 @@ CUpwAUSMPLUSUP2_NEMO::~CUpwAUSMPLUSUP2_NEMO(void) {
6969
delete [] dmRM;
7070
delete [] dpLP;
7171
delete [] dpRM;
72+
delete [] daL;
73+
delete [] daR;
7274
delete [] rhos_i;
7375
delete [] rhos_j;
7476
delete [] u_i;
@@ -254,7 +256,7 @@ CNumerics::ResidualType<> CUpwAUSMPLUSUP2_NEMO::ComputeResidual(const CConfig *c
254256
if (mF >= 0) {
255257

256258
/*--- Jacobian contribution: dFc terms ---*/
257-
for (iVar = 0; iVar < nSpecies+nDim; iVar++) {
259+
for (auto iVar = 0u; iVar < nSpecies+nDim; iVar++) {
258260
for (jVar = 0; jVar < nVar; jVar++) {
259261
Jacobian_i[iVar][jVar] += mF * FcL[iVar]/aF * daL[jVar];
260262
}
@@ -337,7 +339,7 @@ CNumerics::ResidualType<> CUpwAUSMPLUSUP2_NEMO::ComputeResidual(const CConfig *c
337339
if (mF < 0) {
338340

339341
/*--- Jacobian contribution: dFc terms ---*/
340-
for (iVar = 0; iVar < nSpecies+nDim; iVar++) {
342+
for (auto iVar = 0u; iVar < nSpecies+nDim; iVar++) {
341343
for (jVar = 0; jVar < nVar; jVar++) {
342344
Jacobian_j[iVar][jVar] += mF * FcR[iVar]/aF * daR[jVar];
343345
}

‎SU2_CFD/src/numerics/NEMO/convection/lax.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,7 @@ CNumerics::ResidualType<> CCentLax_NEMO::ComputeResidual(const CConfig *config)
139139

140140
cte = Epsilon_0*StretchingFactor*MeanLambda;
141141

142-
for (unsigned short iVar = 0; iVar < nSpecies+nDim; iVar++) {
142+
for (auto iVar = 0u; iVar < nSpecies+nDim; iVar++) {
143143
Jacobian_i[iVar][iVar] += cte;
144144
Jacobian_j[iVar][iVar] -= cte;
145145
}

‎SU2_CFD/src/solvers/CNEMOEulerSolver.cpp‎

Lines changed: 63 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -417,8 +417,6 @@ void CNEMOEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) {
417417
void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
418418
CConfig *config, unsigned short iMesh, unsigned short iRKStep) {
419419
unsigned long iEdge, iPoint, jPoint;
420-
unsigned short iVar, jVar;
421-
bool err;
422420

423421
CNumerics* numerics = numerics_container[CONV_TERM];
424422

@@ -452,17 +450,8 @@ void CNEMOEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_c
452450
/*--- Compute residuals, and Jacobians ---*/
453451
auto residual = numerics->ComputeResidual(config);
454452

455-
/*--- Check for NaNs before applying the residual to the linear system ---*/
456-
err = false;
457-
for (iVar = 0; iVar < nVar; iVar++)
458-
if (residual[iVar] != residual[iVar])
459-
err = true;
460-
if (implicit)
461-
for (iVar = 0; iVar < nVar; iVar++)
462-
for (jVar = 0; jVar < nVar; jVar++)
463-
if ((residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) ||
464-
(residual.jacobian_j[iVar][jVar] != residual.jacobian_j[iVar][jVar]))
465-
err = true;
453+
/*--- Check residuals/Jacobians --*/
454+
bool err = CNumerics::CheckResidualNaNs(implicit, nVar, residual);
466455

467456
/*--- Update the residual and Jacobian ---*/
468457
if (!err) {
@@ -614,15 +603,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
614603
auto residual = numerics->ComputeResidual(config);
615604

616605
/*--- Check for NaNs before applying the residual to the linear system ---*/
617-
bool err = false;
618-
for (iVar = 0; iVar < nVar; iVar++)
619-
if (residual[iVar] != residual[iVar]) err = true;
620-
if (implicit)
621-
for (auto iVar = 0u; iVar < nVar; iVar++)
622-
for (auto jVar = 0u; jVar < nVar; jVar++)
623-
if ((residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) ||
624-
(residual.jacobian_j[iVar][jVar] != residual.jacobian_j[iVar][jVar]) )
625-
err = true;
606+
bool err = CNumerics::CheckResidualNaNs(implicit, nVar, residual);
626607

627608
/*--- Update the residual and Jacobian ---*/
628609
if (!err) {
@@ -761,7 +742,7 @@ bool CNEMOEulerSolver::CheckNonPhys(const su2double *V) const {
761742

762743
void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
763744

764-
unsigned short iVar, jVar;
745+
unsigned short iVar;
765746
unsigned long iPoint;
766747

767748
/*--- Assign booleans ---*/
@@ -817,13 +798,7 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
817798
auto residual = numerics->ComputeChemistry(config);
818799

819800
/*--- Check for errors before applying source to the linear system ---*/
820-
err = false;
821-
for (iVar = 0; iVar < nVar; iVar++)
822-
if (residual[iVar] != residual[iVar]) err = true;
823-
if (implicit)
824-
for (iVar = 0; iVar < nVar; iVar++)
825-
for (jVar = 0; jVar < nVar; jVar++)
826-
if (residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) err = true;
801+
err = CNumerics::CheckResidualNaNs(implicit, nVar, residual);
827802

828803
/*--- Apply the chemical sources to the linear system ---*/
829804
if (!err) {
@@ -842,13 +817,7 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
842817
auto residual = numerics->ComputeVibRelaxation(config);
843818

844819
/*--- Check for errors before applying source to the linear system ---*/
845-
err = false;
846-
for (iVar = 0; iVar < nVar; iVar++)
847-
if (residual[iVar] != residual[iVar]) err = true;
848-
if (implicit)
849-
for (iVar = 0; iVar < nVar; iVar++)
850-
for (jVar = 0; jVar < nVar; jVar++)
851-
if (residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) err = true;
820+
err = CNumerics::CheckResidualNaNs(implicit, nVar, residual);
852821

853822
/*--- Apply the vibrational relaxation terms to the linear system ---*/
854823
if (!err) {
@@ -896,13 +865,7 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
896865
auto residual = numerics->ComputeAxisymmetric(config);
897866

898867
/*--- Check for errors before applying source to the linear system ---*/
899-
err = false;
900-
for (iVar = 0; iVar < nVar; iVar++)
901-
if (residual[iVar] != residual[iVar]) err = true;
902-
if (implicit)
903-
for (iVar = 0; iVar < nVar; iVar++)
904-
for (jVar = 0; jVar < nVar; jVar++)
905-
if (residual.jacobian_i[iVar][jVar] != residual.jacobian_i[iVar][jVar]) err = true;
868+
err = CNumerics::CheckResidualNaNs(implicit, nVar, residual);
906869

907870
/*--- Apply the update to the linear system ---*/
908871
if (!err) {
@@ -965,6 +928,60 @@ void CNEMOEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**,
965928
CompleteImplicitIteration_impl<true>(geometry, config);
966929
}
967930

931+
void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) {
932+
933+
/* Loop over the solution update given by relaxing the linear
934+
system for this nonlinear iteration. */
935+
936+
const su2double allowableRatio = 0.2;
937+
938+
SU2_OMP_FOR_STAT(omp_chunk_size)
939+
for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
940+
su2double localUnderRelaxation = 1.0;
941+
942+
su2double num = 0.0;
943+
su2double denom = 0.0;
944+
945+
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
946+
/* We impose a limit on the maximum percentage that the
947+
density (sum of all species) and energy can change over a nonlinear iteration. */
948+
949+
const unsigned long index = iPoint * nVar + iVar;
950+
if (iVar < config->GetnSpecies()) {
951+
num += fabs(LinSysSol[index]);
952+
denom += fabs(nodes->GetSolution(iPoint, iVar));
953+
954+
/*--- If final density/species, compute Under-relaxation ---*/
955+
if (iVar == (config ->GetnSpecies()-1)){
956+
su2double ratio = (num/(denom+EPS));
957+
if (ratio > allowableRatio) {
958+
localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
959+
}
960+
}
961+
962+
/*--- Energy ---*/
963+
if (iVar == (nVar-2)){
964+
su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
965+
if (ratio > allowableRatio) {
966+
localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
967+
}
968+
}
969+
}
970+
}
971+
972+
/* Threshold the relaxation factor in the event that there is
973+
a very small value. This helps avoid catastrophic crashes due
974+
to non-realizable states by canceling the update. */
975+
976+
if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0;
977+
978+
/* Store the under-relaxation factor for this point. */
979+
980+
nodes->SetUnderRelaxation(iPoint, localUnderRelaxation);
981+
}
982+
END_SU2_OMP_FOR
983+
}
984+
968985
void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMesh) {
969986

970987
su2double
@@ -2327,7 +2344,7 @@ SU2_MPI::Error("BC_SUPERSONIC_INLET: Not operational in NEMO.", CURRENT_FUNCTION
23272344
//
23282345
// /*--- Jacobian contribution for implicit integration ---*/
23292346
// //if (implicit)
2330-
// // Jacobian.AddBlock2Daig(iPoint, residual.jacobian_i);
2347+
// // Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
23312348
//
23322349
// /*--- Viscous contribution ---*/
23332350
// if (viscous) {
@@ -2352,7 +2369,7 @@ SU2_MPI::Error("BC_SUPERSONIC_INLET: Not operational in NEMO.", CURRENT_FUNCTION
23522369
//
23532370
// /*--- Jacobian contribution for implicit integration ---*/
23542371
// //if (implicit)
2355-
// // Jacobian.SubtractBlock(iPoint, residual.jacobian_i);
2372+
// // Jacobian.SubtractBlock2Diag(iPoint, residual.Jacobian_i);
23562373
// }
23572374
//
23582375
// }

0 commit comments

Comments
 (0)