[WIP] Aero-Structural Discrete Adjoint Sensitivities and Python Wrapper Developments

AUTHORS.md

@@ -79,6 +79,7 @@ Giacomo Baldan
 Giulio Gori
 Guillaume Bâty
 Harichand M V
+Harsh Patel
 HL Kline
 HL Zhi
 IndianaStokes

Common/include/CConfig.hpp

@@ -1042,6 +1042,7 @@ class CConfig {
   long ParMETIS_edgeWgt;            /*!< \brief Load balancing weight given to edges. */
   unsigned short DirectDiff;        /*!< \brief Direct Differentation mode. */
   bool DiscreteAdjoint;                /*!< \brief AD-based discrete adjoint mode. */
+  ENUM_DISC_ADJ_TYPE Kind_DiscreteAdjoint; /*!< \brief AD-based discrete adjoint formulation. */
   su2double Const_DES;                 /*!< \brief Detached Eddy Simulation Constant. */
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
@@ -8798,6 +8799,12 @@ class CConfig {
    */
   bool GetDiscrete_Adjoint(void) const { return DiscreteAdjoint; }
 
+  /*!
+   * \brief Get the kind of discrete adjoint solver formulation.
+   * \return the discrete adjoint indicator.
+   */
+  ENUM_DISC_ADJ_TYPE GetKind_DiscreteAdjoint(void) const { return Kind_DiscreteAdjoint; }
+
   /*!
    * \brief Get the number of subiterations while a ramp is applied.
    * \return Number of FSI subiters.

Common/include/option_structure.hpp

@@ -2593,6 +2593,25 @@ static const MapType<std::string, ENUM_STREAMWISE_PERIODIC> Streamwise_Periodic_
   MakePair("MASSFLOW",      ENUM_STREAMWISE_PERIODIC::MASSFLOW)
 };
 
+/*!
+ * \brief Types of discrete adjoint solver formulations.
+ */
+enum class ENUM_DISC_ADJ_TYPE {
+  FIXED_POINT, /*!< \brief Fixed-point discrete-adjoint formulation. */
+  RESIDUALS    /*!< \brief Residual-based discrete-adjoint formulation. */
+};
+static const MapType<std::string, ENUM_DISC_ADJ_TYPE> DiscreteAdjoint_Map = {
+  MakePair("FIXED_POINT", ENUM_DISC_ADJ_TYPE::FIXED_POINT)
+  MakePair("RESIDUALS",   ENUM_DISC_ADJ_TYPE::RESIDUALS)
+};
+
+enum class ENUM_VARIABLE {
+  RESIDUALS,
+  OBJECTIVE,
+  TRACTIONS,
+  COORDINATES
+};
+
 /*!
  * \brief Container to hold Variables for streamwise Periodic flow as they are often used together in places.
  */

Common/src/CConfig.cpp

@@ -1110,6 +1110,8 @@ void CConfig::SetConfig_Options() {
 #endif
   /*!\brief MATH_PROBLEM  \n DESCRIPTION: Mathematical problem \n  Options: DIRECT, ADJOINT \ingroup Config*/
   addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault);
+  /*!\brief DISC_ADJ_TYPE  \n DESCRIPTION: Discrete adjoint problem formulation \n  Options: FIXED_POINT, RESIDUALS \ingroup Config*/
+  addEnumOption("KIND_DISC_ADJ", Kind_DiscreteAdjoint, DiscreteAdjoint_Map, ENUM_DISC_ADJ_TYPE::FIXED_POINT);
   /*!\brief KIND_TURB_MODEL \n DESCRIPTION: Specify turbulence model \n Options: see \link Turb_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("KIND_TURB_MODEL", Kind_Turb_Model, Turb_Model_Map, TURB_MODEL::NONE);
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SST turbulence model options/corrections. \n Options: see \link SST_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/

SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp

@@ -2,7 +2,7 @@
  * \file CDiscAdjSinglezoneDriver.hpp
  * \brief Headers of the main subroutines for driving single or multi-zone problems.
  *        The subroutines and functions are in the <i>driver_structure.cpp</i> file.
- * \author T. Economon, H. Kline, R. Sanchez
+ * \author T. Economon, H. Kline, R. Sanchez, H. Patel, A. Gastaldi
  * \version 8.0.1 "Harrier"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -28,6 +28,10 @@
 
 #pragma once
 #include "CSinglezoneDriver.hpp"
+#include "../../../Common/include/toolboxes/CQuasiNewtonInvLeastSquares.hpp"
+#include "../../../Common/include/linear_algebra/CPreconditioner.hpp"
+#include "../../../Common/include/linear_algebra/CMatrixVectorProduct.hpp"
+#include "../../../Common/include/linear_algebra/CSysSolve.hpp"
 
 /*!
  * \class CDiscAdjSinglezoneDriver
@@ -38,6 +42,11 @@
  */
 class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
 protected:
+#ifdef CODI_FORWARD_TYPE
+    using Scalar = su2double;
+#else
+    using Scalar = su2mixedfloat;
+#endif
 
   unsigned long nAdjoint_Iter;                  /*!< \brief The number of adjoint iterations that are run on the fixed-point solver.*/
   RECORDING RecordingState;                     /*!< \brief The kind of recording the tape currently holds.*/
@@ -55,6 +64,43 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
   COutput *direct_output;
   CNumerics ***numerics;                        /*!< \brief Container vector with all the numerics. */
 
+  /*!< \brief Fixed-Point corrector that can be applied to inner iterations of the residual-based formulation. */
+  CQuasiNewtonInvLeastSquares<passivedouble> Corrector;
+
+  /*!< \brief Members to use GMRES to drive inner iterations (alternative to quasi-Newton) of the residual-based formulation. */
+  static constexpr unsigned long KrylovMinIters = 3;
+  const Scalar KrylovSysTol = 0.00001;
+  const Scalar KrylovPreTol = 0.1;
+  bool KrylovSet = false;
+
+  CSysMatrix<Scalar> CopiedJacobian;
+  CSysSolve<Scalar> AdjSolver;
+  CSysVector<Scalar> AdjRHS;
+  CSysVector<Scalar> AdjSol;
+  CPreconditioner<Scalar>* PrimalPreconditioner = nullptr;
+  CSysMatrixVectorProduct<Scalar>* PrimalJacobian = nullptr;
+
+  class LinOperator : public CMatrixVectorProduct<Scalar> {
+  public:
+      CDiscAdjSinglezoneDriver* const driver;
+      LinOperator(CDiscAdjSinglezoneDriver* d) : driver(d) { }
+
+      inline void operator()(const CSysVector<Scalar> & u, CSysVector<Scalar> & v) const override {
+        driver->ApplyOperator(u, v);
+      }
+  };
+
+  class LinPreconditioner : public CPreconditioner<Scalar> {
+   public:
+    CDiscAdjSinglezoneDriver* const driver;
+
+    LinPreconditioner(CDiscAdjSinglezoneDriver* d): driver(d) { }
+
+    inline void operator()(const CSysVector<Scalar> & u, CSysVector<Scalar> & v) const override {
+      driver->ApplyPreconditioner(u, v);
+    }
+  };
+
   /*!
    * \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)
@@ -65,17 +111,58 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
    * \brief Run one iteration of the solver.
    * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp)
    */
-  void DirectRun(RECORDING kind_recording);
+  void DirectRunFixedPoint(RECORDING kind_recording);
 
   /*!
-   * \brief Set the objective function.
+   * \brief Run one iteration of the solver.
+   * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp)
    */
-  void SetObjFunction(void);
+  void DirectRunResidual(RECORDING kind_recording);
+
+  /*!
+   * \brief Run a single iteration of the main fixed-point discrete adjoint solver with a single zone.
+   */
+  void RunFixedPoint(void);
+
+  /*!
+   * \brief Run the computation of the main residual-based discrete adjoint sensitivities with a single zone.
+   */
+  void RunResidual(void);
+
+  /*!
+   * \brief Run a single iteration of the secondary fixed-point discrete adjoint solver with a single zone.
+   */
+  void SecondaryRunFixedPoint(void);
+
+  /*!
+   * \brief Run the computation of the secondary residual-based discrete adjoint sensitivities with a single zone.
+   */
+  void SecondaryRunResidual(void);
+
+  /*!
+   * \brief Update the fixed-point discrete adjoint solver with a single zone.
+   */
+  void UpdateAdjointsFixedPoint(void);
+
+  /*!
+   * \brief Update the residual-based discrete adjoint solver with a single zone.
+   */
+  void UpdateAdjointsResidual(void);
+
+  /*!
+   * \brief Adjoint problem Jacobian-vector product.
+   */
+  void ApplyOperator(const CSysVector<Scalar>& u, CSysVector<Scalar>& v);
+
+  /*!
+   * \brief Adjoint problem preconditioner (based on the transpose approximate Jacobian of the primal problem).
+   */
+  void ApplyPreconditioner(const CSysVector<Scalar>& u, CSysVector<Scalar>& v);
 
   /*!
    * \brief Initialize the adjoint value of the objective function.
    */
-  void SetAdjObjFunction(void);
+  void SetAdjointObjective(void);
 
   /*!
    * \brief Record the main computational path.
@@ -99,7 +186,6 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
    * \brief Constructor of the class.
    * \param[in] confFile - Configuration file name.
    * \param[in] val_nZone - Total number of zones.
-   * \param[in] val_nDim - Total number of dimensions.
    * \param[in] MPICommunicator - MPI communicator for SU2.
    */
   CDiscAdjSinglezoneDriver(char* confFile,
@@ -122,8 +208,114 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
    */
   void Run(void) override;
 
+  /*!
+   * \brief Update the discrete adjoint solution.
+   */
+  void UpdateAdjoints(void);
+
+  /*!
+   * \brief Update the primal time iteration.
+   */
+  void UpdateTimeIter(void);
+
+  /*!
+   * \brief Update the primal far-field variables.
+   */
+  void UpdateFarfield(void);
+
+  /*!
+   * \brief Update the primal geometry (does not include mesh deformation).
+   */
+  void UpdateGeometry(void);
+
+  /*!
+   * \brief Deform the primal mesh.
+   */
+  void DeformGeometry(void);
+
+  /*!
+   * \brief Update the primal states.
+   */
+  void UpdateStates();
+
+  /*!
+   * \brief Update the primal residuals.
+   */
+  void UpdateResiduals();
+
+  /*!
+   * \brief Update the primal tractions.
+   */
+  void UpdateTractions();
+
+  /*!
+   * \brief Update the primal objective.
+   */
+  void UpdateObjective();
+
+  /*!
+   * \brief Update the primal jacobian.
+   */
+  void UpdateJacobians();
+
   /*!
    * \brief Postprocess the adjoint iteration for ZONE_0.
    */
   void Postprocess(void) override;
+
+  /*!
+   * \brief Get the partial objective sensitivities from all solvers.
+   * \param[out] values - Values object with interface (iPoint, iVar).
+   */
+  template <class Container>
+  void GetAllObjectiveStatesSensitivities(Container& values) const {
+    const auto nPoint = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint();
+
+    /*--- Get all the partial sensitivities (dObjective/dStates) ---*/
+    unsigned short offset = 0;
+
+    for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+      auto solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+
+      if (!solver || !solver->GetAdjoint()) continue;
+
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) {
+          auto value = solver->GetSens_dObjective_dStates(iPoint, iVar);
+
+          values(iPoint, offset + iVar) = -SU2_TYPE::GetValue(value);
+        }
+      }
+
+      offset += solver->GetnVar();
+    }
+  }
+
+  /*!
+   * \brief Get the partial jacobian-adjoint products from all solvers.
+   * \param[out] values - Values object with interface (iPoint, iVar).
+   */
+  template <class Container>
+  void GetAllResidualsStatesSensitivities(Container& values) const {
+    const auto nPoint = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint();
+
+    /*--- Get all the partial jacobian-adjoint products (dResiduals/dStates) ---*/
+    unsigned short offset = 0;
+
+    for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
+      auto solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
+
+      if (!solver || !solver->GetAdjoint()) continue;
+
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) {
+          auto value = solver->GetProd_dResiduals_dStates(iPoint, iVar);
+
+          values(iPoint, offset + iVar) = SU2_TYPE::GetValue(value);
+        }
+      }
+
+      offset += solver->GetnVar();
+    }
+  }
 };