Merge pull request #286 from su2code/feature_unst_discadj

Unsteady Discrete Adjoint

Common/include/config_structure.hpp

@@ -228,7 +228,7 @@ class CConfig {
 	su2double *NRBC_Var1, *NRBC_Var2;    /*!< \brief Specified values for NRBC boundary. */
 	su2double **NRBC_FlowDir;  /*!< \brief Specified flow direction vector (unit vector) for NRBC boundaries. */
 	su2double *Inlet_Ptotal;    /*!< \brief Specified total pressures for inlet boundaries. */
-	su2double **Inlet_FlowDir;  /*!< \brief Specified flow direction vector (unit vector) for inlet boundaries. */
+    su2double **Inlet_FlowDir;  /*!< \brief Specified flow direction vector (unit vector) for inlet boundaries. */
 	su2double *Inlet_Temperature;    /*!< \brief Specified temperatures for a supersonic inlet boundaries. */
 	su2double *Inlet_Pressure;    /*!< \brief Specified static pressures for supersonic inlet boundaries. */
 	su2double **Inlet_Velocity;  /*!< \brief Specified flow velocity vectors for supersonic inlet boundaries. */
@@ -281,6 +281,7 @@ class CConfig {
 	unsigned long Dyn_nIntIter;			/*!< \brief Number of internal iterations (Newton-Raphson Method for nonlinear structural analysis). */
   long Unst_RestartIter;			/*!< \brief Iteration number to restart an unsteady simulation (Dual time Method). */
   long Unst_AdjointIter;			/*!< \brief Iteration number to begin the reverse time integration in the direct solver for the unsteady adjoint. */
+  long Iter_Avg_Objective;			/*!< \brief Iteration the number of time steps to be averaged, counting from the back */
   long Dyn_RestartIter;			/*!< \brief Iteration number to restart a dynamic structural analysis. */
   unsigned short nRKStep;			/*!< \brief Number of steps of the explicit Runge-Kutta method. */
 	su2double *RK_Alpha_Step;			/*!< \brief Runge-Kutta beta coefficients. */
@@ -2195,7 +2196,13 @@ class CConfig {
   long GetUnst_AdjointIter(void);
 
   /*!
-	 * \brief Get the restart iteration number for dynamic structural simulations.
+  * \brief Number of iterations to average (reverse time integration).
+  * \return Starting direct iteration number for the unsteady adjoint.
+  */
+  unsigned long GetIter_Avg_Objective(void);
+
+  /*!
+   * \brief Get the restart iteration number for dynamic structural simulations.
 	 * \return Restart iteration number for dynamic structural simulations.
 	 */
   long GetDyn_RestartIter(void);
@@ -2407,7 +2414,7 @@ class CConfig {
 	 */
 	string GetMarker_EngineExhaust(unsigned short val_marker);
 
-    /*!
+  /*!
 	 * \brief Get the name of the surface defined in the geometry file.
 	 * \param[in] val_marker - Value of the marker in which we are interested.
 	 * \return Name that is in the geometry file for the surface that

Common/include/config_structure.inl

@@ -89,6 +89,8 @@ inline long CConfig::GetUnst_RestartIter(void) { return Unst_RestartIter; }
 
 inline long CConfig::GetUnst_AdjointIter(void) { return Unst_AdjointIter; }
 
+inline unsigned long CConfig::GetIter_Avg_Objective(void) { return Iter_Avg_Objective ; }
+
 inline long CConfig::GetDyn_RestartIter(void) { return Dyn_RestartIter; }
 
 inline string CConfig::GetPlaneTag(unsigned short index) { return PlaneTag[index]; }

Common/src/ad_structure.cpp

@@ -112,6 +112,7 @@ namespace AD {
 
       /* --- Clear local vectors and reset indicator ---*/
 
+
       localInputValues.clear();
       localOutputValues.clear();
 

Common/src/config_structure.cpp

@@ -561,6 +561,8 @@ void CConfig::SetConfig_Options(unsigned short val_iZone, unsigned short val_nZo
   addLongOption("UNST_RESTART_ITER", Unst_RestartIter, 0);
   /* DESCRIPTION: Starting direct solver iteration for the unsteady adjoint */
   addLongOption("UNST_ADJOINT_ITER", Unst_AdjointIter, 0);
+  /* DESCRIPTION: Number of iterations to average the objective */
+  addLongOption("ITER_AVERAGE_OBJ", Iter_Avg_Objective , 0);
   /* DESCRIPTION: Iteration number to begin unsteady restarts (structural analysis) */
   addLongOption("DYN_RESTART_ITER", Dyn_RestartIter, 0);
   /* DESCRIPTION: Time discretization */
@@ -2568,6 +2570,22 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
     /*--- Disable writing of limiters if enabled ---*/
     Wrt_Limiters = false;
 
+    if (Unsteady_Simulation){
+
+      Restart_Flow = false;
+
+      if (Grid_Movement){
+        cout << "Dynamic mesh movement currently not supported for the discrete adjoint solver." << endl;
+        exit(EXIT_FAILURE);
+      }
+
+      /* --- If the averaging interval is not set, we average over all time-steps ---*/
+
+      if (Iter_Avg_Objective == 0.0){
+        Iter_Avg_Objective = nExtIter;
+      }
+    }
+
     switch(Kind_Solver){
       case EULER:
         Kind_Solver = DISC_ADJ_EULER;

Common/src/geometry_structure.cpp

@@ -11811,23 +11811,17 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config){
   bool sa = config->GetKind_Turb_Model() == SA;
   bool grid_movement = config->GetGrid_Movement();
   su2double Sens, dull_val;
-  //su2double delta_T, total_T;
-  unsigned short nExtIter, iDim, iExtIter;
+  unsigned short nExtIter, iDim;
   unsigned long iPoint, index;
-  
+
   Sensitivity = new su2double[nPoint*nDim];
-  
+
   if (config->GetUnsteady_Simulation()){
-    nExtIter = config->GetUnst_AdjointIter();
-    //    delta_T  = config->GetDelta_UnstTimeND();
-    //    delta_T  = 1.0;
-    //total_T  = (su2double)nExtIter*delta_T;
-  } else {
-    //total_T = 1.0;
+    nExtIter = config->GetnExtIter();
+  }else{
     nExtIter = 1;
   }
-
-  int rank = MASTER_NODE;
+    int rank = MASTER_NODE;
 #ifdef HAVE_MPI
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 #endif
@@ -11868,55 +11862,48 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config){
     }
   }
 
-  for (iExtIter = 0; iExtIter < nExtIter; iExtIter++){
-
-    iPoint_Global = 0;
-
-    filename = config->GetSolution_AdjFileName();
-
-    filename = config->GetObjFunc_Extension(filename);
-
-    if (config->GetUnsteady_Simulation()){
-      filename = config->GetUnsteady_FileName(filename, iExtIter);
-    }
-
-    restart_file.open(filename.data(), ios::in);
-    if (restart_file.fail()) {
-      cout << "There is no adjoint restart file!! " << filename.data() << "."<< endl;
-      exit(EXIT_FAILURE);
-    }
-
-    if (rank == MASTER_NODE)
-      cout << "Reading in sensitivity at iteration " << iExtIter << "."<< endl;
-    /*--- The first line is the header ---*/
-    getline (restart_file, text_line);
-
-    while (getline (restart_file, text_line)) {
-      istringstream point_line(text_line);
-
-      /*--- Retrieve local index. If this node from the restart file lives
-       on a different processor, the value of iPoint_Local will be -1.
-       Otherwise, the local index for this node on the current processor
-       will be returned and used to instantiate the vars. ---*/
-      iPoint_Local = Global2Local[iPoint_Global];
-
-      if (iPoint_Local >= 0){
-        point_line >> index;
-        for (iDim = 0; iDim < skipVar; iDim++){ point_line >> dull_val;}
-        for (iDim = 0; iDim < nDim; iDim++){
-          point_line >> Sens;
-          //                	  Sensitivity[iPoint_Local*nDim+iDim] += Sens*delta_T/total_T;
-          Sensitivity[iPoint_Local*nDim+iDim] += Sens;
-
-        }
+
+  iPoint_Global = 0;
+
+  filename = config->GetSolution_AdjFileName();
+
+  filename = config->GetObjFunc_Extension(filename);
+
+  if (config->GetUnsteady_Simulation()){
+    filename = config->GetUnsteady_FileName(filename, nExtIter-1);
+  }
+
+  restart_file.open(filename.data(), ios::in);
+  if (restart_file.fail()) {
+    cout << "There is no adjoint restart file!! " << filename.data() << "."<< endl;
+    exit(EXIT_FAILURE);
+  }
+
+  if (rank == MASTER_NODE)
+    cout << "Reading in sensitivity at iteration " << nExtIter-1 << "."<< endl;
+  /*--- The first line is the header ---*/
+  getline (restart_file, text_line);
+
+  while (getline (restart_file, text_line)) {
+    istringstream point_line(text_line);
+
+    /*--- Retrieve local index. If this node from the restart file lives
+             on a different processor, the value of iPoint_Local will be -1.
+             Otherwise, the local index for this node on the current processor
+             will be returned and used to instantiate the vars. ---*/
+    iPoint_Local = Global2Local[iPoint_Global];
+
+    if (iPoint_Local >= 0){
+      point_line >> index;
+      for (iDim = 0; iDim < skipVar; iDim++){ point_line >> dull_val;}
+      for (iDim = 0; iDim < nDim; iDim++){
+        point_line >> Sens;
+        Sensitivity[iPoint_Local*nDim+iDim] = Sens;
       }
-      iPoint_Global++;
     }
-    restart_file.close();
+    iPoint_Global++;
   }
-
-  delete [] Global2Local;
-
+  restart_file.close();
 }
 
 su2double CPhysicalGeometry::Compute_MaxThickness(su2double *Plane_P0, su2double *Plane_Normal, unsigned short iSection, CConfig *config, vector<su2double> &Xcoord_Airfoil, vector<su2double> &Ycoord_Airfoil, vector<su2double> &Zcoord_Airfoil, bool original_surface) {

SU2_CFD/include/iteration_structure.hpp

@@ -933,6 +933,19 @@ class CDiscAdjMeanFlowIteration : public CIteration {
                        unsigned short iZone,
                        unsigned short kind_recording);
 
+  /*!
+   * \brief load unsteady solution for unsteady problems
+   * \param[in] geometry_container - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config_container - Definition of the particular problem.
+   * \param[in] val_iZone - Index of the zone.
+   * \param[in] val_DirectIter - Direct iteration to load.
+   */
+  void LoadUnsteady_Solution(CGeometry ***geometry_container,
+                      CSolver ****solver_container,
+                      CConfig **config_container,
+                      unsigned short val_iZone,
+                      int val_DirectIter);
 };
 
 

SU2_CFD/include/solver_structure.hpp

@@ -820,15 +820,15 @@ class CSolver {
 	virtual void BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics,
                           CConfig *config, unsigned short val_marker);
 
-	/*!
-	 * \brief A virtual member.
-	 * \param[in] geometry - Geometrical definition of the problem.
-	 * \param[in] solver_container - Container vector with all the solutions.
+    /*!
+     * \brief A virtual member.
+     * \param[in] geometry - Geometrical definition of the problem.
+     * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] conv_numerics - Description of the numerical method.
    * \param[in] visc_numerics - Description of the numerical method.
-	 * \param[in] config - Definition of the particular problem.
-	 * \param[in] val_marker - Surface marker where the boundary condition is applied.
-	 */
+     * \param[in] config - Definition of the particular problem.
+     * \param[in] val_marker - Surface marker where the boundary condition is applied.
+     */
 	virtual void BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_container,
                                      CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker);
 
@@ -2829,6 +2829,12 @@ class CSolver {
    * \param[in] config - Definition of the particular problem.
    */
   virtual void ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config);
+
+  /*!
+   * \brief A virtual member.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void SetFreeStream_Solution(CConfig *config);
 };
 
 /*!
@@ -3548,7 +3554,7 @@ class CEulerSolver : public CSolver {
 	 */
 	void BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics,
                   CConfig *config, unsigned short val_marker);
-    
+
 	/*!
 	 * \brief Impose a supersonic inlet boundary condition.
 	 * \param[in] geometry - Geometrical definition of the problem.
@@ -4630,6 +4636,12 @@ class CEulerSolver : public CSolver {
    * \param[in] Value of freestream temperature.
    */
   void SetTemperature_Inf(su2double t_inf);
+
+  /*!
+   * \brief Set the solution using the Freestream values.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetFreeStream_Solution(CConfig *config);
 };
 
 /*!
@@ -5033,6 +5045,15 @@ class CTurbSolver : public CSolver {
 	void SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config,
                             unsigned short iRKStep, unsigned short iMesh, unsigned short RunTime_EqSystem);
 
+  /*!
+   * \brief Load a solution from a restart file.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver - Container vector with all of the solvers.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_iter - Current external iteration number.
+   */
+  void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter);
+
 };
 
 /*!
@@ -5227,17 +5248,13 @@ class CTurbSASolver: public CTurbSolver {
 	 */
 	void BC_NearField_Boundary(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics,
                                CConfig *config);
-  
+
   /*!
-	 * \brief Load a solution from a restart file.
-	 * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - Container vector with all of the solvers.
-	 * \param[in] config - Definition of the particular problem.
-	 * \param[in] val_iter - Current external iteration number.
-	 */
-	void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter);
+   * \brief Set the solution using the Freestream values.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetFreeStream_Solution(CConfig *config);
 
-
 };
 
 /*!
@@ -5559,6 +5576,12 @@ class CTurbSSTSolver: public CTurbSolver {
 	 * \return A pointer to an array containing a set of constants
 	 */
 	su2double* GetConstants();
