Small discrete adjoint tweaks and other fixes

Common/include/CConfig.hpp

@@ -570,7 +570,7 @@ class CConfig {
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
-  su2double Relaxation_Factor_AdjFlow;  /*!< \brief Relaxation coefficient of the linear solver adjoint mean flow. */
+  su2double Relaxation_Factor_Adjoint;  /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
   su2double Relaxation_Factor_CHT;      /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
   su2double AdjTurb_Linear_Error;       /*!< \brief Min error of the turbulent adjoint linear solver for the implicit formulation. */
   su2double EntropyFix_Coeff;           /*!< \brief Entropy fix coefficient. */
@@ -3998,10 +3998,9 @@ class CConfig {
   su2double GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; }
 
   /*!
-   * \brief Get the relaxation coefficient of the linear solver for the implicit formulation.
-   * \return relaxation coefficient of the linear solver for the implicit formulation.
+   * \brief Get the relaxation factor for solution updates of adjoint solvers.
    */
-  su2double GetRelaxation_Factor_AdjFlow(void) const { return Relaxation_Factor_AdjFlow; }
+  su2double GetRelaxation_Factor_Adjoint(void) const { return Relaxation_Factor_Adjoint; }
 
   /*!
    * \brief Get the relaxation coefficient of the CHT coupling.

Common/src/CConfig.cpp

@@ -1670,8 +1670,8 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0);
   /* DESCRIPTION: Custom number of threads used for additive domain decomposition for ILU and LU_SGS (0 is "auto"). */
   addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0);
-  /* DESCRIPTION: Relaxation of the flow equations solver for the implicit formulation */
-  addDoubleOption("RELAXATION_FACTOR_ADJFLOW", Relaxation_Factor_AdjFlow, 1.0);
+  /* DESCRIPTION: Relaxation factor for updates of adjoint variables. */
+  addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
   /* DESCRIPTION: Relaxation of the CHT coupling */
   addDoubleOption("RELAXATION_FACTOR_CHT", Relaxation_Factor_CHT, 1.0);
   /* DESCRIPTION: Roe coefficient */
