Merge pull request #1178 from su2code/hybrid_parallel_incomp

Hybrid parallel (OpenMP) support for incompressible solvers

Author: pcarruscag

SU2_CFD/include/solvers/CAdjEulerSolver.hpp

@@ -230,13 +230,6 @@ class CAdjEulerSolver : public CSolver {
                        CConfig *config,
                        unsigned short iMesh) final;
 
-  /*!
-   * \brief Compute the undivided laplacian for the adjoint solution.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetUndivided_Laplacian(CGeometry *geometry, CConfig *config);
-
   /*!
    * \brief Value of the characteristic variables at the boundaries.
    * \param[in] val_marker - Surface marker where the coefficient is computed.

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -37,6 +37,8 @@
  */
 class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
 protected:
+  using BaseClass = CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE>;
+
   su2double
   Prandtl_Lam = 0.0,    /*!< \brief Laminar Prandtl number. */
   Prandtl_Turb = 0.0;   /*!< \brief Turbulent Prandtl number. */
@@ -350,20 +352,6 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
                        CConfig *config,
                        unsigned short iMesh) final;
 
-  /*!
-   * \brief Compute the viscous contribution for a particular edge.
-   * \note The convective residual methods include a call to this for each edge,
-   *       this allows convective and viscous loops to be "fused".
-   * \param[in] iEdge - Edge for which the flux and Jacobians are to be computed.
-   * \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.
-   */
-  inline virtual void Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container,
-                                       CNumerics *numerics, CConfig *config) { }
-  using CSolver::Viscous_Residual; /*--- Silence warning ---*/
-
   /*!
    * \brief Recompute the extrapolated quantities, after MUSCL reconstruction,
    *        in a more thermodynamically consistent way.
@@ -1071,20 +1059,6 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
    */
   void UpdateCustomBoundaryConditions(CGeometry **geometry_container, CConfig *config) final;
 
-  /*!
-   * \brief Load a solution from a restart file.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - Container vector with all of the solvers.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] val_iter - Current external iteration number.
-   * \param[in] val_update_geo - Flag for updating coords and grid velocity.
-   */
-  void LoadRestart(CGeometry **geometry,
-                   CSolver ***solver,
-                   CConfig *config,
-                   int val_iter,
-                   bool val_update_geo) final;
-
   /*!
    * \brief Set the initial condition for the Euler Equations.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -1627,7 +1601,7 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
   void PrintVerificationError(const CConfig* config) const final;
 
   /*!
-   * \brief The Euler and NS solvers support MPI+OpenMP (except the BC bits).
+   * \brief The Euler and NS solvers support MPI+OpenMP.
    */
   inline bool GetHasHybridParallel() const final { return true; }
 

SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp

