Merge branch 'develop' into feature_custom_source

Author: bigfooted

Common/include/CConfig.hpp

@@ -700,6 +700,7 @@ class CConfig {
   unsigned long StartConv_Iter;       /*!< \brief Start convergence criteria at iteration. */
   su2double Cauchy_Eps;               /*!< \brief Epsilon used for the convergence. */
   bool Restart,                       /*!< \brief Restart solution (for direct, adjoint, and linearized problems).*/
+  Wrt_Restart_Compact,                /*!< \brief Write compact restart files with minimum nr. of variables. */
   Read_Binary_Restart,                /*!< \brief Read binary SU2 native restart files.*/
   Wrt_Restart_Overwrite,              /*!< \brief Overwrite restart files or append iteration number.*/
   Wrt_Surface_Overwrite,              /*!< \brief Overwrite surface output files or append iteration number.*/
@@ -5503,6 +5504,12 @@ class CConfig {
    */
   bool GetRead_Binary_Restart(void) const { return Read_Binary_Restart; }
 
+  /*!
+   * \brief Flag for whether restart files contain only necessary variables.
+   * \return Flag <code>TRUE</code> then the code will write compact restart files.
+   */
+  bool GetWrt_Restart_Compact(void) const { return Wrt_Restart_Compact; }
+
   /*!
    * \brief Flag for whether restart solution files are overwritten.
    * \return Flag for overwriting. If Flag=false, iteration nr is appended to filename

Common/include/linear_algebra/CSysMatrix.hpp

@@ -805,10 +805,10 @@ class CSysMatrix {
   void EnforceSolutionAtNode(unsigned long node_i, const OtherType* x_i, CSysVector<OtherType>& b);
 
   /*!
-   * \brief Version of EnforceSolutionAtNode for a single degree of freedom.
+   * \brief Similar to EnforceSolutionAtNode, but for 0 projection in a given direction.
    */
   template <class OtherType>
-  void EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, CSysVector<OtherType>& b);
+  void EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector<OtherType>& b);
 
   /*!
    * \brief Sets the diagonal entries of the matrix as the sum of the blocks in the corresponding column.

Common/src/CConfig.cpp

@@ -1172,6 +1172,8 @@ void CConfig::SetConfig_Options() {
 
   /*!\brief RESTART_SOL \n DESCRIPTION: Restart solution from native solution file \n Options: NO, YES \ingroup Config */
   addBoolOption("RESTART_SOL", Restart, false);
+  /*!\brief WRT_RESTART_COMPACT \n DESCRIPTION: Minimize the size of restart files \n Options: NO, YES \ingroup Config */
+  addBoolOption("WRT_RESTART_COMPACT", Wrt_Restart_Compact, true);
   /*!\brief BINARY_RESTART \n DESCRIPTION: Read binary SU2 native restart files. \n Options: YES, NO \ingroup Config */
   addBoolOption("READ_BINARY_RESTART", Read_Binary_Restart, true);
   /*!\brief WRT_RESTART_OVERWRITE \n DESCRIPTION: overwrite restart files or append iteration number. \n Options: YES, NO \ingroup Config */

Common/src/linear_algebra/CSysMatrix.cpp

@@ -1027,36 +1027,59 @@ void CSysMatrix<ScalarType>::EnforceSolutionAtNode(const unsigned long node_i, c
 
 template <class ScalarType>
 template <class OtherType>
-void CSysMatrix<ScalarType>::EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i,
-                                                  CSysVector<OtherType>& b) {
+void CSysMatrix<ScalarType>::EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector<OtherType>& b) {
   for (auto index = row_ptr[node_i]; index < row_ptr[node_i + 1]; ++index) {
     const auto node_j = col_ind[index];
 
-    /*--- Delete row iVar of block j on row i (bij) and ATTEMPT
-     *    to delete column iVar block i on row j (bji). ---*/
+    /*--- Remove product components of block j on row i (bij) and ATTEMPT
+     *    to remove solution components of block i on row j (bji).
+     *    This is identical to symmetry correction applied to gradients
+     *    but extended to the entire matrix. ---*/
 
     auto bij = &matrix[index * nVar * nVar];
     auto bji = GetBlock(node_j, node_i);
 
-    /*--- The "attempt" part. ---*/
+    /*--- Attempt to remove solution components. ---*/
+    ScalarType nbn{};
     if (bji != nullptr) {
-      for (auto jVar = 0ul; jVar < nVar; ++jVar) {
-        /*--- Column product. ---*/
-        b[node_j * nVar + jVar] -= bji[jVar * nVar + iVar] * x_i;
-        /*--- Delete entries. ---*/
-        bji[jVar * nVar + iVar] = 0.0;
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        ScalarType proj{};
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          proj += bji[iVar * nVar + jVar] * PassiveAssign(n[jVar]);
+        }
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          bji[iVar * nVar + jVar] -= proj * PassiveAssign(n[jVar]);
+        }
+        nbn += proj * PassiveAssign(n[iVar]);
       }
     }
 
