Skip to content

Commit 0a219b1

Browse files
pcarruscagpcarruscag
authored
Merge pull request #1224 from su2code/some_adjoint_fixes
Small adjoint fixes
2 parents 048394c + 861285d commit 0a219b1

File tree

10 files changed

+175
-439
lines changed

10 files changed

+175
-439
lines changed

‎Common/src/linear_algebra/CSysSolve.cpp‎

Lines changed: 24 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -835,6 +835,12 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType> & Jacobian, co
835835
#endif
836836
}
837837

838+
SU2_OMP_MASTER
839+
if (KindSolver == RESTARTED_FGMRES) {
840+
xIsZero = false;
841+
tol_type = LinearToleranceType::ABSOLUTE;
842+
}
843+
838844
/*--- Create matrix-vector product, preconditioner, and solve the linear system ---*/
839845

840846
HandleTemporariesIn(LinSysRes, LinSysSol);
@@ -886,9 +892,10 @@ unsigned long CSysSolve<ScalarType>::Solve(CSysMatrix<ScalarType> & Jacobian, co
886892
norm0 = LinSysRes_ptr->norm();
887893
while (IterLinSol < MaxIter) {
888894
/*--- Enforce a hard limit on total number of iterations ---*/
889-
unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol);
890-
IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, IterLimit, residual, ScreenOutput, config);
891-
if ( residual <= SolverTol*norm0 ) break;
895+
auto IterLimit = min(RestartIter, MaxIter-IterLinSol);
896+
auto iter = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, IterLimit, residual, ScreenOutput, config);
897+
IterLinSol += iter;
898+
if (residual <= SolverTol*norm0 || iter < IterLimit) break;
892899
}
893900
break;
894901
case SMOOTHER:
@@ -1018,6 +1025,12 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType> & Jacobian,
10181025

10191026
/*--- Solve the system ---*/
10201027

1028+
SU2_OMP_MASTER
1029+
if (KindSolver == RESTARTED_FGMRES) {
1030+
xIsZero = false;
1031+
tol_type = LinearToleranceType::ABSOLUTE;
1032+
}
1033+
10211034
HandleTemporariesIn(LinSysRes, LinSysSol);
10221035

10231036
switch(KindSolver) {
@@ -1035,19 +1048,23 @@ unsigned long CSysSolve<ScalarType>::Solve_b(CSysMatrix<ScalarType> & Jacobian,
10351048
Norm0 = LinSysRes_ptr->norm();
10361049
while (IterLinSol < MaxIter) {
10371050
/*--- Enforce a hard limit on total number of iterations ---*/
1038-
unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol);
1039-
IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol , IterLimit, Residual, ScreenOutput, config);
1040-
if ( Residual <= SolverTol*Norm0 ) break;
1051+
auto IterLimit = min(RestartIter, MaxIter-IterLinSol);
1052+
auto iter = FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, IterLimit, Residual, ScreenOutput, config);
1053+
IterLinSol += iter;
1054+
if (Residual <= SolverTol*Norm0 || iter < IterLimit) break;
10411055
}
10421056
break;
1057+
case SMOOTHER:
1058+
IterLinSol = Smoother_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, MaxIter, Residual, ScreenOutput, config);
1059+
break;
10431060
case PASTIX_LDLT : case PASTIX_LU:
10441061
Jacobian.BuildPastixPreconditioner(geometry, config, KindSolver, RequiresTranspose);
10451062
Jacobian.ComputePastixPreconditioner(*LinSysRes_ptr, *LinSysSol_ptr, geometry, config);
10461063
IterLinSol = 1;
10471064
Residual = 1e-20;
10481065
break;
10491066
default:
1050-
SU2_MPI::Error("The specified linear solver is not yet implemented for the discrete adjoint method.", CURRENT_FUNCTION);
1067+
SU2_MPI::Error("Unknown type of linear solver.",CURRENT_FUNCTION);
10511068
break;
10521069
}
10531070

