Intermediate class of variable for all flow solver variables + Fix for MUSCL in NEMO solver

Common/include/code_config.hpp

@@ -67,6 +67,17 @@ template<class T, class F> struct su2conditional<false, T, F> { using type = F;
 template<bool B, class T, class F>
 using su2conditional_t = typename su2conditional<B,T,F>::type;
 
+/*! \brief Static cast "In" to "Out", in debug builds a dynamic cast is used. */
+template<class Out, class In>
+FORCEINLINE Out su2staticcast_p(In ptr) {
+  static_assert(std::is_pointer<In>::value, "This expects a pointer");
+#ifndef NDEBUG
+  return static_cast<Out>(ptr);
+#else
+  return dynamic_cast<Out>(ptr);
+#endif
+}
+
 /*--- Detect compilation with OpenMP. ---*/
 #if defined(_OPENMP)
 #define HAVE_OMP

SU2_CFD/include/limiters/CLimiterDetails.hpp

@@ -64,22 +64,28 @@ struct CLimiterDetails
 /*!
  * \brief Common small functions used by limiters.
  */
-namespace LimiterHelpers
+template<class Type = su2double>
+struct LimiterHelpers
 {
-  inline passivedouble epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
+  FORCEINLINE static Type epsilon() {return std::numeric_limits<passivedouble>::epsilon();}
 
-  inline su2double venkatFunction(su2double proj, su2double delta, su2double eps2)
+  FORCEINLINE static Type venkatFunction(const Type& proj, const Type& delta, const Type& eps2)
   {
-    su2double y = delta*(delta+proj) + eps2;
+    Type y = delta*(delta+proj) + eps2;
     return (y + delta*proj) / (y + 2*proj*proj);
   }
 
-  inline su2double raisedSine(su2double dist)
+  FORCEINLINE static Type vanAlbadaFunction(const Type& proj, const Type& delta, const Type& eps)
   {
-    su2double factor = 0.5*(1.0+dist+sin(PI_NUMBER*dist)/PI_NUMBER);
+    return delta * (2*proj + delta) / (4*pow(proj, 2) + pow(delta, 2) + eps);
+  }
+
+  FORCEINLINE static Type raisedSine(const Type& dist)
+  {
+    Type factor = 0.5*(1.0+dist+sin(PI_NUMBER*dist)/PI_NUMBER);
     return max(0.0, min(factor, 1.0));
   }
-}
+};
 
 
 /*!
@@ -94,7 +100,7 @@ struct CLimiterDetails<BARTH_JESPERSEN>
    * \brief Set a small epsilon to avoid divisions by 0.
    */
   template<class... Ts>
-  inline void preprocess(Ts&...) {eps2 = LimiterHelpers::epsilon();}
+  inline void preprocess(Ts&...) {eps2 = LimiterHelpers<>::epsilon();}
 
   /*!
    * \brief No geometric modification for this kind of limiter.
@@ -107,7 +113,7 @@ struct CLimiterDetails<BARTH_JESPERSEN>
    */
   inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
   {
-    return LimiterHelpers::venkatFunction(proj, delta, eps2);
+    return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
   }
 };
 
@@ -130,7 +136,7 @@ struct CLimiterDetails<VENKATAKRISHNAN>
     su2double L = config.GetRefElemLength();
     su2double K = config.GetVenkat_LimiterCoeff();
     su2double eps1 = fabs(L*K);
-    eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
+    eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
   }
 
   /*!
@@ -144,7 +150,7 @@ struct CLimiterDetails<VENKATAKRISHNAN>
    */
   inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
   {
-    return LimiterHelpers::venkatFunction(proj, delta, eps2);
+    return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
   }
 };
 
@@ -229,7 +235,7 @@ struct CLimiterDetails<VENKATAKRISHNAN_WANG>
     for(size_t iVar = varBegin; iVar < varEnd; ++iVar)
     {
       su2double range = sharedMax(iVar) - sharedMin(iVar);
-      eps2(iVar) = max(pow(K*range, 2), LimiterHelpers::epsilon());
+      eps2(iVar) = max(pow(K*range, 2), LimiterHelpers<>::epsilon());
     }
   }
 
