Merge branch 'develop' into improve_mg

Author: pcarruscag

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/drivers/CDiscAdjSinglezoneDriver.hpp

@@ -112,12 +112,6 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
    */
   void SetAdj_ObjFunction(void);
 
-  /*!
-   * \brief Print out the direct residuals.
-   * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp)
-   */
-  void Print_DirectResidual(RECORDING kind_recording);
-
   /*!
    * \brief Record the main computational path.
    */

SU2_CFD/include/drivers/CDriver.hpp

@@ -313,6 +313,12 @@ class CDriver {
    */
   virtual void Update() {}
 
+  /*!
+   * \brief Print out the direct residuals.
+   * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp)
+   */
+  void Print_DirectResidual(RECORDING kind_recording);
+
 public:
 
   /*!

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/numerics/CNumerics.hpp

@@ -130,8 +130,8 @@ class CNumerics {
   *Psi_i,    /*!< \brief Vector of adjoint variables at point i. */
   *Psi_j;    /*!< \brief Vector of adjoint variables at point j. */
   const su2double
-  *TurbVar_i,   /*!< \brief Vector of turbulent variables at point i. */
-  *TurbVar_j;   /*!< \brief Vector of turbulent variables at point j. */
+  *ScalarVar_i,   /*!< \brief Vector of scalar variables at point i. */
+  *ScalarVar_j;   /*!< \brief Vector of scalar variables at point j. */
   const su2double
   *TransVar_i,  /*!< \brief Vector of turbulent variables at point i. */
   *TransVar_j;  /*!< \brief Vector of turbulent variables at point j. */
@@ -146,8 +146,8 @@ class CNumerics {
   PrimVar_Grad_j,  /*!< \brief Gradient of primitive variables at point j. */
   PsiVar_Grad_i,   /*!< \brief Gradient of adjoint variables at point i. */
   PsiVar_Grad_j,   /*!< \brief Gradient of adjoint variables at point j. */
-  TurbVar_Grad_i,  /*!< \brief Gradient of turbulent variables at point i. */
-  TurbVar_Grad_j,  /*!< \brief Gradient of turbulent variables at point j. */
+  ScalarVar_Grad_i,  /*!< \brief Gradient of scalar variables at point i. */
+  ScalarVar_Grad_j,  /*!< \brief Gradient of scalar variables at point j. */
   TransVar_Grad_i, /*!< \brief Gradient of turbulent variables at point i. */
   TransVar_Grad_j, /*!< \brief Gradient of turbulent variables at point j. */
   TurbPsi_Grad_i,  /*!< \brief Gradient of adjoint turbulent variables at point i. */
@@ -347,13 +347,13 @@ class CNumerics {
   }
 
   /*!
-   * \brief Set the value of the turbulent variable.
-   * \param[in] val_turbvar_i - Value of the turbulent variable at point i.
-   * \param[in] val_turbvar_j - Value of the turbulent variable at point j.
+   * \brief Set the value of the scalar variable.
+   * \param[in] val_scalarvar_i - Value of the scalar variable at point i.
+   * \param[in] val_scalarvar_j - Value of the scalar variable at point j.
    */
-  inline void SetTurbVar(const su2double *val_turbvar_i, const su2double *val_turbvar_j) {
-    TurbVar_i = val_turbvar_i;
-    TurbVar_j = val_turbvar_j;
+  inline void SetScalarVar(const su2double *val_scalarvar_i, const su2double *val_scalarvar_j) {
+    ScalarVar_i = val_scalarvar_i;
+    ScalarVar_j = val_scalarvar_j;
   }
 
   /*!
@@ -367,14 +367,14 @@ class CNumerics {
   }
 
   /*!
-   * \brief Set the gradient of the turbulent variables.
-   * \param[in] val_turbvar_grad_i - Gradient of the turbulent variable at point i.
-   * \param[in] val_turbvar_grad_j - Gradient of the turbulent variable at point j.
+   * \brief Set the gradient of the scalar variables.
+   * \param[in] val_scalarvar_grad_i - Gradient of the scalar variable at point i.
+   * \param[in] val_scalarvar_grad_j - Gradient of the scalar variable at point j.
    */
-  inline void SetTurbVarGradient(CMatrixView<const su2double> val_turbvar_grad_i,
-                                 CMatrixView<const su2double> val_turbvar_grad_j) {
-    TurbVar_Grad_i = val_turbvar_grad_i;
-    TurbVar_Grad_j = val_turbvar_grad_j;
+  inline void SetScalarVarGradient(CMatrixView<const su2double> val_scalarvar_grad_i,
+                                   CMatrixView<const su2double> val_scalarvar_grad_j) {
+    ScalarVar_Grad_i = val_scalarvar_grad_i;
+    ScalarVar_Grad_j = val_scalarvar_grad_j;
   }
 
   /*!

SU2_CFD/include/numerics/scalar/scalar_convection.hpp

@@ -0,0 +1,89 @@
+/*!
+ * \file scalar_convection.hpp
+ * \brief Delarations of numerics classes for discretization of
+ *        convective fluxes in scalar problems.
+ * \author F. Palacios, T. Economon
+ * \version 7.2.0 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2021, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "../CNumerics.hpp"
+
+/*!
+ * \class CUpwScalar
+ * \brief Template class for scalar upwind fluxes between nodes i and j.
+ * \details This class serves as a template for the scalar upwinding residual
+ *   classes.  The general structure of a scalar upwinding calculation is the
+ *   same for many different  models, which leads to a lot of repeated code.
+ *   By using the template design pattern, these sections of repeated code are
+ *   moved to this shared base class, and the specifics of each model
+ *   are implemented by derived classes.  In order to add a new residual
+ *   calculation for a convection residual, extend this class and implement
+ *   the pure virtual functions with model-specific behavior.
+ * \ingroup ConvDiscr
+ * \author C. Pederson, A. Bueno., and A. Campos.
+ */
+class CUpwScalar : public CNumerics {
+ protected:
+  su2double a0 = 0.0;               /*!< \brief The maximum of the face-normal velocity and 0 */
+  su2double a1 = 0.0;               /*!< \brief The minimum of the face-normal velocity and 0 */
+  su2double* Flux = nullptr;        /*!< \brief Final result, diffusive flux/residual. */
+  su2double** Jacobian_i = nullptr; /*!< \brief Flux Jacobian w.r.t. node i. */
+  su2double** Jacobian_j = nullptr; /*!< \brief Flux Jacobian w.r.t. node j. */
+
+  const bool implicit = false, incompressible = false, dynamic_grid = false;
+
+  /*!
+   * \brief A pure virtual function; Adds any extra variables to AD
+   */
+  virtual void ExtraADPreaccIn() = 0;
+
+  /*!
+   * \brief Model-specific steps in the ComputeResidual method, derived classes
+   *        compute the Flux and its Jacobians via this method.
+   * \param[in] config - Definition of the particular problem.
+   */
+  virtual void FinishResidualCalc(const CConfig* config) = 0;
+
+ public:
+  /*!
+   * \brief Constructor of the class.
+   * \param[in] val_nDim - Number of dimensions of the problem.
+   * \param[in] val_nVar - Number of variables of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  CUpwScalar(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
+
+  /*!
+   * \brief Destructor of the class.
+   */
+  ~CUpwScalar(void) override;
+
+  /*!
+   * \brief Compute the scalar upwind flux between two nodes i and j.
+   * \param[in] config - Definition of the particular problem.
+   * \return A lightweight const-view (read-only) of the residual/flux and Jacobians.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override;
+};