‎SU2_CFD/include/drivers/CDriver.hpp‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -699,20 +699,30 @@ class CDriver {
699699
* \param[in] iZone - Index of the zone.
700700
* \param[in] adjoint - True to consider adjoint solvers instead of primal.
701701
* \param[in] solution - Solution object with interface (iPoint,iVar).
702+
* \tparam Old - If true set "old solutions" instead.
702703
*/
703-
template<class Container>
704+
template<class Container, bool Old = false>
704705
void SetAllSolutions(unsigned short iZone, bool adjoint, const Container& solution) {
705706
const auto nPoint = geometry_container[iZone][INST_0][MESH_0]->GetnPoint();
706707
for (auto iSol = 0u, offset = 0u; iSol < MAX_SOLS; ++iSol) {
707708
auto solver = solver_container[iZone][INST_0][MESH_0][iSol];
708709
if (!(solver && (solver->GetAdjoint() == adjoint))) continue;
709710
for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint)
710711
for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar)
711-
solver->GetNodes()->SetSolution(iPoint, iVar, solution(iPoint,offset+iVar));
712+
if (!Old) solver->GetNodes()->SetSolution(iPoint, iVar, solution(iPoint,offset+iVar));
713+
else solver->GetNodes()->SetSolution_Old(iPoint, iVar, solution(iPoint,offset+iVar));
712714
offset += solver->GetnVar();
713715
}
714716
}
715717

718+
/*!
719+
* \brief Set the "old solution" of all solvers (adjoint or primal) in a zone.
720+
*/
721+
template<class Container>
722+
void SetAllSolutionsOld(unsigned short iZone, bool adjoint, const Container& solution) {
723+
SetAllSolutions<Container,true>(iZone, adjoint, solution);
724+
}
725+
716726
/*!
717727
* \brief Get the solution of all solvers (adjoint or primal) in a zone.
718728
* \param[in] iZone - Index of the zone.

‎SU2_CFD/include/numerics_simd/flow/convection/common.hpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -79,13 +79,13 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
7979
auto grad_j = gatherVariables<nVar,nDim>(jPoint, gradient);
8080

8181
for (size_t iVar = 0; iVar < nVar; ++iVar) {
82-
auto proj_i = dot(grad_i[iVar], vector_ij);
83-
auto proj_j = dot(grad_j[iVar], vector_ij);
84-
auto delta_ij = V.j.all(iVar) - V.i.all(iVar);
85-
auto delta_ij_2 = pow(delta_ij, 2) + 1e-6;
82+
const Double proj_i = dot(grad_i[iVar], vector_ij);
83+
const Double proj_j = dot(grad_j[iVar], vector_ij);
84+
const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
85+
const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
8686
/// TODO: Customize the limiter function.
87-
auto lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
88-
auto lim_j = (delta_ij_2 + proj_j*delta_ij) / (pow(proj_j,2) + delta_ij_2);
87+
const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
88+
const Double lim_j = (delta_ij_2 + proj_j*delta_ij) / (pow(proj_j,2) + delta_ij_2);
8989
V.i.all(iVar) += lim_i * 0.5 * proj_i;
9090
V.j.all(iVar) -= lim_j * 0.5 * proj_j;
9191
}

‎SU2_CFD/include/solvers/CEulerSolver.hpp‎

Lines changed: 41 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -106,50 +106,45 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
106106
unsigned long Iter_Update_AoA = 0; /*!< \brief Iteration at which AoA was updated last */
107107
su2double dCL_dAlpha; /*!< \brief Value of dCL_dAlpha used to control CL in fixed CL mode */
108108
unsigned long BCThrust_Counter;
109-
unsigned short nSpanWiseSections; /*!< \brief Number of span-wise sections. */
110-
unsigned short nSpanMax; /*!< \brief Max number of maximum span-wise sections for all zones. */
111-
unsigned short nMarkerTurboPerf; /*!< \brief Number of turbo performance. */
112109

113110
vector<CFluidModel*> FluidModel; /*!< \brief fluid model used in the solver. */
114111

115112
/*--- Turbomachinery Solver Variables ---*/
116113

