Merge pull request #1326 from su2code/nemo_update

Fix for axisymmetric terms in NEMO + general NEMO updates

Author: WallyMaier

Common/src/CConfig.cpp

@@ -3601,14 +3601,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("AUSMPW+ is extremely unstable. Feel free to fix me!", CURRENT_FUNCTION);
   }
 
-  if (GetGasModel() == "ARGON" && GetKind_FluidModel() == SU2_NONEQ){
-      SU2_MPI::Error("ARGON is not working with SU2_NONEQ fluid model!", CURRENT_FUNCTION);
-  }
-
-  if (GetKind_FluidModel() == MUTATIONPP && GetFrozen() == true){
-      SU2_MPI::Error("The option of FROZEN_MIXTURE is not yet working with Mutation++ support.", CURRENT_FUNCTION);
-  }
-
   if(GetBoolTurbomachinery()){
     nBlades = new su2double[nZone];
     FreeStreamTurboNormal= new su2double[3];
@@ -5042,7 +5034,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   /*--- Specifying a deforming surface requires a mesh deformation solver. ---*/
   if (GetSurface_Movement(DEFORMING)) Deform_Mesh = true;
 
-  if (GetGasModel() == "ARGON") monoatomic = true;
+  if (GetGasModel() == "ARGON") {monoatomic = true;}
+  else {monoatomic = false;}
 
   // This option is deprecated. After a grace period until 7.2.0 the usage warning should become an error.
   if(OptionIsSet("CONV_CRITERIA") && rank == MASTER_NODE) {

SU2_CFD/include/fluid/CMutationTCLib.hpp

@@ -87,7 +87,7 @@ class CMutationTCLib : public CNEMOGas {
   /*!
    * \brief Compute vector of species V-E energy.
    */
-  vector<su2double>& ComputeSpeciesEve(su2double val_T) final;
+  vector<su2double>& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) final;
 
   /*!
    * \brief Compute species net production rates.

SU2_CFD/include/fluid/CNEMOGas.hpp

@@ -134,7 +134,7 @@ class CNEMOGas : public CFluidModel {
   /*!
    * \brief Compute vector of species V-E energy.
    */
-  virtual vector<su2double>& ComputeSpeciesEve(su2double val_T) = 0;
+  virtual vector<su2double>& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) = 0;
 
   /*!
    * \brief Compute species enthalpies.

SU2_CFD/include/fluid/CSU2TCLib.hpp

@@ -50,7 +50,7 @@ class CSU2TCLib : public CNEMOGas {
   ArrheniusEta,                     /*!< \brief Arrhenius reaction temperature exponent */
   ArrheniusTheta,                   /*!< \brief Arrhenius reaction characteristic temperature */
   CharVibTemp,                      /*!< \brief Characteristic vibrational temperature for e_vib */
-  RotationModes,	          /*!< \brief Rotational modes of energy storage */
+  RotationModes,                  /*!< \brief Rotational modes of energy storage */
   Tcf_a,                          /*!< \brief Rate controlling temperature exponent (fwd) */
   Tcf_b,                          /*!< \brief Rate controlling temperature exponent (fwd) */
   Tcb_a,                          /*!< \brief Rate controlling temperature exponent (bkw) */
@@ -115,7 +115,7 @@ class CSU2TCLib : public CNEMOGas {
   /*!
    * \brief Compute species V-E energy.
    */
-  vector<su2double>& ComputeSpeciesEve(su2double val_T) final;
+  vector<su2double>& ComputeSpeciesEve(su2double val_T, bool vibe_only = false) final;
 
   /*!
    * \brief Compute species net production rates.

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -277,6 +277,11 @@ class CFVMFlowSolverBase : public CSolver {
    */
   void SumEdgeFluxes(const CGeometry* geometry);
 
+  /*!
+   * \brief Computes and sets the required auxilliary vars (and gradients) for axisymmetric flow.
+   */
+  void ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config);
+
   /*!
    * \brief Instantiate a SIMD numerics object.
    */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2923,3 +2923,32 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
 
   return objFun;
 }
