Efficiency improvements for optimization of linear elasticity problems (#2815)

* allow adjoin-only, fix residual and objective reporting in single and multizone solvers

* save main recording for functions that don't need it

* update ref

Author: pcarruscag

Common/src/linear_algebra/CPastixWrapper.cpp

@@ -204,8 +204,7 @@ void CPastixWrapper<ScalarType>::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<ScalarType>::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<ScalarType>::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;
 }

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)

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

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<passivedouble> fixPtCorrector;
   if (config->GetnQuasiNewtonSamples() > 1) {
     fixPtCorrector.resize(config->GetnQuasiNewtonSamples(),
@@ -273,7 +279,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 +315,12 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){
 
   SetObjFunction();
 
+  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;
+  }
+
   if (kind_recording != RECORDING::CLEAR_INDICES && config_container[ZONE_0]->GetWrt_AD_Statistics()) {
     AD::PrintStatistics(SU2_MPI::GetComm(), rank == MASTER_NODE);
   }
@@ -410,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. ---*/
 

SU2_CFD/src/drivers/CDriver.cpp

@@ -2886,8 +2886,6 @@ void CDriver::PrintDirectResidual(RECORDING kind_recording) {
 
   } // for addVals
 
-  cout << "\n-------------------------------------------------------------------------\n" << endl;
-
 }
 
 

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

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<su2mixedfloat, float>;
+
+  /*--- 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]);

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
 

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