117-
su2double ***AverageFlux = nullptr,
118-
***SpanTotalFlux = nullptr,
119-
***AverageVelocity = nullptr,
120-
***AverageTurboVelocity = nullptr,
121-
***OldAverageTurboVelocity = nullptr,
122-
***ExtAverageTurboVelocity = nullptr,
123-
**AveragePressure = nullptr,
124-
**OldAveragePressure = nullptr,
125-
**RadialEquilibriumPressure = nullptr,
126-
**ExtAveragePressure = nullptr,
127-
**AverageDensity = nullptr,
128-
**OldAverageDensity = nullptr,
129-
**ExtAverageDensity = nullptr,
130-
**AverageNu = nullptr,
131-
**AverageKine = nullptr,
132-
**AverageOmega = nullptr,
133-
**ExtAverageNu = nullptr,
134-
**ExtAverageKine = nullptr,
135-
**ExtAverageOmega = nullptr;
136-
137-
su2double **DensityIn = nullptr,
138-
**PressureIn = nullptr,
139-
***TurboVelocityIn = nullptr,
140-
**DensityOut = nullptr,
141-
**PressureOut = nullptr,
142-
***TurboVelocityOut = nullptr,
143-
**KineIn = nullptr,
144-
**OmegaIn = nullptr,
145-
**NuIn = nullptr,
146-
**KineOut = nullptr,
147-
**OmegaOut = nullptr,
148-
**NuOut = nullptr;
149-
150-
complex<su2double> ***CkInflow = nullptr,
151-
***CkOutflow1 = nullptr,
152-
***CkOutflow2 = nullptr;
114+
vector<su2activematrix> AverageFlux;
115+
vector<su2activematrix> SpanTotalFlux;
116+
vector<su2activematrix> AverageVelocity;
117+
vector<su2activematrix> AverageTurboVelocity;
118+
vector<su2activematrix> OldAverageTurboVelocity;
119+
vector<su2activematrix> ExtAverageTurboVelocity;
120+
su2activematrix AveragePressure;
121+
su2activematrix OldAveragePressure;
122+
su2activematrix RadialEquilibriumPressure;
123+
su2activematrix ExtAveragePressure;
124+
su2activematrix AverageDensity;
125+
su2activematrix OldAverageDensity;
126+
su2activematrix ExtAverageDensity;
127+
su2activematrix AverageNu;
128+
su2activematrix AverageKine;
129+
su2activematrix AverageOmega;
130+
su2activematrix ExtAverageNu;
131+
su2activematrix ExtAverageKine;
132+
su2activematrix ExtAverageOmega;
133+
134+
su2activematrix DensityIn;
135+
su2activematrix PressureIn;
136+
vector<su2activematrix> TurboVelocityIn;
137+
su2activematrix DensityOut;
138+
su2activematrix PressureOut;
139+
vector<su2activematrix> TurboVelocityOut;
140+
su2activematrix KineIn;
141+
su2activematrix OmegaIn;
142+
su2activematrix NuIn;
143+
su2activematrix KineOut;
144+
su2activematrix OmegaOut;
145+
su2activematrix NuOut;
146+
147+
vector<su2matrix<complex<su2double> > > CkInflow, CkOutflow1, CkOutflow2;
153148

154149
/*--- End of Turbomachinery Solver Variables ---*/
155150