+
+template <class V, ENUM_REGIME FlowRegime>
+void CFVMFlowSolverBase<V, FlowRegime>::ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config) {
+
+  /*--- Loop through all points to set the auxvargrad --*/
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
+    su2double yCoord          = geometry->nodes->GetCoord(iPoint, 1);
+    su2double yVelocity       = nodes->GetVelocity(iPoint,1);
+    su2double xVelocity       = nodes->GetVelocity(iPoint,0);
+    su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint);
+
+    if (yCoord > EPS){
+      su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord;
+      nodes->SetAuxVar(iPoint, 0, nu_v_on_y);
+      nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity);
+      nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity);
+    }
+  }
+  END_SU2_OMP_FOR
+
+  /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/
+  if (config->GetKind_Gradient_Method() == GREEN_GAUSS) {
+    SetAuxVar_Gradient_GG(geometry, config);
+  }
+  if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) {
+    SetAuxVar_Gradient_LS(geometry, config);
+  }
+}

SU2_CFD/src/fluid/CMutationTCLib.cpp

@@ -48,6 +48,7 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
     transport_model = "Gupta-Yos";
 
   opt.setStateModel("ChemNonEqTTv");
+  if (frozen) opt.setMechanism("none");
   opt.setViscosityAlgorithm(transport_model);
   opt.setThermalConductivityAlgorithm(transport_model);
 
@@ -118,7 +119,7 @@ vector<su2double>& CMutationTCLib::ComputeMixtureEnergies(){
   return energies;
 }
 
