Aero-Structural Discrete Adjoint Sensitivities and Python Wrapper Developments

AUTHORS.md

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

Common/include/CConfig.hpp

@@ -1095,8 +1095,9 @@ class CConfig {
   long ParMETIS_pointWgt;           /*!< \brief Load balancing weight given to points. */
   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. */
-  DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
+  bool DiscreteAdjoint;                /*!< \brief AD-based discrete adjoint mode. */
+  ENUM_DISC_ADJ_TYPE Kind_DiscreteAdjoint; /*!< \brief AD-based discrete adjoint formulation. */
+  bool DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
   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. */
@@ -9188,6 +9189,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 indicator whether a debug run for the discrete adjoint solver will be started.
    * \return the discrete adjoint debug indicator.

Common/include/geometry/CGeometry.hpp

@@ -1400,6 +1400,16 @@ class CGeometry {
    */
   inline virtual void SetSensitivity(unsigned long iPoint, unsigned short iDim, su2double val) {}
 
+  /*!
+   * \brief Get a reference to the sensitivity matrix.
+   * \return Reference to the sensitivity matrix.
+   */
+  inline virtual const su2activematrix& GetSensitivityMatrix() const {
+    SU2_MPI::Error("GetSensitivityMatrix not implemented for this geometry type.", CURRENT_FUNCTION);
+    static su2activematrix sens_matrix;
+    return sens_matrix;
+  }
+
   /*!
    * \brief Get the average normal at a specific span for a given marker in the turbomachinery reference of frame.
    * \param[in] val_marker - marker value.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -779,6 +779,12 @@ class CPhysicalGeometry final : public CGeometry {
     Sensitivity(iPoint, iDim) = val;
   }
 
+  /*!
+   * \brief Get a reference to the sensitivity matrix.
+   * \return Reference to the sensitivity matrix.
+   */
+  inline const su2activematrix& GetSensitivityMatrix() const override { return Sensitivity; }
+
   /*!
    * \brief Check the mesh for periodicity and deactivate multigrid if periodicity is found.
    * \param[in] config - Definition of the particular problem.

Common/include/option_structure.hpp

@@ -2798,6 +2798,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

@@ -1140,6 +1140,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.4.0 "Harrier"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -28,6 +28,19 @@
 
 #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"
+
+#ifdef CODI_FORWARD_TYPE
+  using Scalar = su2double;
+#else
+  using Scalar = su2mixedfloat;
+#endif
+
+class LinOperator;
+class LinPreconditioner;
 
 /*!
  * \class CDiscAdjSinglezoneDriver
@@ -55,6 +68,24 @@ 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;
+  LinOperator* AdjOperator = nullptr;
+  LinPreconditioner* AdjPreconditioner = nullptr;
+
   /*!
    * \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,27 +96,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 Run one iteration of the solver.
+   * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp)
+   */
+  void DirectRunResidual(RECORDING kind_recording);
+
+  /*!
+   * \brief Run a single iteration of the main fixed-point discrete adjoint solver with a single zone.
+   */
+  void RunFixedPoint();
+
+  /*!
+   * \brief Run the computation of the main residual-based discrete adjoint sensitivities with a single zone.
+   */
+  void RunResidual();
+
+  /*!
+   * \brief Run a single iteration of the secondary fixed-point discrete adjoint solver with a single zone.
+   */
+  void SecondaryRunFixedPoint();
+
+  /*!
+   * \brief Run the computation of the secondary residual-based discrete adjoint sensitivities with a single zone.
+   */
+  void SecondaryRunResidual();
+
+  /*!
+   * \brief Update the fixed-point discrete adjoint solver with a single zone.
+   */
+  void UpdateAdjointsFixedPoint();
 
   /*!
-   * \brief Set the objective function.
+   * \brief Update the residual-based discrete adjoint solver with a single zone.
    */
-  void SetObjFunction(void);
+  void UpdateAdjointsResidual();
 
   /*!
    * \brief Initialize the adjoint value of the objective function.
    */
-  void SetAdjObjFunction(void);
+  void SetAdjointObjective();
 
   /*!
    * \brief Record the main computational path.
    */
-  void MainRecording(void);
+  void MainRecording();
 
   /*!
    * \brief Record the secondary computational path.
    */
-  void SecondaryRecording(void);
+  void SecondaryRecording();
 
   /*!
    * \brief gets Convergence on physical time scale, (deactivated in adjoint case)
@@ -99,7 +161,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,
@@ -120,10 +181,144 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
   /*!
    * \brief Run a single iteration of the discrete adjoint solver with a single zone.
    */
-  void Run(void) override;
+  void Run() override;
+
+  /*!
+   * \brief Update the discrete adjoint solution.
+   */
+  void UpdateAdjoints();
+
+  /*!
+   * \brief Update the primal time iteration.
+   */
+  void UpdateTimeIter();
+
+  /*!
+   * \brief Update the primal far-field variables.
+   */
+  void UpdateFarfield();
+
+  /*!
+   * \brief Update the primal geometry (does not include mesh deformation).
+   */
+  void UpdateGeometry();
+
+  /*!
+   * \brief Deform the primal mesh.
+   */
+  void DeformGeometry();
+
+  /*!
+   * \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;
+  void Postprocess() 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;
+
+      const auto& mat = solver->GetPartialMatrix_dObjective_dStates();
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) {
+          values(iPoint, offset + iVar) = -SU2_TYPE::GetValue(mat(iPoint, iVar));
+        }
+      }
+
+      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;
+
+      const auto& mat = solver->GetPartialMatrix_dResiduals_dStates();
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) {
+          values(iPoint, offset + iVar) = SU2_TYPE::GetValue(mat(iPoint, iVar));
+        }
+      }
+
+      offset += solver->GetnVar();
+    }
+  }
+
+  /*!
+   * \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);
+};
+
+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);
+  }
 };