@@ -1219,7 +1214,7 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
12191214
* \param[in] val_marker - bound marker.
12201215
* \return Value of the Average Total Pressure on the surface <i>val_marker</i>.
12211216
*/
1222-
inline su2double* GetAverageTurboVelocity(unsigned short valMarker, unsigned short valSpan) const final {
1217+
inline const su2double* GetAverageTurboVelocity(unsigned short valMarker, unsigned short valSpan) const final {
12231218
return AverageTurboVelocity[valMarker][valSpan];
12241219
}
12251220

@@ -1372,7 +1367,7 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
13721367
* \param[in] inMarkerTP - bound marker.
13731368
* \return Value of the inlet normal velocity.
13741369
*/
1375-
inline su2double* GetTurboVelocityIn(unsigned short inMarkerTP, unsigned short valSpan) const final {
1370+
inline const su2double* GetTurboVelocityIn(unsigned short inMarkerTP, unsigned short valSpan) const final {
13761371
return TurboVelocityIn[inMarkerTP][valSpan];
13771372
}
13781373

@@ -1399,7 +1394,7 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
13991394
* \param[in] inMarkerTP - bound marker.
14001395
* \return Value of the outlet normal velocity.
14011396
*/
1402-
inline su2double* GetTurboVelocityOut(unsigned short inMarkerTP, unsigned short valSpan) const final {
1397+
inline const su2double* GetTurboVelocityOut(unsigned short inMarkerTP, unsigned short valSpan) const final {
14031398
return TurboVelocityOut[inMarkerTP][valSpan];
14041399
}
14051400

@@ -1484,12 +1479,10 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
14841479
* \param[in] value - turboperformance value to set.
14851480
* \param[in] inMarkerTP - turboperformance marker.
14861481
*/
1487-
inline void SetTurboVelocityIn(su2double *value,
1482+
inline void SetTurboVelocityIn(const su2double *value,
14881483
unsigned short inMarkerTP,
14891484
unsigned short valSpan) final {
1490-
unsigned short iDim;
1491-
1492-
for(iDim = 0; iDim < nDim; iDim++)
1485+
for(unsigned short iDim = 0; iDim < nDim; iDim++)
14931486
TurboVelocityIn[inMarkerTP][valSpan][iDim] = value[iDim];
14941487
}
14951488

@@ -1520,12 +1513,10 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
15201513
* \param[in] value - turboperformance value to set.
15211514
* \param[in] inMarkerTP - turboperformance marker.
15221515
*/
1523-
inline void SetTurboVelocityOut(su2double *value,
1516+
inline void SetTurboVelocityOut(const su2double *value,
15241517
unsigned short inMarkerTP,
15251518
unsigned short valSpan) final {
1526-
unsigned short iDim;
1527-
1528-
for(iDim = 0; iDim < nDim; iDim++)
1519+
for(unsigned short iDim = 0; iDim < nDim; iDim++)
15291520
TurboVelocityOut[inMarkerTP][valSpan][iDim] = value[iDim];
15301521
}
15311522

‎SU2_CFD/include/solvers/CSolver.hpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -3983,7 +3983,7 @@ class CSolver {
39833983
* \param[in] val_marker - bound marker.
39843984
* \return Value of the Average Total Pressure on the surface <i>val_marker</i>.
39853985
*/
3986-
inline virtual su2double* GetAverageTurboVelocity(unsigned short valMarker, unsigned short valSpan) const { return nullptr; }
3986+
inline virtual const su2double* GetAverageTurboVelocity(unsigned short valMarker, unsigned short valSpan) const { return nullptr; }
39873987

39883988
/*!
39893989
* \brief A virtual member.
@@ -4101,7 +4101,7 @@ class CSolver {
41014101
* \param[in] inMarkerTP - bound marker.
41024102
* \return Value of the inlet normal velocity.
41034103
*/
4104-
inline virtual su2double* GetTurboVelocityIn(unsigned short inMarkerTP, unsigned short valSpan) const {return nullptr;}
4104+
inline virtual const su2double* GetTurboVelocityIn(unsigned short inMarkerTP, unsigned short valSpan) const {return nullptr;}
41054105

41064106
/*!
41074107
* \brief A virtual member.
@@ -4122,7 +4122,7 @@ class CSolver {
41224122
* \param[in] inMarkerTP - bound marker.
41234123
* \return Value of the outlet normal velocity.
41244124
*/
4125-
inline virtual su2double* GetTurboVelocityOut(unsigned short inMarkerTP, unsigned short valSpan) const {return nullptr;}
4125+
inline virtual const su2double* GetTurboVelocityOut(unsigned short inMarkerTP, unsigned short valSpan) const {return nullptr;}
41264126

41274127
/*!
41284128
* \brief A virtual member.
@@ -4189,7 +4189,7 @@ class CSolver {
41894189
* \param[in] value - turboperformance value to set.
41904190
* \param[in] inMarkerTP - turboperformance marker.
41914191
*/
4192-
inline virtual void SetTurboVelocityIn(su2double *value,
4192+
inline virtual void SetTurboVelocityIn(const su2double *value,
41934193
unsigned short inMarkerTP,
41944194
unsigned short valSpan) { }
41954195

@@ -4216,7 +4216,7 @@ class CSolver {
42164216
* \param[in] value - turboperformance value to set.
42174217
* \param[in] inMarkerTP - turboperformance marker.
42184218
*/
4219-
inline virtual void SetTurboVelocityOut(su2double *value,
4219+
inline virtual void SetTurboVelocityOut(const su2double *value,
42204220
unsigned short inMarkerTP,
42214221
unsigned short valSpan) { }
42224222

‎SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp‎

Lines changed: 24 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -186,14 +186,10 @@ bool CDiscAdjMultizoneDriver::Iterate(unsigned short iZone, unsigned long iInner
186186
if(iInnerIter) SetAllSolutions(iZone, true, FixPtCorrector[iZone]);
187187
}
188188

189-
/*--- Residuals during GMRES iterations have no meaning ---*/
190-
191-
if (KrylovMode && iInnerIter) return false;
192-
193189
/*--- This is done explicitly here for multizone cases, only in inner iterations and not when
194190
* extracting cross terms so that the adjoint residuals in each zone still make sense. ---*/
195191

196-
Set_SolutionOld_To_Solution(iZone);
192+
if (!KrylovMode) Set_SolutionOld_To_Solution(iZone);
197193

198194
/*--- Print out the convergence data to screen and history file. ---*/
199195

@@ -231,7 +227,6 @@ void CDiscAdjMultizoneDriver::Run() {
231227
AdjRHS[iZone].Initialize(nPoint, nPointDomain, nVar, nullptr);
232228
AdjSol[iZone].Initialize(nPoint, nPointDomain, nVar, nullptr);
233229
LinSolver[iZone].SetRecomputeResidual(false);
234-
LinSolver[iZone].SetMonitoringFrequency(config_container[iZone]->GetScreen_Wrt_Freq(2));
235230
}
236231
else if (config_container[iZone]->GetnQuasiNewtonSamples() > 1) {
237232
FixPtCorrector[iZone].resize(config_container[iZone]->GetnQuasiNewtonSamples(), nPoint, nVar, nPointDomain);
@@ -371,20 +366,39 @@ void CDiscAdjMultizoneDriver::Run() {
371366
const bool monitor = config_container[iZone]->GetWrt_ZoneConv();
372367
const auto product = AdjointProduct(this, iZone);
373368

369+
/*--- Manipulate the screen output frequency to avoid printing garbage. ---*/
370+
const auto wrtFreq = config_container[iZone]->GetScreen_Wrt_Freq(2);
371+
config_container[iZone]->SetScreen_Wrt_Freq(2, nInnerIter[iZone]);
372+
LinSolver[iZone].SetMonitoringFrequency(wrtFreq);
373+
374374
Scalar eps = 1.0;
375375
for (auto totalIter = nInnerIter[iZone]; totalIter >= KrylovMinIters && eps > KrylovTol;) {
376376
Scalar eps_l = 0.0;
377377
Scalar tol_l = KrylovTol / eps;
378-
auto iter = min(totalIter-1ul, config_container[iZone]->GetnQuasiNewtonSamples()-1ul);
378+
auto iter = min(totalIter-2ul, config_container[iZone]->GetnQuasiNewtonSamples()-2ul);
379379
iter = LinSolver[iZone].FGMRES_LinSolver(AdjRHS[iZone], AdjSol[iZone], product, Identity(),
380380
tol_l, iter, eps_l, monitor, config_container[iZone]);
381381
totalIter -= iter+1;
382382
eps *= eps_l;
383383
}
384384

385+
/*--- Store the solution and restore user settings. ---*/
385386
SetAllSolutions(iZone, true, AdjSol[iZone]);
387+
config_container[iZone]->SetScreen_Wrt_Freq(2, wrtFreq);
386388

387-
Iterate(iZone, nInnerIter[iZone]-1);
389+
/*--- Set the old solution such that iterating gives meaningful residuals. ---*/
390+
AdjSol[iZone] += AdjRHS[iZone];
391+
SetAllSolutionsOld(iZone, true, AdjSol[iZone]);
392+
393+
/*--- Iterate to evaluate cross terms and residuals, this cannot happen within GMRES
394+
* because the vectors it multiplies by the Jacobian are not the actual solution. ---*/
395+
eval_transfer = true;
396+
Iterate(iZone, product.iInnerIter);
397+
398+
/*--- Set the solution as obtained from GMRES, otherwise it would be GMRES+Iterate once.
399+
* This is set without the "External" (by adding RHS above) so that it can be added
400+
* again in the next outer iteration with new contributions from other zones. ---*/
401+
SetAllSolutions(iZone, true, AdjSol[iZone]);
388402
}
389403

390404
/*--- Off-diagonal (coupling term) BGS update. ---*/
@@ -419,7 +433,8 @@ void CDiscAdjMultizoneDriver::Run() {
419433

420434
/*--- Check for convergence. ---*/
421435

422-
StopCalc = driver_output->GetConvergence() || (iOuterIter == nOuterIter-1);
436+
StopCalc = driver_output->GetConvergence() || (iOuterIter == nOuterIter-1) ||
437+
((nZone==1) && output_container[ZONE_0]->GetConvergence());
423438

424439
/*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/
425440

‎SU2_CFD/src/solvers/CDiscAdjSolver.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -429,8 +429,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
429429
/*--- Relax and store the adjoint solution, compute the residuals. ---*/
430430

431431
for (auto iVar = 0u; iVar < nVar; iVar++) {
432-
su2double residual = relax*(Solution[iVar]-nodes->GetSolution_Old(iPoint,iVar));
433-
nodes->AddSolution(iPoint, iVar, residual);
432+
su2double residual = Solution[iVar]-nodes->GetSolution_Old(iPoint,iVar);
433+
nodes->AddSolution(iPoint, iVar, relax*residual);
434434

435435
residual *= isdomain;
436436
AddRes_RMS(iVar,pow(residual,2));

0 commit comments

Comments
 (0)