Merge pull request #1347 from su2code/nemo_visc_update

Updating some of the NEMO viscous solver routines.

Author: WallyMaier

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -1007,8 +1007,8 @@ class CFVMFlowSolverBase : public CSolver {
    * \brief Evaluate the vorticity and strain rate magnitude.
    * \tparam VelocityOffset - Index in the primitive variables where the velocity starts.
    */
-  template<unsigned long VelocityOffset>
-  void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) {
+  void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh,
+                                    const size_t VelocityOffset = 1) {
 
     auto& StrainMag = nodes->GetStrainMag();
 

SU2_CFD/include/solvers/CNEMONSSolver.hpp

@@ -38,14 +38,10 @@
  * \brief Main class for defining the NEMO Navier-Stokes flow solver.
  * \ingroup Navier_Stokes_Equations
  * \author S. R. Copeland, F. Palacios, W. Maier.
- *
  */
 class CNEMONSSolver final : public CNEMOEulerSolver {
 private:
 
-  su2double Prandtl_Lam,     /*!< \brief Laminar Prandtl number. */
-  Prandtl_Turb;              /*!< \brief Turbulent Prandtl number. */
-
   /*!
    * \brief Compute the velocity^2, SoundSpeed, Pressure, Enthalpy, Viscosity.
    * \param[in] solver_container - Container vector with all the solutions.
@@ -55,6 +51,27 @@ class CNEMONSSolver final : public CNEMOEulerSolver {
    */
   unsigned long SetPrimitive_Variables(CSolver **solver_container,
                                        CConfig *config, bool Output) override;
+
+  /*!
+   * \brief Compute the viscous residuals.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics - Description of the numerical method.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iMesh - Index of the mesh in multigrid computations.
+   * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
+   */
+  void Viscous_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
+                        CConfig *config, unsigned short iMesh, unsigned short iRKStep) override;
+
+  /*!
+   * \brief Computes the wall shear stress (Tau_Wall) on the surface using a wall function.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetTauWall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config);
+
 public:
 
   /*!
@@ -198,17 +215,4 @@ class CNEMONSSolver final : public CNEMOEulerSolver {
                                CNumerics *visc_numerics,
                                CConfig *config,
                                unsigned short val_marker) override;
-
-  /*!
-   * \brief Compute the viscous residuals.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Container vector with all the solutions.
-   * \param[in] numerics - Description of the numerical method.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iMesh - Index of the mesh in multigrid computations.
-   * \param[in] iRKStep - Current step of the Runge-Kutta iteration.
-   */
-  void Viscous_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
-                        CConfig *config, unsigned short iMesh, unsigned short iRKStep) override;
-
 };

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -65,6 +65,10 @@ class CNEMOEulerVariable : public CVariable {
   /*--- New solution container for Classical RK4 ---*/
   MatrixType Solution_New;  /*!< \brief New solution container for Classical RK4. */
 
+  /*--- NS Variables declared here to make it easier to re-use code between compressible and incompressible solvers. ---*/
+  MatrixType Vorticity;       /*!< \brief Vorticity of the fluid. */
+  VectorType StrainMag;       /*!< \brief Magnitude of rate of strain tensor. */
+
   /*--- Other Necessary Variable Definition ---*/
   MatrixType dPdU;   /*!< \brief Partial derivative of pressure w.r.t. conserved variables. */
   MatrixType dTdU;   /*!< \brief Partial derivative of temperature w.r.t. conserved variables. */
@@ -536,6 +540,19 @@ class CNEMOEulerVariable : public CVariable {
       Solution(iPoint, nSpecies+iDim) = Primitive(iPoint,RHO_INDEX) * val_vector[iDim];
   }
 
+  /*!
+   * \brief Get the value of the vorticity.
+   * \return Value of the vorticity.
+   */
+  inline su2double *GetVorticity(unsigned long iPoint) final { return Vorticity[iPoint]; }
+
+  /*!
+   * \brief Get the value of the magnitude of rate of strain.
+   * \return Value of the rate of strain magnitude.
+   */
+  inline su2double GetStrainMag(unsigned long iPoint) const final { return StrainMag(iPoint); }
+  inline su2activevector& GetStrainMag() { return StrainMag; }
+
   /*!
    * \brief Set the momentum part of the truncation error to zero.
    * \param[in] iPoint - Point index.

SU2_CFD/include/variables/CNEMONSVariable.hpp

@@ -51,8 +51,6 @@ class CNEMONSVariable final : public CNEMOEulerVariable {
 
   su2double inv_TimeScale;      /*!< \brief Inverse of the reference time scale. */
 
-  MatrixType Vorticity;         /*!< \brief Vorticity of the fluid. */
-  VectorType StrainMag;         /*!< \brief Magnitude of rate of strain tensor. */
   VectorType Tau_Wall;          /*!< \brief Magnitude of the wall shear stress from a wall function. */
   VectorType DES_LengthScale;   /*!< \brief DES Length Scale. */
   VectorType Roe_Dissipation;   /*!< \brief Roe low dissipation coefficient. */
@@ -155,12 +153,4 @@ class CNEMONSVariable final : public CNEMOEulerVariable {
   inline void SetWallTemperature(unsigned long iPoint, su2double temperature_wall) override {
     Primitive(iPoint,T_INDEX) = temperature_wall;
   }
-
-  /*!
-   * \brief Get the value of the vorticity.
-   * \return Value of the vorticity.
-   */
-  inline su2double *GetVorticity(unsigned long iPoint) override { return Vorticity[iPoint]; }
-
-
 };

SU2_CFD/src/fluid/CNEMOGas.cpp

@@ -240,12 +240,10 @@ void CNEMOGas::ComputedPdU(su2double *V, vector<su2double>& val_eves, su2double
 
 void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
 
-  su2double Vel[3] = {0.0};
-
   /*--- Necessary indexes to assess primitive variables ---*/
-  unsigned long T_INDEX       = nSpecies;
-  unsigned long VEL_INDEX     = nSpecies+2;
-  unsigned long RHOCVTR_INDEX = nSpecies+nDim+6;
+  const unsigned long T_INDEX       = nSpecies;
+  const unsigned long VEL_INDEX     = nSpecies+2;
+  const unsigned long RHOCVTR_INDEX = nSpecies+nDim+6;
 
   /*--- Rename for convenience ---*/
   T                 = V[T_INDEX];
@@ -256,9 +254,7 @@ void CNEMOGas::ComputedTdU(su2double *V, su2double *val_dTdU){
   Ref_Temperature    = GetRefTemperature();
 
   /*--- Calculate supporting quantities ---*/
-  for (iDim = 0; iDim < nDim; iDim++)
-    Vel[iDim] = V[VEL_INDEX+iDim]*V[VEL_INDEX+iDim];
-  su2double v2 = GeometryToolbox::SquaredNorm(nDim,Vel);
+  const su2double v2 = GeometryToolbox::SquaredNorm(nDim, &V[VEL_INDEX]);
 
   /*--- Species density derivatives ---*/
   for (iSpecies = nEl; iSpecies < nSpecies; iSpecies++) {

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -96,7 +96,7 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container
     SetPrimitive_Limiter(geometry, config);
   }
 
-  ComputeVorticityAndStrainMag<1>(*config, iMesh);
+  ComputeVorticityAndStrainMag(*config, iMesh);
 
   /*--- Compute the TauWall from the wall functions ---*/
 

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -46,17 +46,20 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
     description = "Euler";
   }
 
+  const auto nZone = geometry->GetnZone();
+  const bool restart = (config->GetRestart() || config->GetRestart_Flow());
+  const auto direct_diff = config->GetDirectDiff();
+  const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
+                         (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
+  const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
+  const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint();
+
+  int Unst_RestartIter = 0;
   unsigned short iMarker, nLineLets;
-  unsigned short nZone = geometry->GetnZone();
   su2double *Mvec_Inf, Alpha, Beta;
-  bool restart   = (config->GetRestart() || config->GetRestart_Flow());
-  unsigned short direct_diff = config->GetDirectDiff();
-  int Unst_RestartIter = 0;
-  bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
-                    (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
-  bool time_stepping = config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING;
-  bool adjoint = config->GetDiscrete_Adjoint();
-  string filename_ = "flow";
+
+  /*--- A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain ---*/
+  dynamic_grid = config->GetDynamic_Grid();
 
   /*--- Store the multigrid level. ---*/
   MGLevel = iMesh;
@@ -79,6 +82,7 @@ CNEMOEulerSolver::CNEMOEulerSolver(CGeometry *geometry, CConfig *config,
       else Unst_RestartIter = SU2_TYPE::Int(config->GetRestart_Iter())-1;
     }
 
+    string filename_ = "flow";
     filename_ = config->GetFilename(filename_, ".meta", Unst_RestartIter);
 
     /*--- Read and store the restart metadata ---*/
@@ -246,18 +250,24 @@ void CNEMOEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver
   bool center_jst_ke    = (config->GetKind_Centered_Flow() == JST_KE) && (iMesh == MESH_0);
 
   /*--- Set the primitive variables ---*/
-  ErrorCounter  = 0;
+  ompMasterAssignBarrier(ErrorCounter,0);
+  SU2_OMP_ATOMIC
   ErrorCounter += SetPrimitive_Variables(solver_container, config, Output);
+  SU2_OMP_BARRIER
 
-  if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) {
+  SU2_OMP_MASTER { /*--- Ops that are not OpenMP parallel go in this block. ---*/
 
+    if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) {
       unsigned long tmp = ErrorCounter;
       SU2_MPI::Allreduce(&tmp, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       config->SetNonphysical_Points(ErrorCounter);
 
       if ((rank == MASTER_NODE) && (ErrorCounter != 0))
         cout << "Warning. The initial solution contains "<< ErrorCounter << " points that are not physical." << endl;
+    }
   }
+  END_SU2_OMP_MASTER
+  SU2_OMP_BARRIER
 
   /*--- Artificial dissipation ---*/
 
@@ -986,7 +996,7 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short
   Gas_ConstantND         = 0.0, Viscosity_FreeStreamND    = 0.0, sqvel                       = 0.0,
   Tke_FreeStreamND       = 0.0, Total_UnstTimeND          = 0.0, Delta_UnstTimeND            = 0.0,
   soundspeed             = 0.0, GasConstant_Inf           = 0.0, Froude                      = 0.0,
-  Density_FreeStreamND   = 0.0;
+  Density_FreeStreamND   = 0.0, Heat_Flux_Ref             = 0.0;
 
   su2double Velocity_FreeStreamND[3] = {0.0, 0.0, 0.0};
 
@@ -1153,6 +1163,7 @@ void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short
   Time_Ref          = Length_Ref/Velocity_Ref;                                     config->SetTime_Ref(Time_Ref);
   Omega_Ref         = Velocity_Ref/Length_Ref;                                     config->SetOmega_Ref(Omega_Ref);
   Force_Ref         = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref; config->SetForce_Ref(Force_Ref);
+  Heat_Flux_Ref     = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref;           config->SetHeat_Flux_Ref(Heat_Flux_Ref);
   Gas_Constant_Ref  = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref();      config->SetGas_Constant_Ref(Gas_Constant_Ref);
   Viscosity_Ref     = config->GetDensity_Ref()*Velocity_Ref*Length_Ref;            config->SetViscosity_Ref(Viscosity_Ref);
   Conductivity_Ref  = Viscosity_Ref*Gas_Constant_Ref;                              config->SetConductivity_Ref(Conductivity_Ref);
@@ -1640,12 +1651,10 @@ void CNEMOEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contai
 
 void CNEMOEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
                                 CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) {
-
-
-
   SU2_MPI::Error("BC_INLET: Not operational in NEMO.", CURRENT_FUNCTION);
 
   unsigned short iVar, iDim, iSpecies, RHO_INDEX, nSpecies;
+
   unsigned long iVertex, iPoint;
   su2double  T_Total, P_Total, Velocity[3], Velocity2, H_Total, Temperature, Riemann,
   Pressure, Density, Energy, Mach2, SoundSpeed2, SoundSpeed_Total2, Vel_Mag,