+
+  /*!
+   * \brief Set the solution using the Freestream values.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetFreeStream_Solution(CConfig *config);
 
 };
 
@@ -7820,5 +7843,17 @@ class CDiscAdjSolver : public CSolver {
    * \param[in] config - Definition of the particular problem.
    */
   void ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config);
+
+	/*!
+	 * \brief Update the dual-time derivatives.
+	 * \param[in] geometry - Geometrical definition of the problem.
+	 * \param[in] solver_container - Container vector with all the solutions.
+	 * \param[in] config - Definition of the particular problem.
+	 * \param[in] iMesh - Index of the mesh in multigrid computations.
+	 * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
+	* \param[in] RunTime_EqSystem - System of equations which is going to be solved.
+	* \param[in] Output - boolean to determine whether to print output.
+	 */
+  void Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output);
 };
 #include "solver_structure.inl"

SU2_CFD/include/solver_structure.inl

@@ -547,7 +547,7 @@ inline void CSolver::BC_NonReflecting(CGeometry *geometry, CSolver **solver_cont
 
 inline void CSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, 
 										 CConfig *config, unsigned short val_marker) { }
-                     
+
 inline void CSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, 
 										  CConfig *config, unsigned short val_marker) { }
 
@@ -1167,3 +1167,17 @@ inline void CSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo
 
 inline void CSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config){}
 
+inline void CSolver::SetFreeStream_Solution(CConfig *config){}
+
+inline void CTurbSASolver::SetFreeStream_Solution(CConfig *config){
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++)
+    node[iPoint]->SetSolution(0, nu_tilde_Inf);
+}
+
+inline void CTurbSSTSolver::SetFreeStream_Solution(CConfig *config){
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
+    node[iPoint]->SetSolution(0, kine_Inf);
+    node[iPoint]->SetSolution(1, omega_Inf);
+  }
+}
+

SU2_CFD/src/SU2_CFD.cpp

@@ -258,7 +258,7 @@ int main(int argc, char *argv[]) {
           (config_container[iZone]->GetKind_Solver() == DISC_ADJ_RANS))
         geometry_container[iZone][MESH_0]->ComputeWall_Distance(config_container[iZone]);
     }
-    
+
 
   }