Fixed issue with test failue due to missing chemistry source residual term, addressed PR comments (spacing, initialization)

Signed-off-by: jtneedels <jneedels@stanford.edu>

Author: jtneedels

SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp

@@ -342,43 +342,42 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
   residual[nSpecies+2] = yinv*rhov*H*Volume;
   residual[nSpecies+3] = yinv*rhov*U_i[nSpecies+nDim+1]/rho*Volume;
 
-     if (viscous) {
+  if (viscous) {
 
-        if (!rans){ turb_ke_i = 0.0; }
+    if (!rans){ turb_ke_i = 0.0; }
 
-        su2double sumJhs_y  = 0.0;
-        su2double sumJeve_y = 0.0;
-        su2double Mass      = 0.0;
+    su2double sumJhs_y  = 0.0;
+    su2double sumJeve_y = 0.0;
+    su2double Mass      = 0.0;
 
-        for (iSpecies=0; iSpecies<nSpecies; iSpecies++)
-          Mass += V_i[iSpecies]*Ms[iSpecies];
+    for (iSpecies=0; iSpecies<nSpecies; iSpecies++)
+      Mass += V_i[iSpecies]*Ms[iSpecies];
 
-        su2double heat_capacity_cp_i   = V_i[RHOCVTR_INDEX]+Ru/Mass;
-        su2double total_viscosity_i    = Laminar_Viscosity_i + Eddy_Viscosity_i;
-        su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;      
-        su2double u = V_i[VEL_INDEX];
-        su2double v = V_i[VEL_INDEX+1];
-
-        qy_ve = kve*GV[TVE_INDEX][1];
-
-        /*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
-        for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { 
-          sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
-          sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
-        }
-
-        for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
-          residual[iSpecies] -= 0.0;
-        residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0]) 
+    su2double heat_capacity_cp_i   = V_i[RHOCVTR_INDEX]+Ru/Mass;
+    su2double total_viscosity_i    = Laminar_Viscosity_i + Eddy_Viscosity_i;
+    su2double total_conductivity_i = ktr + kve + heat_capacity_cp_i*Eddy_Viscosity_i/Prandtl_Turb;      
+    su2double u                    = V_i[VEL_INDEX];
+    su2double v                    = V_i[VEL_INDEX+1];
+    su2double qy_ve                = kve*GV[TVE_INDEX][1];
+
+    /*--- Enthalpy and vib-el energy transport due to y-direction diffusion---*/
+    for (iSpecies = 0; iSpecies < nHeavy; iSpecies++) { 
+      sumJhs_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * hs[iSpecies];
+      sumJeve_y += (rho*Ds[iSpecies]*GV[RHOS_INDEX+iSpecies][1] - V_i[RHOS_INDEX+iSpecies]*Vector[1]) * eve_i[iSpecies];
+    }
+
+    for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+      residual[iSpecies] -= 0.0;
+    residual[nSpecies] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[nSpecies+2][1]+PrimVar_Grad_i[nSpecies+3][0]) 
                                                                     -TWO3*AuxVar_Grad_i[0][0]);
-        residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
+    residual[nSpecies+1] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[nSpecies+3][1]-v*yinv)
                                                                     -TWO3*AuxVar_Grad_i[0][1]);
-        residual[nSpecies+2] -= Volume*(yinv*(- sumJhs_y + total_viscosity_i*(u*(PrimVar_Grad_i[nSpecies+3][0]+PrimVar_Grad_i[nSpecies+2][1])
+    residual[nSpecies+2] -= Volume*(yinv*(- sumJhs_y + total_viscosity_i*(u*(PrimVar_Grad_i[nSpecies+3][0]+PrimVar_Grad_i[nSpecies+2][1])
                                                      +v*TWO3*(2*PrimVar_Grad_i[nSpecies+2][1]-PrimVar_Grad_i[nSpecies+2][0]
                                                      -v*yinv+rho*turb_ke_i))
                                                      -total_conductivity_i*PrimVar_Grad_i[nSpecies][1])
                                                      -TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1]));
-        residual[nSpecies+3] -= Volume*(yinv*(-sumJeve_y -qy_ve));
+    residual[nSpecies+3] -= Volume*(yinv*(-sumJeve_y -qy_ve));
   }
 
 //  if (implicit) {

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -1196,6 +1196,7 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
         /*--- Set primitive variables for viscous terms and/or generalised source ---*/
         numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));
 
+        /*--- If necessary, set variables needed for viscous computation ---*/       
         if (viscous) {
 
           /*--- Set gradient of primitive variables ---*/
@@ -1229,29 +1230,52 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
           }
         }
 
-          auto residual = numerics->ComputeAxisymmetric(config);
+        auto residual = numerics->ComputeAxisymmetric(config);
 
-          /*--- Check for errors before applying source to the linear system ---*/
-          err = false;
+        /*--- Check for errors before applying source to the linear system ---*/
+        err = false;
+        for (iVar = 0; iVar < nVar; iVar++)
+          if (residual[iVar] != residual[iVar]) err = true;
+        if (implicit)
           for (iVar = 0; iVar < nVar; iVar++)
-            if (residual[iVar] != residual[iVar]) err = true;
+            for (jVar = 0; jVar < nVar; jVar++)
+              if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
+
+        /*--- Apply the update to the linear system ---*/
+        if (!err) {
+          LinSysRes.AddBlock(iPoint, residual);
           if (implicit)
-            for (iVar = 0; iVar < nVar; iVar++)
-              for (jVar = 0; jVar < nVar; jVar++)
-                if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
-
-          /*--- Apply the update to the linear system ---*/
-          if (!err) {
-            LinSysRes.AddBlock(iPoint, residual);
-            if (implicit)
-              Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
-          }
-          else
-            eAxi_local++;
+            Jacobian.AddBlock(iPoint, iPoint, Jacobian_i);
+        }else
+          eAxi_local++;
       }
 
     }
 
+    if(!monoatomic){
+      if(!frozen){
+        /*--- Compute the non-equilibrium chemistry ---*/
+        auto residual = numerics->ComputeChemistry(config);
+
+        /*--- Check for errors before applying source to the linear system ---*/
+        err = false;
+        for (iVar = 0; iVar < nVar; iVar++)
+          if (residual[iVar] != residual[iVar]) err = true;
+        //if (implicit)
+        //  for (iVar = 0; iVar < nVar; iVar++)
+        //    for (jVar = 0; jVar < nVar; jVar++)
+        //      if (Jacobian_i[iVar][jVar] != Jacobian_i[iVar][jVar]) err = true;
+
+        /*--- Apply the chemical sources to the linear system ---*/
+        if (!err) {
+          LinSysRes.SubtractBlock(iPoint, residual);
+          //if (implicit)
+          //  Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
+        } else
+          eChm_local++;
+      }
+    }
+
     /*--- Compute vibrational energy relaxation ---*/
     /// NOTE: Jacobians don't account for relaxation time derivatives
 

SU2_CFD/src/variables/CNEMOEulerVariable.cpp

@@ -78,9 +78,8 @@ CNEMOEulerVariable::CNEMOEulerVariable(su2double val_pressure,
 
   Res_TruncError.resize(nPoint,nVar) = su2double(0.0);
 
-    /*--- Size Grad_AuxVar for axiysmmetric ---*/
-
- if (config->GetAxisymmetric()){
+  /*--- Size Grad_AuxVar for axiysmmetric ---*/
+  if (config->GetAxisymmetric()){
     nAuxVar = 3;
     Grad_AuxVar.resize(nPoint,nAuxVar,nDim,0.0);
     AuxVar.resize(nPoint,nAuxVar) = su2double(0.0);