Merge branch 'develop' into fix/lm-transition-model-adjoint-2606

Author: guptapratykshh

Common/include/CConfig.hpp

@@ -651,6 +651,7 @@ class CConfig {
   unsigned long Linear_Solver_Prec_Threads;      /*!< \brief Number of threads per rank for ILU and LU_SGS preconditioners. */
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
+  su2double MSW_Alpha;                  /*!< \brief Coefficient for blending states in the MSW scheme. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
   su2double Relaxation_Factor_Adjoint;  /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
   su2double Relaxation_Factor_CHT;      /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
@@ -4432,6 +4433,11 @@ class CConfig {
    */
   void SetNewtonKrylovRelaxation(const su2double& relaxation) { NK_Relaxation = relaxation; }
 
+  /*!
+   * \brief Returns the MSW alpha (coefficient of the state blending weight).
+   */
+  su2double GetMSW_Alpha(void) const { return MSW_Alpha; }
+
   /*!
    * \brief Returns the Roe kappa (multipler of the dissipation term).
    */

Common/src/CConfig.cpp

@@ -1916,6 +1916,8 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
   /* DESCRIPTION: Relaxation of the CHT coupling */
   addDoubleOption("RELAXATION_FACTOR_CHT", Relaxation_Factor_CHT, 1.0);
+  /* DESCRIPTION: MSW alpha coefficient */
+  addDoubleOption("MSW_ALPHA", MSW_Alpha, 5.0);
   /* DESCRIPTION: Roe coefficient */
   addDoubleOption("ROE_KAPPA", Roe_Kappa, 0.5);
   /* DESCRIPTION: Roe-Turkel preconditioning for low Mach number flows */

SU2_CFD/include/numerics/CNumerics.hpp

@@ -35,6 +35,7 @@
 
 #include "../../../Common/include/CConfig.hpp"
 #include "../../../Common/include/linear_algebra/blas_structure.hpp"
+#include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
 class CElement;
 class CFluidModel;
@@ -46,7 +47,7 @@ class CFluidModel;
  */
 class CNumerics {
 protected:
-  enum : size_t {MAXNDIM = 3}; /*!< \brief Max number of space dimensions, used in some static arrays. */
+  static constexpr size_t MAXNDIM = 3; /*!< \brief Max number of space dimensions, used in some static arrays. */
 
   unsigned short nDim, nVar;  /*!< \brief Number of dimensions and variables. */
   su2double Gamma;            /*!< \brief Fluid's Gamma constant (ratio of specific heats). */
@@ -1146,15 +1147,77 @@ class CNumerics {
   /*!
    * \brief Computation of the matrix P, this matrix diagonalize the conservative Jacobians in
    *        the form $P^{-1}(A.Normal)P=Lambda$.
-   * \param[in] val_density - Value of the density.
-   * \param[in] val_velocity - Value of the velocity.
-   * \param[in] val_soundspeed - Value of the sound speed.
-   * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
-   * \param[out] val_p_tensor - Pointer to the P matrix.
-   */
-  void GetPMatrix(const su2double *val_density, const su2double *val_velocity,
-                  const su2double *val_soundspeed, const su2double *val_normal,
-                  su2double **val_p_tensor) const;
+   * \param[in] density - Value of the density.
+   * \param[in] velocity - Velocity vector.
+   * \param[in] soundspeed - Value of the sound speed.
+   * \param[in] normal - Normal vector, the norm of the vector is the area of the face.
+   * \param[out] p_tensor - P matrix.
+   */
+  template <typename Matrix>
+  void GetPMatrix(const su2double& density, const su2double* velocity,
+                  const su2double& soundspeed, const su2double* normal,
+                  Matrix& p_tensor) const {
+    const su2double rhooc = density / soundspeed;
+    const su2double rhoxc = density * soundspeed;
+
+    if (nDim == 2) {
+      const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity);
+      const su2double projvel = GeometryToolbox::DotProduct(2, velocity, normal);
+
+      p_tensor[0][0] = 1.0;
+      p_tensor[0][1] = 0.0;
+      p_tensor[0][2] = 0.5 * rhooc;
+      p_tensor[0][3] = 0.5 * rhooc;
+
+      p_tensor[1][0] = velocity[0];
+      p_tensor[1][1] = density * normal[1];
+      p_tensor[1][2] = 0.5 * (velocity[0] * rhooc + normal[0] * density);
+      p_tensor[1][3] = 0.5 * (velocity[0] * rhooc - normal[0] * density);
+
+      p_tensor[2][0] = velocity[1];
+      p_tensor[2][1] = -density * normal[0];
+      p_tensor[2][2] = 0.5 * (velocity[1] * rhooc + normal[1] * density);
+      p_tensor[2][3] = 0.5 * (velocity[1] * rhooc - normal[1] * density);
+
+      p_tensor[3][0] = ke;
+      p_tensor[3][1] = density * (velocity[0] * normal[1] - velocity[1] * normal[0]);
+      p_tensor[3][2] = 0.5 * (ke * rhooc + density * projvel + rhoxc / Gamma_Minus_One);
+      p_tensor[3][3] = 0.5 * (ke * rhooc - density * projvel + rhoxc / Gamma_Minus_One);
+    } else {
+      const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity);
+      const su2double projvel = GeometryToolbox::DotProduct(3, velocity, normal);
+
+      p_tensor[0][0] = normal[0];
+      p_tensor[0][1] = normal[1];
+      p_tensor[0][2] = normal[2];
+      p_tensor[0][3] = 0.5 * rhooc;
+      p_tensor[0][4] = 0.5 * rhooc;
+
+      p_tensor[1][0] = velocity[0] * normal[0];
+      p_tensor[1][1] = velocity[0] * normal[1] - density * normal[2];
+      p_tensor[1][2] = velocity[0] * normal[2] + density * normal[1];
+      p_tensor[1][3] = 0.5 * (velocity[0] * rhooc + density * normal[0]);
+      p_tensor[1][4] = 0.5 * (velocity[0] * rhooc - density * normal[0]);
+
+      p_tensor[2][0] = velocity[1] * normal[0] + density * normal[2];
+      p_tensor[2][1] = velocity[1] * normal[1];
+      p_tensor[2][2] = velocity[1] * normal[2] - density * normal[0];
+      p_tensor[2][3] = 0.5 * (velocity[1] * rhooc + density * normal[1]);
+      p_tensor[2][4] = 0.5 * (velocity[1] * rhooc - density * normal[1]);
+
+      p_tensor[3][0] = velocity[2] * normal[0] - density * normal[1];
+      p_tensor[3][1] = velocity[2] * normal[1] + density * normal[0];
+      p_tensor[3][2] = velocity[2] * normal[2];
+      p_tensor[3][3] = 0.5 * (velocity[2] * rhooc + density * normal[2]);
+      p_tensor[3][4] = 0.5 * (velocity[2] * rhooc - density * normal[2]);
+
+      p_tensor[4][0] = ke * normal[0] + density * (velocity[1] * normal[2] - velocity[2] * normal[1]);
+      p_tensor[4][1] = ke * normal[1] + density * (velocity[2] * normal[0] - velocity[0] * normal[2]);
+      p_tensor[4][2] = ke * normal[2] + density * (velocity[0] * normal[1] - velocity[1] * normal[0]);
+      p_tensor[4][3] = 0.5 * (ke * rhooc + density * projvel + rhoxc / Gamma_Minus_One);
+      p_tensor[4][4] = 0.5 * (ke * rhooc - density * projvel + rhoxc / Gamma_Minus_One);
+    }
+  }
 
   /*!
    * \brief Computation of the matrix Rinv*Pe.
@@ -1261,15 +1324,83 @@ class CNumerics {
   /*!
    * \brief Computation of the matrix P^{-1}, this matrix diagonalize the conservative Jacobians
    *        in the form $P^{-1}(A.Normal)P=Lambda$.
-   * \param[in] val_density - Value of the density.
-   * \param[in] val_velocity - Value of the velocity.
-   * \param[in] val_soundspeed - Value of the sound speed.
-   * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
-   * \param[out] val_invp_tensor - Pointer to inverse of the P matrix.
-   */
-  void GetPMatrix_inv(const su2double *val_density, const su2double *val_velocity,
-                      const su2double *val_soundspeed, const su2double *val_normal,
-                      su2double **val_invp_tensor) const;
+   * \param[in] density - Value of the density.
+   * \param[in] velocity - Velocity vector.
+   * \param[in] soundspeed - Value of the sound speed.
+   * \param[in] normal - Normal vector, the norm of the vector is the area of the face.
+   * \param[out] inv_p_tensor - Pointer to inverse of the P matrix.
+   */
+  template <typename Matrix>
+  void GetPMatrix_inv(const su2double& density, const su2double* velocity,
+                      const su2double& soundspeed, const su2double* normal,
+                      Matrix& inv_p_tensor) const {
+    const su2double rhoxc = density * soundspeed;
+    const su2double c2 = pow(soundspeed, 2);
+    const su2double gm1 = Gamma_Minus_One;
+    const su2double k0orho = normal[0] / density;
+    const su2double k1orho = normal[1] / density;
+    const su2double gm1_o_c2 = gm1 / c2;
+    const su2double gm1_o_rhoxc = gm1 / rhoxc;
+
+    if (nDim == 3) {
+      const su2double k2orho = normal[2] / density;
+      const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity);
+      const su2double projvel_o_rho = GeometryToolbox::DotProduct(3, velocity, normal) / density;
+
+      inv_p_tensor[0][0] = normal[0] + k1orho * velocity[2] - k2orho * velocity[1] - normal[0] * gm1_o_c2 * ke;
+      inv_p_tensor[0][1] = normal[0] * gm1_o_c2 * velocity[0];
+      inv_p_tensor[0][2] = k2orho + normal[0] * gm1_o_c2 * velocity[1];
+      inv_p_tensor[0][3] = -k1orho + normal[0] * gm1_o_c2 * velocity[2];
+      inv_p_tensor[0][4] = -normal[0] * gm1_o_c2;
+
+      inv_p_tensor[1][0] = normal[1] + k2orho * velocity[0] - k0orho * velocity[2] - normal[1] * gm1_o_c2 * ke;
+      inv_p_tensor[1][1] = -k2orho + normal[1] * gm1_o_c2 * velocity[0];
+      inv_p_tensor[1][2] = normal[1] * gm1_o_c2 * velocity[1];
+      inv_p_tensor[1][3] = k0orho + normal[1] * gm1_o_c2 * velocity[2];
+      inv_p_tensor[1][4] = -normal[1] * gm1_o_c2;
+
+      inv_p_tensor[2][0] = normal[2] + k0orho * velocity[1] - k1orho * velocity[0] - normal[2] * gm1_o_c2 * ke;
+      inv_p_tensor[2][1] = k1orho + normal[2] * gm1_o_c2 * velocity[0];
+      inv_p_tensor[2][2] = -k0orho + normal[2] * gm1_o_c2 * velocity[1];
+      inv_p_tensor[2][3] = normal[2] * gm1_o_c2 * velocity[2];
+      inv_p_tensor[2][4] = -normal[2] * gm1_o_c2;
+
+      inv_p_tensor[3][0] = -projvel_o_rho + gm1_o_rhoxc * ke;
+      inv_p_tensor[3][1] = k0orho - gm1_o_rhoxc * velocity[0];
+      inv_p_tensor[3][2] = k1orho - gm1_o_rhoxc * velocity[1];
+      inv_p_tensor[3][3] = k2orho - gm1_o_rhoxc * velocity[2];
+      inv_p_tensor[3][4] = gm1_o_rhoxc;
+
+      inv_p_tensor[4][0] = projvel_o_rho + gm1_o_rhoxc * ke;
+      inv_p_tensor[4][1] = -k0orho - gm1_o_rhoxc * velocity[0];
+      inv_p_tensor[4][2] = -k1orho - gm1_o_rhoxc * velocity[1];
+      inv_p_tensor[4][3] = -k2orho - gm1_o_rhoxc * velocity[2];
+      inv_p_tensor[4][4] = gm1_o_rhoxc;
+    } else {
+      const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity);
+      const su2double projvel_o_rho = GeometryToolbox::DotProduct(2, velocity, normal) / density;
+
+      inv_p_tensor[0][0] = 1 - gm1_o_c2 * ke;
+      inv_p_tensor[0][1] = gm1_o_c2 * velocity[0];
+      inv_p_tensor[0][2] = gm1_o_c2 * velocity[1];
+      inv_p_tensor[0][3] = -gm1_o_c2;
+
+      inv_p_tensor[1][0] = -k1orho * velocity[0] + k0orho * velocity[1];
+      inv_p_tensor[1][1] = k1orho;
+      inv_p_tensor[1][2] = -k0orho;
+      inv_p_tensor[1][3] = 0;
+
+      inv_p_tensor[2][0] = -projvel_o_rho + gm1_o_rhoxc * ke;
+      inv_p_tensor[2][1] = k0orho - gm1_o_rhoxc * velocity[0];
+      inv_p_tensor[2][2] = k1orho - gm1_o_rhoxc * velocity[1];
+      inv_p_tensor[2][3] = gm1_o_rhoxc;
+
+      inv_p_tensor[3][0] = projvel_o_rho + gm1_o_rhoxc * ke;
+      inv_p_tensor[3][1] = -k0orho - gm1_o_rhoxc * velocity[0];
+      inv_p_tensor[3][2] = -k1orho - gm1_o_rhoxc * velocity[1];
+      inv_p_tensor[3][3] = gm1_o_rhoxc;
+    }
+  }
 
   /*!
    * \brief Compute viscous residual and jacobian.

SU2_CFD/include/numerics/NEMO/convection/msw.hpp

@@ -38,6 +38,8 @@
  */
 class CUpwMSW_NEMO : public CNEMONumerics {
 private:
+    const su2double alpha;
+
     su2double *Diff_U;
     su2double *ust_i, *ust_j;
     su2double *Fc_i, *Fc_j;

SU2_CFD/include/numerics/flow/convection/fvs.hpp

@@ -38,15 +38,16 @@
  */
 class CUpwMSW_Flow final : public CNumerics {
 private:
-  bool implicit;
-  su2double *u_i, *u_j, *ust_i, *ust_j;
-  su2double *Fc_i, *Fc_j;
-  su2double *Lambda_i, *Lambda_j;
-  su2double *Vst_i, *Vst_j, *Velst_i, *Velst_j;
-  su2double **P_Tensor, **invP_Tensor;
+  static constexpr auto MAXNVAR = MAXNDIM + 2;
 
-  su2double** Jacobian_i; /*!< \brief The Jacobian w.r.t. point i after computation. */
-  su2double** Jacobian_j; /*!< \brief The Jacobian w.r.t. point j after computation. */
+  const su2double alpha;
+  const bool dynamic_grid;
+
+  su2double buf_Jacobian_i[MAXNVAR * MAXNVAR], buf_Jacobian_j[MAXNVAR * MAXNVAR];
+  su2double* Jacobian_i[MAXNVAR]; /*!< \brief The Jacobian w.r.t. point i after computation. */
+  su2double* Jacobian_j[MAXNVAR]; /*!< \brief The Jacobian w.r.t. point j after computation. */
+
+  su2double Fc[MAXNVAR];
 
 public:
   /*!
@@ -57,11 +58,6 @@ class CUpwMSW_Flow final : public CNumerics {
    */
   CUpwMSW_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
 
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CUpwMSW_Flow(void) override;
-
   /*!
    * \brief Compute the Roe's flux between two nodes i and j.
    * \param[in] config - Definition of the particular problem.

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -340,8 +340,11 @@ void CFluidIteration::UpdateRamp(CGeometry**** geometry_container, CConfig** con
     const long unsigned updateFreq = RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::UPDATE_FREQ];
     const long unsigned rampLength = RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::FINAL_ITER];
     auto iterFrac = (static_cast<double>(iter - startIter)/static_cast<double>((rampLength + startIter) - startIter));
-    if (iter < startIter) return;
-    if ((iter == startIter) && (rank == MASTER_NODE)) cout << "Beginning to ramp MUSCL scheme..." << endl;
+    if (iter < startIter) {
+      config->SetMUSCLRampValue(0);
+      return;
+    }
+    if (iter == startIter && rank == MASTER_NODE) cout << "Beginning to ramp MUSCL scheme..." << endl;
     if ((iter % updateFreq == 0 && iter < (rampLength + startIter)) || (iter == (rampLength + startIter))) {
       switch (RampMUSCLParam.Kind_MUSCLRamp) {
         case MUSCL_RAMP_TYPE::ITERATION: