add protection for zero electron massfraction for GY

Author: WallyMaier

Common/src/CConfig.cpp

@@ -3936,8 +3936,10 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Transport model not available for NEMO solver using SU2TCLIB. Please use the WILKE, SUTHERLAND or GUPTAYOS transport model instead.", CURRENT_FUNCTION);
     }
 
-    if (Kind_FluidModel == SU2_NONEQ && GasModel == "AIR-7" && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) {
-      SU2_MPI::Error("Only Gupta-Yos transport model available for ionized flows using SU2TCLIB.", CURRENT_FUNCTION);
+    if (Kind_Solver == NEMO_NAVIER_STOKES) {
+      if (Kind_FluidModel == SU2_NONEQ && GasModel == "AIR-7" && Kind_TransCoeffModel != TRANSCOEFFMODEL::GUPTAYOS) {
+        SU2_MPI::Error("Only Gupta-Yos transport model available for ionized flows using SU2TCLIB.", CURRENT_FUNCTION);
+      }
     }
 
     if (Kind_FluidModel == MUTATIONPP &&

SU2_CFD/src/fluid/CMutationTCLib.cpp

@@ -58,7 +58,9 @@ CMutationTCLib::CMutationTCLib(const CConfig* config, unsigned short val_nDim):
   /* Initialize mixture object */
   mix.reset(new Mutation::Mixture(opt));
 
-  for(iSpecies = 0; iSpecies < nSpecies; iSpecies++) MolarMass[iSpecies] = 1000* mix->speciesMw(iSpecies); // x1000 to have Molar Mass in kg/kmol
+  // x1000 to have Molar Mass in kg/kmol
+  for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+    MolarMass[iSpecies] = 1000* mix->speciesMw(iSpecies);
 
   if (mix->hasElectrons()) { nHeavy = nSpecies-1; nEl = 1; }
   else { nHeavy = nSpecies; nEl = 0; }

SU2_CFD/src/fluid/CSU2TCLib.cpp

@@ -1952,21 +1952,23 @@ void CSU2TCLib::DiffusionCoeffGY(){
     DiffusionCoeff[iSpecies] = 0.0;
 
     /*--- Calculate molar concentration ---*/
-    const su2double Mi    = MolarMass[iSpecies];
+    const su2double Mi    = (MolarMass[iSpecies] + EPS);
     const su2double gam_i = rhos[iSpecies] / (Density*Mi);
     su2double denom = 0.0;
 
     for (jSpecies = 0; jSpecies < nSpecies; jSpecies++) {
       if (jSpecies != iSpecies) { 
 
-        const su2double Mj    = MolarMass[jSpecies];
+        const su2double Mj    = (MolarMass[jSpecies] + EPS);
         const su2double gam_j = rhos[jSpecies] / (Density*Mj);
 
         const su2double kb = BOLTZMANN_CONSTANT;
 
         const su2double T_col = (iSpecies == 0 && ionization) ? Tve : T; 
 
-        const su2double d1_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, true);
+        su2double d1_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, true);
+        if (d1_ij > 1E16) d1_ij = 1E16;
+
         const su2double D_ij = kb*T_col/(Pressure*d1_ij);
         denom += gam_j/D_ij;
       }
@@ -1987,16 +1989,18 @@ void CSU2TCLib::ViscosityGY(){
     su2double denom = 0.0;
 
     /*--- Calculate molar concentration ---*/
-    const su2double Mi    = MolarMass[iSpecies];
+    const su2double Mi    = (MolarMass[iSpecies] + EPS);
     const su2double gam_i = rhos[iSpecies] / (Density*Mi);
 
     for (jSpecies = 0; jSpecies < nSpecies; jSpecies++) {
-      const su2double Mj    = MolarMass[jSpecies];
+      const su2double Mj    = (MolarMass[jSpecies] + EPS);
       const su2double gam_j = rhos[jSpecies] / (Density*Mj);
 
       const su2double T_col = (iSpecies == 0 && ionization) ? Tve : T; 
 
-      const su2double d2_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, false);
+      su2double d2_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, false);
+      if (d2_ij > 1E16) d2_ij = 1E16;
+
       denom += gam_j*d2_ij;
     }
     /*--- Calculate species laminar viscosity ---*/
@@ -2028,33 +2032,38 @@ void CSU2TCLib::ThermalConductivitiesGY(){
   for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 
     /*--- Calculate molar concentration ---*/
-    const su2double Mi      = MolarMass[iSpecies];
-    const su2double mi      = Mi/Na;
-    const su2double gam_i   = rhos[iSpecies] / (Density*Mi);
+    const su2double Mi    = (MolarMass[iSpecies] + EPS);
+    const su2double mi    = Mi/Na;
+    const su2double gam_i = rhos[iSpecies] / (Density*Mi);
     su2double denom_t = 0.0;
     su2double denom_r = 0.0;
     su2double denom_re = 0.0;
 
     for (jSpecies = 0; jSpecies < nSpecies; jSpecies++) {
-      const su2double Mj    = MolarMass[jSpecies];
+      const su2double Mj    = (MolarMass[jSpecies] + EPS);
       const su2double mj    = Mj/Na;
       const su2double gam_j = rhos[iSpecies] / (Density*Mj);
       const su2double a_ij = 1.0 + (1.0 - mi/mj)*(0.45 - 2.54*mi/mj) / ((1.0 + mi/mj)*(1.0 + mi/mj));
 
       const su2double T_col = ((iSpecies == 0 && ionization) || (jSpecies == 0 && ionization)) ? Tve : T; 
 
-      const su2double d1_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, true);
-      const su2double d2_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, false);
-      
-      if (jSpecies == 0 && ionization) {
-        denom_t += 3.54*gam_j*d2_ij;
-      } else {
-        denom_t += a_ij*gam_j*d2_ij;
-      }
+      su2double d1_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, true);
+      su2double d2_ij = ComputeCollisionDelta(iSpecies, jSpecies, Mi, Mj, T_col, false);
+      if (d1_ij >1E16) d1_ij = 1E16;
+      if (d2_ij >1E16) d2_ij = 1E16;      
+
+      if (jSpecies == 0 && ionization) { denom_t += 3.54*gam_j*d2_ij; }
+      else { denom_t += a_ij*gam_j*d2_ij; }
+
       denom_r += gam_j*d1_ij;
       denom_re += gam_j*d2_ij;
     }
 
