Merge pull request #1170 from su2code/feature_nemo_axi_comp

Restructure source residual computation to fix axisymmetric chemsitry/vib source computation

Author: pcarruscag

SU2_CFD/src/numerics/NEMO/NEMO_sources.cpp

@@ -372,12 +372,12 @@ CNumerics::ResidualType<> CSource_NEMO::ComputeAxisymmetric(const CConfig *confi
                                                                     -TWO3*AuxVar_Grad_i[0][0]);
     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

@@ -1147,6 +1147,57 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
     numerics->SetCoord(geometry->nodes->GetCoord(iPoint),
                        geometry->nodes->GetCoord(iPoint) );
 
+    /*--- Compute finite rate chemistry ---*/
+
+    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
+
+    if (!monoatomic){
+      auto residual = numerics->ComputeVibRelaxation(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 vibrational relaxation terms to the linear system ---*/
+      if (!err) {
+        LinSysRes.SubtractBlock(iPoint, residual);
+        //if (implicit)
+        //  Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
+      } else
+        eVib_local++;
+    }
+
+  }
     /*--- Compute axisymmetric source terms (if needed) ---*/
     if (config->GetAxisymmetric()) {
 
@@ -1180,18 +1231,6 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
       SU2_OMP_FOR_DYN(omp_chunk_size)
       for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
-        /*--- Set solution  ---*/
-        numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));
-
-        /*--- Set control volume ---*/
-        numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
-
-        /*--- Set y coordinate ---*/
-        numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint));
-
-        /*--- 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) {
 
@@ -1216,9 +1255,6 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
           /*--- Vib-el. thermal conductivity ---*/
           numerics->SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), nodes->GetThermalConductivity_ve(iPoint));
 
-          /*--- Vib-el energy ---*/
-          numerics->SetEve(nodes->GetEve(iPoint), nodes->GetEve(iPoint));
-
           /*--- Set turbulence kinetic energy ---*/
           if (rans){
             CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes();
@@ -1247,55 +1283,6 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con
       }
     }
 
-    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
-
-    if (!monoatomic){
-      auto residual = numerics->ComputeVibRelaxation(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 vibrational relaxation terms to the linear system ---*/
-      if (!err) {
-        LinSysRes.SubtractBlock(iPoint, residual);
-        //if (implicit)
-        //  Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i);
-      } else
-        eVib_local++;
-    }
-  }
-
   /*--- Checking for NaN ---*/
   eAxi_global = eAxi_local;
   eChm_global = eChm_local;