New Discrete Adjoint Multizone Capabilities (e.g. Adjoint CHT)

Common/include/ad_structure.hpp

@@ -65,17 +65,42 @@ namespace AD{
   bool TapeActive();
 
   /*!
-   * \brief Registers the variable as an input. I.e. as a leaf of the computational graph.
+   * \brief Prints out tape statistics.
+   */
+  void PrintStatistics();
+
+  /*!
+   * \brief Registers the variable as an input and saves internal data (indices). I.e. as a leaf of the computational graph.
    * \param[in] data - The variable to be registered as input.
    */
   void RegisterInput(su2double &data);
 
+  /*!
+   * \brief Registers the variable as an input. I.e. as a leaf of the computational graph.
+   * \param[in] data - The variable to be registered as input.
+   */
+  void RegisterInput_intIndexBased(su2double &data);
+
   /*!
    * \brief Registers the variable as an output. I.e. as the root of the computational graph.
    * \param[in] data - The variable to be registered as output.
    */
   void RegisterOutput(su2double &data);
 
+  /*!
+   * \brief Sets the adjoint value at index to val
+   * \param[in] index - Position in the adjoint vector.
+   * \param[in] val - adjoint value to be set.
+   */
+  void SetDerivative(int index, const double val);
+
+  /*!
+   * \brief Extracts the adjoint value at index
+   * \param[in] index - position in the adjoint vector where the derivative will be extracted.
+   * \return Derivative value.
+   */
+  double GetDerivative(int index);
+
   /*!
    * \brief Clears the currently stored adjoints but keeps the computational graph.
    */
@@ -84,7 +109,14 @@ namespace AD{
   /*!
    * \brief Computes the adjoints, i.e. the derivatives of the output with respect to the input variables.
    */
-  void ComputeAdjoint();
+  void ComputeAdjoint();  
+
+  /*!
+   * \brief Computes the adjoints, i.e. the derivatives of the output with respect to the input variables.
+   * \param[in] enter - Position where we start evaluating the tape.
+   * \param[in] leave - Position where we stop evaluating the tape.
+   */
+  void ComputeAdjoint(unsigned short enter, unsigned short leave);
 
   /*!
    * \brief Reset the tape structure to be ready for a new recording.
@@ -209,6 +241,17 @@ namespace AD{
    */
   void EndExtFunc();
 
+  /*!
+   * \brief Evaluates and saves gradient data from a variable.
+   * \param[in] data - variable whose gradient information will be extracted.
+   * \param[in] index - where obtained gradient information will be stored.
+   */
+  void SetAdjIndex(int &index, const su2double &data);
+
+  /*!
+   * \brief Pushes back the current tape position to the tape position's vector.
+   */
+  void Push_TapePosition();
 }
 
 /*--- Macro to begin and end sections with a passive tape ---*/

Common/include/ad_structure.inl

@@ -68,6 +68,8 @@ namespace AD{
 
   extern su2double::TapeType::Position StartPosition, EndPosition;
 
+  extern std::vector<su2double::TapeType::Position> TapePositions;
+
   extern std::vector<su2double::GradientData> localInputValues;
 
   extern std::vector<su2double*> localOutputValues;
@@ -77,6 +79,8 @@ namespace AD{
   inline void RegisterInput(su2double &data) {AD::globalTape.registerInput(data);
                                              inputValues.push_back(data.getGradientData());}
 
+  inline void RegisterInput_intIndexBased(su2double &data) {AD::globalTape.registerInput(data);}
+
   inline void RegisterOutput(su2double& data) {AD::globalTape.registerOutput(data);}
 
   inline void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();}
@@ -85,19 +89,43 @@ namespace AD{
 
   inline void StopRecording() {AD::globalTape.setPassive();}
 
-  inline bool TapeActive() {return AD::globalTape.isActive();}
+  inline bool TapeActive() { return AD::globalTape.isActive(); }
+
+  inline void PrintStatistics() {AD::globalTape.printStatistics();}
 
   inline void ClearAdjoints() {AD::globalTape.clearAdjoints(); }
 
   inline void ComputeAdjoint() {AD::globalTape.evaluate();
                                adjointVectorPosition = 0;}
 
+  inline void ComputeAdjoint(unsigned short enter, unsigned short leave) {
+    AD::globalTape.evaluate(TapePositions[enter], TapePositions[leave]);
+    if (leave == 0) {
+      adjointVectorPosition = 0;
+    }
+  }
+
   inline void Reset() {
+    globalTape.reset();
     if (inputValues.size() != 0) {
-      globalTape.reset();
       adjointVectorPosition = 0;
       inputValues.clear();
     }
+    if (TapePositions.size() != 0) {
+      TapePositions.clear();
+    }    
+  }
+
+  inline void SetAdjIndex(int &index, const su2double &data) {
+    index = data.getGradientData();
+  }
+
+  inline void SetDerivative(int index, const double val) {
+    AD::globalTape.setGradient(index, val);
+  }
+
+  inline double GetDerivative(int index) {
+    return AD::globalTape.getGradient(index);
   }
 
   inline void SetPreaccIn(const su2double &data) {
@@ -167,6 +195,10 @@ namespace AD{
     }
   }
 
+  inline void Push_TapePosition() {
+    TapePositions.push_back(AD::globalTape.getPosition());
+  }
+
   inline void EndPreacc(){
     if (PreaccActive) {
       PreaccHelper.finish(false);
@@ -239,18 +271,30 @@ namespace AD{
 
   inline void RegisterInput(su2double &data) {}
 
+  inline void RegisterInput_intIndexBased(su2double &data) {}
+
   inline void RegisterOutput(su2double& data) {}
 
   inline void StartRecording() {}
 
   inline void StopRecording() {}
 
-  inline bool TapeActive() {return false;}
+  inline bool TapeActive() { return false; }
+
+  inline void PrintStatistics() {}
 
   inline void ClearAdjoints() {}
 
   inline void ComputeAdjoint() {}
 
+  inline void ComputeAdjoint(unsigned short enter, unsigned short leave) {}
+
+  inline void SetAdjIndex(int &index, const su2double &data) {}
+
+  inline void SetDerivative(int index, const double val) {}
+
+  inline double GetDerivative(int position) { return 0.0; }
+
   inline void Reset() {}
 
   inline void ResetInput(su2double &data) {}
@@ -270,7 +314,9 @@ namespace AD{
   inline void StartPreacc() {}
 
   inline void EndPreacc() {}
-
+
+  inline void Push_TapePosition() {}
+
   inline void StartExtFunc(bool storePrimalInput, bool storePrimalOutput){}
 
   inline void SetExtFuncIn(const su2double &data) {}

Common/include/config_structure.hpp

@@ -417,7 +417,6 @@ class CConfig {
   CFLRedCoeff_AdjFlow,	/*!< \brief CFL reduction coefficient for the adjoint problem. */
   CFLRedCoeff_AdjTurb,	/*!< \brief CFL reduction coefficient for the adjoint problem. */
   CFLFineGrid,		/*!< \brief CFL of the finest grid. */
-  CFLSolid,       /*!< \brief CFL in (heat) solid solvers. */
   Max_DeltaTime,  		/*!< \brief Max delta time. */
   Unst_CFL;		/*!< \brief Unsteady CFL number. */
   bool ReorientElements;		/*!< \brief Flag for enabling element reorientation. */
@@ -758,6 +757,7 @@ class CConfig {
   Wrt_SharpEdges,              /*!< \brief Write residuals to solution file */
   Wrt_Halo,                   /*!< \brief Write rind layers in solution files */
   Wrt_Performance,            /*!< \brief Write the performance summary at the end of a calculation.  */
+  Wrt_AD_Statistics,          /*!< \brief Write the tape statistics (discrete adjoint).  */
   Wrt_MeshQuality,            /*!< \brief Write the mesh quality statistics to the visualization files.  */
   Wrt_InletFile,                   /*!< \brief Write a template inlet profile file */
   Wrt_Slice,                   /*!< \brief Write 1D slice of a 2D cartesian solution */
@@ -874,6 +874,7 @@ class CConfig {
   bool Sine_Load;             /*!< \brief option for sine load */
   su2double *SineLoad_Coeff;  /*!< \brief Stores the load coefficient */
   su2double Thermal_Diffusivity;			/*!< \brief Thermal diffusivity used in the heat solver. */
+  bool CHT_Robin;             /*!< \brief Option for boundary condition method at CHT interfaces. */
   su2double Cyclic_Pitch,     /*!< \brief Cyclic pitch for rotorcraft simulations. */
   Collective_Pitch;           /*!< \brief Collective pitch for rotorcraft simulations. */
   su2double Mach_Motion;			/*!< \brief Mach number based on mesh velocity and freestream quantities. */
@@ -975,7 +976,8 @@ class CConfig {
   Nonphys_Reconstr;      /*!< \brief Current number of non-physical reconstructions for 2nd-order upwinding. */
   bool ParMETIS;      /*!< \brief Boolean for activating ParMETIS mode (while testing). */
   unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */
-  bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */
+  bool DiscreteAdjoint,       /*!< \brief AD-based discrete adjoint mode. */
+  FullTape;                     /*!< \brief Full tape mode for coupled discrete adjoints. */
   unsigned long Wrt_Surf_Freq_DualTime;	/*!< \brief Writing surface solution frequency for Dual Time. */
   su2double Const_DES;   /*!< \brief Detached Eddy Simulation Constant. */
   unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */
@@ -1404,7 +1406,7 @@ class CConfig {
   /*!
    * \brief Constructor of the class which reads the input file.
    */
-  CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, unsigned short val_nZone, bool verb_high);
+  CConfig(char case_filename[MAX_STRING_SIZE], unsigned short val_software, bool verb_high);
 
   /*!
    * \brief Constructor of the class which reads the input file and uses default options from another config.
@@ -2757,13 +2759,6 @@ class CConfig {
    * \return CFL number for each grid.
    */
   su2double GetCFL(unsigned short val_mesh);
-
-  /*!
-   * \brief Get the Courant Friedrich Levi number for solid solvers.
-   * \param[in] val_mesh - Index of the mesh were the CFL is applied.
-   * \return CFL number for each grid.
-   */
-  su2double GetCFL_Solid(void);
 
   /*!
    * \brief Get the Courant Friedrich Levi number for each grid.
@@ -3252,6 +3247,12 @@ class CConfig {
    * \return <code>TRUE</code> means that the performance summary will be written at the end of a calculation.
    */
   bool GetWrt_Performance(void);
+
+  /*!
+   * \brief Get information about the computational graph (e.g. memory usage) when using AD in reverse mode.
+   * \return <code>TRUE</code> means that the tape statistics will be written after each recording.
+   */
+  bool GetWrt_AD_Statistics(void);
 
   /*!
    * \brief Get information about writing the mesh quality metrics to the visualization files.
@@ -8434,6 +8435,12 @@ class CConfig {
    * \return the discrete adjoint indicator.
    */
   bool GetDiscrete_Adjoint(void);
+
+  /*!
+  * \brief Get the indicator whether we want to use full (coupled) tapes.
+  * \return the full tape indicator.
+  */
+  bool GetFull_Tape(void);
 
   /*!
    * \brief Get the indicator whether we want to benchmark the MPI performance of FSI problems
@@ -8940,6 +8947,12 @@ class CConfig {
    */
   bool GetWeakly_Coupled_Heat(void);
 
+  /*!
+   * \brief Get the boundary condition method for CHT.
+   * \return YES if Robin BC is used.
+   */
+  bool GetCHT_Robin(void);
+
   /*!
    * \brief Check if values passed to the BC_HeatFlux-Routine are already integrated.
    * \return YES if the passed values is the integrated heat flux over the marker's surface.

Common/include/config_structure.inl

@@ -753,8 +753,6 @@ inline unsigned short CConfig::GetGeometryMode(void) { return GeometryMode; }
 
 inline su2double CConfig::GetCFL(unsigned short val_mesh) {	return CFL[val_mesh]; }
 
-inline su2double CConfig::GetCFL_Solid(void) {	return CFLSolid; }
-
 inline void CConfig::SetCFL(unsigned short val_mesh, su2double val_cfl) { CFL[val_mesh] = val_cfl; }
 
 inline su2double CConfig::GetUnst_CFL(void) {	return Unst_CFL; }
@@ -1669,6 +1667,8 @@ inline bool CConfig::GetWrt_Halo(void) { return Wrt_Halo; }
 
 inline bool CConfig::GetWrt_Performance(void) { return Wrt_Performance; }
 
+inline bool CConfig::GetWrt_AD_Statistics(void) { return Wrt_AD_Statistics; }
+
 inline bool CConfig::GetWrt_MeshQuality(void) { return Wrt_MeshQuality; }
 
 inline bool CConfig::GetWrt_InletFile(void) { return Wrt_InletFile; }
@@ -1905,6 +1905,8 @@ inline unsigned short CConfig::GetDirectDiff() { return DirectDiff;}
 
 inline bool CConfig::GetDiscrete_Adjoint() { return DiscreteAdjoint;}
 
+inline bool CConfig::GetFull_Tape() { return FullTape; }
+
 inline unsigned short CConfig::GetRiemann_Solver_FEM(void) {return Riemann_Solver_FEM;}
 
 inline su2double CConfig::GetQuadrature_Factor_Straight(void) {return Quadrature_Factor_Straight;}
@@ -1925,6 +1927,8 @@ inline bool CConfig::GetJacobian_Spatial_Discretization_Only(void) {return Jacob
 
 inline bool CConfig::GetWeakly_Coupled_Heat(void) { return Weakly_Coupled_Heat; }
 
+inline bool CConfig::GetCHT_Robin(void) { return CHT_Robin; }
+
 inline bool CConfig::GetIntegrated_HeatFlux(void) { return Integrated_HeatFlux; }
 
 inline bool CConfig::GetAD_Mode(void) { return AD_Mode;}

Common/include/dual_grid_structure.hpp

@@ -176,6 +176,8 @@ class CPoint : public CDualGrid {
   unsigned short nNeighbor;           /*!< \brief Number of neighbors. */
   bool Flip_Orientation;              /*!< \brief Flip the orientation of the normal. */
   su2double MaxLength;                /*!< \brief The maximum cell-center to cell-center length. */
+  int *Input_AdjIndices,              /*!< \brief Indices of Coord variables in the adjoint vector. */
+  *Output_AdjIndices;                 /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */
 
 public:
 
@@ -289,7 +291,25 @@ class CPoint : public CDualGrid {
 	 * \param[in] val_dim - Position to store the coordinate.		 
 	 * \param[in] val_coord - Coordinate for val_dim.			 
 	 */
-	void SetCoord(unsigned short val_dim, su2double val_coord);
+  void SetCoord(unsigned short val_dim, su2double val_coord);
+
+  /*!
+   * \brief Set the adjoint vector indices of Coord vector.
+   * \param[in] input - Save them to the input or output indices vector.
+   */
+  void SetAdjIndices(bool input);
+
+  /*!
+   * \brief Set the adjoint values of the (geometric) coordinates.
+   * \param[in] adj_sol - Adjoint values of the Coord variables.
+   */
+  void SetAdjointSolution(const su2double *adj_sol);
+
+  /*!
+   * \brief Get the adjoint values of the (geometric) coordinates.
+   * \param[in] adj_sol - Adjoint values of the Coord variables.
+   */
+  su2double GetAdjointSolution(unsigned short iDim);
 
   /*!
 	 * \brief Get the coordinates of the control volume.

Common/src/ad_structure.cpp

@@ -49,6 +49,7 @@ namespace AD {
 
   su2double::TapeType& globalTape = su2double::getGlobalTape();
   su2double::TapeType::Position StartPosition, EndPosition;
+  std::vector<su2double::TapeType::Position> TapePositions;
 
   bool Status = false;
   bool PreaccActive = false;