Move MUSCL reconstruction to CSolver

Author: bmunguia

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -223,14 +223,12 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
           }
 
           for (auto iVar = 0u; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) {
-            su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
-            su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
+            su2double V_ij = 0.0;
+            if (umusclFlow)
+              V_ij = V_j[iVar] - V_i[iVar];
 
-            if (umusclFlow) {
-              const su2double V_ij = V_j[iVar] - V_i[iVar];
-              Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappaFlow);
-              Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappaFlow);
-            }
+            su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umusclFlow, kappaFlow);
+            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umusclFlow, kappaFlow);
 
             if (limiterFlow) {
               Project_Grad_i *= Limiter_i[iVar];
@@ -256,14 +254,12 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
           }
 
           for (auto iVar = 0u; iVar < nVar; iVar++) {
-            su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
-            su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
+            su2double U_ij = 0.0;
+            if (umuscl)
+              U_ij = Scalar_j[iVar] - Scalar_i[iVar];
 
-            if (umuscl) {
-              const su2double U_ij = Scalar_j[iVar] - Scalar_i[iVar];
-              Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, U_ij, kappa);
-              Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, U_ij, kappa);
-            }
+            su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], U_ij, Vector_ij, umuscl, kappa);
+            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], U_ij, Vector_ij, umuscl, kappa);
 
             if (limiter) {
               Project_Grad_i *= Limiter_i[iVar];

SU2_CFD/include/solvers/CSolver.hpp

@@ -578,6 +578,17 @@ class CSolver {
    */
   inline virtual void SetPrimitive_Limiter(CGeometry *geometry, const CConfig *config) { }
 
+  /*!
+   * \brief Compute the projection of a variable for MUSCL reconstruction.
+   * \param[in] grad - Gradient vector.
+   * \param[in] delta_ij - Centered difference.
+   * \param[in] vector_ij - Distance vector.
+   * \param[in] umuscl - Whether to use a U-MUSCL (blended difference) reconstruction.
+   * \param[in] kappa - Blending coefficient for U-MUSCL reconstruction.
+   * \return - Projected variable.
+   */
+  su2double MUSCL_Reconstruction(su2double* grad, su2double delta_ij, su2double* vector_ij, bool umuscl, su2double kappa);
+
   /*!
    * \brief Set the old solution variables to the current solution value for Runge-Kutta iteration.
    *        It is a virtual function, because for the DG-FEM solver a different version is needed.

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -1881,21 +1881,15 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
       auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint);
 
       for (auto iVar = 0u; iVar < nPrimVarGrad; iVar++) {
-        su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
-        su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
-
         su2double V_ij = 0.0;
         if (umuscl || van_albada)
           V_ij = V_j[iVar] - V_i[iVar];
 
-        if (umuscl) {
-          Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappa);
-          Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappa);
-        }
+        const su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
+        const su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
 
         su2double lim_i = 1.0;
         su2double lim_j = 1.0;
-
         if (van_albada) {
           lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
           lim_j = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_j, V_ij, EPS);

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1279,21 +1279,15 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
       auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint);
 
       for (auto iVar = 0u; iVar < nPrimVarGrad; iVar++) {
-        su2double Project_Grad_i = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
-        su2double Project_Grad_j = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
-
         su2double V_ij = 0.0;
         if (umuscl || van_albada)
           V_ij = V_j[iVar] - V_i[iVar];
 
-        if (umuscl) {
-          Project_Grad_i = LimiterHelpers<>::umusclProjection(Project_Grad_i, V_ij, kappa);
-          Project_Grad_j = LimiterHelpers<>::umusclProjection(Project_Grad_j, V_ij, kappa);
-        }
+        const su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
+        const su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
 
         su2double lim_i = 1.0;
         su2double lim_j = 1.0;
-
         if (van_albada) {
           lim_i = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_i, V_ij, EPS);
           lim_j = LimiterHelpers<>::vanAlbadaFunction(Project_Grad_j, V_ij, EPS);

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -542,17 +542,12 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
       su2double lim_j = 2.0;
 
       for (auto iVar = 0ul; iVar < nPrimVarGrad; iVar++) {
-        Project_Grad_i[iVar] = GeometryToolbox::DotProduct(nDim, Gradient_i[iVar], Vector_ij);
-        Project_Grad_j[iVar] = GeometryToolbox::DotProduct(nDim, Gradient_j[iVar], Vector_ij);
-
         su2double V_ij = 0.0;
         if (umuscl || van_albada)
           V_ij = V_j[iVar] - V_i[iVar];
 
-        if (umuscl) {
-          Project_Grad_i[iVar] = LimiterHelpers<>::umusclProjection(Project_Grad_i[iVar], V_ij, kappa);
-          Project_Grad_j[iVar] = LimiterHelpers<>::umusclProjection(Project_Grad_j[iVar], V_ij, kappa);
-        }
+        Project_Grad_i[iVar] = MUSCL_Reconstruction(Gradient_i[iVar], V_ij, Vector_ij, umuscl, kappa);
+        Project_Grad_j[iVar] = MUSCL_Reconstruction(Gradient_j[iVar], V_ij, Vector_ij, umuscl, kappa);
 
         if (limiter) {
           if (van_albada) {

SU2_CFD/src/solvers/CSolver.cpp

@@ -2279,6 +2279,14 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) {
                   *geometry, *config, 0, nVar, umuscl, umusclKappa, solution, gradient, solMin, solMax, limiter);
 }
 
+su2double CSolver::MUSCL_Reconstruction(su2double* grad, su2double delta_ij, su2double* vector_ij, bool umuscl, su2double kappa) {
+  su2double project_grad = GeometryToolbox::DotProduct(nDim, grad, vector_ij);
+  if (umuscl) {
+    return LimiterHelpers<>::umusclProjection(project_grad, delta_ij, kappa);
+  }
+  return project_grad;
+}
+
 void CSolver::Gauss_Elimination(su2double** A, su2double* rhs, unsigned short nVar) {
 
   short iVar, jVar, kVar;