@@ -245,7 +251,7 @@ struct CLimiterDetails<VENKATAKRISHNAN_WANG>
   inline su2double limiterFunction(size_t iVar, su2double proj, su2double delta) const
   {
     AD::SetPreaccIn(eps2(iVar));
-    return LimiterHelpers::venkatFunction(proj, delta, eps2(iVar));
+    return LimiterHelpers<>::venkatFunction(proj, delta, eps2(iVar));
   }
 };
 
@@ -268,7 +274,7 @@ struct CLimiterDetails<SHARP_EDGES>
     su2double L = config.GetRefElemLength();
     su2double K = config.GetVenkat_LimiterCoeff();
     eps1 = fabs(L*K);
-    eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
+    eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
   }
 
   /*!
@@ -278,15 +284,15 @@ struct CLimiterDetails<SHARP_EDGES>
   {
     AD::SetPreaccIn(geometry.nodes->GetSharpEdge_Distance(iPoint));
     su2double dist = geometry.nodes->GetSharpEdge_Distance(iPoint)/(sharpCoeff*eps1)-1.0;
-    return LimiterHelpers::raisedSine(dist);
+    return LimiterHelpers<>::raisedSine(dist);
   }
 
   /*!
    * \brief Smooth function that disables limiting in smooth regions.
    */
   inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
   {
-    return LimiterHelpers::venkatFunction(proj, delta, eps2);
+    return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
   }
 };
 
@@ -309,7 +315,7 @@ struct CLimiterDetails<WALL_DISTANCE>
     su2double L = config.GetRefElemLength();
     su2double K = config.GetVenkat_LimiterCoeff();
     eps1 = fabs(L*K);
-    eps2 = max(eps1*eps1*eps1, LimiterHelpers::epsilon());
+    eps2 = max(eps1*eps1*eps1, LimiterHelpers<>::epsilon());
   }
 
   /*!
@@ -319,14 +325,14 @@ struct CLimiterDetails<WALL_DISTANCE>
   {
     AD::SetPreaccIn(geometry.nodes->GetWall_Distance(iPoint));
     su2double dist = geometry.nodes->GetWall_Distance(iPoint)/(sharpCoeff*eps1)-1.0;
-    return LimiterHelpers::raisedSine(dist);
+    return LimiterHelpers<>::raisedSine(dist);
   }
 
   /*!
    * \brief Smooth function that disables limiting in smooth regions.
    */
   inline su2double limiterFunction(size_t, su2double proj, su2double delta) const
   {
-    return LimiterHelpers::venkatFunction(proj, delta, eps2);
+    return LimiterHelpers<>::venkatFunction(proj, delta, eps2);
   }
 };

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -684,12 +684,6 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
                          CConfig *config,
                          unsigned short val_marker) final;
 
-  /*!
-   * \brief Set the new solution variables to the current solution value for classical RK.
-   * \param[in] geometry - Geometrical definition of the problem.
-   */
-  inline void Set_NewSolution() final { nodes->SetSolution_New(); }
-
   /*!
    * \brief Update the solution using a Runge-Kutta scheme.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -1102,6 +1102,11 @@ class CFVMFlowSolverBase : public CSolver {
 
  public:
 
+  /*!
+   * \brief Set the new solution variables to the current solution value for classical RK.
+   */
+  inline void Set_NewSolution() final { nodes->SetSolution_New(); }
+
   /*!
    * \brief Load a solution from a restart file.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -27,6 +27,7 @@
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../../include/solvers/CScalarSolver.hpp"
+#include "../../include/variables/CFlowVariable.hpp"
 
 template <class VariableType>
 CScalarSolver<VariableType>::CScalarSolver(bool conservative) : CSolver(), Conservative(conservative) {}
@@ -92,10 +93,10 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
   const bool limiterFlow =
       (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (config->GetKind_SlopeLimit_Flow() != VAN_ALBADA_EDGE);
 
-  CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
+  auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
   /*--- Pick one numerics object per thread. ---*/
-  CNumerics* numerics = numerics_container[CONV_TERM + omp_get_thread_num() * MAX_TERMS];
+  auto* numerics = numerics_container[CONV_TERM + omp_get_thread_num() * MAX_TERMS];
 
   /*--- Static arrays of MUSCL-reconstructed flow primitives and turbulence variables (thread safety). ---*/
   su2double solution_i[MAXNVAR] = {0.0}, flowPrimVar_i[MAXNVARFLOW] = {0.0};