Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hybrid parallel (OpenMP) support for incompressible solvers #1178

Merged
merged 17 commits into from
Jan 30, 2021
Merged

Hybrid parallel (OpenMP) support for incompressible solvers #1178

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
9c958b7
Merge branch 'stress_penalty' into hybrid_parallel_incomp
pcarruscag Jan 24, 2021
203369d
BC's, fluid model, restart
pcarruscag Jan 24, 2021
f6bc631
Vorticity/StrainMag and UndivLaplacian
pcarruscag Jan 24, 2021
d858f70
fluxes and sources
pcarruscag Jan 24, 2021
bbfd9c0
BetaInc2
pcarruscag Jan 24, 2021
fc4cfd1
add incompressible hybrid regressions
pcarruscag Jan 24, 2021
cc2946e
Merge remote-tracking branch 'upstream/develop' into hybrid_parallel_…
pcarruscag Jan 25, 2021
ce04458
Merge branch 'stress_penalty' into hybrid_parallel_incomp
pcarruscag Jan 26, 2021
5bad047
fix #1175
pcarruscag Jan 26, 2021
dec8abf
Merge branch 'develop' into hybrid_parallel_incomp
pcarruscag Jan 27, 2021
396e8d5
Merge remote-tracking branch 'upstream/develop' into hybrid_parallel_…
pcarruscag Jan 27, 2021
7bd274a
fix leak set incomp hybrid regression residuals
pcarruscag Jan 27, 2021
6118dff
fix the Roe kappa thing and some compiler warnings
pcarruscag Jan 27, 2021
30184d3
fix directdiff compilation with OpenMP
pcarruscag Jan 27, 2021
726f4e7
one LoadRestart to rule them all
pcarruscag Jan 28, 2021
afe4928
simplify SetInitialCondition
pcarruscag Jan 28, 2021
a7e6695
even cleaner
pcarruscag Jan 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 0 additions & 7 deletions SU2_CFD/include/solvers/CAdjEulerSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
32 changes: 3 additions & 29 deletions SU2_CFD/include/solvers/CEulerSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@
*/
class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
protected:
using BaseClass = CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE>;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need this to call an implementation of the base class later like BaseClass::SetInitialCondition in this very CEulerSolver. But SetInitialCondition is also public in the FVMBase class so it should be callable anyway because we inherit that via class CEulerSolver : public CFVMFlowSolverBase ... where am I wrong?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you call SetInitialCondition in SetInitialCondition you get an infinite loop, I just wanted to make more expressive I was calling the implementation of the parent first.


su2double
Prandtl_Lam = 0.0, /*!< \brief Laminar Prandtl number. */
Prandtl_Turb = 0.0; /*!< \brief Turbulent Prandtl number. */
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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; }

Expand Down
6 changes: 6 additions & 0 deletions SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down Expand Up @@ -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).
Expand Down
142 changes: 126 additions & 16 deletions SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down Expand Up @@ -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".
Comment on lines +256 to +257
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a speed-up due to this? Or is this just for simplification?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fast code goes brrrrrrrrrrrrrrrrrrrr

* \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 ---*/
Comment on lines +264 to +268
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea behind this to be able to 'fake' override Viscous_Residual for viscous flows with just calling Viscous_Residual_impl because for inviscid flow the empty implementation from here is taken. This trick has to be used because FVMBase -> Euler -> NS is the inheritance and Viscous_Residual itself cannot be overloaded because Euler is the middleman?
And what happens to the call in CIntegration to Viscous_Residual for NS flow. Is there a doubled contribution then because it is called there and in the centered/upwind residual?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, in Cintegration the empty impl of CSolver is used and in the centered/upwind residual the correct Viscous_Residual impl with the Viscous_Residual_impl init ...
Seems that i have to go back and read a couple of pages in 'Inheritance for complete morons'. fml

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep that's all correct.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Except the self flagellation part, no need for that xD


/*!
* \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.
Expand Down Expand Up @@ -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);
}

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose you do the reduction over the threads yourself because of efficiency, when the if statement is false, but wouldn't it be more readable to have a reduction clause in the OMP for loop?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for having a look at the code Edwin, the issues with these kind of reductions were:

  • The microsoft compilers do not support max/min reductions;
  • For the code to compile with AD (forward) we would need to define custom reduction operators, which might not be very readable;
  • Earlier discussions about OpenMP+AD (reverse) suggested that avoiding reduction clauses would make it easier to develop a thread-safe AD tool.

For performance the reduction clause would be more efficient, as the compiler can make min/max atomic, however I do not know of a nice simple way to implement atomic min/max.

SU2_OMP_BARRIER
SU2_OMP_MASTER
{
SU2_OMP_MASTER {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For efficiency reasons then the Allreduce could be done in a single call with two arguments, as for both Omega and the Strain you compute the maximum.

su2double MyOmega_Max = Omega_Max;
su2double MyStrainMag_Max = StrainMag_Max;

Expand All @@ -988,7 +1074,6 @@ class CFVMFlowSolverBase : public CSolver {
SU2_OMP_BARRIER
}

AD::EndPassive(wasActive);
}

/*!
Expand All @@ -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.
Expand Down Expand Up @@ -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).
Expand Down
Loading