-vector<su2double>& CMutationTCLib::ComputeSpeciesEve(su2double val_T){
+vector<su2double>& CMutationTCLib::ComputeSpeciesEve(su2double val_T, bool vibe_only){
 
   SetTDStateRhosTTv(rhos, T, val_T);
 

SU2_CFD/src/fluid/CNEMOGas.cpp

@@ -40,7 +40,7 @@ CNEMOGas::CNEMOGas(const CConfig* config, unsigned short val_nDim): CFluidModel(
   Cvtrs.resize(nSpecies,0.0);
   Cvves.resize(nSpecies,0.0);
   eves.resize(nSpecies,0.0);
-  hs.resize(nSpecies,0.0);
+  hs.resize(2*nSpecies,0.0);
   ws.resize(nSpecies,0.0);
   DiffusionCoeff.resize(nSpecies,0.0);
   Enthalpy_Formation.resize(nSpecies,0.0);
@@ -59,23 +59,21 @@ CNEMOGas::CNEMOGas(const CConfig* config, unsigned short val_nDim): CFluidModel(
 void CNEMOGas::SetTDStatePTTv(su2double val_pressure, const su2double *val_massfrac,
                               su2double val_temperature, su2double val_temperature_ve){
 
-  su2double denom;
-
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++)
+  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     MassFrac[iSpecies] = val_massfrac[iSpecies];
   Pressure = val_pressure;
   T        = val_temperature;
   Tve      = val_temperature_ve;
 
-  denom   = 0.0;
+  su2double denom = 0.0;
 
   /*--- Calculate mixture density from supplied primitive quantities ---*/
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++)
-    denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * T;
   for (iSpecies = 0; iSpecies < nEl; iSpecies++)
-    denom += MassFrac[nSpecies-1] * (Ru/MolarMass[nSpecies-1]) * Tve;
+    denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * Tve;
+  for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++)
+    denom += MassFrac[iSpecies] * (Ru/MolarMass[iSpecies]) * T;
   Density = Pressure / denom;
-
+  
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++){
     rhos[iSpecies]     = MassFrac[iSpecies]*Density;
     MassFrac[iSpecies] = rhos[iSpecies]/Density;
@@ -84,18 +82,16 @@ void CNEMOGas::SetTDStatePTTv(su2double val_pressure, const su2double *val_massf
 
 su2double CNEMOGas::ComputeSoundSpeed(){
 
-  su2double conc, rhoCvtr;
-
-  conc    = 0.0;
-  rhoCvtr = 0.0;
+  su2double conc    = 0.0;
+  su2double rhoCvtr = 0.0;
   Density = 0.0;
 
   auto& Cvtrs = GetSpeciesCvTraRot();
 
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     Density+=rhos[iSpecies];
 
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++){
+  for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++){
     conc += rhos[iSpecies]/MolarMass[iSpecies];
     rhoCvtr += rhos[iSpecies] * Cvtrs[iSpecies];
   }
@@ -108,10 +104,10 @@ su2double CNEMOGas::ComputePressure(){
 
   su2double P = 0.0;
 
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++)
-    P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * T;
   for (iSpecies = 0; iSpecies < nEl; iSpecies++)
-    P += rhos[nSpecies-1] * Ru/MolarMass[nSpecies-1] * Tve;
+    P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * Tve;
+  for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++)
+    P += rhos[iSpecies] * Ru/MolarMass[iSpecies] * T;
 
   Pressure = P;
 
@@ -124,7 +120,7 @@ su2double CNEMOGas::ComputeGasConstant(){
   su2double Mass = 0.0;
 
   // TODO - extend for ionization
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++)
+  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
     Mass += MassFrac[iSpecies] * MolarMass[iSpecies];
   GasConstant = Ru / Mass;
 
@@ -163,15 +159,13 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
 
   // Note: Electron energy not included properly.
 
-  su2double CvtrBAR, rhoCvtr, rhoCvve, rho_el, sqvel, conc, ef;
-
   if (val_dPdU == NULL) {
     SU2_MPI::Error("Array dPdU not allocated!", CURRENT_FUNCTION);
   }
 
   /*--- Determine the electron density (if ionized) ---*/
-  if (ionization) { rho_el = rhos[nSpecies-1]; }
-  else {            rho_el = 0.0;              }
+  su2double rho_el = 0.0;
+  if (ionization) { rho_el = rhos[0]; }
 
   /*--- Necessary indexes to assess primitive variables ---*/
   unsigned long RHOS_INDEX    = 0;
@@ -189,14 +183,14 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
   Ref_Temperature    = GetRefTemperature();
 
   /*--- Rename for convenience ---*/
-  rhoCvtr = V[RHOCVTR_INDEX];
-  rhoCvve = V[RHOCVVE_INDEX];
+  su2double rhoCvtr = V[RHOCVTR_INDEX];
+  su2double rhoCvve = V[RHOCVVE_INDEX];
   T       = V[T_INDEX];
 
   /*--- Pre-compute useful quantities ---*/
-  CvtrBAR = 0.0;
-  sqvel   = 0.0;
-  conc    = 0.0;
+  su2double CvtrBAR = 0.0;
+  su2double sqvel   = 0.0;
+  su2double conc    = 0.0;
   for (iDim = 0; iDim < nDim; iDim++)
     sqvel += V[VEL_INDEX+iDim] * V[VEL_INDEX+iDim];
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
@@ -205,15 +199,16 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
   }
 
   /*--- Species density derivatives ---*/
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
-    ef                 = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies];
+  su2double ef = 0.0;
+  for (iSpecies = nEl; iSpecies < nHeavy; iSpecies++) {
+    ef = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies];
     val_dPdU[iSpecies] = T*Ru/MolarMass[iSpecies] + Ru*conc/rhoCvtr *
                          (-Cvtrs[iSpecies]*(T-Ref_Temperature[iSpecies]) -
                          ef + 0.5*sqvel);
 
   }
   if (ionization) {
-    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
+    for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) {
       //      evibs = Ru/MolarMass[iSpecies] * thetav[iSpecies]/(exp(thetav[iSpecies]/Tve)-1.0);
       //      num = 0.0;
       //      denom = g[iSpecies][0] * exp(-thetae[iSpecies][0]/Tve);
@@ -223,12 +218,11 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
       //      }
       //      eels = Ru/MolarMass[iSpecies] * (num/denom);
 
-      val_dPdU[iSpecies] -= rho_el * Ru/MolarMass[nSpecies-1] * (val_eves[iSpecies])/rhoCvve;
+      val_dPdU[iSpecies] -= rho_el * Ru/MolarMass[0] * (val_eves[iSpecies])/rhoCvve;
     }