@@ -2832,20 +2832,9 @@ void CConfig::SetConfig_Parsing(char case_filename[MAX_STRING_SIZE]) {
           newString.append(": invalid option name");
           newString.append(". Check current SU2 options in config_template.cfg.");
           newString.append("\n");
-          if (!option_name.compare("EXT_ITER")) newString.append("Option EXT_ITER is deprecated as of v7.0. Please use TIME_ITER, OUTER_ITER or ITER \n"
-                                                                 "to specify the number of time iterations, outer multizone iterations or iterations, respectively.");
-          if (!option_name.compare("UNST_TIMESTEP")) newString.append("UNST_TIMESTEP is now TIME_STEP.\n");
-          if (!option_name.compare("UNST_TIME")) newString.append("UNST_TIME is now MAX_TIME.\n");
-          if (!option_name.compare("UNST_INT_ITER")) newString.append("UNST_INT_ITER is now INNER_ITER.\n");
-          if (!option_name.compare("RESIDUAL_MINVAL")) newString.append("RESIDUAL_MINVAL is now CONV_RESIDUAL_MINVAL.\n");
-          if (!option_name.compare("STARTCONV_ITER")) newString.append("STARTCONV_ITER is now CONV_STARTITER.\n");
-          if (!option_name.compare("CAUCHY_ELEMS")) newString.append("CAUCHY_ELEMS is now CONV_CAUCHY_ELEMS.\n");
-          if (!option_name.compare("CAUCHY_EPS")) newString.append("CAUCHY_EPS is now CONV_CAUCHY_EPS.\n");
-          if (!option_name.compare("OUTPUT_FORMAT")) newString.append("OUTPUT_FORMAT is now TABULAR_FORMAT.\n");
-          if (!option_name.compare("PHYSICAL_PROBLEM")) newString.append("PHYSICAL_PROBLEM is now SOLVER.\n");
-          if (!option_name.compare("REGIME_TYPE")) newString.append("REGIME_TYPE has been removed.\n "
-                                                                    "If you want use the incompressible solver, \n"
-                                                                    "use INC_EULER, INC_NAVIER_STOKES or INC_RANS as value of the SOLVER option.");
+          if (!option_name.compare("RELAXATION_FACTOR_ADJFLOW"))
+            newString.append("Option RELAXATION_FACTOR_ADJFLOW is now RELAXATION_FACTOR_ADJOINT, "
+                             "and it also applies to discrete adjoint problems\n.");
           errorString.append(newString);
           err_count++;
         continue;

SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp

@@ -180,6 +180,12 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
    */
   void Add_External_To_Solution(unsigned short iZone);
 
+  /*!
+   * \brief Puts Solution into SolutionOld.
+   * \param[in] iZone - Zone index.
+   */
+  void Set_SolutionOld_To_Solution(unsigned short iZone);
+
   /*!
    * \brief Extract contribution of iZone to jZone with BGS relaxation.
    * \param[in] iZone - Source zone (the one that was initialized).

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -55,8 +55,8 @@ class CFEASolver : public CSolver {
   su2double Conv_Check[3];          /*!< \brief Current values for convergence check: UTOL, RTOL, ETOL. */
   su2double FSI_Conv[2];            /*!< \brief Values to check the convergence of the FSI problem. */
 
-  unsigned long idxIncrement;       /*!< \brief Index of the current load increment */
-  su2double loadIncrement;          /*!< \brief Coefficient that determines the amount of load which is applied. */
+  unsigned long idxIncrement = 0;   /*!< \brief Index of the current load increment */
+  su2double loadIncrement = 1.0;    /*!< \brief Coefficient that determines the amount of load which is applied. */
 
   su2double WAitken_Dyn;            /*!< \brief Aitken's dynamic coefficient. */
   su2double WAitken_Dyn_tn1;        /*!< \brief Aitken's dynamic coefficient in the previous iteration. */

SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp

@@ -256,7 +256,7 @@ void CDiscAdjMultizoneDriver::Run() {
           Add_External_To_Solution(iZone);
         }
         else {
-          /*--- If we restarted, Solution already has all contribution,
+          /*--- If we restarted, Solution already has all contributions,
            *    we run only one inner iter to compute the cross terms. ---*/
           eval_transfer = true;
         }
@@ -274,7 +274,12 @@ void CDiscAdjMultizoneDriver::Run() {
                                                     solver_container, numerics_container, config_container,
                                                     surface_movement, grid_movement, FFDBox, iZone, INST_0);
 
-        /*--- Print out the convergence data to screen and history file ---*/
+        /*--- This is done explicitly here for multizone cases, only in inner iterations and not when
+         *    extracting cross terms so that the adjoint residuals in each zone still make sense. ---*/
+
+        Set_SolutionOld_To_Solution(iZone);
+
+        /*--- Print out the convergence data to screen and history file. ---*/
 
         bool converged = iteration_container[iZone][INST_0]->Monitor(output_container[iZone], integration_container,
                                                     geometry_container, solver_container, numerics_container,
@@ -293,6 +298,7 @@ void CDiscAdjMultizoneDriver::Run() {
           /*--- Extracting adjoints for solvers in jZone w.r.t. to the output of all solvers in iZone,
            *    that is, for the cases iZone != jZone we are evaluating cross derivatives between zones. ---*/
 
+          config_container[jZone]->SetInnerIter(0);
           iteration_container[jZone][INST_0]->Iterate(output_container[jZone], integration_container, geometry_container,
                                                       solver_container, numerics_container, config_container,
                                                       surface_movement, grid_movement, FFDBox, jZone, INST_0);
@@ -898,6 +904,15 @@ void CDiscAdjMultizoneDriver::Add_External_To_Solution(unsigned short iZone) {
   }
 }
 
+void CDiscAdjMultizoneDriver::Set_SolutionOld_To_Solution(unsigned short iZone) {
+
+  for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {
+    auto solver = solver_container[iZone][INST_0][MESH_0][iSol];
+    if (solver && solver->GetAdjoint())
+      solver->GetNodes()->Set_OldSolution();
+  }
+}
+
 void CDiscAdjMultizoneDriver::Update_Cross_Term(unsigned short iZone, unsigned short jZone) {
 
   for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {

SU2_CFD/src/iteration_structure.cpp

@@ -473,11 +473,12 @@ void CFluidIteration::Iterate(COutput *output,
                                                                      RUNTIME_RADIATION_SYS, val_iZone, val_iInst);
   }
 
-  /*--- Adapt the CFL number using an exponential progression
-   with under-relaxation approach. ---*/
+  /*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/
 
   if (config[val_iZone]->GetCFL_Adapt() == YES) {
-    solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone]);
+    SU2_OMP_PARALLEL
+    solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst],
+                                                                   solver[val_iZone][val_iInst], config[val_iZone]);
   }
 
   /*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/

SU2_CFD/src/solvers/CAdjEulerSolver.cpp

@@ -2301,7 +2301,7 @@ void CAdjEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **sol
 
   for (iPoint = 0; iPoint < nPointDomain; iPoint++)
     for (iVar = 0; iVar < nVar; iVar++) {
-      nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_AdjFlow()*LinSysSol[iPoint*nVar+iVar]);
+      nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_Adjoint()*LinSysSol[iPoint*nVar+iVar]);
     }
 
   /*--- MPI solution ---*/

SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp

@@ -667,8 +667,6 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co
     }
   }
 
-  if(multizone) nodes->Set_OldSolution();
-
   SetResidual_RMS(geometry, config);
 
 }

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -476,26 +476,26 @@ void CDiscAdjSolver::SetAdj_ObjFunc(CGeometry *geometry, CConfig *config) {
 
 void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config){
 
-  bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
-  bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
-  bool multizone = config->GetMultizone_Problem();
+  const bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
+  const bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
+  const bool multizone = config->GetMultizone_Problem();
 
-  unsigned short iVar;
-  unsigned long iPoint;
-  su2double residual;
+  const su2double relax = (config->GetInnerIter()==0)? 1.0 : config->GetRelaxation_Factor_Adjoint();
 
   /*--- Set Residuals to zero ---*/
 
-  for (iVar = 0; iVar < nVar; iVar++) {
-      SetRes_RMS(iVar,0.0);
-      SetRes_Max(iVar,0.0,0);
+  for (auto iVar = 0u; iVar < nVar; iVar++) {
+    SetRes_RMS(iVar,0.0);
+    SetRes_Max(iVar,0.0,0);
   }
 
-  /*--- Set the old solution ---*/
+  /*--- Set the old solution and compute residuals. ---*/
 
   if(!multizone) nodes->Set_OldSolution();
 
-  for (iPoint = 0; iPoint < nPoint; iPoint++) {
+  for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
+
+    const su2double isdomain = (iPoint < nPointDomain)? 1.0 : 0.0;
 
     /*--- Extract the adjoint solution ---*/
 
@@ -506,13 +506,22 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
       direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);
     }
 
-    /*--- Store the adjoint solution ---*/
+    /*--- Relax and store the adjoint solution, compute the residuals. ---*/
 
-    nodes->SetSolution(iPoint,Solution);
+    for (auto iVar = 0u; iVar < nVar; iVar++) {
+      su2double residual = relax*(Solution[iVar]-nodes->GetSolution_Old(iPoint,iVar));
+      nodes->AddSolution(iPoint, iVar, residual);
+
+      residual *= isdomain;
+      AddRes_RMS(iVar,pow(residual,2));
+      AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
+    }
   }
 
+  SetResidual_RMS(geometry, config);
+
   if (time_n_needed) {
-    for (iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
 
       /*--- Extract the adjoint solution at time n ---*/
 
@@ -523,8 +532,9 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
       nodes->Set_Solution_time_n(iPoint,Solution);
     }
   }
+
   if (time_n1_needed) {
-    for (iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {
 
       /*--- Extract the adjoint solution at time n-1 ---*/
 
@@ -536,21 +546,6 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
     }
   }
 
-  /*--- Set the residuals ---*/
-
-  for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
-    for (iVar = 0; iVar < nVar; iVar++) {
-      residual = nodes->GetSolution(iPoint,iVar) - nodes->GetSolution_Old(iPoint,iVar);
-
-      AddRes_RMS(iVar,residual*residual);
-      AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
-    }
-  }
-
-  if(multizone) nodes->Set_OldSolution();
-
-  SetResidual_RMS(geometry, config);
-
 }
 
 void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) {

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -48,8 +48,6 @@ CFEASolver::CFEASolver(bool mesh_deform_mode) : CSolver(mesh_deform_mode) {
   Total_CFEA = 0.0;
   WAitken_Dyn = 0.0;
   WAitken_Dyn_tn1 = 0.0;
-  idxIncrement = 0;
-  loadIncrement = 1.0;
 
   element_container = new CElement** [MAX_TERMS]();
   for (unsigned short iTerm = 0; iTerm < MAX_TERMS; iTerm++)
@@ -124,8 +122,6 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() {
   Total_CFEA        = 0.0;
   WAitken_Dyn       = 0.0;
   WAitken_Dyn_tn1   = 0.0;
-  idxIncrement      = 0;
-  loadIncrement     = 0.0;
 
   SetFSI_ConvValue(0,0.0);
   SetFSI_ConvValue(1,0.0);

SU2_CFD/src/solvers/CSolver.cpp

@@ -2320,8 +2320,6 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
                              CSolver   ***solver_container,
                              CConfig   *config) {
 
-  /// TODO: Add OpenMP stuff to this method.
-
   /* Adapt the CFL number on all multigrid levels using an
    exponential progression with under-relaxation approach. */
 
@@ -2330,6 +2328,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
   const su2double CFLFactorIncrease = config->GetCFL_AdaptParam(1);
   const su2double CFLMin            = config->GetCFL_AdaptParam(2);
   const su2double CFLMax            = config->GetCFL_AdaptParam(3);
+  const bool fullComms              = (config->GetComm_Level() == COMM_FULL);
 
   for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
 
@@ -2365,6 +2364,9 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
     /* Check that we are meeting our nonlinear residual reduction target
      over time so that we do not get stuck in limit cycles. */
 
+    SU2_OMP_MASTER
+    { /* Only the master thread updates the shared variables. */
+
     Old_Func = New_Func;
     unsigned short Res_Count = 100;
     if (NonLinRes_Series.size() == 0) NonLinRes_Series.resize(Res_Count,0.0);
@@ -2408,17 +2410,30 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
      Reset the array so that we delay the next decrease for some iterations. */
 
     if (fabs(NonLinRes_Value) < 0.1*New_Func) {
-      reduceCFL = true;
       NonLinRes_Counter = 0;
       for (unsigned short iCounter = 0; iCounter < Res_Count; iCounter++)
         NonLinRes_Series[iCounter] = New_Func;
     }
 
+    } /* End SU2_OMP_MASTER, now all threads update the CFL number. */
+    SU2_OMP_BARRIER
+
+    if (fabs(NonLinRes_Value) < 0.1*New_Func) {
+      reduceCFL = true;
+    }
+
     /* Loop over all points on this grid and apply CFL adaption. */
 
-    su2double myCFLMin = 1e30;
-    su2double myCFLMax = 0.0;
-    su2double myCFLSum = 0.0;
+    su2double myCFLMin = 1e30, myCFLMax = 0.0, myCFLSum = 0.0;
+
+    SU2_OMP_MASTER
+    if ((iMesh == MESH_0) && fullComms) {
+      Min_CFL_Local = 1e30;
+      Max_CFL_Local = 0.0;
+      Avg_CFL_Local = 0.0;
+    }
+
+    SU2_OMP_FOR_STAT(roundUpDiv(geometry[iMesh]->GetnPointDomain(),omp_get_max_threads()))
     for (unsigned long iPoint = 0; iPoint < geometry[iMesh]->GetnPointDomain(); iPoint++) {
 
       /* Get the current local flow CFL number at this point. */
@@ -2476,20 +2491,37 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
         solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL);
       }
 
-      /* Store min and max CFL for reporting on fine grid. */
+      /* Store min and max CFL for reporting on the fine grid. */
 
-      myCFLMin = min(CFL,myCFLMin);
-      myCFLMax = max(CFL,myCFLMax);
-      myCFLSum += CFL;
+      if ((iMesh == MESH_0) && fullComms) {
+        myCFLMin = min(CFL,myCFLMin);
+        myCFLMax = max(CFL,myCFLMax);
+        myCFLSum += CFL;
+      }
 
     }
 
     /* Reduce the min/max/avg local CFL numbers. */
 
-    SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
-    SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
-    SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-    Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
+    if ((iMesh == MESH_0) && fullComms) {
+      SU2_OMP_CRITICAL
+      { /* OpenMP reduction. */
+        Min_CFL_Local = min(Min_CFL_Local,myCFLMin);
+        Max_CFL_Local = max(Max_CFL_Local,myCFLMax);
+        Avg_CFL_Local += myCFLSum;
+      }
+      SU2_OMP_BARRIER
+
+      SU2_OMP_MASTER
+      { /* MPI reduction. */
+        myCFLMin = Min_CFL_Local; myCFLMax = Max_CFL_Local; myCFLSum = Avg_CFL_Local;
+        SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+        SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+        SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
+      }
+      SU2_OMP_BARRIER
+    }
 
   }
 

TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg

@@ -174,7 +174,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
 %
 % Relaxation coefficient
-RELAXATION_FACTOR_ADJFLOW= 1.0
+RELAXATION_FACTOR_ADJOINT= 1.0
 %
 % Reduction factor of the CFL coefficient in the adjoint problem
 CFL_REDUCTION_ADJFLOW= 0.8