-    /*--- Delete row. ---*/
-    for (auto jVar = 0ul; jVar < nVar; ++jVar) bij[iVar * nVar + jVar] = 0.0;
+    /*--- Product components. ---*/
+    for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+      ScalarType proj{};
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        proj += bij[iVar * nVar + jVar] * PassiveAssign(n[iVar]);
+      }
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        bij[iVar * nVar + jVar] -= proj * PassiveAssign(n[iVar]);
+      }
+    }
 
-    /*--- Set the diagonal entry of the block to 1. ---*/
-    if (node_j == node_i) bij[iVar * (nVar + 1)] = 1.0;
+    /*--- This part doesn't have the "*2" factor because the product components
+     *    were removed from the result of removing the solution components
+     *    instead of from the original block (bji == bij). ---*/
+    if (node_i == node_j) {
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          bij[iVar * nVar + jVar] += PassiveAssign(n[iVar]) * nbn * PassiveAssign(n[jVar]);
+        }
+      }
+    }
   }
 
-  /*--- Set known solution in rhs vector. ---*/
-  b(node_i, iVar) = x_i;
+  OtherType proj{};
+  for (auto iVar = 0ul; iVar < nVar; ++iVar) proj += b(node_i, iVar) * n[iVar];
+  for (auto iVar = 0ul; iVar < nVar; ++iVar) b(node_i, iVar) -= proj * n[iVar];
 }
 
 template <class ScalarType>
@@ -1203,8 +1226,7 @@ void CSysMatrix<ScalarType>::ComputePastixPreconditioner(const CSysVector<Scalar
 #define INSTANTIATE_MATRIX(TYPE)                                                                                  \
   template class CSysMatrix<TYPE>;                                                                                \
   template void CSysMatrix<TYPE>::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector<su2double>&); \
-  template void CSysMatrix<TYPE>::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double,                   \
-                                                       CSysVector<su2double>&);                                   \
+  template void CSysMatrix<TYPE>::EnforceZeroProjection(unsigned long, const su2double*, CSysVector<su2double>&); \
   INSTANTIATE_COMMS(TYPE)
 
 #ifdef CODI_FORWARD_TYPE

SU2_CFD/include/output/COutput.hpp

@@ -241,15 +241,16 @@ class COutput {
 
   /*----------------------------- Volume output ----------------------------*/
 
-  CParallelDataSorter* volumeDataSorter;    //!< Volume data sorter
-  CParallelDataSorter* surfaceDataSorter;   //!< Surface data sorter
+  CParallelDataSorter* volumeDataSorter;        //!< Volume data sorter.
+  CParallelDataSorter* volumeDataSorterCompact; //!< Volume data sorter for compact files.
+  CParallelDataSorter* surfaceDataSorter;       //!< Surface data sorter.
 
-  vector<string> volumeFieldNames;     //!< Vector containing the volume field names
-  unsigned short nVolumeFields;        //!< Number of fields in the volume output
+  vector<string> volumeFieldNames;          //!< Vector containing the volume field names.
+  vector<string> requiredVolumeFieldNames;  //!< Vector containing the minimum required volume field names.
 
-  string volumeFilename,               //!< Volume output filename
-  surfaceFilename,                     //!< Surface output filename
-  restartFilename;                     //!< Restart output filename
+  string volumeFilename,               //!< Volume output filename.
+  surfaceFilename,                     //!< Surface output filename.
+  restartFilename;                     //!< Restart output filename.
 
   /** \brief Structure to store information for a volume output field.
    *
@@ -259,40 +260,48 @@ class COutput {
     /*! \brief The name of the field, i.e. the name that is printed in the file header.*/
     string fieldName;
     /*! \brief This value identifies the position of the values of this field at each node in the ::Local_Data array. */
-    short  offset;
+    short offset = -1;
+    /*! \brief This offset is used for the compact formulation. */
+    short offsetCompact = -1;
     /*! \brief The group this field belongs to. */
     string outputGroup;
-    /*! \brief String containing the description of the field */
+    /*! \brief String containing the description of the field. */
     string description;
     /*! \brief Default constructor. */
-    VolumeOutputField () {}
+    VolumeOutputField() = default;
     /*! \brief Constructor to initialize all members. */
-    VolumeOutputField(string fieldName_, int offset_, string volumeOutputGroup_, string description_):
-      fieldName(std::move(fieldName_)), offset(std::move(offset_)),
-      outputGroup(std::move(volumeOutputGroup_)), description(std::move(description_)){}
+    VolumeOutputField(string fieldName_, string volumeOutputGroup_, string description_):
+      fieldName(std::move(fieldName_)),
+      outputGroup(std::move(volumeOutputGroup_)),
+      description(std::move(description_)) {}
   };
 
   /*! \brief Associative map to access data stored in the volume output fields by a string identifier. */
-  std::map<string, VolumeOutputField >          volumeOutput_Map;
+  std::map<string, VolumeOutputField > volumeOutput_Map;
   /*! \brief Vector that contains the keys of the ::volumeOutput_Map in the order of their insertion. */
-  std::vector<string>                           volumeOutput_List;
-
-  /*! \brief Vector to cache the positions of the field in the data array */
-  std::vector<short>                            fieldIndexCache;
-  /*! \brief Current value of the cache index */
-  unsigned short                                cachePosition;
-  /*! \brief Boolean to store whether the field index cache should be build. */
-  bool                                          buildFieldIndexCache;
-  /*! \brief Vector to cache the positions of the field in the data array */
-  std::vector<short>                            fieldGetIndexCache;
-  /*! \brief Current value of the cache index */
-  unsigned short                                curGetFieldIndex;
+  std::vector<string> volumeOutput_List;
+
+  /*! \brief Whether the field index caches should be build. */
+  bool buildFieldIndexCache;
+
+  /*! \brief Vectors to cache the positions of the fields in the data array. */
+  std::vector<short> fieldIndexCache, fieldIndexCacheCompact;
+  /*! \brief Current value of the cache indices. */
+  unsigned short cachePosition;
+
+  /*! \brief Vector to cache the positions of the field in the data array. */
+  std::vector<short> fieldGetIndexCache;
+  /*! \brief Current value of the cache index. */
+  unsigned short curGetFieldIndex;
 
   /*! \brief Requested volume field names in the config file. */
   std::vector<string> requestedVolumeFields;
   /*! \brief Number of requested volume field names in the config file. */
   unsigned short nRequestedVolumeFields;
 
+  /*! \brief Minimum required volume fields for restart file. */
+  const std::vector<string> restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY"};
+
   /*----------------------------- Convergence monitoring ----------------------------*/
 
   su2double cauchyValue,         /*!< \brief Summed value of the convergence indicator. */