-    ef = Enthalpy_Formation[nSpecies-1] - Ru/MolarMass[nSpecies-1]*Ref_Temperature[nSpecies-1];
-    val_dPdU[nSpecies-1] = Ru*conc/rhoCvtr * (-ef + 0.5*sqvel)
-        + Ru/MolarMass[nSpecies-1]*Tve
-        - rho_el*Ru/MolarMass[nSpecies-1] * (-3.0/2.0*Ru/MolarMass[nSpecies-1]*Tve)/rhoCvve;
+    ef = Enthalpy_Formation[0] - Ru/MolarMass[0]*Ref_Temperature[0];
+    val_dPdU[0] = Ru*conc/rhoCvtr * (-ef + 0.5*sqvel) + Ru/MolarMass[0]*Tve
+                - rho_el*Ru/MolarMass[0] * (-3.0/2.0*Ru/MolarMass[0]*Tve)/rhoCvve;
   }
 
   /*--- Momentum derivatives ---*/
@@ -240,13 +234,12 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
 
   /*--- Vib.-el energy derivative ---*/
   val_dPdU[nSpecies+nDim+1] = -val_dPdU[nSpecies+nDim] +
-                               rho_el*Ru/MolarMass[nSpecies-1]*1.0/rhoCvve;
+                               rho_el*Ru/MolarMass[0]*1.0/rhoCvve;
 
 }
 
 void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
 
-  su2double v2, ef, rhoCvtr;
   su2double Vel[3] = {0.0};
 
   /*--- Necessary indexes to assess primitive variables ---*/
@@ -255,8 +248,8 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
   unsigned long RHOCVTR_INDEX = nSpecies+nDim+6;
 
   /*--- Rename for convenience ---*/
-  T       = V[T_INDEX];
-  rhoCvtr = V[RHOCVTR_INDEX];
+  T                 = V[T_INDEX];
+  su2double rhoCvtr = V[RHOCVTR_INDEX];
 
   Cvtrs              = GetSpeciesCvTraRot();
   Enthalpy_Formation = GetSpeciesFormationEnthalpy();
@@ -265,11 +258,11 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
   /*--- Calculate supporting quantities ---*/
   for (iDim = 0; iDim < nDim; iDim++)
     Vel[iDim] = V[VEL_INDEX+iDim]*V[VEL_INDEX+iDim];
-  v2 = GeometryToolbox::SquaredNorm(nDim,Vel);
+  su2double v2 = GeometryToolbox::SquaredNorm(nDim,Vel);
 
   /*--- Species density derivatives ---*/
-  for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) {
-    ef    = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies];
+  for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) {
+    su2double ef    = Enthalpy_Formation[iSpecies] - Ru/MolarMass[iSpecies]*Ref_Temperature[iSpecies];
     val_dTdU[iSpecies]   = (-ef + 0.5*v2 + Cvtrs[iSpecies]*(Ref_Temperature[iSpecies]-T)) / rhoCvtr;
   }
 
@@ -289,13 +282,11 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
 
 void CNEMOGas::ComputedTvedU(su2double *V, vector<su2double>& val_eves, su2double *val_dTvedU){
 
-  su2double rhoCvve;
-
   /*--- Necessary indexes to assess primitive variables ---*/
   unsigned long RHOCVVE_INDEX = nSpecies+nDim+7;
 
   /*--- Rename for convenience ---*/
-  rhoCvve = V[RHOCVVE_INDEX];
+  su2double rhoCvve = V[RHOCVVE_INDEX];
 
   /*--- Species density derivatives ---*/
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {