Merge branch 'develop' into fix_multigrid_agglomeration

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.*/
@@ -1004,7 +1005,7 @@ class CConfig {
   bool ExtraOutput;           /*!< \brief Check if extra output need. */
   bool Wall_Functions;           /*!< \brief Use wall functions with the turbulence model */
   long ExtraHeatOutputZone;      /*!< \brief Heat solver zone with extra screen output */
-  bool DeadLoad;                 /*!< \brief Application of dead loads to the FE analysis */
+  bool CentrifugalForce;         /*!< \brief Application of centrifugal forces to the FE analysis */
   bool PseudoStatic;             /*!< \brief Application of dead loads to the FE analysis */
   bool SteadyRestart;            /*!< \brief Restart from a steady state for FSI problems. */
   su2double Newmark_beta,        /*!< \brief Parameter alpha for Newmark method. */
@@ -1013,10 +1014,13 @@ class CConfig {
   su2double *Int_Coeffs;         /*!< \brief Time integration coefficients for structural method. */
   unsigned short nElasticityMod, /*!< \brief Number of different values for the elasticity modulus. */
   nPoissonRatio,                    /*!< \brief Number of different values for the Poisson ratio modulus. */
-  nMaterialDensity;                 /*!< \brief Number of different values for the Material density. */
+  nMaterialDensity,                 /*!< \brief Number of different values for the Material density. */
+  nMaterialThermalExpansion;        /*!< \brief Number of different values for thermal expansion coefficient. */
   su2double *ElasticityMod,         /*!< \brief Value of the elasticity moduli. */
   *PoissonRatio,                    /*!< \brief Value of the Poisson ratios. */
-  *MaterialDensity;                 /*!< \brief Value of the Material densities. */
+  *MaterialDensity,                 /*!< \brief Value of the Material densities. */
+  *MaterialThermalExpansion,        /*!< \brief Value of the thermal expansion coefficients. */
+  MaterialReferenceTemperature;     /*!< \brief Value of the reference temperature for thermal expansion. */
   unsigned short nElectric_Field,   /*!< \brief Number of different values for the electric field in the membrane. */
   nDim_Electric_Field;              /*!< \brief Dimensionality of the problem. */
   unsigned short nDim_RefNode;      /*!< \brief Dimensionality of the vector . */
@@ -2389,6 +2393,16 @@ class CConfig {
    */
   su2double GetMaterialDensity(unsigned short id_val) const { return MaterialDensity[id_val]; }
 
+  /*!
+   * \brief Get the thermal expansion coefficient.
+   */
+  su2double GetMaterialThermalExpansion(unsigned short id_val) const { return MaterialThermalExpansion[id_val]; }
+
+  /*!
+   * \brief Temperature at which there is no stress from thermal expansion.
+   */
+  su2double GetMaterialReferenceTemperature() const { return MaterialReferenceTemperature; }
+
   /*!
    * \brief Compressibility/incompressibility of the solids analysed using the structural solver.
    * \return Compressible or incompressible.
@@ -5490,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
@@ -8929,10 +8949,9 @@ class CConfig {
   su2double GetAitkenDynMinInit(void) const { return AitkenDynMinInit; }
 
   /*!
-   * \brief Decide whether to apply dead loads to the model.
-   * \return <code>TRUE</code> if the dead loads are to be applied, <code>FALSE</code> otherwise.
+   * \brief Decide whether to apply centrifugal forces to the model.
    */
-  bool GetDeadLoad(void) const { return DeadLoad; }
+  bool GetCentrifugalForce(void) const { return CentrifugalForce; }
 
   /*!
    * \brief Identifies if the mesh is matching or not (temporary, while implementing interpolation procedures).

Common/include/geometry/elements/CElement.hpp

@@ -67,6 +67,8 @@ class CElement {
   su2activematrix NodalExtrap;  /*!< \brief Coordinates of the nodal points for Gaussian extrapolation. */
   su2activematrix NodalStress;  /*!< \brief Stress at the nodes. */
 
+  su2activevector NodalTemperature; /*!< \brief Temperature at the nodes. */
+
   /*--- Stiffness and load matrices. ---*/
   std::vector<su2activematrix> Kab; /*!< \brief Structure for the constitutive component of the tangent matrix. */
   su2activematrix Mab;              /*!< \brief Structure for the nodal components of the mass matrix. */
@@ -151,7 +153,7 @@ class CElement {
    * \param[in] iDim - Dimension.
    * \param[in] val_CoordRef - Value of the coordinate.
    */
-  inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordRef) {
+  inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordRef) {
     RefCoord(iNode, iDim) = val_CoordRef;
   }
 
@@ -161,10 +163,22 @@ class CElement {
    * \param[in] iDim - Dimension.
    * \param[in] val_CoordRef - Value of the coordinate.
    */
-  inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordCurr) {
+  inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordCurr) {
     CurrentCoord(iNode, iDim) = val_CoordCurr;
   }
 
+  /*!
+   * \brief Set the value of the temperature of a node.
+   */
+  inline void SetTemperature(unsigned short iNode, const su2double& val_Temperature) {
+    NodalTemperature[iNode] = val_Temperature;
+  }
+
+  /*!
+   * \brief Set the value of the temperature of all nodes.
+   */
+  inline void SetTemperature(const su2double& val_Temperature) { NodalTemperature = val_Temperature; }
+
   /*!
    * \brief Get the value of the coordinate of the nodes in the reference configuration.
    * \param[in] iNode - Index of node.
@@ -181,6 +195,17 @@ class CElement {
    */
   inline su2double GetCurr_Coord(unsigned short iNode, unsigned short iDim) const { return CurrentCoord(iNode, iDim); }
 
+  /*!
+   * \brief Get the value of the temperature at a Gaussian integration point.
+   */
+  inline su2double GetTemperature(unsigned short iGauss) const {
+    su2double Temperature = 0;
+    for (auto iNode = 0u; iNode < nNodes; ++iNode) {
+      Temperature += GetNi(iNode, iGauss) * NodalTemperature[iNode];
+    }
+    return Temperature;
+  }
+
   /*!
    * \brief Get the weight of the corresponding Gaussian Point.
    * \param[in] iGauss - index of the Gaussian point.
@@ -214,7 +239,9 @@ class CElement {
    * \param[in] nodeB - index of Node b.
    * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution.
    */
-  inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, su2double val_Mab) { Mab(nodeA, nodeB) += val_Mab; }
+  inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Mab) {
+    Mab(nodeA, nodeB) += val_Mab;
+  }
 
   /*!
    * \brief Add the value of a submatrix K relating nodes a and b, for the constitutive term.
@@ -243,7 +270,7 @@ class CElement {
    * \param[in] nodeB - index of Node b.
    * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution.
    */
-  inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, su2double val_Ks_ab) {
+  inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Ks_ab) {
     Ks_ab(nodeA, nodeB) += val_Ks_ab;
   }
 
@@ -354,7 +381,7 @@ class CElement {
    * \param[in] iVar - Variable index.
    * \param[in] val_Stress - Value of the stress added.
    */
-  inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, su2double val_Stress) {
+  inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, const su2double& val_Stress) {
     NodalStress(iNode, iVar) += val_Stress;
   }
 
@@ -789,7 +816,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> {
   /*!
    * \brief Shape functions (Ni) evaluated at point Xi,Eta.
    */
-  inline static void ShapeFunctions(su2double Xi, su2double Eta, su2double* Ni) {
+  inline static void ShapeFunctions(const su2double& Xi, const su2double& Eta, su2double* Ni) {
     Ni[0] = 0.25 * (1.0 - Xi) * (1.0 - Eta);
     Ni[1] = 0.25 * (1.0 + Xi) * (1.0 - Eta);
     Ni[2] = 0.25 * (1.0 + Xi) * (1.0 + Eta);
@@ -799,7 +826,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> {
   /*!
    * \brief Shape function Jacobian (dNi) evaluated at point Xi,Eta.
    */
-  inline static void ShapeFunctionJacobian(su2double Xi, su2double Eta, su2double dNi[][2]) {
+  inline static void ShapeFunctionJacobian(const su2double& Xi, const su2double& Eta, su2double dNi[][2]) {
     dNi[0][0] = -0.25 * (1.0 - Eta);
     dNi[0][1] = -0.25 * (1.0 - Xi);
     dNi[1][0] = 0.25 * (1.0 - Eta);

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

@@ -914,7 +914,10 @@ void CConfig::SetPointersNull() {
   Inlet_Velocity              = nullptr;   Inflow_Mach             = nullptr;   Inflow_Pressure       = nullptr;
   Outlet_Pressure             = nullptr;   Isothermal_Temperature  = nullptr;
 
-  ElasticityMod             = nullptr;     PoissonRatio                = nullptr;     MaterialDensity       = nullptr;
+  ElasticityMod = nullptr;
+  PoissonRatio = nullptr;
+  MaterialDensity = nullptr;
+  MaterialThermalExpansion = nullptr;
 
   Load_Dir = nullptr;            Load_Dir_Value = nullptr;          Load_Dir_Multiplier = nullptr;
   Disp_Dir = nullptr;            Disp_Dir_Value = nullptr;          Disp_Dir_Multiplier = nullptr;
@@ -1169,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 */
@@ -2440,6 +2445,10 @@ void CConfig::SetConfig_Options() {
   addDoubleListOption("POISSON_RATIO", nPoissonRatio, PoissonRatio);
   /* DESCRIPTION: Material density */
   addDoubleListOption("MATERIAL_DENSITY", nMaterialDensity, MaterialDensity);
+  /* DESCRIPTION: Material thermal expansion coefficient */
+  addDoubleListOption("MATERIAL_THERMAL_EXPANSION_COEFF", nMaterialThermalExpansion, MaterialThermalExpansion);
+  /* DESCRIPTION: Temperature at which there is no stress from thermal expansion */
+  addDoubleOption("MATERIAL_REFERENCE_TEMPERATURE", MaterialReferenceTemperature, 288.15);
   /* DESCRIPTION: Knowles B constant */
   addDoubleOption("KNOWLES_B", Knowles_B, 1.0);
   /* DESCRIPTION: Knowles N constant */
@@ -2500,11 +2509,11 @@ void CConfig::SetConfig_Options() {
   addEnumOption("NONLINEAR_FEM_SOLUTION_METHOD", Kind_SpaceIteScheme_FEA, Space_Ite_Map_FEA, STRUCT_SPACE_ITE::NEWTON);
   /* DESCRIPTION: Formulation for bidimensional elasticity solver */
   addEnumOption("FORMULATION_ELASTICITY_2D", Kind_2DElasForm, ElasForm_2D, STRUCT_2DFORM::PLANE_STRAIN);
-  /*  DESCRIPTION: Apply dead loads
-  *  Options: NO, YES \ingroup Config */
-  addBoolOption("DEAD_LOAD", DeadLoad, false);
+  /*  DESCRIPTION: Apply centrifugal forces
+   *  Options: NO, YES \ingroup Config */
+  addBoolOption("CENTRIFUGAL_FORCE", CentrifugalForce, false);
   /*  DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis)
-  *  Options: NO, YES \ingroup Config */
+   *  Options: NO, YES \ingroup Config */
   addBoolOption("PSEUDO_STATIC", PseudoStatic, false);
   /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */
   addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25);
@@ -3063,6 +3072,8 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){
             newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN instead.\n\n");
           else if (!option_name.compare("SPECIES_USE_STRONG_BC"))
             newString.append("SPECIES_USE_STRONG_BC is deprecated. Use MARKER_SPECIES_STRONG_BC= (marker1, ...) instead.\n\n");
+          else if (!option_name.compare("DEAD_LOAD"))
+            newString.append("DEAD_LOAD is deprecated. Use GRAVITY_FORCE or BODY_FORCE instead.\n\n");
           else {
             /*--- Find the most likely candidate for the unrecognized option, based on the length
              of start and end character sequences shared by candidates and the option. ---*/
@@ -4834,9 +4845,15 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854;
   }
 
-  if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity) {
-    SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, and MATERIAL_DENSITY need to have the same number "
-                   "of entries (the number of materials).", CURRENT_FUNCTION);
+  if (nMaterialThermalExpansion == 0) {
+    nMaterialThermalExpansion = 1;
+    MaterialThermalExpansion = new su2double[1]();
+  }
+
+  if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity ||
+      nElasticityMod != nMaterialThermalExpansion) {
+    SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, MATERIAL_DENSITY, and THERMAL_EXPANSION_COEFF need "
+                   "to have the same number of entries (the number of materials).", CURRENT_FUNCTION);
   }
 
   if (nElectric_Constant == 0) {

Common/src/geometry/elements/CElement.cpp

@@ -47,6 +47,8 @@ CElement::CElement(unsigned short ngauss, unsigned short nnodes, unsigned short
 
   NodalStress.resize(nNodes, 6) = su2double(0.0);
 
+  NodalTemperature.resize(nNodes) = su2double(0.0);
+
   Mab.resize(nNodes, nNodes);
   Ks_ab.resize(nNodes, nNodes);
   Kab.resize(nNodes);

Common/src/geometry/elements/CQUAD4.cpp

@@ -67,7 +67,6 @@ CQUAD4::CQUAD4() : CElementWithKnownSizes<NGAUSS, NNODE, NDIM>() {
   /*--- Store the extrapolation functions (used to compute nodal stresses) ---*/
 
   su2double ExtrapCoord[4][2], sqrt3 = 1.732050807568877;
-  ;
 
   ExtrapCoord[0][0] = -sqrt3;
   ExtrapCoord[0][1] = -sqrt3;

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