@@ -736,8 +745,9 @@ class COutput {
    * \param[in] groupname - The name of the group this field belongs to.
    * \param[in] description - Description of the volume field.
    */
-  inline void AddVolumeOutput(string name, string field_name, string groupname, string description){
-    volumeOutput_Map[name] = VolumeOutputField(field_name, -1, groupname, description);
+  inline void AddVolumeOutput(const string& name, const string& field_name,
+                              const string& group_name, const string& description) {
+    volumeOutput_Map[name] = VolumeOutputField(field_name, group_name, description);
     volumeOutput_List.push_back(name);
   }
 
@@ -959,14 +969,14 @@ class COutput {
 
   /*!
    * \brief Sets the turboperformance screen output
-   * \param[in] TurboPerf - Turboperformance class 
+   * \param[in] TurboPerf - Turboperformance class
    * \param[in] config - Definition of the particular problem
    * \param[in] TimeIter - Index of the current time-step
    * \param[in] OuterIter - Index of current outer iteration
    * \param[in] InnerIter - Index of current inner iteration
    */
   inline virtual void SetTurboPerformance_Output(std::shared_ptr<CTurboOutput> TurboPerf, CConfig *config, unsigned long TimeIter, unsigned long OuterIter, unsigned long InnerIter) {}
-  
+
   /*!
    * \brief Sets the multizone turboperformacne screen output
    * \param[in] TurboStagePerf - Stage turboperformance class
@@ -982,7 +992,7 @@ class COutput {
    * \param[in] config - Definition of the particular problem
    */
   inline virtual void LoadTurboHistoryData(std::shared_ptr<CTurbomachineryStagePerformance> TurboStagePerf, std::shared_ptr<CTurboOutput> TurboPerf, CConfig *config) {}
-                                  
+
   /*!
    * \brief Write the kinematic and thermodynamic variables at each spanwise division
    * \param[in] solver - The container hold all solution data

SU2_CFD/include/output/filewriter/CParallelDataSorter.hpp

@@ -108,7 +108,8 @@ class CParallelDataSorter{
   int nSends,                          //!< Number of sends
   nRecvs;                              //!< Number of receives
 
-  vector<string> fieldNames;           //!< Vector with names of the output fields
+  vector<string> fieldNames;           //!< Vector with names of all the output fields
+  vector<string> requiredFieldNames;   //!< Vector with names of the required output fields that we write to file
 
   unsigned short nDim;                 //!< Spatial dimension of the data
 
@@ -340,6 +341,22 @@ class CParallelDataSorter{
     return fieldNames;
   }
 
+ /*!
+   * \brief Get the vector containing the names of the required output fields
+   * \return Vector of strings containing the required field names
+   */
+  const vector<string>& GetRequiredFieldNames() const{
+    return requiredFieldNames;
+  }
+
+  /*!
+   * \brief Set the vector of required output fields.
+   * \return None.
+   */
+  void SetRequiredFieldNames(const vector<string>& req_field_names) {
+    requiredFieldNames = req_field_names;
+  }
+
   /*!
    * \brief Get the spatial dimension
    * \return The spatial dimension

SU2_CFD/src/output/CAdjFlowCompOutput.cpp

@@ -61,7 +61,7 @@ CAdjFlowCompOutput::CAdjFlowCompOutput(CConfig *config, unsigned short nDim) : C
 
   if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), string("SENSITIVITY")) == requestedVolumeFields.end()) {
     requestedVolumeFields.emplace_back("SENSITIVITY");
-    nRequestedVolumeFields ++;
+    nRequestedVolumeFields++;
   }
 
   stringstream ss;

SU2_CFD/src/output/CAdjFlowIncOutput.cpp

@@ -67,7 +67,7 @@ CAdjFlowIncOutput::CAdjFlowIncOutput(CConfig *config, unsigned short nDim) : CAd
 
   if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), string("SENSITIVITY")) == requestedVolumeFields.end()) {
     requestedVolumeFields.emplace_back("SENSITIVITY");
-    nRequestedVolumeFields ++;
+    nRequestedVolumeFields++;
   }
 
   stringstream ss;

SU2_CFD/src/output/CAdjHeatOutput.cpp

@@ -58,7 +58,7 @@ CAdjHeatOutput::CAdjHeatOutput(CConfig *config, unsigned short nDim) : COutput(c
 
   if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), string("SENSITIVITY")) == requestedVolumeFields.end()) {
     requestedVolumeFields.emplace_back("SENSITIVITY");
-    nRequestedVolumeFields ++;
+    nRequestedVolumeFields++;
   }
 
   stringstream ss;

SU2_CFD/src/output/CElasticityOutput.cpp

@@ -235,21 +235,20 @@ void CElasticityOutput::SetVolumeOutputFields(CConfig *config){
   // Grid coordinates
   AddVolumeOutput("COORD-X", "x", "COORDINATES", "x-component of the coordinate vector");
   AddVolumeOutput("COORD-Y", "y", "COORDINATES", "y-component of the coordinate vector");
-  if (nDim == 3)
-    AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
+  if (nDim == 3) AddVolumeOutput("COORD-Z", "z", "COORDINATES", "z-component of the coordinate vector");
 
   AddVolumeOutput("DISPLACEMENT-X",    "Displacement_x", "SOLUTION", "x-component of the displacement vector");
   AddVolumeOutput("DISPLACEMENT-Y",    "Displacement_y", "SOLUTION", "y-component of the displacement vector");
   if (nDim == 3) AddVolumeOutput("DISPLACEMENT-Z", "Displacement_z", "SOLUTION", "z-component of the displacement vector");
 
-  if(dynamic){
-    AddVolumeOutput("VELOCITY-X",    "Velocity_x", "VELOCITY", "x-component of the velocity vector");
-    AddVolumeOutput("VELOCITY-Y",    "Velocity_y", "VELOCITY", "y-component of the velocity vector");
-    if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "VELOCITY", "z-component of the velocity vector");
+  if (dynamic) {
+    AddVolumeOutput("VELOCITY-X",    "Velocity_x", "SOLUTION", "x-component of the velocity vector");
+    AddVolumeOutput("VELOCITY-Y",    "Velocity_y", "SOLUTION", "y-component of the velocity vector");
+    if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector");
 
-    AddVolumeOutput("ACCELERATION-X",    "Acceleration_x", "ACCELERATION", "x-component of the acceleration vector");
-    AddVolumeOutput("ACCELERATION-Y",    "Acceleration_y", "ACCELERATION", "y-component of the acceleration vector");
-    if (nDim == 3) AddVolumeOutput("ACCELERATION-Z", "Acceleration_z", "ACCELERATION", "z-component of the acceleration vector");
+    AddVolumeOutput("ACCELERATION-X",    "Acceleration_x", "SOLUTION", "x-component of the acceleration vector");
+    AddVolumeOutput("ACCELERATION-Y",    "Acceleration_y", "SOLUTION", "y-component of the acceleration vector");
+    if (nDim == 3) AddVolumeOutput("ACCELERATION-Z", "Acceleration_z", "SOLUTION", "z-component of the acceleration vector");
   }
 
   AddVolumeOutput("STRESS-XX",    "Sxx", "STRESS", "x-component of the normal stress vector");