+    /*--- Prevent divide by 0 ---*/
+    if (denom_t == 0.0) denom_t = EPS;
+    if (denom_r == 0.0) denom_r = EPS;
+    if (denom_re == 0.0) denom_re = EPS;
+
     /*--- Translational contribution to thermal conductivity ---*/
     if (!(ionization && iSpecies == 0)) ThermalCond_tr += ((15.0/4.0)*kb*gam_i/denom_t);
 

TestCases/parallel_regression.py

@@ -188,7 +188,7 @@ def main():
     ion_gy.cfg_dir = "nonequilibrium/visc_cylinder"
     ion_gy.cfg_file = "cyl_ion_gy.cfg"
     ion_gy.test_iter = 10
-    ion_gy.test_vals = [-11.629871, -4.156742, -4.692113, -4.949447, -5.146196, -4.991664, -6.893329, 5.990109, 5.990004, -0.014849, 0.000000, 90090.000000]
+    ion_gy.test_vals = [-11.629873, -4.165563, -4.702662, -4.950351, -5.146155, -4.993878, -6.893332, 5.990109, 5.990004, -0.014849, 0.000000, 90090.000000]
     ion_gy.su2_exec = "mpirun -n 2 SU2_CFD"
     ion_gy.timeout = 1600
     ion_gy.new_output = True