Merge branch 'develop' into fix_wallfunctions

Author: TobiKattmann

Common/include/linear_algebra/CSysMatrix.hpp

@@ -111,7 +111,7 @@ class CSysMatrix {
   const int rank;     /*!< \brief MPI Rank. */
   const int size;     /*!< \brief MPI Size. */
 
-  enum : size_t { MAXNVAR = 8 };    /*!< \brief Maximum number of variables the matrix can handle. The static
+  enum : size_t { MAXNVAR = 20 };   /*!< \brief Maximum number of variables the matrix can handle. The static
                                                 size is needed for fast, per-thread, static memory allocation. */
 
   enum { OMP_MAX_SIZE_L = 8192 };   /*!< \brief Max. chunk size used in light parallel for loops. */

Common/src/CConfig.cpp

@@ -6471,8 +6471,8 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
           break;
         case EULER_IMPLICIT:
           cout << "Euler implicit method for the flow equations." << endl;
-          if((Kind_Solver == NEMO_EULER) || (Kind_Solver == NEMO_NAVIER_STOKES))
-            SU2_MPI::Error("Implicit time scheme is not working with NEMO. Use EULER_EXPLICIT.", CURRENT_FUNCTION);
+          if (Kind_Solver == NEMO_NAVIER_STOKES)
+            SU2_MPI::Error("Implicit time scheme is not working with NEMO for viscous problems. Use EULER_EXPLICIT.", CURRENT_FUNCTION);
           switch (Kind_Linear_Solver) {
             case BCGSTAB:
             case FGMRES:

SU2_CFD/include/fluid/CMutationTCLib.hpp

@@ -77,7 +77,7 @@ class CMutationTCLib : public CNEMOGas {
   /*!
    * \brief Compute species V-E specific heats at constant volume.
    */
-  vector<su2double>& ComputeSpeciesCvVibEle() final;
+  vector<su2double>& ComputeSpeciesCvVibEle(su2double val_T) final;
 
   /*!
    * \brief Compute mixture energies (total internal energy and vibrational energy).
@@ -92,7 +92,9 @@ class CMutationTCLib : public CNEMOGas {
   /*!
    * \brief Compute species net production rates.
    */
-  vector<su2double>& ComputeNetProductionRates() final;
+  vector<su2double>& ComputeNetProductionRates(bool implicit, const su2double *V, const su2double* eve,
+                                               const su2double* cvve, const su2double* dTdU, const su2double* dTvedU,
+                                               su2double **val_jacobian) final;
 
   /*!
    * \brief Compute vibrational energy source term.

SU2_CFD/include/fluid/CNEMOGas.hpp

@@ -79,6 +79,7 @@ class CNEMOGas : public CFluidModel {
   hs,                                    /*!< \brief Species enthalpies */
   MolarFractions,                        /*!< \brief Species molar fractions */
   ws,                                    /*!< \brief Species net production rates */
+  taus,                                  /*!< \brief Relaxtion time scales */
   DiffusionCoeff,                        /*!< \brief Species diffusion coefficients*/
   Enthalpy_Formation,                    /*!< \brief Enthalpy of formation */
   Ref_Temperature;                       /*!< \brief Reference temperature for thermodynamic relations */
@@ -101,6 +102,7 @@ class CNEMOGas : public CFluidModel {
   /*!
    * \brief Set mixture thermodynamic state.
    * \param[in] P    - Pressure.
+   * \param[in] Ms   - Mass fractions of the gas.
    * \param[in] T    - Translational/Rotational temperature.
    * \param[in] Tve  - Vibrational/Electronic temperature.
    */
@@ -114,7 +116,7 @@ class CNEMOGas : public CFluidModel {
   /*!
    * \brief Compute species V-E specific heats at constant volume.
    */
-  virtual vector<su2double>& ComputeSpeciesCvVibEle() = 0;
+  virtual vector<su2double>& ComputeSpeciesCvVibEle(su2double val_T) = 0;
 
   /*!
    * \brief Compute mixture energies (total internal energy and vibrational energy).
@@ -124,13 +126,29 @@ class CNEMOGas : public CFluidModel {
   /*!
    * \brief Compute species net production rates.
    */
-  virtual vector<su2double>& ComputeNetProductionRates() = 0;
+  virtual vector<su2double>& ComputeNetProductionRates(bool implicit, const su2double *V, const su2double* eve,
+                                                       const su2double* cvve, const su2double* dTdU, const su2double* dTvedU,
+                                                       su2double **val_jacobian) = 0;
+
+  /*!
+   * \brief Populate chemical source term jacobian. 
+   */
+  virtual void ChemistryJacobian(unsigned short iReaction, const su2double *V, const su2double* eve,
+                                 const su2double* cvve, const su2double* dTdU, const su2double* dTvedU,
+                                 su2double **val_jacobian){};
 
   /*!
    * \brief Compute vibrational energy source term.
    */
   virtual su2double ComputeEveSourceTerm() { return 0; }
 
+  /*!
+   * \brief Compute vibration enery source term jacobian.
+   */
+  virtual void GetEveSourceTermJacobian(const su2double *V, const su2double *eve, const su2double *cvve,
+                                        const su2double *dTdU, const su2double* dTvedU,
+                                        su2double **val_jacobian){};
+
   /*!
    * \brief Compute vector of species V-E energy.
    */
@@ -184,17 +202,17 @@ class CNEMOGas : public CFluidModel {
   /*!
    * \brief Compute derivative of pressure w.r.t. conservative variables.
    */
-  void ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double *val_dPdU);
+  void ComputedPdU(const su2double *V, const vector<su2double>& val_eves, su2double *val_dPdU);
 
   /*!
    * \brief Compute derivative of temperature w.r.t. conservative variables.
    */
-  void ComputedTdU(su2double *V, su2double *val_dTdU);
+  void ComputedTdU(const su2double *V, su2double *val_dTdU);
 
   /*!
    * \brief Compute derivative of vibrational temperature w.r.t. conservative variables.
    */
-  void ComputedTvedU(su2double *V, vector<su2double>& val_eves, su2double *val_dTvedU);
+  void ComputedTvedU(const su2double *V, const vector<su2double>& val_eves, su2double *val_dTvedU);
 
   /*!
    * \brief Set the translational temperature.

SU2_CFD/include/fluid/CSU2TCLib.hpp

@@ -55,6 +55,7 @@ class CSU2TCLib : public CNEMOGas {
   Tcf_b,                          /*!< \brief Rate controlling temperature exponent (fwd) */
   Tcb_a,                          /*!< \brief Rate controlling temperature exponent (bkw) */
   Tcb_b,                          /*!< \brief Rate controlling temperature exponent (bkw) */
+  taus,                           /*!< \brief Relaxtion time scales */
   Diss,                           /*!< \brief Dissociation potential. */
   MassFrac_FreeStream,            /*!< \brief Mixture mass fractions of the fluid. */
   Wall_Catalycity,                /*!< \brief Specified wall species mass-fractions for catalytic boundaries. */
@@ -72,6 +73,25 @@ class CSU2TCLib : public CNEMOGas {
   C3DDoubleMatrix Omega00,       /*!< \brief Collision integrals (Omega(0,0)) */
   Omega11;                       /*!< \brief Collision integrals (Omega(1,1)) */
 
+  /*--- Implicit variables ---*/
+  su2double                     /*!< \brief Derivatives w.r.t. conservative variables */
+  *dPdU, *dTdU, *dTvedU;
+
+  su2double fwdRxn, bkwRxn,
+  kf,kfb,kb,
+  coeff, eta, epsilon, T_min,
+  Trxnf, Trxnb,
+  Thf, Thb, dThf, dThb,
+  theta, af, bf, ab, bb;
+
+  vector<su2double>
+  dkf, dkb,
+  dRfok, dRbok,
+  eve, eve_eq, cvve, cvve_eq;
+
+  vector<int>
+  alphak, betak;
+
 public:
 
   /*!
@@ -105,7 +125,7 @@ class CSU2TCLib : public CNEMOGas {
   /*!
    * \brief Compute species V-E specific heats at constant volume.
    */
-  vector<su2double>& ComputeSpeciesCvVibEle() final;
+  vector<su2double>& ComputeSpeciesCvVibEle(su2double val_T) final;
 
   /*!
    * \brief Compute mixture energies (total internal energy and vibrational energy).
@@ -120,13 +140,29 @@ class CSU2TCLib : public CNEMOGas {
   /*!
    * \brief Compute species net production rates.
    */
-  vector<su2double>& ComputeNetProductionRates() final;
+  vector<su2double>& ComputeNetProductionRates(bool implicit, const su2double *V, const su2double* eve,
+                                               const su2double* cvve, const su2double* dTdU, const su2double* dTvedU,
+                                               su2double **val_jacobian) final;
+
+  /*!
+   * \brief Compute chemical source term jacobian. 
+   */
+  void ChemistryJacobian(unsigned short iReaction, const su2double *V, const su2double* eve,
+                         const su2double* cvve, const su2double* dTdU, const su2double* dTvedU,
+                         su2double **val_jacobian) final;
 
   /*!
    * \brief Compute vibrational energy source term.
    */
   su2double ComputeEveSourceTerm() final;
 
+  /*!
+   * \brief Compute relaxation source term jacobian.
+   */
+  void GetEveSourceTermJacobian(const su2double *V, const su2double *eve, const su2double *cvve,
+                                const su2double *dTdU, const su2double* dTvedU,
+                                su2double **val_jacobian) final;
+
   /*!
    * \brief Compute species enthalpies.
    */

SU2_CFD/include/numerics/CNumerics.hpp

@@ -1400,6 +1400,31 @@ class CNumerics {
    * \param[out] val_Jacobian_i - Jacobian of the source terms
    */
   inline virtual ResidualType<> ComputeChemistry(const CConfig* config) { return ResidualType<>(nullptr,nullptr,nullptr); }
+
+  /*!
+   * \brief Check if residual constains a NaN value
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] val_residual - residual of the numeric function.
+   * \param[out] ERR - Presencse of NaN in vector
+   */
+  static bool CheckResidualNaNs(bool implicit, int nVar, const ResidualType<> residual) {
+
+    bool ERR = false;
+    const bool jac_j = residual.jacobian_j != nullptr;
+
+    for (auto iVar = 0; iVar<nVar; iVar++){
+      if (std::isnan(SU2_TYPE::GetValue(residual[iVar]))) ERR = true;
+
+      if (implicit) {
+        for (auto jVar = 0; jVar < nVar; jVar++){
+          if (std::isnan(SU2_TYPE::GetValue(residual.jacobian_i[iVar][jVar]))) ERR = true;
+          if ((jac_j) && (std::isnan(SU2_TYPE::GetValue(residual.jacobian_j[iVar][jVar])))) ERR = true;
+        }
+      }
+    }
+    return ERR;
+  }
+
   /*!
    * \brief Set intermittency for numerics (used in SA with LM transition model)
    */

SU2_CFD/include/numerics/NEMO/CNEMONumerics.hpp

@@ -44,19 +44,28 @@ class CNEMONumerics : public CNumerics {
   su2double *rhos_j, *u_j;
   su2double a_i, P_i, h_i;
   su2double a_j, P_j, h_j;
+  su2double rho_i, rho_j;
+  su2double e_ve_i, e_ve_j;
+  su2double rhoCvtr_i, rhoCvtr_j;
+  su2double rhoCvve_i, rhoCvve_j;
   unsigned short nPrimVar, nPrimVarGrad;
 
   su2double* Flux = nullptr;            /*!< \brief The flux / residual across the edge. */
+  su2double** Jacobian_i = nullptr;
+  su2double** Jacobian_j = nullptr;
 
   unsigned short nSpecies, nHeavy, nEl; /*!< \brief Number of species present in plasma */
 
+  /*--- Graidents w.r.t. conservative variables. ---*/
   su2double *dPdU_i, *dPdU_j;
   su2double *dTdU_i, *dTdU_j;
   su2double *dTvedU_i, *dTvedU_j;
   su2double Gamma_i, Gamma_j;
 
   vector<su2double> hs;
+  vector<su2double> Cvtr;
   su2double *eve_i, *eve_j, *Cvve_i, *Cvve_j;
+  su2double *Ys, *Ys_i, *Ys_j;
 
   unsigned short RHOS_INDEX, T_INDEX, TVE_INDEX, VEL_INDEX, P_INDEX,
   RHO_INDEX, H_INDEX, A_INDEX, RHOCVTR_INDEX, RHOCVVE_INDEX,

SU2_CFD/include/numerics/NEMO/NEMO_diffusion.hpp

@@ -54,6 +54,7 @@ class CAvgGrad_NEMO : public CNEMONumerics {
   Mean_Eddy_Viscosity,          /*!< \brief Mean value of the eddy viscosity. */
   Mean_Thermal_Conductivity,    /*!< \brief Mean value of the thermal conductivity. */
   Mean_Thermal_Conductivity_ve, /*!< \brief Mean value of the vib-el. thermal conductivity. */
+  *ProjFlux,                    /*!< \brief Projection of the viscous fluxes. */
   dist_ij;                      /*!< \brief Length of the edge and face. */
 
 public:
@@ -108,6 +109,7 @@ class CAvgGradCorrected_NEMO : public CNEMONumerics {
   Mean_Eddy_Viscosity,          /*!< \brief Mean value of the eddy viscosity. */
   Mean_Thermal_Conductivity,    /*!< \brief Mean value of the thermal conductivity. */
   Mean_Thermal_Conductivity_ve, /*!< \brief Mean value of the vib-el. thermal conductivity. */
+  *ProjFlux,                    /*!< \brief Projection of the viscous fluxes. */
   dist_ij;                      /*!< \brief Length of the edge and face. */
   bool implicit;                /*!< \brief Implicit calculus. */
 

SU2_CFD/include/numerics/NEMO/NEMO_sources.hpp

@@ -39,14 +39,10 @@
 class CSource_NEMO : public CNEMONumerics {
 private:
 
-  int    *alphak, *betak;
-  su2double *X; // Mole fraction
-  su2double *Y, **dYdr; // Mass fraction
-  su2double *dkf, *dkb, *dRfok, *dRbok;
-  vector<su2double> Cvvsst, ws;
-
-  su2double* residual = nullptr;        /*!< \brief The source residual. */
+  su2double *Y, **dYdr;                  // Mass fraction
 
+  su2double*  residual = nullptr;        /*!< \brief The source residual. */
+  su2double** jacobian = nullptr;
 public:
 
   /*!

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

@@ -37,20 +37,19 @@
  */
 class CUpwAUSM_NEMO : public CNEMONumerics {
 private:
-  su2double *FcL, *FcR;
-  //su2double *FcLR;
-  //su2double *dmLP, *dmRM, *dpLP, *dpRM;
-  //su2double *daL, *daR;
+  su2double *FcL, *FcR, *FcLR;
+  su2double *dmLP, *dmRM, *dpLP, *dpRM;
+  su2double *daL, *daR;
   su2double ProjVel_i, ProjVel_j;
-  su2double sq_vel;
-  //su2double  Proj_ModJac_Tensor_ij;
 
 public:
 
   /*!
    * \brief Constructor of the class.
    * \param[in] val_nDim - Number of dimensions of the problem.
    * \param[in] val_nVar - Number of variables of the problem.
+   * \param[in] val_nPrimVar - Number of primitive variables of the problem
+   * \param[in] val_nPrimVarGrad - Number of grad primitive variables of the problem
    * \param[in] config - Definition of the particular problem.
    */
   CUpwAUSM_NEMO(unsigned short val_nDim, unsigned short val_nVar, unsigned short val_nPrimVar, unsigned short val_nPrimVarGrad, CConfig *config);