@@ -92,6 +92,7 @@ class CFEM_DG_EulerSolver : public CSolver {
   AllBound_CEff_Inv;  /*!< \brief Total efficiency (Cl/Cd) (inviscid contribution) for all the boundaries. */
 
   su2double
+  AeroCoeffForceRef,   /*!< \brief Reference force for coefficients */
   Total_CL,      /*!< \brief Total lift coefficient for all the boundaries. */
   Total_CD,        /*!< \brief Total drag coefficient for all the boundaries. */
   Total_CSF,       /*!< \brief Total sideforce coefficient for all the boundaries. */
@@ -1127,6 +1128,11 @@ class CFEM_DG_EulerSolver : public CSolver {
    */
   inline void SetTotal_CL(su2double val_Total_CL) final { Total_CL = val_Total_CL; }
 
+  /*!
+   * \brief Get the reference force used to compute CL, CD, etc.
+   */
+  inline su2double GetAeroCoeffsReferenceForce() const final { return AeroCoeffForceRef; }
+
   /*!
    * \brief Provide the total (inviscid + viscous) non dimensional lift coefficient.
    * \return Value of the lift coefficient (inviscid + viscous contribution).

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -125,6 +125,8 @@ class CFVMFlowSolverBase : public CSolver {
   AeroCoeffsArray SurfaceCoeff; /*!< \brief Totals for each monitoring surface. */
   AeroCoeffs TotalCoeff;        /*!< \brief Totals for all boundaries. */
 
+  su2double AeroCoeffForceRef = 1.0;    /*!< \brief Reference force for aerodynamic coefficients. */
+
   su2double InverseDesign = 0.0;        /*!< \brief Inverse design functional for each boundary. */
   su2double Total_ComboObj = 0.0;       /*!< \brief Total 'combo' objective for all monitored boundaries */
   su2double Total_Custom_ObjFunc = 0.0; /*!< \brief Total custom objective function for all the boundaries. */
@@ -249,6 +251,36 @@ class CFVMFlowSolverBase : public CSolver {
    */
   inline virtual void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) {}
 
+  /*!
+   * \brief Compute the viscous contribution for a particular edge.
+   * \note The convective residual methods include a call to this for each edge,
+   *       this allows convective and viscous loops to be "fused".
+   * \param[in] iEdge - Edge for which the flux and Jacobians are to be computed.
+   * \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.
+   */
+  inline virtual void Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container,
+                                       CNumerics *numerics, CConfig *config) { }
+  void Viscous_Residual_impl(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container,
+                             CNumerics *numerics, CConfig *config);
+  using CSolver::Viscous_Residual; /*--- Silence warning ---*/
+
+  /*!
+   * \brief General implementation to load a flow solution from a restart file.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver - Container vector with all of the solvers.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iter - Current external iteration number.
+   * \param[in] update_geo - Flag for updating coords and grid velocity.
+   * \param[in] RestartSolution - Optional buffer to load restart vars into,
+   *            this allows default values to be given when nVar > nVar_Restart.
+   * \param[in] nVar_Restart - Number of restart variables, if 0 defaults to nVar.
+   */
+  void LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo,
+                        su2double* RestartSolution = nullptr, unsigned short nVar_Restart = 0);
+
   /*!
    * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -949,36 +981,90 @@ class CFVMFlowSolverBase : public CSolver {
 
   /*!
    * \brief Evaluate the vorticity and strain rate magnitude.
+   * \tparam VelocityOffset: Index in the primitive variables where the velocity starts.
    */
