From c1ecb6055f324ab5a4538c9dfe46638f3152c160 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 10 May 2026 19:10:02 -0700 Subject: [PATCH 1/3] allow adjoin-only, fix residual and objective reporting in single and multizone solvers --- Common/src/linear_algebra/CPastixWrapper.cpp | 12 ++++------ .../src/drivers/CDiscAdjMultizoneDriver.cpp | 11 ++++++--- .../src/drivers/CDiscAdjSinglezoneDriver.cpp | 7 +++++- SU2_CFD/src/drivers/CDriver.cpp | 2 -- SU2_CFD/src/output/CElasticityOutput.cpp | 6 ++--- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 23 +++++++++++++++---- SU2_CFD/src/solvers/CFEASolver.cpp | 15 ++++++++---- 7 files changed, 50 insertions(+), 26 deletions(-) diff --git a/Common/src/linear_algebra/CPastixWrapper.cpp b/Common/src/linear_algebra/CPastixWrapper.cpp index 1633a88bbde..bcdf38ed9cc 100644 --- a/Common/src/linear_algebra/CPastixWrapper.cpp +++ b/Common/src/linear_algebra/CPastixWrapper.cpp @@ -204,8 +204,7 @@ void CPastixWrapper::Initialize(CGeometry* geometry, const CConfig* SU2_MPI::Error("Error analyzing matrix: " + std::to_string(rc), CURRENT_FUNCTION); } - if (mpi_rank == MASTER_NODE && verb > 0) - cout << " +--------------------------------------------------------------------+" << endl; + if (mpi_rank == MASTER_NODE && verb > 0) cout << "+-------------------------------------------------+" << endl; isinitialized = true; } @@ -250,10 +249,8 @@ void CPastixWrapper::Factorize(CGeometry* geometry, const CConfig* c /*--- Yes ---*/ if (mpi_rank == MASTER_NODE && verb > 0) { - cout << endl; - cout << " +--------------------------------------------------------------------+" << endl; - cout << " + PaStiX : Parallel Sparse matriX package +" << endl; - cout << " +--------------------------------------------------------------------+" << endl; + cout << "\n+-------------------------------------------------+"; + cout << "\n+ PaStiX : Parallel Sparse matriX package +" << endl; } const unsigned long szBlk = matrix.nVar * matrix.nVar, nNonZero = values.size(); @@ -297,8 +294,7 @@ void CPastixWrapper::Factorize(CGeometry* geometry, const CConfig* c SU2_MPI::Error("Error factorizing matrix: " + std::to_string(rc), CURRENT_FUNCTION); } - if (mpi_rank == MASTER_NODE && verb > 0) - cout << " +--------------------------------------------------------------------+" << endl << endl; + if (mpi_rank == MASTER_NODE && verb > 0) cout << "+-------------------------------------------------+\n" << endl; isfactorized = true; } diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 61673bdad1b..c9529050da0 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -719,7 +719,9 @@ void CDiscAdjMultizoneDriver::SetRecording(RECORDING kind_recording, Kind_Tape t switch(kind_recording) { case RECORDING::CLEAR_INDICES: cout << "Clearing the computational graph." << endl; break; case RECORDING::MESH_COORDS: cout << "Storing computational graph wrt MESH COORDINATES." << endl; break; - case RECORDING::SOLUTION_VARIABLES: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << endl; break; + case RECORDING::SOLUTION_VARIABLES: + cout << "Storing computational graph wrt CONSERVATIVE VARIABLES.\n"; + cout << "Computing residuals to check the convergence of the direct problem." << endl; break; case RECORDING::TAG_INIT_SOLVER_VARIABLES: cout << "Simulating recording with tag 1 on conservative variables." << endl; AD::SetTag(1); break; case RECORDING::TAG_CHECK_SOLVER_VARIABLES: cout << "Checking first recording with tag 2 on conservative variables." << endl; AD::SetTag(2); break; case RECORDING::TAG_INIT_SOLVER_AND_MESH: cout << "Simulating recording with tag 1 on conservative variables and mesh coordinates." << endl; AD::SetTag(1); break; @@ -850,7 +852,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { solvers[FLOW_SOL]->Momentum_Forces(geometry, config); solvers[FLOW_SOL]->Friction_Forces(geometry, config); - if(config->GetWeakly_Coupled_Heat()) { + if (config->GetWeakly_Coupled_Heat()) { solvers[HEAT_SOL]->Heat_Fluxes(geometry, solvers, config); } @@ -865,6 +867,9 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { break; case MAIN_SOLVER::DISC_ADJ_FEM: + if (config->GetWeakly_Coupled_Heat()) { + solvers[HEAT_SOL]->Heat_Fluxes(geometry, solvers, config); + } solvers[FEA_SOL]->Postprocessing(geometry, config, numerics_container[iZone][INST_0][MESH_0][FEA_SOL], true); direct_output[iZone]->SetHistoryOutput(geometry, solvers, config); ObjFunc += solvers[FEA_SOL]->GetTotal_ComboObj(); @@ -883,7 +888,7 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { kind_recording == RECORDING::TAG_CHECK_SOLVER_VARIABLES || kind_recording == RECORDING::TAG_INIT_SOLVER_AND_MESH || kind_recording == RECORDING::TAG_CHECK_SOLVER_AND_MESH) { - cout << " Objective function : " << ObjFunc << endl; + cout << "Objective function value: " << std::setprecision(driver_config->GetOutput_Precision()) << ObjFunc << endl; } } } diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index ce6f7d93b10..b2ef750d7ea 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -273,7 +273,7 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ case RECORDING::CLEAR_INDICES: cout << "Clearing the computational graph." << endl; break; case RECORDING::MESH_COORDS: cout << "Storing computational graph wrt MESH COORDINATES." << endl; break; case RECORDING::SOLUTION_VARIABLES: - cout << "Direct iteration to store the primal computational graph." << endl; + cout << "Direct iteration to store the primal computational graph.\n"; cout << "Computing residuals to check the convergence of the direct problem." << endl; break; default: break; } @@ -309,6 +309,11 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ SetObjFunction(); + if (rank == MASTER_NODE && kind_recording == RECORDING::SOLUTION_VARIABLES) { + cout << "\nObjective function value: " << std::setprecision(config_container[ZONE_0]->GetOutput_Precision()) << ObjFunc; + cout << "\n-------------------------------------------------------------------------\n" << endl; + } + if (kind_recording != RECORDING::CLEAR_INDICES && config_container[ZONE_0]->GetWrt_AD_Statistics()) { AD::PrintStatistics(SU2_MPI::GetComm(), rank == MASTER_NODE); } diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index daf5037a538..8689f42ab9f 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2886,8 +2886,6 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) { } // for addVals - cout << "\n-------------------------------------------------------------------------\n" << endl; - } diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 2694b034674..1a78dbd25e6 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -107,10 +107,10 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS /*--- Nonlinear analysis: UTOL, RTOL and DTOL (defined in the Postprocessing function) ---*/ if (linear_analysis){ - SetHistoryOutputValue("RMS_DISP_X", log10(fea_solver->GetRes_RMS(0))); - SetHistoryOutputValue("RMS_DISP_Y", log10(fea_solver->GetRes_RMS(1))); + SetHistoryOutputValue("RMS_DISP_X", log10(fea_solver->GetRes_FEM(0))); + SetHistoryOutputValue("RMS_DISP_Y", log10(fea_solver->GetRes_FEM(1))); if (nDim == 3){ - SetHistoryOutputValue("RMS_DISP_Z", log10(fea_solver->GetRes_RMS(2))); + SetHistoryOutputValue("RMS_DISP_Z", log10(fea_solver->GetRes_FEM(2))); } } else if (nonlinear_analysis){ SetHistoryOutputValue("RMS_UTOL", log10(fea_solver->GetRes_FEM(0))); diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 8e1a6499db5..f8e2844d3fd 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -126,16 +126,31 @@ CDiscAdjFEASolver::~CDiscAdjFEASolver() { delete nodes; } void CDiscAdjFEASolver::SetRecording(CGeometry* geometry, CConfig *config){ SU2_ZONE_SCOPED - /*--- Reset the solution to the initial (converged) solution ---*/ + /*--- Under some conditions, linear elasticity problems converge in one iteration. + * This means that the "clear indices" step of the discrete adjoint solver produces a + * converged solution regardless of the primal solution given to the adjoint solver. + * We can take advantage of this for optimization to skip the primal solver. ---*/ + + const bool linear = config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL; + const bool heat = config->GetWeakly_Coupled_Heat(); + const bool time_domain = config->GetTime_Domain(); + const bool keep_solution = linear && !heat && !time_domain && !std::is_same_v; + + /*--- Restore the solution to the initial (converged) solution or reset AD indices. ---*/ for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - for (auto iVar = 0u; iVar < nVar; iVar++) - direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]); + if (keep_solution) { + for (auto iVar = 0u; iVar < nVar; iVar++) + AD::ResetInput(direct_solver->GetNodes()->GetSolution(iPoint)[iVar]); + } else { + for (auto iVar = 0u; iVar < nVar; iVar++) + direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]); + } } /*--- Reset the input for time n ---*/ - if (config->GetTime_Domain()) { + if (time_domain) { for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) for (auto iVar = 0u; iVar < nVar; iVar++) AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n(iPoint)[iVar]); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index e12f92140a4..31cb41c52d9 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1833,14 +1833,13 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CConfig *config, CNumerics /*--- RTOL = norm(Residual(k): ABSOLUTE, norm of the residual (T-F) ---*/ /*--- ETOL = Delta_U(k) * Residual(k): ABSOLUTE, energy norm ---*/ - SU2_OMP_PARALLEL - { + SU2_OMP_PARALLEL { + su2double utol = LinSysSol.norm(); su2double rtol = LinSysRes.norm(); su2double etol = fabs(LinSysSol.dot(LinSysRes)); - SU2_OMP_MASTER - { + SU2_OMP_MASTER { Conv_Check[0] = utol; Conv_Check[1] = rtol; Conv_Check[2] = etol; @@ -1872,8 +1871,14 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CConfig *config, CNumerics END_SU2_OMP_FOR /*--- "Add" residuals from all threads to global residual variables. ---*/ - ResidualReductions_FromAllThreads(geometry, config, resRMS,resMax,idxMax); + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); + SU2_OMP_MASTER { + Conv_Check[0] = Residual_RMS[0]; + Conv_Check[1] = Residual_RMS[1]; + Conv_Check[2] = nDim == 3 ? Residual_RMS[2] : 0; + } + END_SU2_OMP_MASTER } END_SU2_OMP_PARALLEL From a2be4c9766e3d5875aff28682be38acd44ef8aa1 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 16 May 2026 11:24:28 -0700 Subject: [PATCH 2/3] save main recording for functions that don't need it --- .../drivers/CDiscAdjSinglezoneDriver.hpp | 8 ++++++++ .../src/drivers/CDiscAdjSinglezoneDriver.cpp | 19 ++++++++++++++++++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp index 137c961f148..6a2f00932ef 100644 --- a/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp @@ -55,6 +55,14 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver { COutput *direct_output; CNumerics ***numerics; /*!< \brief Container vector with all the numerics. */ + /*! + * \brief Returns true if the objective function does not depend on the main variables. In which case, + * the adjoint variables are 0 and the sensitivities can be computed just with the secondary recording. + */ + bool TrivialFunction() const { + return config_container[ZONE_0]->GetnObj() == 1 && config_container[ZONE_0]->GetKind_ObjFunc() == VOLUME_FRACTION; + } + /*! * \brief Record one iteration of a flow iteration in within multiple zones. * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp) diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index b2ef750d7ea..5fcd203b84d 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -156,6 +156,12 @@ void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) { void CDiscAdjSinglezoneDriver::Run() { SU2_ZONE_SCOPED + /*--- No need to solve anything, the tape for the main recording is empty. ---*/ + if (TrivialFunction()) { + SetAllSolutions(ZONE_0, true, [](auto, auto) { return 0; }); + return; + } + CQuasiNewtonInvLeastSquares fixPtCorrector; if (config->GetnQuasiNewtonSamples() > 1) { fixPtCorrector.resize(config->GetnQuasiNewtonSamples(), @@ -309,7 +315,8 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ SetObjFunction(); - if (rank == MASTER_NODE && kind_recording == RECORDING::SOLUTION_VARIABLES) { + if (rank == MASTER_NODE && + (kind_recording == RECORDING::SOLUTION_VARIABLES || (TrivialFunction() && kind_recording == RECORDING::MESH_COORDS))) { cout << "\nObjective function value: " << std::setprecision(config_container[ZONE_0]->GetOutput_Precision()) << ObjFunc; cout << "\n-------------------------------------------------------------------------\n" << endl; } @@ -415,6 +422,16 @@ void CDiscAdjSinglezoneDriver::DirectRun(RECORDING kind_recording){ void CDiscAdjSinglezoneDriver::MainRecording(){ SU2_ZONE_SCOPED + + /*--- We know this function only depends on the secondary variables, hence skip the main recording. ---*/ + + if (TrivialFunction()) { + if (rank == MASTER_NODE) { + cout << "Trivial objective function, skipping the solution of the adjoint equations." << endl; + } + return; + } + /*--- SetRecording stores the computational graph on one iteration of the direct problem. Calling it with * RECORDING::CLEAR_INDICES as argument ensures that all information from a previous recording is removed. ---*/ From 8cb37589eb62d9182ff4d6e44f61e1d40308740e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 23 May 2026 15:32:54 -0700 Subject: [PATCH 3/3] update ref --- TestCases/fea_topology/grad_ref_node.dat.ref | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/fea_topology/grad_ref_node.dat.ref b/TestCases/fea_topology/grad_ref_node.dat.ref index 6f3a4d4a685..15afff813c1 100644 --- a/TestCases/fea_topology/grad_ref_node.dat.ref +++ b/TestCases/fea_topology/grad_ref_node.dat.ref @@ -6157,7 +6157,7 @@ -9.57383e-11 -2.60763e-10 -2.54078e-10 --1.57598e-10 +-1.57599e-10 -8.68931e-12 -4.26927e-13 0