Merge pull request #1139 from su2code/feature_NEMO_primitive_limiter

NEMO - transition to primitive limiter

Author: WallyMaier

SU2_CFD/include/numerics/NEMO/NEMO_diffusion.hpp

@@ -34,7 +34,7 @@
  * \brief Class for computing viscous term using the average of gradients.
  * \ingroup ViscDiscr
  * \author S.R. Copeland, W. Maier, C. Garbacz
- * \version 7.0.6 "falcon"
+ * \version 7.0.8 "Blackbird"
  */
 class CAvgGrad_NEMO : public CNEMONumerics {
 private:

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

@@ -34,15 +34,17 @@
  * \brief Class for solving an approximate Riemann AUSM.
  * \ingroup ConvDiscr
  * \author F. Palacios, S.R. Copeland, W. Maier, C. Garbacz
- * \version 7.0.6 'Falcon'
+ * \version 7.0.8 'Blackbird'
  */
 class CUpwAUSM_NEMO : public CNEMONumerics {
 private:
-  su2double *FcL, *FcR, *FcLR;
-  su2double *dmLP, *dmRM, *dpLP, *dpRM;
-  su2double *daL, *daR;
+  su2double *FcL, *FcR;
+  //su2double *FcLR;
+  //su2double *dmLP, *dmRM, *dpLP, *dpRM;
+  //su2double *daL, *daR;
   su2double ProjVel_i, ProjVel_j;
-  su2double sq_vel, Proj_ModJac_Tensor_ij;
+  su2double sq_vel;
+  //su2double  Proj_ModJac_Tensor_ij;
 
 public:
 
@@ -64,4 +66,4 @@ class CUpwAUSM_NEMO : public CNEMONumerics {
    * \param[in] config - Definition of the particular problem.
    */
   ResidualType<> ComputeResidual(const CConfig* config) final;
-};
+};

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

@@ -38,26 +38,14 @@
 class CUpwAUSMPLUSUP2_NEMO : public CNEMONumerics {
 private:
   su2double *FcL, *FcR, *FcLR;
-  su2double *dmLP, *dmRM, *dpLP, *dpRM;
-  su2double *daL, *daR;
+  //su2double *dmLP, *dmRM, *dpLP, *dpRM;
+  //su2double *daL, *daR;
   su2double ProjVel_i, ProjVel_j;
-  su2double sq_vel, Proj_ModJac_Tensor_ij;
   su2double mL, mR, mLP, mRM, mF, pLP, pRM, pFi, pF, Phi;
   su2double CstarL, CstarR, ChatL, ChatR, aF, rhoF, MFsq, Mrefsq, Mp, fa;
   su2double Kp, sigma, alpha, beta, param1, mfP, mfM;
   su2double rhoCvtr_i, rhoCvtr_j, rhoCvve_i, rhoCvve_j;
 
-  /*--- Roe Only ---*/
-  su2double *Diff_U;
-  su2double *RoeU, *RoeV, *RoeEve;
-  su2double *ProjFlux_i, *ProjFlux_j;
-  su2double *Lambda, *Epsilon;
-  su2double **P_Tensor, **invP_Tensor;
-  su2double RoeSoundSpeed;
-  su2double ProjVelocity, ProjVelocity_i, ProjVelocity_j;
-  su2double R;
-  su2double *RoedPdU;
-
 public:
 
   /*!

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

@@ -34,21 +34,18 @@
  * \brief Class for computing the Lax-Friedrich centered scheme.
  * \ingroup ConvDiscr
  * \author F. Palacios, S.R. Copeland, W. Maier, C. Garbacz
- * \version 2.0.6
+ * \version 7.0.8
  */
 class CCentLax_NEMO : public CNEMONumerics {
 private:
-  unsigned short iDim, iVar, jVar; /*!< \brief Iteration on dimension and variables. */
+  unsigned short iDim, iVar; /*!< \brief Iteration on dimension and variables. */
   su2double *Diff_U; /*!< \brief Difference of conservative variables. */
   su2double *MeanU, *MeanV;
-  vector<su2double> mean_eves;
-  su2double *MeandPdU;
   su2double *ProjFlux;  /*!< \brief Projected inviscid flux tensor. */
   su2double Param_p, Param_Kappa_0; /*!< \brief Artificial dissipation parameters. */
   su2double Local_Lambda_i, Local_Lambda_j, MeanLambda; /*!< \brief Local eigenvalues. */
   su2double Phi_i, Phi_j, sc0, StretchingFactor; /*!< \brief Streching parameters. */
-  su2double Epsilon_0, cte; /*!< \brief Artificial dissipation values. */
-  bool stretching;
+  su2double Epsilon_0; /*!< \brief Artificial dissipation values. */
 
 public:
 

SU2_CFD/include/solvers/CNEMOEulerSolver.hpp

@@ -56,6 +56,9 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, COMPRESSI
 
   unsigned long ErrorCounter = 0; /*!< \brief Counter for number of un-physical states. */
 
+  su2double Global_Delta_Time = 0.0, /*!< \brief Time-step for TIME_STEPPING time marching strategy. */
+  Global_Delta_UnstTimeND = 0.0;     /*!< \brief Unsteady time step for the dual time strategy. */
+
   CNEMOGas  *FluidModel;          /*!< \brief fluid model used in the solver */
 
   CNEMOEulerVariable* node_infty = nullptr;
@@ -75,30 +78,27 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, COMPRESSI
   CNEMOEulerSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, const bool navier_stokes = false);
 
   /*!
-     * \brief Destructor of the class.
-     */
+   * \brief Destructor of the class.
+   */
   ~CNEMOEulerSolver(void) override;
 
   /*!
-     * \brief Set the maximum value of the eigenvalue.
-     * \param[in] geometry - Geometrical definition of the problem.
-     * \param[in] config - Definition of the particular problem.
-     */
+   * \brief Set the maximum value of the eigenvalue.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
   void SetMax_Eigenvalue(CGeometry *geometry, CConfig *config);
 
-   /*!
-  * \brief Compute the time step for solving the Euler equations.
-  * \param[in] geometry - Geometrical definition of the problem.
-  * \param[in] solver_container - Container vector with all the solutions.
-  * \param[in] config - Definition of the particular problem.
-  * \param[in] iMesh - Index of the mesh in multigrid computations.
-  * \param[in] Iteration - Value of the current iteration.
-    */
-  void SetTime_Step(CGeometry *geometry,
-                    CSolver **solver_container,
-                    CConfig *config,
-                    unsigned short iMesh,
-                    unsigned long Iteration) final;
+  /*!
+   * \brief Compute the time step for solving the Euler equations.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] Iteration - Value of the current iteration.
+   */
+  void SetTime_Step(CGeometry *geometry, CSolver **solver_container,
+                    CConfig *config, unsigned short iMesh, unsigned long Iteration) final;
 
   /*!
    * \brief Set the initial condition for the Euler Equations.
@@ -138,12 +138,21 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, COMPRESSI
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
-  void Upwind_Residual(CGeometry *geometry,
-                       CSolver **solver_container,
-                       CNumerics **numerics_container,
-                       CConfig *config,
-                       unsigned short iMesh) final;
+  void Upwind_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
+                       CConfig *config, unsigned short iMesh) final;
 
+  /*!
+   * \brief Recompute the extrapolated quantities, after MUSCL reconstruction,
+   *        in a more thermodynamically consistent way.
+   * \param[in] V - primitve variables.
+   * \param[out] d*dU - reconstructed secondaryvariables.
+   * \param[out] val_eves - reconstructed eve per species.
+   * \param[out] val_cvves - reconstructed cvve per species.
+   * \param[out] Gamma - reconstructed gamma.
+   */
+  static su2double ComputeConsistentExtrapolation(CNEMOGas *fluidmodel, unsigned short nSpecies, su2double *V,
+                                                  su2double* dPdU, su2double* dTdU, su2double* dTvedU,
+                                                  su2double* val_eves, su2double* val_cvves);
   /*!
    * \brief Source term integration.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -153,11 +162,8 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, COMPRESSI
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
-  void Source_Residual(CGeometry *geometry,
-                       CSolver **solver_container,
-                       CNumerics **numerics_container,
-                       CConfig *config,
-                       unsigned short iMesh) final;
+  void Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
+                       CConfig *config, unsigned short iMesh) final;
 
   /*!
    * \brief Preprocessing actions common to the Euler and NS solvers.
@@ -199,6 +205,23 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, COMPRESSI
    */
   void SetNondimensionalization(CConfig *config, unsigned short iMesh) final;
 
+ /*!
+  * \brief Set all the conserved variables from the primitive vector..
+  */
+  void RecomputeConservativeVector(su2double *U, const su2double *V) const;
+
+ /*!
+  * \brief Check for unphysical points.
+  * \return Boolean value of physical point
+  */
+  bool CheckNonPhys(const su2double *V) const;
+
+    /*!
+   * \brief Compute the pressure at the infinity.
+   * \return Value of the pressure at the infinity.
+   */
+  inline CNEMOGas* GetFluidModel(void) const final { return FluidModel;}
+
   /*!
    * \brief Impose the far-field boundary condition using characteristics.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CNEMONSSolver.hpp

@@ -38,7 +38,8 @@
  * \brief Main class for defining the NEMO Navier-Stokes flow solver.
  * \ingroup Navier_Stokes_Equations
  * \author S. R. Copeland, F. Palacios, W. Maier.
- * \version 6.1
+ * \version 7.0.8
+ *
  */
 class CNEMONSSolver final : public CNEMOEulerSolver {
 private:

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -35,7 +35,7 @@
  * \brief Main class for defining the variables of the NEMO Euler's solver.
  * \ingroup Euler_Equations
  * \author S. R. Copeland, F. Palacios, W. Maier, C. Garbacz
- * \version 7.0.6
+ * \version 7.0.8
  */
 class CNEMOEulerVariable : public CVariable {
 public:
@@ -157,12 +157,29 @@ class CNEMOEulerVariable : public CVariable {
    * \brief Get the primitive variables limiter.
    * \return Primitive variables limiter for the entire domain.
    */
-  inline MatrixType& GetLimiter_Primitive(void) {
+  inline MatrixType& GetLimiter_Primitive(void) {return Limiter_Primitive; }
 
-    SU2_MPI::Error(string("Limiters (associated to MUSCL) are computed for conserved variables in the NEMO solver.") +
-                   string("Limiters for primitive variables are not allocated/computed."),
-                   CURRENT_FUNCTION);
-    return Primitive;
+  /*!
+   * \brief Set the gradient of the primitive variables.
+   * \param[in] iVar - Index of the variable.
+   * \param[in] iDim - Index of the dimension.
+   * \param[in] value - Value of the gradient.
+   */
+  inline su2double GetLimiter_Primitive(unsigned long iPoint, unsigned long iVar) const final {return Limiter_Primitive(iPoint,iVar); }
+
+  /*!
+   * \brief Get the value of the primitive variables gradient.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double *GetLimiter_Primitive(unsigned long iPoint) final { return Limiter_Primitive[iPoint]; }
+
+  /*!
+   * \brief Set the gradient of the primitive variables.
+   * \param[in] iVar - Index of the variable.
+   * \param[in] value - Value of the gradient.
+   */
+  inline void SetLimiter_Primitive(unsigned long iPoint, unsigned long iVar, su2double value) final {
+    Limiter_Primitive(iPoint,iVar) = value;
   }
 
   /*!
@@ -332,21 +349,13 @@ class CNEMOEulerVariable : public CVariable {
    */
   bool SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) override;
 
- /*!
-  * \brief Set all the conserved variables.
+   /*!
+  * \brief Set all the primitive and secondary variables from the conserved vector.
   */
   bool Cons2PrimVar(su2double *U, su2double *V, su2double *dPdU,
                     su2double *dTdU, su2double *dTvedU, su2double *val_eves,
                     su2double *val_Cvves);
 
- /*!
-  * \brief Check for unphysical points.
-  * \return Boolean value of physical point 
-  */
-  bool CheckNonPhys(su2double *U, su2double *V, su2double *dPdU,
-                    su2double *dTdU, su2double *dTvedU, su2double *val_eves,
-                    su2double *val_Cvves);
-
   /*---------------------------------------*/
   /*---   Specific variable routines    ---*/
   /*---------------------------------------*/

SU2_CFD/include/variables/CNEMONSVariable.hpp

@@ -34,7 +34,7 @@
  * \brief Main class for defining the variables of the NEMO Navier-Stokes' solver.
  * \ingroup Navier_Stokes_Equations
  * \author C. Garbacz, W. Maier, S.R. Copeland.
- * \version 7.0.6
+ * \version 7.0.8
  */
 class CNEMONSVariable final : public CNEMOEulerVariable {
 private:

SU2_CFD/include/variables/CVariable.hpp

@@ -1426,16 +1426,6 @@ class CVariable {
    */
   inline virtual void SetSecondaryVar(unsigned long iPoint, CFluidModel *FluidModel) {}
 
-  /*!
-   * \brief A virtual member.
-   */
-  inline virtual bool Cons2PrimVar(CConfig *config, unsigned long iPoint, su2double *U, su2double *V, su2double *dPdU,
-                                   su2double *dTdU, su2double *dTvedU) { return false; }
-  /*!
-   * \brief A virtual member.
-   */
-  inline virtual void Prim2ConsVar(CConfig *config, unsigned long iPoint, su2double *V, su2double *U) { }
-
   /*!
    * \brief A virtual member.
    */
@@ -1590,45 +1580,6 @@ class CVariable {
    */
   inline virtual bool SetPressure(unsigned long iPoint, su2double Gamma, su2double turb_ke) { return false; }
 
-  /*!
-   * \brief Calculates vib.-el. energy per mass, \f$e^{vib-el}_s\f$, for input species (not including KE)
-   */
-  inline virtual su2double CalcEve(unsigned long iPoint, su2double *V, CConfig *config, unsigned long val_Species) { return 0.0; }
-
-  /*!
-   * \brief Calculates enthalpy per mass, \f$h_s\f$, for input species (not including KE)
-   */
-  inline virtual su2double CalcHs(unsigned long iPoint, su2double *V, CConfig *config, unsigned long val_Species) { return 0.0; }
-
-  /*!
-   * \brief Calculates enthalpy per mass, \f$Cv_s\f$, for input species (not including KE)
-   */
-  inline virtual su2double CalcCvve(unsigned long iPoint, su2double val_Tve, CConfig *config, unsigned long val_Species) { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \param[in] V
-   * \param[in] config - Configuration settings
-   * \param[in] dPdU
-   */
-  inline virtual void CalcdPdU(unsigned long iPoint, su2double *V, CConfig *config, su2double *dPdU) {}
-
-  /*!
-   * \brief Set partial derivative of temperature w.r.t. density \f$\frac{\partial P}{\partial \rho_s}\f$
-   * \param[in] V
-   * \param[in] config - Configuration settings
-   * \param[in] dTdU
-   */
-  inline virtual void CalcdTdU(unsigned long iPoint, su2double *V, CConfig *config, su2double *dTdU) {}
-
-  /*!
-   * \brief Set partial derivative of temperature w.r.t. density \f$\frac{\partial P}{\partial \rho_s}\f$
-   * \param[in] V
-   * \param[in] config - Configuration settings
-   * \param[in] dTdU
-   */
-  inline virtual void CalcdTvedU(unsigned long iPoint, su2double *V, CConfig *config, su2double *dTdU) {}
-
   /*!
    * \brief A virtual member.
    */

SU2_CFD/src/fluid/CNEMOGas.cpp

@@ -176,6 +176,7 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
 
   /*--- Necessary indexes to assess primitive variables ---*/
   unsigned long RHOS_INDEX    = 0;
+  unsigned long T_INDEX       = nSpecies;
   unsigned long VEL_INDEX     = nSpecies+2;
   unsigned long RHOCVTR_INDEX = nSpecies+nDim+6;
   unsigned long RHOCVVE_INDEX = nSpecies+nDim+7;
@@ -191,6 +192,7 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
   /*--- Rename for convenience ---*/
   rhoCvtr = V[RHOCVTR_INDEX];
   rhoCvve = V[RHOCVVE_INDEX];
+  T       = V[T_INDEX];
 
   /*--- Pre-compute useful quantities ---*/
   CvtrBAR = 0.0;
@@ -249,10 +251,12 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
   su2double Vel[3] = {0.0};
 
   /*--- Necessary indexes to assess primitive variables ---*/
+  unsigned long T_INDEX       = nSpecies;
   unsigned long VEL_INDEX     = nSpecies+2;
   unsigned long RHOCVTR_INDEX = nSpecies+nDim+6;
 
   /*--- Rename for convenience ---*/
+  T       = V[T_INDEX];
   rhoCvtr = V[RHOCVTR_INDEX];
 
   Cvtrs              = GetSpeciesCvTraRot();