-  inline void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) {
+  template<size_t VelocityOffset>
+  void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) {
+
+    const auto& Gradient_Primitive = nodes->GetGradient_Primitive();
+    auto& StrainMag = nodes->GetStrainMag();
 
     SU2_OMP_MASTER {
       StrainMag_Max = 0.0;
       Omega_Max = 0.0;
     }
     SU2_OMP_BARRIER
 
-    nodes->SetVorticity_StrainMag();
+    su2double strainMax = 0.0, omegaMax = 0.0;
 
-    /*--- Min and Max are not really differentiable ---*/
-    const bool wasActive = AD::BeginPassive();
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
 
-    su2double strainMax = 0.0, omegaMax = 0.0;
+      constexpr size_t u = VelocityOffset;
+      constexpr size_t v = VelocityOffset+1;
+      constexpr size_t w = VelocityOffset+2;
 
-    SU2_OMP(for schedule(static,omp_chunk_size) nowait)
-    for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
-      strainMax = max(strainMax, nodes->GetStrainMag(iPoint));
-      omegaMax = max(omegaMax, GeometryToolbox::Norm(3, nodes->GetVorticity(iPoint)));
-    }
-    SU2_OMP_CRITICAL {
-      StrainMag_Max = max(StrainMag_Max, strainMax);
-      Omega_Max = max(Omega_Max, omegaMax);
+      /*--- Vorticity ---*/
+
+      su2double* Vorticity = nodes->GetVorticity(iPoint);
+
+      Vorticity[0] = 0.0; Vorticity[1] = 0.0;
+
+      Vorticity[2] = Gradient_Primitive(iPoint,v,0)-Gradient_Primitive(iPoint,u,1);
+
+      if (nDim == 3) {
+        Vorticity[0] = Gradient_Primitive(iPoint,w,1)-Gradient_Primitive(iPoint,v,2);
+        Vorticity[1] = -(Gradient_Primitive(iPoint,w,0)-Gradient_Primitive(iPoint,u,2));
+      }
+
+      /*--- Strain Magnitude ---*/
+
+      AD::StartPreacc();
+      AD::SetPreaccIn(&Gradient_Primitive[iPoint][VelocityOffset], nDim, nDim);
+
+      su2double Div = 0.0;
+      for (unsigned long iDim = 0; iDim < nDim; iDim++)
+        Div += Gradient_Primitive(iPoint, iDim+VelocityOffset, iDim);
+      Div /= 3.0;
+
+      StrainMag(iPoint) = 0.0;
+
+      /*--- Add diagonal part ---*/
+
+      for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+        StrainMag(iPoint) += pow(Gradient_Primitive(iPoint, iDim+VelocityOffset, iDim) - Div, 2);
+      }
+      if (nDim == 2) {
+        StrainMag(iPoint) += pow(Div, 2);
+      }
+
+      /*--- Add off diagonals ---*/
+
+      StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,u,1) + Gradient_Primitive(iPoint,v,0)), 2);
+
+      if (nDim == 3) {
+        StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,u,2) + Gradient_Primitive(iPoint,w,0)), 2);
+        StrainMag(iPoint) += 2.0*pow(0.5*(Gradient_Primitive(iPoint,v,2) + Gradient_Primitive(iPoint,w,1)), 2);
+      }
+
+      StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint));
+      AD::SetPreaccOut(StrainMag(iPoint));
+
+      /*--- Max is not differentiable, so we not register them for preacc. ---*/
+      strainMax = max(strainMax, StrainMag(iPoint));
+      omegaMax = max(omegaMax, GeometryToolbox::Norm(3, Vorticity));
+
+      AD::EndPreacc();
     }
 
     if ((iMesh == MESH_0) && (config.GetComm_Level() == COMM_FULL)) {
+      SU2_OMP_CRITICAL {
+        StrainMag_Max = max(StrainMag_Max, strainMax);
+        Omega_Max = max(Omega_Max, omegaMax);
+      }
+
       SU2_OMP_BARRIER
-      SU2_OMP_MASTER
-      {
+      SU2_OMP_MASTER {
         su2double MyOmega_Max = Omega_Max;
         su2double MyStrainMag_Max = StrainMag_Max;
 
@@ -988,7 +1074,6 @@ class CFVMFlowSolverBase : public CSolver {
       SU2_OMP_BARRIER
     }
 
-    AD::EndPassive(wasActive);
   }
 
   /*!
@@ -997,6 +1082,26 @@ class CFVMFlowSolverBase : public CSolver {
   ~CFVMFlowSolverBase();
 
  public:
+
+  /*!
+   * \brief Load a solution from a restart file.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver - Container vector with all of the solvers.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iter - Current external iteration number.
+   * \param[in] update_geo - Flag for updating coords and grid velocity.
+   */
+  void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo) override;
+
+  /*!
+   * \brief Set the initial condition for the Euler Equations.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] ExtIter - External iteration.
+   */
+  void SetInitialCondition(CGeometry **geometry, CSolver ***solver_container, CConfig *config, unsigned long ExtIter) override;
+
   /*!
    * \brief Compute the gradient of the primitive variables using Green-Gauss method,
    *        and stores the result in the <i>Gradient_Primitive</i> variable.
@@ -1490,6 +1595,11 @@ class CFVMFlowSolverBase : public CSolver {
    */
   inline su2double GetTotal_CEff() const final { return TotalCoeff.CEff; }
 
+  /*!
+   * \brief Get the reference force used to compute CL, CD, etc.
+   */
+  inline su2double GetAeroCoeffsReferenceForce() const final { return AeroCoeffForceRef; }
+
   /*!
    * \brief Provide the total (inviscid + viscous) non dimensional lift coefficient.
    * \return Value of the lift coefficient (inviscid + viscous contribution).