From 148d948a4f4b7fe318813dc0b662310defcd9ce8 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Mon, 2 Mar 2020 15:19:17 +0000 Subject: [PATCH 01/28] prevent restarted FGMRES from going into infinite loop when RHS is zero --- Common/src/linear_algebra/CSysSolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 9bafc45d8620..462f269dbd16 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -919,7 +919,7 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, co /*--- Enforce a hard limit on total number of iterations ---*/ unsigned long IterLimit = min(RestartIter, MaxIter-IterLinSol); IterLinSol += FGMRES_LinSolver(*LinSysRes_ptr, *LinSysSol_ptr, mat_vec, *precond, SolverTol, IterLimit, residual, ScreenOutput, config); - if ( residual < SolverTol*norm0 ) break; + if ( residual <= SolverTol*norm0 ) break; } break; case SMOOTHER: From 4dddf6dfa8d63e716b79e6d6589787999a387ffe Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Mon, 2 Mar 2020 15:27:59 +0000 Subject: [PATCH 02/28] fix for "old compiler" compatibility and legacy build system --- Common/include/toolboxes/C1DInterpolation.hpp | 10 ---------- SU2_CFD/obj/Makefile.am | 4 ++-- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/Common/include/toolboxes/C1DInterpolation.hpp b/Common/include/toolboxes/C1DInterpolation.hpp index a336b33a2448..c5653de37f27 100644 --- a/Common/include/toolboxes/C1DInterpolation.hpp +++ b/Common/include/toolboxes/C1DInterpolation.hpp @@ -84,11 +84,6 @@ class CAkimaInterpolation final: public C1DInterpolation{ SetSpline(X,Data); } - /*! - * \brief Destructor of the CAkimaInterpolation class. - */ - ~CAkimaInterpolation(){} - /*! * \brief for setting the cofficients for the Akima spline. * \param[in] X - the x values. @@ -119,11 +114,6 @@ class CLinearInterpolation final: public C1DInterpolation{ SetSpline(X,Data); } - /*! - * \brief Destructor of the CInletInterpolation class. - */ - ~CLinearInterpolation(){} - /*! * \brief for setting the cofficients for Linear 'spline'. * \param[in] X - the x values. diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 70d94546b765..d7ec34d0ddd7 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -137,7 +137,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/solvers/CFEASolver.cpp \ ../src/solvers/CFEM_DG_EulerSolver.cpp \ ../src/solvers/CFEM_DG_NSSolver.cpp \ - ../src/solvers/CHeatSolverFVM.cpp \ + ../src/solvers/CHeatSolver.cpp \ ../src/solvers/CIncEulerSolver.cpp \ ../src/solvers/CIncNSSolver.cpp \ ../src/solvers/CMeshSolver.cpp \ @@ -170,7 +170,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/variables/CMeshBoundVariable.cpp \ ../src/variables/CMeshElement.cpp \ ../src/variables/CMeshVariable.cpp \ - ../src/variables/CHeatFVMVariable.cpp \ + ../src/variables/CHeatVariable.cpp \ ../src/variables/CVariable.cpp \ ../src/variables/CAdjNSVariable.cpp \ ../src/variables/CTurbSSTVariable.cpp \ From 9344f886bb522ae50278722015db8351c7c50a26 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Mon, 2 Mar 2020 15:55:40 +0000 Subject: [PATCH 03/28] update old build system --- SU2_CFD/obj/Makefile.am | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index d7ec34d0ddd7..a46d07ef0b79 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -59,6 +59,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/integration/CMultiGridIntegration.cpp \ ../src/integration/CStructuralIntegration.cpp \ ../src/integration/CFEM_DG_Integration.cpp \ + ../src/integration/CIntegrationFactory.cpp \ ../src/drivers/CMultizoneDriver.cpp \ ../src/drivers/CSinglezoneDriver.cpp \ ../src/drivers/CDiscAdjSinglezoneDriver.cpp \ @@ -123,6 +124,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/output/CAdjFlowCompOutput.cpp \ ../src/output/CAdjFlowIncOutput.cpp \ ../src/output/CMultizoneOutput.cpp \ + ../src/output/COutputFactory.cpp \ ../src/output/output_structure_legacy.cpp \ ../src/python_wrapper_structure.cpp \ ../src/solvers/CAdjEulerSolver.cpp \ @@ -150,6 +152,7 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/solvers/CTurbSolver.cpp \ ../src/solvers/CTurbSASolver.cpp \ ../src/solvers/CTurbSSTSolver.cpp \ + ../src/solvers/CSolverFactory.cpp \ ../src/limiters/CLimiterDetails.cpp \ ../src/CMarkerProfileReaderFVM.cpp \ ../src/interfaces/CInterface.cpp \ From f4e8185980cc431aa3c4a8ef79a8a08c80a89f5e Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Mon, 2 Mar 2020 16:23:32 +0000 Subject: [PATCH 04/28] unnecessary initialization of stiffness matrix in CMeshSolver --- SU2_CFD/src/solvers/CMeshSolver.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index c0ba9e1abd0b..bf0ab73da6d6 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -436,15 +436,15 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (multizone) nodes->Set_BGSSolution_k(); - /*--- Initialize sparse matrix ---*/ - Jacobian.SetValZero(); - /*--- Compute the stiffness matrix. ---*/ Compute_StiffMatrix(geometry[MESH_0], numerics, config); - /*--- Initialize vectors and clean residual ---*/ - LinSysSol.SetValZero(); - LinSysRes.SetValZero(); + /*--- Initialize vectors and clean residual. ---*/ + SU2_OMP_PARALLEL + { + LinSysSol.SetValZero(); + LinSysRes.SetValZero(); + } /*--- LinSysSol contains the non-transformed displacements in the periodic halo cells. Hence we still need a communication of the transformed coordinates, otherwise periodicity From 101b52d926258fe603fc1e2d987eb37c5ef92475 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Mon, 2 Mar 2020 17:56:38 +0000 Subject: [PATCH 05/28] fix OpenMP bug in SetMesh_Stiffness, allow upper bound on element stiffness, cleanup unused vars --- .../elasticity/CFEALinearElasticity.hpp | 2 +- SU2_CFD/include/solvers/CMeshSolver.hpp | 6 --- .../elasticity/CFEALinearElasticity.cpp | 4 +- SU2_CFD/src/solvers/CMeshSolver.cpp | 42 +++++++++++-------- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp index acda3944b9b2..dd564505bdc4 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp @@ -109,7 +109,7 @@ class CFEAMeshElasticity final : public CFEALinearElasticity { * \param[in] val_nVar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CFEAMeshElasticity(unsigned short val_nDim, unsigned short val_nVar, unsigned long val_nElem, CConfig *config); + CFEAMeshElasticity(unsigned short val_nDim, unsigned short val_nVar, unsigned long val_nElem, const CConfig *config); /*! * \brief Destructor of the class. diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index a39c571bd3e1..0c33951be053 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -48,12 +48,6 @@ class CMeshSolver final : public CFEASolver { su2double MinDistance; su2double MaxDistance; - su2double E; /*!< \brief Young's modulus of elasticity. */ - su2double Nu; /*!< \brief Poisson's ratio. */ - - su2double Mu; /*!< \brief Lame's coeficient. */ - su2double Lambda; /*!< \brief Lame's coeficient. */ - vector element; /*!< \brief Vector which stores element information for each problem. */ /*! diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index c7e71511aaa9..7e09c1410b90 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -345,7 +345,7 @@ void CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CFEAMeshElasticity::CFEAMeshElasticity(unsigned short val_nDim, unsigned short val_nVar, - unsigned long val_nElem, CConfig *config) : + unsigned long val_nElem, const CConfig *config) : CFEALinearElasticity() { DV_Val = NULL; FAux_Dead_Load = NULL; @@ -358,7 +358,7 @@ CFEAMeshElasticity::CFEAMeshElasticity(unsigned short val_nDim, unsigned short v unsigned long iVar; - E = config->GetDeform_ElasticityMod(); + E = 1.0; Nu = config->GetDeform_PoissonRatio(); Compute_Lame_Parameters(); diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index bf0ab73da6d6..da137653ae18 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -120,15 +120,7 @@ CMeshSolver::CMeshSolver(CGeometry *geometry, CConfig *config) : CFEASolver(true ElemColoring[0] = DummyGridColor<>(nElement); #endif - /*--- Structural parameters ---*/ - - E = config->GetDeform_ElasticityMod(); - Nu = config->GetDeform_PoissonRatio(); - - Mu = E / (2.0*(1.0 + Nu)); - Lambda = Nu*E/((1.0+Nu)*(1.0-2.0*Nu)); - - /*--- Element container structure ---*/ + /*--- Element container structure. ---*/ if (nDim == 2) { for(int thread = 0; thread < omp_get_max_threads(); ++thread) { @@ -410,23 +402,37 @@ void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) { void CMeshSolver::SetMesh_Stiffness(CGeometry **geometry, CNumerics **numerics, CConfig *config){ - unsigned long iElem; + /*--- Use the config option as an upper bound on elasticity modulus. + * For RANS meshes the range of element volume or wall distance is + * very large and leads to an ill-conditioned stiffness matrix. + * Absolute values of elasticity modulus are not important for + * mesh deformation, since linear elasticity is used and all + * boundary conditions are essential (Dirichlet). ---*/ + const su2double maxE = config->GetDeform_ElasticityMod(); if (!stiffness_set) { - for (iElem = 0; iElem < nElement; iElem++) { + /*--- All threads must execute the entire loop (no worksharing), + * each sets the stiffnesses for its numerics instance. ---*/ + SU2_OMP_PARALLEL + { + CNumerics* myNumerics = numerics[FEA_TERM + omp_get_thread_num()*MAX_TERMS]; + + for (unsigned long iElem = 0; iElem < nElement; iElem++) { + + su2double E = 1.0; switch (config->GetDeform_Stiffness_Type()) { - /*--- Stiffness inverse of the volume of the element ---*/ - case INVERSE_VOLUME: E = 1.0 / element[iElem].GetRef_Volume(); break; - /*--- Stiffness inverse of the distance of the element to the closest wall ---*/ - case SOLID_WALL_DISTANCE: E = 1.0 / element[iElem].GetWallDistance(); break; + /*--- Stiffness inverse of the volume of the element ---*/ + case INVERSE_VOLUME: E = 1.0 / element[iElem].GetRef_Volume(); break; + + /*--- Stiffness inverse of the distance of the element to the closest wall ---*/ + case SOLID_WALL_DISTANCE: E = 1.0 / element[iElem].GetWallDistance(); break; } /*--- Set the element elastic properties in the numerics container ---*/ - numerics[FEA_TERM]->SetMeshElasticProperties(iElem, E); - + myNumerics->SetMeshElasticProperties(iElem, min(E,maxE)); + } } - stiffness_set = true; } From aba265fda0df5208b136bb98b98244e2de8bd0c1 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Tue, 3 Mar 2020 14:27:49 +0000 Subject: [PATCH 06/28] potential fix for potential cause of observed deadlock --- Common/src/linear_algebra/CSysSolve.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index 462f269dbd16..6f041e89f8bb 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -206,6 +206,7 @@ unsigned long CSysSolve::CG_LinSolver(const CSysVector & * do this since the working vectors are shared. ---*/ if (!cg_ready) { + SU2_OMP_BARRIER SU2_OMP_MASTER { auto nVar = b.GetNVar(); @@ -348,6 +349,7 @@ unsigned long CSysSolve::FGMRES_LinSolver(const CSysVector::BCGSTAB_LinSolver(const CSysVector::Smoother_LinSolver(const CSysVector Date: Wed, 4 Mar 2020 15:27:44 +0000 Subject: [PATCH 07/28] add dummy locks and functions to omp_structure --- Common/include/omp_structure.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/Common/include/omp_structure.hpp b/Common/include/omp_structure.hpp index 0a4ac692ae92..c352483aa4bd 100644 --- a/Common/include/omp_structure.hpp +++ b/Common/include/omp_structure.hpp @@ -80,6 +80,19 @@ inline void omp_set_num_threads(int) { } */ inline constexpr int omp_get_thread_num(void) {return 0;} +/*! + * \brief Dummy lock type and associated functions. + */ +struct omp_lock_t {}; +struct DummyVectorOfLocks { + omp_lock_t l; + inline omp_lock_t& operator[](int) {return l;} +}; +inline void omp_init_lock(omp_lock_t*){} +inline void omp_set_lock(omp_lock_t*){} +inline void omp_unset_lock(omp_lock_t*){} +inline void omp_destroy_lock(omp_lock_t*){} + #endif /*--- Convenience macros (do not use excessive nesting of macros). ---*/ From 7e25a93387354495e38506db909f34da3aa9f0e5 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 15:32:42 +0000 Subject: [PATCH 08/28] bad coloring fallback strategy for CFEASolver --- Common/include/geometry/CGeometry.hpp | 15 ++- Common/include/option_structure.hpp | 2 + Common/include/toolboxes/graph_toolbox.hpp | 54 +++++++++- Common/src/geometry/CGeometry.cpp | 66 ++++++++++--- SU2_CFD/include/solvers/CFEASolver.hpp | 19 +++- SU2_CFD/src/solvers/CFEASolver.cpp | 110 ++++++++++++++++----- SU2_CFD/src/solvers/CMeshSolver.cpp | 22 +---- 7 files changed, 217 insertions(+), 71 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 3066dcd97fd6..b8a23510b5f8 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1621,7 +1621,7 @@ class CGeometry { * \param[in] fillLvl - Level of fill of the pattern. * \return Reference to the sparse pattern. */ - const CCompressedSparsePatternUL& GetSparsePattern(ConnectivityType type, unsigned long fillLvl); + const CCompressedSparsePatternUL& GetSparsePattern(ConnectivityType type, unsigned long fillLvl = 0); /*! * \brief Get the edge to sparse pattern map. @@ -1630,12 +1630,20 @@ class CGeometry { */ const CEdgeToNonZeroMapUL& GetEdgeToSparsePatternMap(void); + /*! + * \brief Get the transpose of the (main, i.e 0 fill) sparse pattern (e.g. CSR becomes CSC). + * \param[in] type - Finite volume or finite element. + * \return Reference to the map. + */ + const su2vector& GetTransposeSparsePatternMap(ConnectivityType type); + /*! * \brief Get the edge coloring. * \note This method computes the coloring if that has not been done yet. + * \param[out] efficiency - optional output of the coloring efficiency. * \return Reference to the coloring. */ - const CCompressedSparsePatternUL& GetEdgeColoring(void); + const CCompressedSparsePatternUL& GetEdgeColoring(su2double* efficiency = nullptr); /*! * \brief Get the group size used in edge coloring. @@ -1646,9 +1654,10 @@ class CGeometry { /*! * \brief Get the element coloring. * \note This method computes the coloring if that has not been done yet. + * \param[out] efficiency - optional output of the coloring efficiency. * \return Reference to the coloring. */ - const CCompressedSparsePatternUL& GetElementColoring(void); + const CCompressedSparsePatternUL& GetElementColoring(su2double* efficiency = nullptr); /*! * \brief Get the group size used in element coloring. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index f2a115236b8c..1a9430e06ab4 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -133,6 +133,8 @@ const int SU2_CONN_SIZE = 10; /*!< \brief Size of the connectivity array that that we read from a mesh file in the format [[globalID vtkType n0 n1 n2 n3 n4 n5 n6 n7 n8]. */ const int SU2_CONN_SKIP = 2; /*!< \brief Offset to skip the globalID and VTK type at the start of the element connectivity list for each CGNS element. */ +const su2double COLORING_EFF_THRESH = 0.8; /*!< \brief Below this value the reducer strategy is used instead. */ + /*! * \brief Boolean answers */ diff --git a/Common/include/toolboxes/graph_toolbox.hpp b/Common/include/toolboxes/graph_toolbox.hpp index 5f5ea4cdf5bc..d5806b135169 100644 --- a/Common/include/toolboxes/graph_toolbox.hpp +++ b/Common/include/toolboxes/graph_toolbox.hpp @@ -59,6 +59,7 @@ class CCompressedSparsePattern { su2vector m_outerPtr; /*!< \brief Start positions of the inner indices for each outer index. */ su2vector m_innerIdx; /*!< \brief Inner indices of the non zero entries. */ su2vector m_diagPtr; /*!< \brief Position of the diagonal entry. */ + su2vector m_innerIdxTransp; /*!< \brief Position of the transpose non zero entries, requires symmetry. */ public: using IndexType = Index_t; @@ -111,6 +112,23 @@ class CCompressedSparsePattern { m_diagPtr(k) = findInnerIdx(k,k); } + /*! + * \brief Build a list of pointers to the transpose entries of the pattern, requires symmetry. + */ + void buildTransposePtr() { + if(!m_innerIdxTransp.empty()) return; + + m_innerIdxTransp.resize(getNumNonZeros()); + + for(Index_t i = 0; i < getOuterSize(); ++i) { + for(Index_t k = m_outerPtr(i); k < m_outerPtr(i+1); ++k) { + auto j = m_innerIdx(k); + m_innerIdxTransp(k) = findInnerIdx(j,i); + assert(m_innerIdxTransp(k) != m_innerIdx.size() && "The pattern is not symmetric."); + } + } + } + /*! * \return True if the pattern is empty, i.e. has not been built yet. */ @@ -224,6 +242,14 @@ class CCompressedSparsePattern { return m_diagPtr.data(); } + /*! + * \return Raw pointer to the transpose pointer vector. + */ + inline const su2vector& transposePtr() const { + assert(!m_innerIdxTransp.empty() && "Transpose map has not been built."); + return m_innerIdxTransp; + } + /*! * \return The minimum inner index. */ @@ -404,7 +430,7 @@ CEdgeToNonZeroMap mapEdgesToSparsePattern(Geometry_t& geometry, * \param[out] indexColor - Optional, vector with colors given to the outer indices. * \return Coloring in the same type of the input pattern. */ -template +template T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors = false, std::vector* indexColor = nullptr) { @@ -553,9 +579,11 @@ struct GridColor static_assert(std::is_integral::value,""); const T size; + T groupSize; const T* const indices; - GridColor(const T* idx = nullptr, T sz = 0) : size(sz), indices(idx) { } + GridColor(const T* idx = nullptr, T sz = 0, T grp = 0) : + size(sz), groupSize(grp), indices(idx) { } inline const T* begin() const {return indices;} inline const T* end() const {return indices+size;} @@ -592,3 +620,25 @@ struct DummyGridColor inline IteratorLikeInt begin() const {return IteratorLikeInt(0);} inline IteratorLikeInt end() const {return IteratorLikeInt(size);} }; + + +/*! + * \brief Computes the efficiency of a grid coloring for given number of threads and chunk size. + */ +template +su2double coloringEfficiency(const SparsePattern& coloring, int numThreads, int chunkSize) +{ + using Index_t = typename SparsePattern::IndexType; + + /*--- Ideally compute time is proportional to total work over number of threads. ---*/ + su2double ideal = coloring.getNumNonZeros() / su2double(numThreads); + + /*--- In practice the total work is quantized first by colors and then by chunks. ---*/ + auto roundUpDiv = [](int n, int d){return (n+d-1)/d;}; + + Index_t real = 0; + for(Index_t color = 0; color < coloring.getOuterSize(); ++color) + real += chunkSize * roundUpDiv(roundUpDiv(coloring.getNumNonZeros(color), chunkSize), numThreads); + + return ideal / real; +} diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 05af0f1c15c0..5b1710f7af9d 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -4036,13 +4036,25 @@ const CEdgeToNonZeroMapUL& CGeometry::GetEdgeToSparsePatternMap(void) return edgeToCSRMap; } -const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(void) +const su2vector& CGeometry::GetTransposeSparsePatternMap(ConnectivityType type) { - if (edgeColoring.empty() && nEdge) { + /*--- Yes the const cast is weird but it is still better than repeating code. ---*/ + auto& pattern = const_cast(GetSparsePattern(type)); + pattern.buildTransposePtr(); + return pattern.transposePtr(); +} + +const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficiency) +{ + if (nEdge==0) return edgeColoring; + + /*--- Build if required. ---*/ + if (edgeColoring.empty()) { - /*--- When not using threading, and for coarse grids, use the natural coloring. ---*/ - if ((omp_get_max_threads() == 1) || (MGLevel != MESH_0)) { + /*--- When not using threading use the natural coloring. ---*/ + if (omp_get_max_threads() == 1) { edgeColoring = createNaturalColoring(nEdge); + if (efficiency != nullptr) *efficiency = 1.0; // by definition return edgeColoring; } @@ -4051,7 +4063,7 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(void) su2vector outerPtr(nEdge+1); su2vector innerIdx(nEdge*2); - for(unsigned long iEdge = 0; iEdge < nEdge; ++iEdge) { + for (unsigned long iEdge = 0; iEdge < nEdge; ++iEdge) { outerPtr(iEdge) = 2*iEdge; innerIdx(iEdge*2+0) = edge[iEdge]->GetNode(0); innerIdx(iEdge*2+1) = edge[iEdge]->GetNode(1); @@ -4060,23 +4072,37 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(void) CCompressedSparsePatternUL pattern(move(outerPtr), move(innerIdx)); - /*--- Color the edges, only balance sizes on coarse levels. ---*/ - bool balanceColors = (MGLevel != MESH_0); + /*--- Color the edges. ---*/ + const bool balanceColors = true; edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, balanceColors); - if(edgeColoring.empty()) - SU2_MPI::Error("Edge coloring failed.", CURRENT_FUNCTION); + /*--- If the coloring fails use the natural coloring and set the group size + * to nEdge to prevent client code from looping in parallel. This is a + * "soft" failure as the "bad" coloring should be detected downstream + * and a fallback strategy put in place. ---*/ + if (edgeColoring.empty()) { + edgeColoring = createNaturalColoring(nEdge); + edgeColorGroupSize = nEdge; + } + } + + if (efficiency != nullptr) { + *efficiency = coloringEfficiency(edgeColoring, omp_get_max_threads(), edgeColorGroupSize); } return edgeColoring; } -const CCompressedSparsePatternUL& CGeometry::GetElementColoring(void) +const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* efficiency) { - if (elemColoring.empty() && nElem) { + if (nElem==0) return elemColoring; + + /*--- Build if required. ---*/ + if (elemColoring.empty()) { /*--- When not using threading use the natural coloring. ---*/ if (omp_get_max_threads() == 1) { elemColoring = createNaturalColoring(nElem); + if (efficiency != nullptr) *efficiency = 1.0; // by definition return elemColoring; } @@ -4085,10 +4111,10 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(void) vector outerPtr(nElem+1); vector innerIdx; innerIdx.reserve(nElem); - for(unsigned long iElem = 0; iElem < nElem; ++iElem) { + for (unsigned long iElem = 0; iElem < nElem; ++iElem) { outerPtr[iElem] = innerIdx.size(); - for(unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); ++iNode) { + for (unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); ++iNode) { innerIdx.push_back(elem[iElem]->GetNode(iNode)); } } @@ -4097,10 +4123,18 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(void) CCompressedSparsePatternUL pattern(outerPtr, innerIdx); /*--- Color the elements. ---*/ - elemColoring = colorSparsePattern(pattern, elemColorGroupSize); + const bool balanceColors = true; + elemColoring = colorSparsePattern(pattern, elemColorGroupSize, balanceColors); + + /*--- Same as for the edge coloring. ---*/ + if (elemColoring.empty()) { + elemColoring = createNaturalColoring(nElem); + elemColorGroupSize = nElem; + } + } - if(elemColoring.empty()) - SU2_MPI::Error("Element coloring failed.", CURRENT_FUNCTION); + if (efficiency != nullptr) { + *efficiency = coloringEfficiency(elemColoring, omp_get_max_threads(), elemColorGroupSize); } return elemColoring; } diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index eac3ff18a3bb..b4a20c377d97 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -77,15 +77,18 @@ class CFEASolver : public CSolver { #ifdef HAVE_OMP vector > ElemColoring; /*!< \brief Element colors. */ + bool LockStrategy = false; /*!< \brief Whether to use an OpenMP lock to guard updates of the Jacobian. */ + vector UpdateLocks; /*!< \brief Locks that may be used to protect accesses to CSysMatrix/Vector in element loops. */ #else - array,1> ElemColoring; + array,1> ElemColoring; /*--- Behaves like a normal integer type. ---*/ + static constexpr bool LockStrategy = false; /*--- Lock strategy is never needed for MPI-only. ---*/ + DummyVectorOfLocks UpdateLocks; #endif - unsigned long ColorGroupSize; /*!< \brief Group size used for coloring, chunk size must be a multiple of this. */ - bool element_based; /*!< \brief Bool to determine if an element-based file is used. */ - bool topol_filter_applied; /*!< \brief True if density filtering has been performed. */ + bool element_based; /*!< \brief Bool to determine if an element-based file is used. */ + bool topol_filter_applied; /*!< \brief True if density filtering has been performed. */ - unsigned long nElement; /*!< \brief Number of elements. */ + unsigned long nElement; /*!< \brief Number of elements. */ /*! * \brief The highest level in the variable hierarchy this solver can safely use, @@ -116,6 +119,12 @@ class CFEASolver : public CSolver { } } + /*! + * \brief Actions required to initialize the supporting variables for hybrid parallel execution. + * \param[in] geometry - Geometrical definition of the problem. + */ + void HybridParallelInitialization(CGeometry* geometry); + /*! * \brief Set container of element properties. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 8f0066180225..2e0559ebabbc 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -240,25 +240,8 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() { LinSysReact.Initialize(nPoint, nPointDomain, nVar, 0.0); -#ifdef HAVE_OMP - /*--- Get the element coloring. ---*/ - - const auto& coloring = geometry->GetElementColoring(); - - if (!coloring.empty()) { - auto nColor = coloring.getOuterSize(); - ElemColoring.reserve(nColor); - - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - ElemColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor)); - } - } - ColorGroupSize = geometry->GetElementColorGroupSize(); - - omp_chunk_size = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_SIZE); -#else - ElemColoring[0] = DummyGridColor<>(nElement); -#endif + /*--- Initialize structures for hybrid-parallel mode. ---*/ + HybridParallelInitialization(geometry); iElem_iDe = nullptr; @@ -338,6 +321,51 @@ CFEASolver::~CFEASolver(void) { if (iElem_iDe != nullptr) delete [] iElem_iDe; if (nodes != nullptr) delete nodes; + + if (LockStrategy) { + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + omp_destroy_lock(&UpdateLocks[iPoint]); + } +} + +void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { +#ifdef HAVE_OMP + /*--- Get the element coloring. ---*/ + + su2double parallelEff = 1.0; + const auto& coloring = geometry->GetElementColoring(¶llelEff); + + if (!coloring.empty()) { + auto nColor = coloring.getOuterSize(); + ElemColoring.reserve(nColor); + + for(auto iColor = 0ul; iColor < nColor; ++iColor) { + ElemColoring.emplace_back(coloring.innerIdx(iColor), + coloring.getNumNonZeros(iColor), + geometry->GetElementColorGroupSize()); + } + } + + su2double minEff = 1.0; + SU2_MPI::Reduce(¶llelEff, &minEff, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); + + if (minEff < COLORING_EFF_THRESH) { + cout << "WARNING: The element coloring efficiency was " << minEff << ", a fallback strategy is in use.\n" + << " Better performance may be possible by reducing the number of threads per rank." << endl; + } + + if (LockStrategy) { + UpdateLocks.resize(nPoint); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + omp_init_lock(&UpdateLocks[iPoint]); + /*--- We are no longer constrained by the color group size. ---*/ + for(auto& color : ElemColoring) color.groupSize = OMP_MIN_SIZE; + } + + omp_chunk_size = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_SIZE); +#else + ElemColoring[0] = DummyGridColor<>(nElement); +#endif } void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { @@ -829,7 +857,7 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -879,6 +907,8 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, /*--- Update residual and stiffness matrix with contributions from the element. ---*/ for (iNode = 0; iNode < nNodes; iNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); + auto Ta = element->Get_Kt_a(iNode); for (iVar = 0; iVar < nVar; iVar++) LinSysRes(indexNode[iNode], iVar) -= simp_penalty*Ta[iVar]; @@ -887,6 +917,8 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, auto Kab = element->Get_Kab(iNode, jNode); Jacobian.AddBlock(indexNode[iNode], indexNode[jNode], simp_penalty, Kab); } + + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); } } // end iElem loop @@ -917,7 +949,7 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -986,6 +1018,8 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri /*--- Update residual and stiffness matrix with contributions from the element. ---*/ for (iNode = 0; iNode < nNodes; iNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); + auto Ta = fea_elem->Get_Kt_a(iNode); for (iVar = 0; iVar < nVar; iVar++) LinSysRes(indexNode[iNode], iVar) -= simp_penalty*Ta[iVar]; @@ -1023,6 +1057,8 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri Kij[iVar*(nVar+1)] += SU2_TYPE::GetValue(simp_penalty*Ks_ab_DE); } } + + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); } } // end iElem loop @@ -1048,7 +1084,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1091,6 +1127,8 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C /*--- Add contributions of this element to the mass matrix. ---*/ for (iNode = 0; iNode < nNodes; iNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); + for (jNode = 0; jNode < nNodes; jNode++) { auto Mij = MassMatrix.GetBlock(indexNode[iNode], indexNode[jNode]); @@ -1099,6 +1137,8 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C for (iVar = 0; iVar < nVar; iVar++) Mij[iVar*(nVar+1)] += simp_penalty*Mab; } + + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); } } // end iElem loop @@ -1125,7 +1165,7 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CCon for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1165,7 +1205,10 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CCon numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Mass_Matrix(element, config); - /*--- Add contributions of this element to the mass matrix. ---*/ + /*--- Add contributions of this element to the mass matrix. + * In case we need to use locks we guard the entire update. ---*/ + if (LockStrategy) omp_set_lock(&UpdateLocks[0]); + for (iNode = 0; iNode < nNodes; iNode++) { for (jNode = 0; jNode < nNodes; jNode++) { @@ -1177,6 +1220,7 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CCon } } } + if (LockStrategy) omp_unset_lock(&UpdateLocks[0]); } // end iElem loop @@ -1205,7 +1249,7 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1260,9 +1304,13 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric numerics[NUM_TERM]->Compute_NodalStress_Term(element, config); for (iNode = 0; iNode < nNodes; iNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); + auto Ta = element->Get_Kt_a(iNode); for (iVar = 0; iVar < nVar; iVar++) LinSysRes(indexNode[iNode], iVar) -= simp_penalty*Ta[iVar]; + + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); } } // end iElem loop @@ -1302,7 +1350,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1359,6 +1407,8 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, auto iPoint = indexNode[iNode]; + if (LockStrategy) omp_set_lock(&UpdateLocks[iPoint]); + auto Ta = element->Get_Kt_a(iNode); for (iVar = 0; iVar < nVar; iVar++) LinSysReact(iPoint,iVar) += simp_penalty*Ta[iVar]; @@ -1368,6 +1418,8 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, for (iStress = 0; iStress < nStress; iStress++) nodes->AddStress_FEM(iPoint,iStress, weight*element->Get_NodalStress(iNode,iStress)); + + if (LockStrategy) omp_unset_lock(&UpdateLocks[iPoint]); } } // end iElem loop @@ -1564,7 +1616,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1604,6 +1656,8 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo /*--- Add contributions of this element to the mass matrix. ---*/ for (iNode = 0; iNode < nNodes; iNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); + auto Dead_Load = element->Get_FDL_a(iNode); su2double Aux_Dead_Load[MAXNVAR]; @@ -1611,6 +1665,8 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo Aux_Dead_Load[iVar] = simp_penalty*Dead_Load[iVar]; nodes->Add_BodyForces_Res(indexNode[iNode], Aux_Dead_Load); + + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); } } // end iElem loop diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index da137653ae18..25e4b711fc9f 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -96,29 +96,15 @@ CMeshSolver::CMeshSolver(CGeometry *geometry, CConfig *config) : CFEASolver(true /*--- Initialize matrix, solution, and r.h.s. structures for the linear solver. ---*/ + if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (Mesh Deformation)." << endl; + LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, false, geometry, config); -#ifdef HAVE_OMP - /*--- Get the element coloring. ---*/ - - const auto& coloring = geometry->GetElementColoring(); - - if (!coloring.empty()) { - auto nColor = coloring.getOuterSize(); - ElemColoring.reserve(nColor); - - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - ElemColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor)); - } - } - ColorGroupSize = geometry->GetElementColorGroupSize(); + /*--- Initialize structures for hybrid-parallel mode. ---*/ - omp_chunk_size = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_SIZE); -#else - ElemColoring[0] = DummyGridColor<>(nElement); -#endif + HybridParallelInitialization(geometry); /*--- Element container structure. ---*/ From 32ef1941f92a38d3d455ffff091f2228855cfc9e Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 15:33:26 +0000 Subject: [PATCH 09/28] use the "reducer strategy" as a fallback for when grid coloring is bad also on the fine grid --- Common/include/linear_algebra/CSysMatrix.hpp | 6 +- Common/src/linear_algebra/CSysMatrix.cpp | 30 ++----- SU2_CFD/include/solvers/CEulerSolver.hpp | 14 +-- SU2_CFD/include/solvers/CTurbSolver.hpp | 15 +++- SU2_CFD/src/solvers/CEulerSolver.cpp | 89 ++++++++++++-------- SU2_CFD/src/solvers/CTurbSASolver.cpp | 2 +- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 2 +- SU2_CFD/src/solvers/CTurbSolver.cpp | 22 +++-- 8 files changed, 104 insertions(+), 76 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index b01be62fe13f..a992549ce590 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -95,7 +95,7 @@ class CSysMatrix { const unsigned long *row_ptr; /*!< \brief Pointers to the first element in each row. */ const unsigned long *dia_ptr; /*!< \brief Pointers to the diagonal element in each row. */ const unsigned long *col_ind; /*!< \brief Column index for each of the elements in val(). */ - vector col_ptr;/*!< \brief The transpose of col_ind, pointer to blocks with the same column index. */ + const unsigned long *col_ptr; /*!< \brief The transpose of col_ind, pointer to blocks with the same column index. */ ScalarType *ILU_matrix; /*!< \brief Entries of the ILU sparse matrix. */ unsigned long nnz_ilu; /*!< \brief Number of possible nonzero entries in the matrix (ILU). */ @@ -349,10 +349,12 @@ class CSysMatrix { * \param[in] neqn - Number of equations. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. + * \param[in] needTranspPtr - If "col_ptr" should be created. */ void Initialize(unsigned long npoint, unsigned long npointdomain, unsigned short nvar, unsigned short neqn, - bool EdgeConnect, CGeometry *geometry, CConfig *config); + bool EdgeConnect, CGeometry *geometry, + CConfig *config, bool needTranspPtr = false); /*! * \brief Sets to zero all the entries of the sparse matrix. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 577caff65fc7..6d97a3518e98 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -89,8 +89,9 @@ CSysMatrix::~CSysMatrix(void) { template void CSysMatrix::Initialize(unsigned long npoint, unsigned long npointdomain, - unsigned short nvar, unsigned short neqn, - bool EdgeConnect, CGeometry *geometry, CConfig *config) { + unsigned short nvar, unsigned short neqn, + bool EdgeConnect, CGeometry *geometry, + CConfig *config, bool needTranspPtr) { assert(omp_get_thread_num()==0 && "Only the master thread is allowed to initialize the matrix."); @@ -126,10 +127,13 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi const auto& csr = geometry->GetSparsePattern(type,0); + nnz = csr.getNumNonZeros(); row_ptr = csr.outerPtr(); col_ind = csr.innerIdx(); dia_ptr = csr.diagPtr(); - nnz = csr.getNumNonZeros(); + + if (needTranspPtr) + col_ptr = geometry->GetTransposeSparsePatternMap(type).data(); if (type == ConnectivityType::FiniteVolume) edge_ptr.ptr = geometry->GetEdgeToSparsePatternMap().data(); @@ -189,24 +193,6 @@ void CSysMatrix::Initialize(unsigned long npoint, unsigned long npoi omp_partitions[part] = part * pts_per_part; omp_partitions[omp_num_parts] = nPointDomain; - /*--- For coarse grid levels setup a structure that allows doing a column sum efficiently, - * essentially the transpose of the col_ind, this allows populating the matrix by setting - * the off-diagonal entries and then setting the diagonal ones as the sum of column - * (excluding the diagonal itself). We use the fact that the pattern is symmetric. ---*/ - - if ((geometry->GetMGLevel() != MESH_0) && (omp_get_max_threads() > 1)) { - - col_ptr.resize(nnz, nullptr); - - SU2_OMP_PARALLEL_(for schedule(static,omp_heavy_size)) - for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { - for (auto k = row_ptr[iPoint]; k < row_ptr[iPoint+1]; ++k) { - auto jPoint = col_ind[k]; - col_ptr[k] = GetBlock(jPoint, iPoint); - } - } - } - /*--- Generate MKL Kernels ---*/ #ifdef USE_MKL @@ -1334,7 +1320,7 @@ void CSysMatrix::SetDiagonalAsColumnSum() { auto block_ii = &matrix[dia_ptr[iPoint]*nVar*nEqn]; for (auto k = row_ptr[iPoint]; k < row_ptr[iPoint+1]; ++k) { - auto block_ji = col_ptr[k]; + auto block_ji = &matrix[col_ptr[k]*nVar*nEqn]; if (block_ji != block_ii) MatrixSubtraction(block_ii, block_ji, block_ii); } } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index df1944c6129a..820a17ffb368 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -264,16 +264,18 @@ class CEulerSolver : public CSolver { #ifdef HAVE_OMP vector > EdgeColoring; /*!< \brief Edge colors. */ + bool ReducerStrategy = false; /*!< \brief If the reducer strategy is in use. */ #else array,1> EdgeColoring; + /*--- Never use the reducer strategy if compiling for MPI-only. ---*/ + static constexpr bool ReducerStrategy = false; #endif - unsigned long ColorGroupSize; /*!< \brief Group size used for coloring, chunk size in edge loops must be a multiple of this. */ - /*--- Edge fluxes, for OpenMP parallelization on coarse grids. As it is difficult to - * color them, we first store the fluxes and then compute the sum for each cell. - * This strategy is thread-safe but lower performance than writting to both end - * points of each edge, so we only use it when necessary, i.e. coarse grids and - * with more than one thread per MPI rank. ---*/ + /*--- Edge fluxes, for OpenMP parallelization off difficult-to-color grids. + * We first store the fluxes and then compute the sum for each cell. + * This strategy is thread-safe but lower performance than writting to both + * end points of each edge, so we only use it when necessary, i.e. when the + * coloring does not allow "enough" parallelism. ---*/ CSysVector EdgeFluxes; /*!< \brief Flux across each edge. */ diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index cfbc0aa44bae..1b9b311fa1c9 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -55,7 +55,7 @@ class CTurbSolver : public CSolver { Gamma_Minus_One, /*!< \brief Fluids's Gamma - 1.0 . */ ***Inlet_TurbVars = nullptr; /*!< \brief Turbulence variables at inlet profiles */ - /* Sliding meshes variables */ + /*--- Sliding meshes variables. ---*/ su2double ****SlidingState = nullptr; int **SlidingStateNodes = nullptr; @@ -64,10 +64,15 @@ class CTurbSolver : public CSolver { #ifdef HAVE_OMP vector > EdgeColoring; /*!< \brief Edge colors. */ + bool ReducerStrategy = false; /*!< \brief If the reducer strategy is in use. */ #else array,1> EdgeColoring; + /*--- Never use the reducer strategy if compiling for MPI-only. ---*/ + static constexpr bool ReducerStrategy = false; #endif - unsigned long ColorGroupSize; /*!< \brief Group size used for coloring, chunk size in edge loops must be a multiple of this. */ + + /*--- Edge fluxes for reducer strategy (see the notes in CEulerSolver.hpp). ---*/ + CSysVector EdgeFluxes; /*!< \brief Flux across each edge. */ /*! * \brief The highest level in the variable hierarchy this solver can safely use. @@ -79,6 +84,12 @@ class CTurbSolver : public CSolver { */ inline CVariable* GetBaseClassPointerToNodes() final { return nodes; } + /*! + * \brief Sum the edge fluxes for each cell to populate the residual vector, only used on coarse grids. + * \param[in] geometry - Geometrical definition of the problem. + */ + void SumEdgeFluxes(CGeometry* geometry); + public: /*! diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index addb18101fe3..36996595b935 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -208,6 +208,50 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); +#ifdef HAVE_OMP + /*--- Get the edge coloring. If the expected parallel efficiency becomes too low setup the + * reducer strategy. Where one loop is performed over edges followed by a point loop to + * sum the fluxes for each cell and set the diagonal of the system matrix. ---*/ + + su2double parallelEff = 1.0; + const auto& coloring = geometry->GetEdgeColoring(¶llelEff); + + if (!coloring.empty()) { + auto nColor = coloring.getOuterSize(); + EdgeColoring.reserve(nColor); + + for(auto iColor = 0ul; iColor < nColor; ++iColor) { + EdgeColoring.emplace_back(coloring.innerIdx(iColor), + coloring.getNumNonZeros(iColor), + geometry->GetEdgeColorGroupSize()); + } + } + + /*--- The decision to use the strategy is local to each rank. ---*/ + ReducerStrategy = parallelEff < COLORING_EFF_THRESH; + + /*--- The warning messages are global. ---*/ + su2double minEff = 1.0; + SU2_MPI::Reduce(¶llelEff, &minEff, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); + + int tmp = ReducerStrategy, numRanksUsingReducer = 0; + SU2_MPI::Reduce(&tmp, &numRanksUsingReducer, 1, MPI_INT, MPI_SUM, MASTER_NODE, MPI_COMM_WORLD); + + if (minEff < COLORING_EFF_THRESH) { + cout << "WARNING: On " << numRanksUsingReducer << " MPI ranks the coloring efficiency was less than " + << COLORING_EFF_THRESH << " (min value was " << minEff << ").\n" + << " Those ranks will now use a fallback strategy, better performance may be possible\n" + << " with different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl; + } + + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + + omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); +#else + EdgeColoring[0] = DummyGridColor<>(geometry->GetnEdge()); +#endif + /*--- Jacobians and vector structures for implicit computations ---*/ if (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT) { @@ -222,7 +266,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (" << description << "). MG level: " << iMesh <<"." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); @@ -420,32 +464,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, } SetBaseClassPointerToNodes(); -#ifdef HAVE_OMP - /*--- Get the edge coloring, on coarse grids (which are difficult to color) setup the reducer strategy, - * if required (i.e. in parallel). One loop is performed over edges followed by a point loop to sum - * the fluxes for each cell and set the diagonal of the system matrix. ---*/ - - const auto& coloring = geometry->GetEdgeColoring(); - - if (!coloring.empty()) { - auto nColor = coloring.getOuterSize(); - EdgeColoring.reserve(nColor); - - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - EdgeColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor)); - } - } - ColorGroupSize = geometry->GetEdgeColorGroupSize(); - - if ((MGLevel != MESH_0) && (omp_get_max_threads() > 1)) { - EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); - } - - omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); -#else - EdgeColoring[0] = DummyGridColor<>(geometry->GetnEdge()); -#endif - /*--- Check that the initial solution is physical, report any non-physical nodes ---*/ counter_local = 0; @@ -485,7 +503,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, SU2_MPI::Reduce(&counter_local, &counter_global, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, MPI_COMM_WORLD); if ((rank == MASTER_NODE) && (counter_global != 0)) - cout << "Warning. The original solution contains "<< counter_global << " points that are not physical." << endl; + cout << "Warning. The original solution contains " << counter_global << " points that are not physical." << endl; } /*--- Initialize the BGS residuals in FSI problems. ---*/ @@ -538,6 +556,10 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, /*--- Add the solver name (max 8 characters) ---*/ SolverName = "C.FLOW"; + /*--- Finally, check that the static arrays will be large enough (keep this + * check at the bottom to make sure we consider the "final" values). ---*/ + if((nDim > MAXNDIM) || (nPrimVar > MAXNVAR) || (nSecondaryVar > MAXNVAR)) + SU2_MPI::Error("Oops! The CEulerSolver static array sizes are not large enough.",CURRENT_FUNCTION); } CEulerSolver::~CEulerSolver(void) { @@ -2917,7 +2939,7 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -2991,9 +3013,6 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { - assert(nDim <= MAXNDIM && nPrimVar <= MAXNVAR && nSecondaryVar <= MAXNVAR && - "Oops! The CEulerSolver static array sizes are not large enough."); - const auto InnerIter = config->GetInnerIter(); const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || @@ -3027,7 +3046,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -3640,7 +3659,7 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -3714,7 +3733,7 @@ void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig * for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 284edf7ad6ee..fae8274a64d9 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -96,7 +96,7 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor /*--- Initialization of the structure of the whole Jacobian ---*/ if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SA model)." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 148462ef1f77..3a87faa570b2 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -98,7 +98,7 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Initialization of the structure of the whole Jacobian ---*/ if (rank == MASTER_NODE) cout << "Initialize Jacobian structure (SST model)." << endl; - Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); + Jacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config, ReducerStrategy); if (config->GetKind_Linear_Solver_Prec() == LINELET) { nLineLets = Jacobian.BuildLineletPreconditioner(geometry, config); diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index ca1024b209c3..8f32da44a089 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -50,17 +50,25 @@ CTurbSolver::CTurbSolver(CGeometry* geometry, CConfig *config) : CSolver() { #ifdef HAVE_OMP /*--- Get the edge coloring. ---*/ - const auto& coloring = geometry->GetEdgeColoring(); + su2double parallelEff = 1.0; + const auto& coloring = geometry->GetEdgeColoring(¶llelEff); if (!coloring.empty()) { auto nColor = coloring.getOuterSize(); - EdgeColoring.resize(nColor); + EdgeColoring.reserve(nColor); for(auto iColor = 0ul; iColor < nColor; ++iColor) { - EdgeColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor)); + EdgeColoring.emplace_back(coloring.innerIdx(iColor), + coloring.getNumNonZeros(iColor), + geometry->GetEdgeColorGroupSize()); } } - ColorGroupSize = geometry->GetEdgeColorGroupSize(); + + /*--- Local (to the rank) decision to use the reducer strategy, any + * warnings have been, or will be, printed by the flow solver. ---*/ + ReducerStrategy = parallelEff < COLORING_EFF_THRESH; + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); nPoint = geometry->GetnPoint(); omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); @@ -105,7 +113,7 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -236,7 +244,7 @@ void CTurbSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -685,7 +693,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, ColorGroupSize)*ColorGroupSize) + SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; From 14ed268f6b6811618a7103e21d315dcc63b1b380 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 15:39:37 +0000 Subject: [PATCH 10/28] build diagonal and transpose map in parallel --- Common/include/toolboxes/graph_toolbox.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Common/include/toolboxes/graph_toolbox.hpp b/Common/include/toolboxes/graph_toolbox.hpp index d5806b135169..b0a26d579383 100644 --- a/Common/include/toolboxes/graph_toolbox.hpp +++ b/Common/include/toolboxes/graph_toolbox.hpp @@ -28,6 +28,7 @@ #pragma once #include "C2DContainer.hpp" +#include "../omp_structure.hpp" #include #include @@ -108,6 +109,8 @@ class CCompressedSparsePattern { if(!m_diagPtr.empty()) return; m_diagPtr.resize(getOuterSize()); + + SU2_OMP_PARALLEL_(for schedule(static,roundUpDiv(getOuterSize(),omp_get_max_threads()))) for(Index_t k = 0; k < getOuterSize(); ++k) m_diagPtr(k) = findInnerIdx(k,k); } @@ -120,6 +123,7 @@ class CCompressedSparsePattern { m_innerIdxTransp.resize(getNumNonZeros()); + SU2_OMP_PARALLEL_(for schedule(static,roundUpDiv(getOuterSize(),omp_get_max_threads()))) for(Index_t i = 0; i < getOuterSize(); ++i) { for(Index_t k = m_outerPtr(i); k < m_outerPtr(i+1); ++k) { auto j = m_innerIdx(k); @@ -634,8 +638,6 @@ su2double coloringEfficiency(const SparsePattern& coloring, int numThreads, int su2double ideal = coloring.getNumNonZeros() / su2double(numThreads); /*--- In practice the total work is quantized first by colors and then by chunks. ---*/ - auto roundUpDiv = [](int n, int d){return (n+d-1)/d;}; - Index_t real = 0; for(Index_t color = 0; color < coloring.getOuterSize(); ++color) real += chunkSize * roundUpDiv(roundUpDiv(coloring.getNumNonZeros(color), chunkSize), numThreads); From d7bb8418de188353d7d89862371e2db2622490c7 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 16:11:46 +0000 Subject: [PATCH 11/28] use a more expressive function to round up to next multiple --- Common/include/omp_structure.hpp | 8 ++++++++ SU2_CFD/src/solvers/CEulerSolver.cpp | 26 ++++++++++---------------- SU2_CFD/src/solvers/CFEASolver.cpp | 14 +++++++------- 3 files changed, 25 insertions(+), 23 deletions(-) diff --git a/Common/include/omp_structure.hpp b/Common/include/omp_structure.hpp index c352483aa4bd..ad6e1c29d57b 100644 --- a/Common/include/omp_structure.hpp +++ b/Common/include/omp_structure.hpp @@ -121,6 +121,14 @@ inline constexpr size_t roundUpDiv(size_t numerator, size_t denominator) return (numerator+denominator-1)/denominator; } +/*! + * \brief Round up to next multiple. + */ +inline constexpr size_t nextMultiple(size_t argument, size_t multiple) +{ + return roundUpDiv(argument, multiple) * multiple; +} + /*! * \brief Compute a chunk size based on totalWork and number of threads such that * all threads get the same number of chunks (with limited size). diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 36996595b935..cd70602a06f9 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2932,14 +2932,11 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Pick one numerics object per thread. ---*/ CNumerics* numerics = numerics_container[CONV_TERM + omp_get_thread_num()*MAX_TERMS]; - /*--- Determine if using the reducer strategy is necessary, see CEulerSolver::SumEdgeFluxes(). ---*/ - const bool reducer_strategy = (MGLevel != MESH_0) && (omp_get_num_threads() > 1); - /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -2981,7 +2978,7 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Update convective and artificial dissipation residuals. ---*/ - if (reducer_strategy) { + if (ReducerStrategy) { EdgeFluxes.SetBlock(iEdge, residual); if (implicit) Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); @@ -3002,7 +2999,7 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta } } // end color loop - if (reducer_strategy) { + if (ReducerStrategy) { SumEdgeFluxes(geometry); if (implicit) Jacobian.SetDiagonalAsColumnSum(); @@ -3039,14 +3036,11 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain su2double Primitive_i[MAXNVAR] = {0.0}, Primitive_j[MAXNVAR] = {0.0}; su2double Secondary_i[MAXNVAR] = {0.0}, Secondary_j[MAXNVAR] = {0.0}; - /*--- Determine if using the reducer strategy is necessary, see CEulerSolver::SumEdgeFluxes(). ---*/ - const bool reducer_strategy = (MGLevel != MESH_0) && (omp_get_num_threads() > 1); - /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -3212,7 +3206,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Update residual value ---*/ - if (reducer_strategy) { + if (ReducerStrategy) { EdgeFluxes.SetBlock(iEdge, residual); if (implicit) Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); @@ -3233,7 +3227,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } } // end color loop - if (reducer_strategy) { + if (ReducerStrategy) { SumEdgeFluxes(geometry); if (implicit) Jacobian.SetDiagonalAsColumnSum(); @@ -3659,7 +3653,7 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -3733,7 +3727,7 @@ void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig * for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 2e0559ebabbc..81a1990a19b6 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -857,7 +857,7 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -949,7 +949,7 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1084,7 +1084,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1165,7 +1165,7 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CCon for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1249,7 +1249,7 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1350,7 +1350,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; @@ -1616,7 +1616,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo for(auto color : ElemColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iElem = color.indices[k]; From 2fc6f30d870802ad18ff6a30c5fd8f2cd46d5a07 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 16:23:22 +0000 Subject: [PATCH 12/28] use reducer strategy for turbulence solvers --- SU2_CFD/include/solvers/CTurbSolver.hpp | 16 ++-- SU2_CFD/src/solvers/CNSSolver.cpp | 2 +- SU2_CFD/src/solvers/CTurbSolver.cpp | 110 ++++++++++++++---------- 3 files changed, 71 insertions(+), 57 deletions(-) diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index 1b9b311fa1c9..5056e0dd9ec8 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -125,20 +125,18 @@ class CTurbSolver : public CSolver { unsigned short iMesh) override; /*! - * \brief Compute the viscous residuals for the turbulent equation. + * \brief Compute the viscous flux for the turbulent equation at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux * \param[in] geometry - Geometrical definition of the problem. * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics_container - Description of the numerical method. + * \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, + void Viscous_Residual(unsigned long iEdge, + CGeometry *geometry, CSolver **solver_container, - CNumerics **numerics_container, - CConfig *config, - unsigned short iMesh, - unsigned short iRKStep) override; + CNumerics *numerics, + CConfig *config); /*! * \brief Impose the Symmetry Plane boundary condition. diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 2484ce2ed0b7..9d0d57494a24 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -440,7 +440,7 @@ void CNSSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolv auto residual = numerics->ComputeResidual(config); - if ((MGLevel != MESH_0) && (omp_get_num_threads() > 1)) { + if (ReducerStrategy) { EdgeFluxes.SubtractBlock(iEdge, residual); if (implicit) Jacobian.UpdateBlocksSub(iEdge, residual.jacobian_i, residual.jacobian_j); diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 8f32da44a089..6cff4b71cd12 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -112,8 +112,8 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -216,80 +216,96 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe } - /*--- Add and subtract residual ---*/ + /*--- Update convective residual value ---*/ auto residual = numerics->ComputeResidual(config); - LinSysRes.AddBlock(iPoint, residual); - LinSysRes.SubtractBlock(jPoint, residual); - - /*--- Implicit part ---*/ + if (ReducerStrategy) { + EdgeFluxes.SetBlock(iEdge, residual); + Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + } + else { + LinSysRes.AddBlock(iPoint, residual); + LinSysRes.SubtractBlock(jPoint, residual); + Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + } - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + /*--- Viscous contribution. ---*/ + Viscous_Residual(iEdge, geometry, solver_container, + numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS], config); } } // end color loop + if (ReducerStrategy) { + SumEdgeFluxes(geometry); + Jacobian.SetDiagonalAsColumnSum(); + } } -void CTurbSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, - CConfig *config, unsigned short iMesh, unsigned short iRKStep) { +void CTurbSolver::Viscous_Residual(unsigned long iEdge, CGeometry *geometry, CSolver **solver_container, + CNumerics *numerics, CConfig *config) { CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); - /*--- Pick one numerics object per thread. ---*/ - CNumerics* numerics = numerics_container[VISC_TERM + omp_get_thread_num()*MAX_TERMS]; - - /*--- Loop over edge colors. ---*/ - for (auto color : EdgeColoring) - { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) - for(auto k = 0ul; k < color.size; ++k) { + /*--- Points in edge ---*/ - auto iEdge = color.indices[k]; + auto iPoint = geometry->edge[iEdge]->GetNode(0); + auto jPoint = geometry->edge[iEdge]->GetNode(1); - /*--- Points in edge ---*/ + /*--- Points coordinates, and normal vector ---*/ - auto iPoint = geometry->edge[iEdge]->GetNode(0); - auto jPoint = geometry->edge[iEdge]->GetNode(1); + numerics->SetCoord(geometry->node[iPoint]->GetCoord(), + geometry->node[jPoint]->GetCoord()); + numerics->SetNormal(geometry->edge[iEdge]->GetNormal()); - /*--- Points coordinates, and normal vector ---*/ + /*--- Conservative variables w/o reconstruction ---*/ - numerics->SetCoord(geometry->node[iPoint]->GetCoord(), - geometry->node[jPoint]->GetCoord()); - numerics->SetNormal(geometry->edge[iEdge]->GetNormal()); + numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), + flowNodes->GetPrimitive(jPoint)); - /*--- Conservative variables w/o reconstruction ---*/ + /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - numerics->SetPrimitive(flowNodes->GetPrimitive(iPoint), - flowNodes->GetPrimitive(jPoint)); + numerics->SetTurbVar(nodes->GetSolution(iPoint), + nodes->GetSolution(jPoint)); + numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), + nodes->GetGradient(jPoint)); - /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ + /*--- Menter's first blending function (only SST)---*/ + if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) + numerics->SetF1blending(nodes->GetF1blending(iPoint), + nodes->GetF1blending(jPoint)); - numerics->SetTurbVar(nodes->GetSolution(iPoint), - nodes->GetSolution(jPoint)); - numerics->SetTurbVarGradient(nodes->GetGradient(iPoint), - nodes->GetGradient(jPoint)); + /*--- Compute residual, and Jacobians ---*/ - /*--- Menter's first blending function (only SST)---*/ - if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) - numerics->SetF1blending(nodes->GetF1blending(iPoint), - nodes->GetF1blending(jPoint)); + auto residual = numerics->ComputeResidual(config); - /*--- Compute residual, and Jacobians ---*/ + if (ReducerStrategy) { + EdgeFluxes.SubtractBlock(iEdge, residual); + Jacobian.UpdateBlocksSub(iEdge, residual.jacobian_i, residual.jacobian_j); + } + else { + LinSysRes.SubtractBlock(iPoint, residual); + LinSysRes.AddBlock(jPoint, residual); + Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + } +} - auto residual = numerics->ComputeResidual(config); +void CTurbSolver::SumEdgeFluxes(CGeometry* geometry) { - /*--- Add and subtract residual, and update Jacobians ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { - LinSysRes.SubtractBlock(iPoint, residual); - LinSysRes.AddBlock(jPoint, residual); + for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) { - Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + auto iEdge = geometry->node[iPoint]->GetEdge(iNeigh); + if (iPoint == geometry->edge[iEdge]->GetNode(0)) + LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); + else + LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge)); + } } - } // end color loop } @@ -693,7 +709,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con for (auto color : EdgeColoring) { /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(roundUpDiv(OMP_MIN_SIZE, color.groupSize)*color.groupSize) + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; From c88ea7a9a3585016af0f5d0085ede0efa39580f8 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 17:56:38 +0000 Subject: [PATCH 13/28] polish up the reducer strategy, avoid unnecessary resets of CSysMatrix --- Common/include/linear_algebra/CSysMatrix.hpp | 16 ++++-- Common/include/toolboxes/graph_toolbox.hpp | 52 +++++++++++--------- Common/src/geometry/CMultiGridGeometry.cpp | 7 +-- Common/src/geometry/CPhysicalGeometry.cpp | 6 +-- Common/src/linear_algebra/CSysMatrix.cpp | 2 + SU2_CFD/src/solvers/CEulerSolver.cpp | 24 ++++----- SU2_CFD/src/solvers/CNSSolver.cpp | 13 ++--- SU2_CFD/src/solvers/CTurbSASolver.cpp | 12 +++-- SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 12 +++-- SU2_CFD/src/solvers/CTurbSolver.cpp | 6 +-- 10 files changed, 86 insertions(+), 64 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index a992549ce590..73d803a4fbff 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -586,11 +586,13 @@ class CSysMatrix { * \brief Update 2 blocks ij and ji (add to i* sub from j*). * \note The template parameter Sign, can be used create a "subtractive" * update i.e. subtract from row i and add to row j instead. + * The parameter Overwrite allows completely writing over the + * current values held by the matrix. * \param[in] edge - Index of edge that connects iPoint and jPoint. * \param[in] block_i - Subs from ji. * \param[in] block_j - Adds to ij. */ - template + template inline void UpdateBlocks(unsigned long iEdge, const OtherType* const* block_i, const OtherType* const* block_j) { ScalarType *bij = &matrix[edge_ptr(iEdge,0)*nVar*nEqn]; @@ -600,8 +602,8 @@ class CSysMatrix { for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nEqn; jVar++) { - bij[offset] += PassiveAssign(block_j[iVar][jVar]) * Sign; - bji[offset] -= PassiveAssign(block_i[iVar][jVar]) * Sign; + bij[offset] = (Overwrite? ScalarType(0) : bij[offset]) + PassiveAssign(block_j[iVar][jVar]) * Sign; + bji[offset] = (Overwrite? ScalarType(0) : bji[offset]) - PassiveAssign(block_i[iVar][jVar]) * Sign; ++offset; } } @@ -615,6 +617,14 @@ class CSysMatrix { UpdateBlocks(iEdge, block_i, block_j); } + /*! + * \brief Short-hand for the "additive overwrite" version of UpdateBlocks. + */ + template + inline void SetBlocks(unsigned long iEdge, const OtherType* const* block_i, const OtherType* const* block_j) { + UpdateBlocks(iEdge, block_i, block_j); + } + /*! * \brief Adds the specified block to the (i, i) subblock of the matrix-by-blocks structure. * \param[in] block_i - Diagonal index. diff --git a/Common/include/toolboxes/graph_toolbox.hpp b/Common/include/toolboxes/graph_toolbox.hpp index b0a26d579383..e87d48402f9a 100644 --- a/Common/include/toolboxes/graph_toolbox.hpp +++ b/Common/include/toolboxes/graph_toolbox.hpp @@ -414,6 +414,30 @@ CEdgeToNonZeroMap mapEdgesToSparsePattern(Geometry_t& geometry, } +/*! + * \brief Create the natural coloring (equivalent to the normal sequential loop + * order) for a given number of inner indexes. + * \note This is to reduce overhead in "OpenMP-ready" code when only 1 thread is used. + * \param[in] numInnerIndexes - Number of indexes that are to be colored. + * \return Natural (sequential) coloring of the inner indices. + */ +template +T createNaturalColoring(Index_t numInnerIndexes) +{ + /*--- One color. ---*/ + su2vector outerPtr(2); + outerPtr(0) = 0; + outerPtr(1) = numInnerIndexes; + + /*--- Containing all indexes in ascending order. ---*/ + su2vector innerIdx(numInnerIndexes); + std::iota(innerIdx.data(), innerIdx.data()+numInnerIndexes, 0); + + return T(std::move(outerPtr), std::move(innerIdx)); +} + + /*! * \brief Color contiguous groups of outer indices of a sparse pattern such that * within each color, any two groups do not have inner indices in common. @@ -445,6 +469,10 @@ T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors const Index_t grpSz = groupSize; const Index_t nOuter = pattern.getOuterSize(); + + /*--- Trivial case. ---*/ + if(groupSize >= nOuter) return createNaturalColoring(nOuter); + const Index_t minIdx = pattern.getMinInnerIdx(); const Index_t nInner = pattern.getMaxInnerIdx()+1-minIdx; @@ -550,30 +578,6 @@ T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors } -/*! - * \brief Create the natural coloring (equivalent to the normal sequential loop - * order) for a given number of inner indexes. - * \note This is to reduce overhead in "OpenMP-ready" code when only 1 thread is used. - * \param[in] numInnerIndexes - Number of indexes that are to be colored. - * \return Natural (sequential) coloring of the inner indices. - */ -template -T createNaturalColoring(Index_t numInnerIndexes) -{ - /*--- One color. ---*/ - su2vector outerPtr(2); - outerPtr(0) = 0; - outerPtr(1) = numInnerIndexes; - - /*--- Containing all indexes in ascending order. ---*/ - su2vector innerIdx(numInnerIndexes); - std::iota(innerIdx.data(), innerIdx.data()+numInnerIndexes, 0); - - return T(std::move(outerPtr), std::move(innerIdx)); -} - - /*! * \brief A way to represent one grid color that allows range-for syntax. */ diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index 05746f165dc4..830d49b4fdcd 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -624,13 +624,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry **geometry, CConfig *config_con Local_nPointCoarse = nPoint; Local_nPointFine = fine_grid->GetnPoint(); -#ifdef HAVE_MPI SU2_MPI::Allreduce(&Local_nPointCoarse, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); SU2_MPI::Allreduce(&Local_nPointFine, &Global_nPointFine, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); -#else - Global_nPointCoarse = Local_nPointCoarse; - Global_nPointFine = Local_nPointFine; -#endif su2double Coeff = 1.0, CFL = 0.0, factor = 1.5; @@ -670,6 +665,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry **geometry, CConfig *config_con } } + edgeColorGroupSize = config->GetEdgeColoringGroupSize(); + delete [] copy_marker; } diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index cf4375463fb8..8f8ea74fe17a 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -468,15 +468,11 @@ CPhysicalGeometry::CPhysicalGeometry(CGeometry *geometry, nLocal_Pris + nLocal_Pyra); nLocal_Bound_Elem = nLocal_Line + nLocal_BoundTria + nLocal_BoundQuad; -#ifndef HAVE_MPI - nGlobal_Elem = nLocal_Elem; - nGlobal_Bound_Elem = nLocal_Bound_Elem; -#else + SU2_MPI::Allreduce(&nLocal_Elem, &nGlobal_Elem, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); SU2_MPI::Allreduce(&nLocal_Bound_Elem, &nGlobal_Bound_Elem, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); -#endif /*--- With the distribution of all points, elements, and markers based on the ParMETIS coloring complete, as a final step, load this data into diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 6d97a3518e98..2a47d4785a60 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1319,6 +1319,8 @@ void CSysMatrix::SetDiagonalAsColumnSum() { auto block_ii = &matrix[dia_ptr[iPoint]*nVar*nEqn]; + for (auto k = 0ul; k < nVar*nEqn; ++k) block_ii[k] = 0.0; + for (auto k = row_ptr[iPoint]; k < row_ptr[iPoint+1]; ++k) { auto block_ji = &matrix[col_ptr[k]*nVar*nEqn]; if (block_ji != block_ii) MatrixSubtraction(block_ii, block_ji, block_ii); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index cd70602a06f9..a1316b3db04d 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -241,7 +241,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, cout << "WARNING: On " << numRanksUsingReducer << " MPI ranks the coloring efficiency was less than " << COLORING_EFF_THRESH << " (min value was " << minEff << ").\n" << " Those ranks will now use a fallback strategy, better performance may be possible\n" - << " with different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl; + << " with a different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl; } if (ReducerStrategy) @@ -2627,9 +2627,14 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container } } - /*--- Initialize the Jacobian matrices ---*/ + /*--- Initialize the Jacobian matrix and residual, not needed for the reducer strategy + * as we set blocks (including diagonal ones) and completely overwrite. ---*/ - if (implicit && !disc_adjoint) Jacobian.SetValZero(); + if(!ReducerStrategy && !Output) { + LinSysRes.SetValZero(); + if (implicit && !config->GetDiscrete_Adjoint()) Jacobian.SetValZero(); + else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes + } } @@ -2650,11 +2655,6 @@ unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, C /* Check for non-realizable states for reporting. */ if (!physical) nonPhysicalPoints++; - - /*--- Initialize the convective, source and viscous residual vector ---*/ - - if (!Output) LinSysRes.SetBlock_Zero(iPoint); - } return nonPhysicalPoints; @@ -2981,13 +2981,11 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta if (ReducerStrategy) { EdgeFluxes.SetBlock(iEdge, residual); if (implicit) - Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); } else { LinSysRes.AddBlock(iPoint, residual); LinSysRes.SubtractBlock(jPoint, residual); - - /*--- Set implicit computation ---*/ if (implicit) Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } @@ -3209,7 +3207,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain if (ReducerStrategy) { EdgeFluxes.SetBlock(iEdge, residual); if (implicit) - Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); } else { LinSysRes.AddBlock(iPoint, residual); @@ -3258,6 +3256,8 @@ void CEulerSolver::SumEdgeFluxes(CGeometry* geometry) { SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + LinSysRes.SetBlock_Zero(iPoint); + for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) { auto iEdge = geometry->node[iPoint]->GetEdge(iNeigh); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 9d0d57494a24..0b1185e6ac73 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -342,9 +342,14 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C SU2_OMP_BARRIER } - /*--- Initialize the Jacobian matrices ---*/ + /*--- Initialize the Jacobian matrix and residual, not needed for the reducer strategy + * as we set blocks (including diagonal ones) and completely overwrite. ---*/ - if (implicit && !config->GetDiscrete_Adjoint()) Jacobian.SetValZero(); + if(!ReducerStrategy && !Output) { + LinSysRes.SetValZero(); + if (implicit && !config->GetDiscrete_Adjoint()) Jacobian.SetValZero(); + else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes + } } @@ -383,10 +388,6 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, CCon nonPhysicalPoints += !physical; - /*--- Initialize the convective, source and viscous residual vector ---*/ - - if (!Output) LinSysRes.SetBlock_Zero(iPoint); - } return nonPhysicalPoints; diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index fae8274a64d9..9d78a28e2b13 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -106,6 +106,9 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + if (config->GetExtraOutput()) { if (nDim == 2) { nOutputVariables = 13; } else if (nDim == 3) { nOutputVariables = 19; } @@ -277,9 +280,12 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe const su2double* Vorticity = nullptr; su2double Laminar_Viscosity = 0.0; - /*--- Clear residual and system matrix. ---*/ - LinSysRes.SetValZero(); - Jacobian.SetValZero(); + /*--- Clear residual and system matrix, not needed for + * reducer strategy as we write over the entire matrix. ---*/ + if (!ReducerStrategy) { + LinSysRes.SetValZero(); + Jacobian.SetValZero(); + } /*--- Upwind second order reconstruction and gradients ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 3a87faa570b2..41fc01aa4890 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -108,6 +108,9 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + if (ReducerStrategy) + EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); + /*--- Initialize the BGS residuals in multizone problems. ---*/ if (multizone){ Residual_BGS = new su2double[nVar](); @@ -274,9 +277,12 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain const bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (config->GetInnerIter() <= config->GetLimiterIter()); - /*--- Clear residual and system matrix. ---*/ - LinSysRes.SetValZero(); - Jacobian.SetValZero(); + /*--- Clear residual and system matrix, not needed for + * reducer strategy as we write over the entire matrix. ---*/ + if (!ReducerStrategy) { + LinSysRes.SetValZero(); + Jacobian.SetValZero(); + } /*--- Upwind second order reconstruction and gradients ---*/ diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 6cff4b71cd12..024912da8d25 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -67,8 +67,6 @@ CTurbSolver::CTurbSolver(CGeometry* geometry, CConfig *config) : CSolver() { /*--- Local (to the rank) decision to use the reducer strategy, any * warnings have been, or will be, printed by the flow solver. ---*/ ReducerStrategy = parallelEff < COLORING_EFF_THRESH; - if (ReducerStrategy) - EdgeFluxes.Initialize(geometry->GetnEdge(), geometry->GetnEdge(), nVar, nullptr); nPoint = geometry->GetnPoint(); omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); @@ -222,7 +220,7 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe if (ReducerStrategy) { EdgeFluxes.SetBlock(iEdge, residual); - Jacobian.UpdateBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); + Jacobian.SetBlocks(iEdge, residual.jacobian_i, residual.jacobian_j); } else { LinSysRes.AddBlock(iPoint, residual); @@ -296,6 +294,8 @@ void CTurbSolver::SumEdgeFluxes(CGeometry* geometry) { SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + LinSysRes.SetBlock_Zero(iPoint); + for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) { auto iEdge = geometry->node[iPoint]->GetEdge(iNeigh); From bcf6af496082eb27c6acf5b3bf35b88213ff48bc Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 19:03:15 +0000 Subject: [PATCH 14/28] small regression changes --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6d15e0560ace..075fdae3d93b 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -243,7 +243,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-12.078361, -16.147829, 1.064326, 0.019770] #last 4 columns + turb_naca0012_sa.test_vals = [-12.078401, -16.147829, 1.064326, 0.019770] #last 4 columns turb_naca0012_sa.su2_exec = "parallel_computation.py -f" turb_naca0012_sa.timeout = 3200 turb_naca0012_sa.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 09657a95b35c..c53f1bfdc723 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -281,7 +281,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 10 - turb_naca0012_sa.test_vals = [-12.075893, -16.146770, 1.064326, 0.019770] #last 4 columns + turb_naca0012_sa.test_vals = [-12.075861, -16.146770, 1.064326, 0.019770] #last 4 columns turb_naca0012_sa.su2_exec = "SU2_CFD" turb_naca0012_sa.new_output = True turb_naca0012_sa.timeout = 3200 @@ -790,7 +790,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.003327, 1.312032, 6.085201, 2.413460] #last 4 columns + turb_naca0012_p1c1.test_vals = [-5.003335, 1.312021, 6.085201, 2.413460] #last 4 columns turb_naca0012_p1c1.su2_exec = "SU2_CFD" turb_naca0012_p1c1.new_output = True turb_naca0012_p1c1.timeout = 1600 @@ -802,7 +802,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.263992, 1.251332, 6.085705, 2.413434] #last 4 columns + turb_naca0012_p1c2.test_vals = [-5.264009, 1.251324, 6.085705, 2.413434] #last 4 columns turb_naca0012_p1c2.su2_exec = "SU2_CFD" turb_naca0012_p1c2.new_output = True turb_naca0012_p1c2.timeout = 1600 From 19f2fe5818d12fa88792ceed35ace060cf4aee18 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 19:04:12 +0000 Subject: [PATCH 15/28] write Undivided Laplacian as point loop --- SU2_CFD/include/variables/CVariable.hpp | 16 ++------ SU2_CFD/src/solvers/CEulerSolver.cpp | 52 ++++++++----------------- SU2_CFD/src/solvers/CSolver.cpp | 2 +- 3 files changed, 21 insertions(+), 49 deletions(-) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index c9da45da7762..5e28e972abc8 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1076,16 +1076,6 @@ class CVariable { */ inline virtual su2double GetSensor(unsigned long iPoint, unsigned long iSpecies) const { return 0.0; } - /*! - * \brief Set the value of the undivided laplacian of the solution. - * \param[in] iPoint - Point index. - * \param[in] iVar - Index of the variable. - * \param[in] val_undivided_laplacian - Value of the undivided solution for the index iVar. - */ - inline void SetUndivided_Laplacian(unsigned long iPoint, unsigned long iVar, su2double val_undivided_laplacian) { - Undivided_Laplacian(iPoint,iVar) = val_undivided_laplacian; - } - /*! * \brief Add the value of the undivided laplacian of the solution. * \param[in] iPoint - Point index. @@ -1107,13 +1097,13 @@ class CVariable { } /*! - * \brief Subtract the value of the undivided laplacian of the solution. + * \brief Increment the value of the undivided laplacian of the solution. * \param[in] iPoint - Point index. * \param[in] iVar - Variable of the undivided laplacian. * \param[in] val_und_lapl - Value of the undivided solution. */ - inline void SubtractUnd_Lapl(unsigned long iPoint, unsigned long iVar, su2double val_und_lapl) { - Undivided_Laplacian(iPoint, iVar) -= val_und_lapl; + inline void AddUnd_Lapl(unsigned long iPoint, unsigned long iVar, su2double val_und_lapl) { + Undivided_Laplacian(iPoint, iVar) += val_und_lapl; } /*! diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index a1316b3db04d..c3da31843ce3 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3647,53 +3647,35 @@ void CEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) { void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) { - nodes->SetUnd_LaplZero(); - - /*--- Loop over edge colors. ---*/ - for (auto color : EdgeColoring) - { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) - for(auto k = 0ul; k < color.size; ++k) { - - auto iEdge = color.indices[k]; - - auto iPoint = geometry->edge[iEdge]->GetNode(0); - auto jPoint = geometry->edge[iEdge]->GetNode(1); + /*--- Loop domain points. ---*/ - /*--- Solution differences ---*/ + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - su2double Diff[MAXNVAR] = {0.0}; + const bool boundary_i = geometry->node[iPoint]->GetPhysicalBoundary(); + /*--- Initialize. ---*/ for (unsigned short iVar = 0; iVar < nVar; iVar++) - Diff[iVar] = nodes->GetSolution(iPoint,iVar) - nodes->GetSolution(jPoint,iVar); + nodes->SetUnd_Lapl(iPoint, iVar, 0.0); - /*--- Correction for compressible flows which use the enthalpy ---*/ + /*--- Loop over the neighbors of point i. ---*/ + for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) + { + auto jPoint = geometry->node[iPoint]->GetPoint(iNeigh); + bool boundary_j = geometry->node[jPoint]->GetPhysicalBoundary(); - Diff[nVar-1] += nodes->GetPressure(iPoint) - nodes->GetPressure(jPoint); + /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ + if (boundary_i && !boundary_j) continue; - bool boundary_i = geometry->node[iPoint]->GetPhysicalBoundary(); - bool boundary_j = geometry->node[jPoint]->GetPhysicalBoundary(); + /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ - /*--- Both points inside the domain, or both in the boundary ---*/ + for (unsigned short iVar = 0; iVar < nVar; iVar++) + nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)); - if ((!boundary_i && !boundary_j) || (boundary_i && boundary_j)) { - if (geometry->node[iPoint]->GetDomain()) nodes->SubtractUnd_Lapl(iPoint, Diff); - if (geometry->node[jPoint]->GetDomain()) nodes->AddUnd_Lapl(jPoint, Diff); + nodes->AddUnd_Lapl(iPoint, nVar-1, nodes->GetPressure(jPoint)-nodes->GetPressure(iPoint)); } - /*--- iPoint inside the domain, jPoint on the boundary ---*/ - - if (!boundary_i && boundary_j) - if (geometry->node[iPoint]->GetDomain()) nodes->SubtractUnd_Lapl(iPoint, Diff); - - /*--- jPoint inside the domain, iPoint on the boundary ---*/ - - if (boundary_i && !boundary_j) - if (geometry->node[jPoint]->GetDomain()) nodes->AddUnd_Lapl(jPoint, Diff); - } - } // end color loop SU2_OMP_MASTER { diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index a4735cec3c17..112a912fddcb 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -2219,7 +2219,7 @@ void CSolver::CompleteComms(CGeometry *geometry, break; case UNDIVIDED_LAPLACIAN: for (iVar = 0; iVar < nVar; iVar++) - base_nodes->SetUndivided_Laplacian(iPoint, iVar, bufDRecv[buf_offset+iVar]); + base_nodes->SetUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; case SOLUTION_LIMITER: for (iVar = 0; iVar < nVar; iVar++) From 13e9572005e4a3cc542214755acfe7bdb6674545 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 19:17:53 +0000 Subject: [PATCH 16/28] write Centered_Dissipation_Sensor as a point loop --- SU2_CFD/src/solvers/CEulerSolver.cpp | 82 +++++++++++++--------------- 1 file changed, 37 insertions(+), 45 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index c3da31843ce3..9780ecb0716e 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3697,67 +3697,59 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig *config) { - /*--- Reset variables to store the undivided pressure ---*/ + /*--- We can access memory more efficiently if there are no periodic boundaries. ---*/ - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - iPoint_UndLapl[iPoint] = 0.0; - jPoint_UndLapl[iPoint] = 0.0; - } + const bool isPeriodic = (config->GetnMarker_Periodic() > 0); - /*--- Loop over edge colors. ---*/ - for (auto color : EdgeColoring) - { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) - for(auto k = 0ul; k < color.size; ++k) { + /*--- Loop domain points. ---*/ - auto iEdge = color.indices[k]; + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - auto iPoint = geometry->edge[iEdge]->GetNode(0); - auto jPoint = geometry->edge[iEdge]->GetNode(1); + const bool boundary_i = geometry->node[iPoint]->GetPhysicalBoundary(); + const su2double Pressure_i = nodes->GetPressure(iPoint); - su2double Pressure_i = nodes->GetPressure(iPoint); - su2double Pressure_j = nodes->GetPressure(jPoint); + /*--- Initialize. ---*/ + iPoint_UndLapl[iPoint] = 0.0; + jPoint_UndLapl[iPoint] = 0.0; - bool boundary_i = geometry->node[iPoint]->GetPhysicalBoundary(); - bool boundary_j = geometry->node[jPoint]->GetPhysicalBoundary(); + /*--- Loop over the neighbors of point i. ---*/ + for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) + { + auto jPoint = geometry->node[iPoint]->GetPoint(iNeigh); + bool boundary_j = geometry->node[jPoint]->GetPhysicalBoundary(); - /*--- Both points inside the domain, or both on the boundary ---*/ + /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ + if (boundary_i && !boundary_j) continue; - if ((!boundary_i && !boundary_j) || (boundary_i && boundary_j)) { - if (geometry->node[iPoint]->GetDomain()) { iPoint_UndLapl[iPoint] += (Pressure_j - Pressure_i); jPoint_UndLapl[iPoint] += (Pressure_i + Pressure_j); } - if (geometry->node[jPoint]->GetDomain()) { iPoint_UndLapl[jPoint] += (Pressure_i - Pressure_j); jPoint_UndLapl[jPoint] += (Pressure_i + Pressure_j); } + /*--- Add pressure difference and pressure sum. ---*/ + su2double Pressure_j = nodes->GetPressure(jPoint); + iPoint_UndLapl[iPoint] += Pressure_j - Pressure_i; + jPoint_UndLapl[iPoint] += Pressure_j + Pressure_i; } - /*--- iPoint inside the domain, jPoint on the boundary ---*/ - - if (!boundary_i && boundary_j) - if (geometry->node[iPoint]->GetDomain()) { iPoint_UndLapl[iPoint] += (Pressure_j - Pressure_i); jPoint_UndLapl[iPoint] += (Pressure_i + Pressure_j); } - - /*--- jPoint inside the domain, iPoint on the boundary ---*/ - - if (boundary_i && !boundary_j) - if (geometry->node[jPoint]->GetDomain()) { iPoint_UndLapl[jPoint] += (Pressure_i - Pressure_j); jPoint_UndLapl[jPoint] += (Pressure_i + Pressure_j); } + if (!isPeriodic) + nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); } - } // end color loop - /*--- Correct the sensor values across any periodic boundaries. ---*/ + if (isPeriodic) { + /*--- Correct the sensor values across any periodic boundaries. ---*/ - SU2_OMP_MASTER - { - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); + SU2_OMP_MASTER + { + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); + } } - } - SU2_OMP_BARRIER + SU2_OMP_BARRIER - /*--- Set pressure switch for each point ---*/ + /*--- Set pressure switch for each point ---*/ - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) - nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); + } /*--- MPI parallelization ---*/ From 733aef15e2c372f4394bd793f9703f62f8b7feb8 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 4 Mar 2020 19:48:37 +0000 Subject: [PATCH 17/28] allow forcing of "reducer strategy" without warnings, fix some indentation in CConfig.cpp --- Common/include/option_structure.hpp | 2 +- Common/src/CConfig.cpp | 83 ++++++++++++++-------------- SU2_CFD/src/solvers/CEulerSolver.cpp | 27 +++++---- 3 files changed, 59 insertions(+), 53 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1a9430e06ab4..1387b3227064 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -133,7 +133,7 @@ const int SU2_CONN_SIZE = 10; /*!< \brief Size of the connectivity array that that we read from a mesh file in the format [[globalID vtkType n0 n1 n2 n3 n4 n5 n6 n7 n8]. */ const int SU2_CONN_SKIP = 2; /*!< \brief Offset to skip the globalID and VTK type at the start of the element connectivity list for each CGNS element. */ -const su2double COLORING_EFF_THRESH = 0.8; /*!< \brief Below this value the reducer strategy is used instead. */ +const su2double COLORING_EFF_THRESH = 0.875; /*!< \brief Below this value fallback strategies are used instead. */ /*! * \brief Boolean answers diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 283c0768abd5..84dd8cef7eae 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2807,7 +2807,7 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Level of fill for PaStiX incomplete LU factorization. */ addUnsignedShortOption("PASTIX_FILL_LEVEL", pastix_fill_lvl, 1); - /* DESCRIPTION: Size of the edge groups colored for OpenMP parallelization of edge loops. */ + /* DESCRIPTION: Size of the edge groups colored for thread parallel edge loops (0 forces the reducer strategy). */ addUnsignedLongOption("EDGE_COLORING_GROUP_SIZE", edgeColorGroupSize, 512); /* END_CONFIG_OPTIONS */ @@ -3677,96 +3677,96 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerMotion_Origin[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerMotion_Origin/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_MOTION_ORIGIN must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerTranslation == 0){ nMarkerTranslation = 3*nMarker_Moving; MarkerTranslation_Rate = new su2double[nMarkerTranslation]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerTranslation_Rate[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerTranslation/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_TRANSLATION_RATE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerRotation_Rate == 0){ nMarkerRotation_Rate = 3*nMarker_Moving; MarkerRotation_Rate = new su2double[nMarkerRotation_Rate]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerRotation_Rate[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerRotation_Rate/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_ROTATION_RATE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerPlunging_Ampl == 0){ nMarkerPlunging_Ampl = 3*nMarker_Moving; MarkerPlunging_Ampl = new su2double[nMarkerPlunging_Ampl]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerPlunging_Ampl[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerPlunging_Ampl/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PLUNGING_AMPL must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerPlunging_Omega == 0){ nMarkerPlunging_Omega = 3*nMarker_Moving; MarkerPlunging_Omega = new su2double[nMarkerPlunging_Omega]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerPlunging_Omega[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerPlunging_Omega/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PLUNGING_OMEGA must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerPitching_Ampl == 0){ nMarkerPitching_Ampl = 3*nMarker_Moving; MarkerPitching_Ampl = new su2double[nMarkerPitching_Ampl]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerPitching_Ampl[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerPitching_Ampl/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_AMPL must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerPitching_Omega == 0){ nMarkerPitching_Omega = 3*nMarker_Moving; MarkerPitching_Omega = new su2double[nMarkerPitching_Omega]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerPitching_Omega[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerPitching_Omega/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_OMEGA must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMarkerPitching_Phase == 0){ nMarkerPitching_Phase = 3*nMarker_Moving; MarkerPitching_Phase = new su2double[nMarkerPitching_Phase]; for (iMarker = 0; iMarker < nMarker_Moving; iMarker++){ for (iDim = 0; iDim < 3; iDim++){ MarkerPitching_Phase[3*iMarker+iDim] = 0.0; - } - } + } + } } if (nMarkerPitching_Phase/3 != nMarker_Moving){ SU2_MPI::Error("Number of SURFACE_PITCHING_PHASE must be three times the number of MARKER_MOVING, (x,y,z) per marker.", CURRENT_FUNCTION); - } + } if (nMoveMotion_Origin == 0){ nMoveMotion_Origin = nMarker_Moving; @@ -4289,39 +4289,39 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } if (nElasticityMod == 0) { - nElasticityMod = 1; - ElasticityMod = new su2double[1]; ElasticityMod[0] = 2E11; + nElasticityMod = 1; + ElasticityMod = new su2double[1]; ElasticityMod[0] = 2E11; } if (nPoissonRatio == 0) { - nPoissonRatio = 1; - PoissonRatio = new su2double[1]; PoissonRatio[0] = 0.30; + nPoissonRatio = 1; + PoissonRatio = new su2double[1]; PoissonRatio[0] = 0.30; } if (nMaterialDensity == 0) { - nMaterialDensity = 1; - MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854; + nMaterialDensity = 1; + MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854; } if (nElectric_Constant == 0) { - nElectric_Constant = 1; - Electric_Constant = new su2double[1]; Electric_Constant[0] = 0.0; + nElectric_Constant = 1; + Electric_Constant = new su2double[1]; Electric_Constant[0] = 0.0; } if (nElectric_Field == 0) { - nElectric_Field = 1; - Electric_Field_Mod = new su2double[1]; Electric_Field_Mod[0] = 0.0; + nElectric_Field = 1; + Electric_Field_Mod = new su2double[1]; Electric_Field_Mod[0] = 0.0; } if (nDim_RefNode == 0) { - nDim_RefNode = 3; - RefNode_Displacement = new su2double[3]; - RefNode_Displacement[0] = 0.0; RefNode_Displacement[1] = 0.0; RefNode_Displacement[2] = 0.0; + nDim_RefNode = 3; + RefNode_Displacement = new su2double[3]; + RefNode_Displacement[0] = 0.0; RefNode_Displacement[1] = 0.0; RefNode_Displacement[2] = 0.0; } if (nDim_Electric_Field == 0) { - nDim_Electric_Field = 2; - Electric_Field_Dir = new su2double[2]; Electric_Field_Dir[0] = 0.0; Electric_Field_Dir[1] = 1.0; + nDim_Electric_Field = 2; + Electric_Field_Dir = new su2double[2]; Electric_Field_Dir[0] = 0.0; Electric_Field_Dir[1] = 1.0; } if ((Kind_SU2 == SU2_CFD) && (Kind_Solver == NO_SOLVER)) { @@ -4920,6 +4920,9 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ CURRENT_FUNCTION); } + /*--- 0 in the config file means "disable" which can be done using a very large group. ---*/ + if (edgeColorGroupSize==0) edgeColorGroupSize = 1<<30; + } void CConfig::SetMarkers(unsigned short val_software) { diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 9780ecb0716e..74bdbe570652 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -230,18 +230,21 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, /*--- The decision to use the strategy is local to each rank. ---*/ ReducerStrategy = parallelEff < COLORING_EFF_THRESH; - /*--- The warning messages are global. ---*/ - su2double minEff = 1.0; - SU2_MPI::Reduce(¶llelEff, &minEff, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - - int tmp = ReducerStrategy, numRanksUsingReducer = 0; - SU2_MPI::Reduce(&tmp, &numRanksUsingReducer, 1, MPI_INT, MPI_SUM, MASTER_NODE, MPI_COMM_WORLD); - - if (minEff < COLORING_EFF_THRESH) { - cout << "WARNING: On " << numRanksUsingReducer << " MPI ranks the coloring efficiency was less than " - << COLORING_EFF_THRESH << " (min value was " << minEff << ").\n" - << " Those ranks will now use a fallback strategy, better performance may be possible\n" - << " with a different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl; + /*--- If the reducer strategy is not being forced (by EDGE_COLORING_GROUP_SIZE=0) print some messages. ---*/ + if (config->GetEdgeColoringGroupSize() != 1<<30) { + + su2double minEff = 1.0; + SU2_MPI::Reduce(¶llelEff, &minEff, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); + + int tmp = ReducerStrategy, numRanksUsingReducer = 0; + SU2_MPI::Reduce(&tmp, &numRanksUsingReducer, 1, MPI_INT, MPI_SUM, MASTER_NODE, MPI_COMM_WORLD); + + if (minEff < COLORING_EFF_THRESH) { + cout << "WARNING: On " << numRanksUsingReducer << " MPI ranks the coloring efficiency was less than " + << COLORING_EFF_THRESH << " (min value was " << minEff << ").\n" + << " Those ranks will now use a fallback strategy, better performance may be possible\n" + << " with a different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl; + } } if (ReducerStrategy) From debb952a95f6653c9a752d22b166175f0fea7c1c Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 5 Mar 2020 08:55:35 +0000 Subject: [PATCH 18/28] update fixedcl regression --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 075fdae3d93b..82b58fbf9076 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -91,7 +91,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-12.137879, -6.705803, 0.300000, 0.019470] #last 4 columns + fixedCL_naca0012.test_vals = [-12.137437, -6.705109, 0.300000, 0.019470] #last 4 columns fixedCL_naca0012.su2_exec = "parallel_computation.py -f" fixedCL_naca0012.timeout = 1600 fixedCL_naca0012.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c53f1bfdc723..2b2ee1dd9ae8 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -101,7 +101,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-12.129044, -6.702294, 0.300000, 0.019470] #last 4 columns + fixedCL_naca0012.test_vals = [-12.128275, -6.700329, 0.300000, 0.019470] #last 4 columns fixedCL_naca0012.su2_exec = "SU2_CFD" fixedCL_naca0012.new_output = True fixedCL_naca0012.timeout = 1600 From b406fa3d55f9f526dd91329962d316c561d36a17 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 5 Mar 2020 10:16:57 +0000 Subject: [PATCH 19/28] methods to set natural colorings --- Common/include/geometry/CGeometry.hpp | 10 ++++++ Common/src/geometry/CGeometry.cpp | 45 +++++++++++++++++---------- 2 files changed, 38 insertions(+), 17 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index b8a23510b5f8..4473ac72b72a 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1645,6 +1645,11 @@ class CGeometry { */ const CCompressedSparsePatternUL& GetEdgeColoring(su2double* efficiency = nullptr); + /*! + * \brief Force the natural (sequential) edge coloring. + */ + void SetNaturalEdgeColoring(); + /*! * \brief Get the group size used in edge coloring. * \return Group size. @@ -1659,6 +1664,11 @@ class CGeometry { */ const CCompressedSparsePatternUL& GetElementColoring(su2double* efficiency = nullptr); + /*! + * \brief Force the natural (sequential) element coloring. + */ + void SetNaturalElementColoring(); + /*! * \brief Get the group size used in element coloring. * \return Group size. diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 5b1710f7af9d..a7a502b30b63 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -4046,14 +4046,15 @@ const su2vector& CGeometry::GetTransposeSparsePatternMap(Connecti const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficiency) { + /*--- Check for dry run mode with dummy geometry. ---*/ if (nEdge==0) return edgeColoring; /*--- Build if required. ---*/ if (edgeColoring.empty()) { - /*--- When not using threading use the natural coloring. ---*/ + /*--- When not using threading use the natural coloring to reduce overhead. ---*/ if (omp_get_max_threads() == 1) { - edgeColoring = createNaturalColoring(nEdge); + SetNaturalEdgeColoring(); if (efficiency != nullptr) *efficiency = 1.0; // by definition return edgeColoring; } @@ -4073,17 +4074,13 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficien CCompressedSparsePatternUL pattern(move(outerPtr), move(innerIdx)); /*--- Color the edges. ---*/ - const bool balanceColors = true; + constexpr bool balanceColors = true; edgeColoring = colorSparsePattern(pattern, edgeColorGroupSize, balanceColors); - /*--- If the coloring fails use the natural coloring and set the group size - * to nEdge to prevent client code from looping in parallel. This is a - * "soft" failure as the "bad" coloring should be detected downstream - * and a fallback strategy put in place. ---*/ - if (edgeColoring.empty()) { - edgeColoring = createNaturalColoring(nEdge); - edgeColorGroupSize = nEdge; - } + /*--- If the coloring fails use the natural coloring. This is a + * "soft" failure as this "bad" coloring should be detected + * downstream and a fallback strategy put in place. ---*/ + if (edgeColoring.empty()) SetNaturalEdgeColoring(); } if (efficiency != nullptr) { @@ -4092,8 +4089,17 @@ const CCompressedSparsePatternUL& CGeometry::GetEdgeColoring(su2double* efficien return edgeColoring; } +void CGeometry::SetNaturalEdgeColoring() +{ + if (nEdge == 0) return; + edgeColoring = createNaturalColoring(nEdge); + /*--- In parallel, set the group size to nEdge to protect client code. ---*/ + if (omp_get_max_threads() > 1) edgeColorGroupSize = nEdge; +} + const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* efficiency) { + /*--- Check for dry run mode with dummy geometry. ---*/ if (nElem==0) return elemColoring; /*--- Build if required. ---*/ @@ -4101,7 +4107,7 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* effic /*--- When not using threading use the natural coloring. ---*/ if (omp_get_max_threads() == 1) { - elemColoring = createNaturalColoring(nElem); + SetNaturalElementColoring(); if (efficiency != nullptr) *efficiency = 1.0; // by definition return elemColoring; } @@ -4123,14 +4129,11 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* effic CCompressedSparsePatternUL pattern(outerPtr, innerIdx); /*--- Color the elements. ---*/ - const bool balanceColors = true; + constexpr bool balanceColors = true; elemColoring = colorSparsePattern(pattern, elemColorGroupSize, balanceColors); /*--- Same as for the edge coloring. ---*/ - if (elemColoring.empty()) { - elemColoring = createNaturalColoring(nElem); - elemColorGroupSize = nElem; - } + if (elemColoring.empty()) SetNaturalElementColoring(); } if (efficiency != nullptr) { @@ -4138,3 +4141,11 @@ const CCompressedSparsePatternUL& CGeometry::GetElementColoring(su2double* effic } return elemColoring; } + +void CGeometry::SetNaturalElementColoring() +{ + if (nElem == 0) return; + elemColoring = createNaturalColoring(nElem); + /*--- In parallel, set the group size to nElem to protect client code. ---*/ + if (omp_get_max_threads() > 1) elemColorGroupSize = nElem; +} From fb500f226f9fa90fabc54cea775ba1a9840c9d25 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 5 Mar 2020 10:17:53 +0000 Subject: [PATCH 20/28] fix FEA solver lock strategy --- SU2_CFD/src/solvers/CFEASolver.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 81a1990a19b6..f1a7104a235f 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -335,15 +335,22 @@ void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { su2double parallelEff = 1.0; const auto& coloring = geometry->GetElementColoring(¶llelEff); + /*--- If the coloring is too bad use lock-guarded accesses + * to CSysMatrix/Vector in element loops instead. ---*/ + LockStrategy = parallelEff < COLORING_EFF_THRESH; + + /*--- When using locks force a single color to reduce the color loop overhead. ---*/ + if (LockStrategy && (coloring.getOuterSize()>1)) + geometry->SetNaturalElementColoring(); + if (!coloring.empty()) { + /*--- We are not constrained by the color group size when using locks. ---*/ + auto groupSize = LockStrategy? 1ul : geometry->GetElementColorGroupSize(); auto nColor = coloring.getOuterSize(); ElemColoring.reserve(nColor); - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - ElemColoring.emplace_back(coloring.innerIdx(iColor), - coloring.getNumNonZeros(iColor), - geometry->GetElementColorGroupSize()); - } + for(auto iColor = 0ul; iColor < nColor; ++iColor) + ElemColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor), groupSize); } su2double minEff = 1.0; @@ -358,8 +365,6 @@ void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { UpdateLocks.resize(nPoint); for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) omp_init_lock(&UpdateLocks[iPoint]); - /*--- We are no longer constrained by the color group size. ---*/ - for(auto& color : ElemColoring) color.groupSize = OMP_MIN_SIZE; } omp_chunk_size = computeStaticChunkSize(nPointDomain, omp_get_max_threads(), OMP_MAX_SIZE); From db8d771398207e81a2acf84aaba040f880321dae Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 5 Mar 2020 10:18:19 +0000 Subject: [PATCH 21/28] make overhead of reducer strategy due to bad coloring same as when coloring fails --- SU2_CFD/include/solvers/CTurbSolver.hpp | 31 ++++++++++---------- SU2_CFD/src/solvers/CEulerSolver.cpp | 28 ++++++++++-------- SU2_CFD/src/solvers/CTurbSolver.cpp | 38 ++++++++++--------------- 3 files changed, 47 insertions(+), 50 deletions(-) diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index 5056e0dd9ec8..ce4adf446c8f 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -84,6 +84,22 @@ class CTurbSolver : public CSolver { */ inline CVariable* GetBaseClassPointerToNodes() final { return nodes; } +private: + + /*! + * \brief Compute the viscous flux for the turbulent equation at a particular edge. + * \param[in] iEdge - Edge for which we want to compute the flux + * \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. + */ + void Viscous_Residual(unsigned long iEdge, + CGeometry *geometry, + CSolver **solver_container, + CNumerics *numerics, + CConfig *config); + /*! * \brief Sum the edge fluxes for each cell to populate the residual vector, only used on coarse grids. * \param[in] geometry - Geometrical definition of the problem. @@ -117,27 +133,12 @@ class CTurbSolver : public CSolver { * \param[in] config - Definition of the particular problem. * \param[in] iMesh - Index of the mesh in multigrid computations. */ - void Upwind_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) override; - /*! - * \brief Compute the viscous flux for the turbulent equation at a particular edge. - * \param[in] iEdge - Edge for which we want to compute the flux - * \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. - */ - void Viscous_Residual(unsigned long iEdge, - CGeometry *geometry, - CSolver **solver_container, - CNumerics *numerics, - CConfig *config); - /*! * \brief Impose the Symmetry Plane boundary condition. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 74bdbe570652..93c2a47c5719 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -216,20 +216,24 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, su2double parallelEff = 1.0; const auto& coloring = geometry->GetEdgeColoring(¶llelEff); + /*--- The decision to use the strategy is local to each rank. ---*/ + ReducerStrategy = parallelEff < COLORING_EFF_THRESH; + + /*--- When using the reducer force a single color to reduce the color loop overhead. ---*/ + if (ReducerStrategy && (coloring.getOuterSize()>1)) + geometry->SetNaturalEdgeColoring(); + if (!coloring.empty()) { + /*--- If the reducer strategy is used we are not constrained by group + * size as we have no other edge loops in the Euler/NS solvers. ---*/ + auto groupSize = ReducerStrategy? 1ul : geometry->GetEdgeColorGroupSize(); auto nColor = coloring.getOuterSize(); EdgeColoring.reserve(nColor); - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - EdgeColoring.emplace_back(coloring.innerIdx(iColor), - coloring.getNumNonZeros(iColor), - geometry->GetEdgeColorGroupSize()); - } + for(auto iColor = 0ul; iColor < nColor; ++iColor) + EdgeColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor), groupSize); } - /*--- The decision to use the strategy is local to each rank. ---*/ - ReducerStrategy = parallelEff < COLORING_EFF_THRESH; - /*--- If the reducer strategy is not being forced (by EDGE_COLORING_GROUP_SIZE=0) print some messages. ---*/ if (config->GetEdgeColoringGroupSize() != 1<<30) { @@ -2938,8 +2942,8 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -3040,8 +3044,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 024912da8d25..d814efdb7348 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -48,26 +48,24 @@ CTurbSolver::CTurbSolver(CGeometry* geometry, CConfig *config) : CSolver() { dynamic_grid = config->GetDynamic_Grid(); #ifdef HAVE_OMP - /*--- Get the edge coloring. ---*/ - + /*--- Get the edge coloring, see notes in CEulerSolver's constructor. ---*/ su2double parallelEff = 1.0; const auto& coloring = geometry->GetEdgeColoring(¶llelEff); + ReducerStrategy = parallelEff < COLORING_EFF_THRESH; + + if (ReducerStrategy && (coloring.getOuterSize()>1)) + geometry->SetNaturalEdgeColoring(); + if (!coloring.empty()) { + auto groupSize = ReducerStrategy? 1ul : geometry->GetEdgeColorGroupSize(); auto nColor = coloring.getOuterSize(); EdgeColoring.reserve(nColor); - for(auto iColor = 0ul; iColor < nColor; ++iColor) { - EdgeColoring.emplace_back(coloring.innerIdx(iColor), - coloring.getNumNonZeros(iColor), - geometry->GetEdgeColorGroupSize()); - } + for(auto iColor = 0ul; iColor < nColor; ++iColor) + EdgeColoring.emplace_back(coloring.innerIdx(iColor), coloring.getNumNonZeros(iColor), groupSize); } - /*--- Local (to the rank) decision to use the reducer strategy, any - * warnings have been, or will be, printed by the flow solver. ---*/ - ReducerStrategy = parallelEff < COLORING_EFF_THRESH; - nPoint = geometry->GetnPoint(); omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); #else @@ -110,8 +108,8 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe /*--- Loop over edge colors. ---*/ for (auto color : EdgeColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size (unless we use the reducer). ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, ReducerStrategy? 1 : color.groupSize)) + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) for(auto k = 0ul; k < color.size; ++k) { auto iEdge = color.indices[k]; @@ -705,14 +703,9 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con we will loop over the edges and boundaries to compute the GCL component of the dual time source term that depends on grid velocities. ---*/ - /*--- Loop over edge colors. ---*/ - for (auto color : EdgeColoring) - { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) - for(auto k = 0ul; k < color.size; ++k) { - - auto iEdge = color.indices[k]; + /// TODO: Make this a point loop. + SU2_OMP_MASTER + for(auto iEdge = 0ul; iEdge < geometry->GetnEdge(); ++iEdge) { /*--- Get indices for nodes i & j plus the face normal ---*/ @@ -771,9 +764,8 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con for (iVar = 0; iVar < nVar; iVar++) LinSysRes(jPoint,iVar) -= U_time_n[iVar]*Residual_GCL; } - } - } // end color loop + SU2_OMP_BARRIER /*--- Loop over the boundary edges ---*/ From f4dd41aae782b5b5cbd7714c3eed9b21397ad69d Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Thu, 5 Mar 2020 13:28:46 +0000 Subject: [PATCH 22/28] fuse "JST dissipation" loops, cleanup Euler/NS preprocessing --- SU2_CFD/include/solvers/CEulerSolver.hpp | 23 ++-- SU2_CFD/src/solvers/CEulerSolver.cpp | 151 ++++++++++------------- SU2_CFD/src/solvers/CNSSolver.cpp | 117 +++--------------- 3 files changed, 95 insertions(+), 196 deletions(-) diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 820a17ffb368..30a7c455fa73 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -302,23 +302,27 @@ class CEulerSolver : public CSolver { void SumEdgeFluxes(CGeometry* geometry); /*! - * \brief Update the AoA and freestream velocity at the farfield. + * \brief Preprocessing actions common to the Euler and NS solvers. * \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. - * \param[in] iMesh - current mesh level for the multigrid. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. * \param[in] Output - boolean to determine whether to print output. */ - void SetFarfield_AoA(CGeometry *geometry, CSolver **solver_container, - CConfig *config, unsigned short iMesh, bool Output); + void CommonPreprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output); /*! - * \brief Compute a pressure sensor switch. + * \brief Update the AoA and freestream velocity at the farfield. * \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. + * \param[in] iMesh - current mesh level for the multigrid. + * \param[in] Output - boolean to determine whether to print output. */ - void SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig *config); + void SetFarfield_AoA(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iMesh, bool Output); /*! * \brief Compute Ducros Sensor for Roe Dissipation. @@ -362,11 +366,13 @@ class CEulerSolver : public CSolver { void SetMax_Eigenvalue(CGeometry *geometry, CConfig *config); /*! - * \brief Compute the undivided laplacian for the solution, except the energy equation. + * \brief Compute the undivided laplacian for the solution and the + * dissipation sensor for centered schemes. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SetUndivided_Laplacian(CGeometry *geometry, CConfig *config); + void SetUndivided_Laplacian_And_Centered_Dissipation_Sensor(CGeometry *geometry, + CConfig *config); /*! * \brief A virtual member. @@ -440,7 +446,6 @@ class CEulerSolver : public CSolver { return sqrt(Vel2); } - /*! * \brief Compute the density multiply by energy at the infinity. * \return Value of the density multiply by energy at the infinity. diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 93c2a47c5719..c45bdc5ba616 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2518,22 +2518,18 @@ void CEulerSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_c } -void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, - unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { +void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - unsigned long InnerIter = config->GetInnerIter(); bool cont_adjoint = config->GetContinuous_Adjoint(); bool disc_adjoint = config->GetDiscrete_Adjoint(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); - bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); - bool center_jst = center && (config->GetKind_Centered_Flow() == JST); + bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); bool engine = ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)); bool actuator_disk = ((config->GetnMarker_ActDiskInlet() != 0) || (config->GetnMarker_ActDiskOutlet() != 0)); bool nearfield = (config->GetnMarker_NearFieldBound() != 0); bool fixed_cl = config->GetFixed_CL_Mode(); - bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; unsigned short kind_row_dissipation = config->GetKind_RoeLowDiss(); bool roe_low_dissipation = (kind_row_dissipation != NO_ROELOWDISS) && (config->GetKind_Upwind_Flow() == ROE || @@ -2596,33 +2592,12 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container SU2_OMP_BARRIER } - /*--- Upwind second order reconstruction ---*/ - - if ((muscl && !center) && (iMesh == MESH_0) && !Output) { - - /*--- Gradient computation for MUSCL reconstruction. ---*/ - - if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) - SetPrimitive_Gradient_GG(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - - /*--- Limiter computation ---*/ - - if (limiter && (iMesh == MESH_0) && !Output && !van_albada) - SetPrimitive_Limiter(geometry, config); - } - /*--- Artificial dissipation ---*/ if (center && !Output) { SetMax_Eigenvalue(geometry, config); - if ((center_jst) && (iMesh == MESH_0)) { - SetCentered_Dissipation_Sensor(geometry, config); - SetUndivided_Laplacian(geometry, config); - } + if (center_jst) + SetUndivided_Laplacian_And_Centered_Dissipation_Sensor(geometry, config); } /*--- Roe Low Dissipation Sensor ---*/ @@ -2645,6 +2620,44 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container } +void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { + + unsigned long InnerIter = config->GetInnerIter(); + bool cont_adjoint = config->GetContinuous_Adjoint(); + bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); + bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); + bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); + bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; + + /*--- Common preprocessing steps. ---*/ + + CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); + + /*--- Upwind second order reconstruction ---*/ + + if ((muscl && !center) && (iMesh == MESH_0) && !Output) { + + /*--- Gradient computation for MUSCL reconstruction. ---*/ + + switch (config->GetKind_Gradient_Method_Recon()) { + case GREEN_GAUSS: + SetPrimitive_Gradient_GG(geometry, config, true); break; + case LEAST_SQUARES: + case WEIGHTED_LEAST_SQUARES: + SetPrimitive_Gradient_LS(geometry, config, true); break; + default: + break; + } + + /*--- Limiter computation ---*/ + + if (limiter && (iMesh == MESH_0) && !Output && !van_albada) + SetPrimitive_Limiter(geometry, config); + } + +} + unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { /*--- Number of non-physical points, local to the thread, needs @@ -3652,57 +3665,7 @@ void CEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) { } -void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) { - - /*--- Loop domain points. ---*/ - - SU2_OMP_FOR_DYN(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - const bool boundary_i = geometry->node[iPoint]->GetPhysicalBoundary(); - - /*--- Initialize. ---*/ - for (unsigned short iVar = 0; iVar < nVar; iVar++) - nodes->SetUnd_Lapl(iPoint, iVar, 0.0); - - /*--- Loop over the neighbors of point i. ---*/ - for (unsigned short iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnPoint(); ++iNeigh) - { - auto jPoint = geometry->node[iPoint]->GetPoint(iNeigh); - bool boundary_j = geometry->node[jPoint]->GetPhysicalBoundary(); - - /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ - if (boundary_i && !boundary_j) continue; - - /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ - - for (unsigned short iVar = 0; iVar < nVar; iVar++) - nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)); - - nodes->AddUnd_Lapl(iPoint, nVar-1, nodes->GetPressure(jPoint)-nodes->GetPressure(iPoint)); - } - - } - - SU2_OMP_MASTER - { - /*--- Correct the Laplacian values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); - } - - /*--- MPI parallelization ---*/ - - InitiateComms(geometry, config, UNDIVIDED_LAPLACIAN); - CompleteComms(geometry, config, UNDIVIDED_LAPLACIAN); - } - SU2_OMP_BARRIER - -} - -void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig *config) { +void CEulerSolver::SetUndivided_Laplacian_And_Centered_Dissipation_Sensor(CGeometry *geometry, CConfig *config) { /*--- We can access memory more efficiently if there are no periodic boundaries. ---*/ @@ -3717,6 +3680,9 @@ void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig * const su2double Pressure_i = nodes->GetPressure(iPoint); /*--- Initialize. ---*/ + for (unsigned short iVar = 0; iVar < nVar; iVar++) + nodes->SetUnd_Lapl(iPoint, iVar, 0.0); + iPoint_UndLapl[iPoint] = 0.0; jPoint_UndLapl[iPoint] = 0.0; @@ -3729,8 +3695,15 @@ void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig * /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ if (boundary_i && !boundary_j) continue; - /*--- Add pressure difference and pressure sum. ---*/ + /*--- Add solution differences, with correction for compressible flows which use the enthalpy. ---*/ + + for (unsigned short iVar = 0; iVar < nVar; iVar++) + nodes->AddUnd_Lapl(iPoint, iVar, nodes->GetSolution(jPoint,iVar)-nodes->GetSolution(iPoint,iVar)); + su2double Pressure_j = nodes->GetPressure(jPoint); + nodes->AddUnd_Lapl(iPoint, nVar-1, Pressure_j-Pressure_i); + + /*--- Dissipation sensor, add pressure difference and pressure sum. ---*/ iPoint_UndLapl[iPoint] += Pressure_j - Pressure_i; jPoint_UndLapl[iPoint] += Pressure_j + Pressure_i; } @@ -3758,10 +3731,20 @@ void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig * nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); } - /*--- MPI parallelization ---*/ - SU2_OMP_MASTER { + /*--- Correct the Laplacian values across any periodic boundaries. ---*/ + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); + } + + /*--- MPI parallelization ---*/ + + InitiateComms(geometry, config, UNDIVIDED_LAPLACIAN); + CompleteComms(geometry, config, UNDIVIDED_LAPLACIAN); + InitiateComms(geometry, config, SENSOR); CompleteComms(geometry, config, SENSOR); } diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 0b1185e6ac73..d2031ca454ec 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -171,108 +171,28 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C unsigned long InnerIter = config->GetInnerIter(); bool cont_adjoint = config->GetContinuous_Adjoint(); - bool disc_adjoint = config->GetDiscrete_Adjoint(); - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); - bool center_jst = center && config->GetKind_Centered_Flow() == JST; bool limiter_flow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); bool limiter_adjflow = (cont_adjoint && (config->GetKind_SlopeLimit_AdjFlow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter())); - bool fixed_cl = config->GetFixed_CL_Mode(); - bool engine = ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)); - bool actuator_disk = ((config->GetnMarker_ActDiskInlet() != 0) || (config->GetnMarker_ActDiskOutlet() != 0)); - bool nearfield = (config->GetnMarker_NearFieldBound() != 0); bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; - unsigned short kind_row_dissipation = config->GetKind_RoeLowDiss(); - bool roe_low_dissipation = (kind_row_dissipation != NO_ROELOWDISS) && - (config->GetKind_Upwind_Flow() == ROE || - config->GetKind_Upwind_Flow() == SLAU || - config->GetKind_Upwind_Flow() == SLAU2); bool wall_functions = config->GetWall_Functions(); - /*--- Update the angle of attack at the far-field for fixed CL calculations (only direct problem). ---*/ + /*--- Common preprocessing steps (implemented by CEulerSolver) ---*/ - if (fixed_cl && !disc_adjoint && !cont_adjoint) { - SU2_OMP_MASTER - SetFarfield_AoA(geometry, solver_container, config, iMesh, Output); - SU2_OMP_BARRIER - } - - /*--- Set the primitive variables ---*/ - - SU2_OMP_MASTER - ErrorCounter = 0; - SU2_OMP_BARRIER - - SU2_OMP_ATOMIC - ErrorCounter += SetPrimitive_Variables(solver_container, config, Output); - - if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) { - SU2_OMP_BARRIER - SU2_OMP_MASTER - { - unsigned long tmp = ErrorCounter; - SU2_MPI::Allreduce(&tmp, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); - config->SetNonphysical_Points(ErrorCounter); - } - SU2_OMP_BARRIER - } - - /*--- Compute the engine properties ---*/ - - if (engine) { - SU2_OMP_MASTER - GetPower_Properties(geometry, config, iMesh, Output); - SU2_OMP_BARRIER - } - - /*--- Compute the actuator disk properties and distortion levels ---*/ - - if (actuator_disk) { - SU2_OMP_MASTER - { - Set_MPI_ActDisk(solver_container, geometry, config); - SetActDisk_BCThrust(geometry, solver_container, config, iMesh, Output); - } - SU2_OMP_BARRIER - } - - /*--- Compute NearField MPI ---*/ - - if (nearfield) { - SU2_OMP_MASTER - Set_MPI_Nearfield(geometry, config); - SU2_OMP_BARRIER - } - - /*--- Artificial dissipation ---*/ - - if (center && !Output) { - SetMax_Eigenvalue(geometry, config); - if ((center_jst) && (iMesh == MESH_0)) { - SetCentered_Dissipation_Sensor(geometry, config); - SetUndivided_Laplacian(geometry, config); - } - } - - /*--- Roe Low Dissipation Sensor ---*/ - - if (roe_low_dissipation){ - SetRoe_Dissipation(geometry, config); - if (kind_row_dissipation == FD_DUCROS || kind_row_dissipation == NTS_DUCROS){ - SetUpwind_Ducros_Sensor(geometry, config); - } - } + CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); /*--- Compute gradient for MUSCL reconstruction. ---*/ if (config->GetReconstructionGradientRequired() && (iMesh == MESH_0)) { - if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) - SetPrimitive_Gradient_GG(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); + switch (config->GetKind_Gradient_Method_Recon()) { + case GREEN_GAUSS: + SetPrimitive_Gradient_GG(geometry, config, true); break; + case LEAST_SQUARES: + case WEIGHTED_LEAST_SQUARES: + SetPrimitive_Gradient_LS(geometry, config, true); break; + default: + break; + } } /*--- Compute gradient of the primitive variables ---*/ @@ -280,12 +200,12 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { SetPrimitive_Gradient_GG(geometry, config); } - if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + else if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { SetPrimitive_Gradient_LS(geometry, config); } - /*--- Compute the limiter in case we need it in the turbulence model - or to limit the viscous terms (check this logic with JST and 2nd order turbulence model) ---*/ + /*--- Compute the limiter in case we need it in the turbulence model or to limit the + * viscous terms (check this logic with JST and 2nd order turbulence model) ---*/ if ((iMesh == MESH_0) && (limiter_flow || limiter_turb || limiter_adjflow) && !Output && !van_albada) { SetPrimitive_Limiter(geometry, config); @@ -342,15 +262,6 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C SU2_OMP_BARRIER } - /*--- Initialize the Jacobian matrix and residual, not needed for the reducer strategy - * as we set blocks (including diagonal ones) and completely overwrite. ---*/ - - if(!ReducerStrategy && !Output) { - LinSysRes.SetValZero(); - if (implicit && !config->GetDiscrete_Adjoint()) Jacobian.SetValZero(); - else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes - } - } unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { From 775980e4d0531d13fc1191f293c88270e1f163d2 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Thu, 5 Mar 2020 14:10:05 +0000 Subject: [PATCH 23/28] fix virtual bug in SetPrimitiveVariables --- SU2_CFD/include/solvers/CEulerSolver.hpp | 6 ++---- SU2_CFD/include/solvers/CNSSolver.hpp | 4 +--- SU2_CFD/src/solvers/CEulerSolver.cpp | 19 +++++++------------ SU2_CFD/src/solvers/CNSSolver.cpp | 3 +-- 4 files changed, 11 insertions(+), 21 deletions(-) diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 30a7c455fa73..8b335d7b3384 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -381,7 +381,6 @@ class CEulerSolver : public CSolver { */ inline virtual void SetRoe_Dissipation(CGeometry *geometry, CConfig *config) { } -private: /*! * \brief Compute the velocity^2, SoundSpeed, Pressure, Enthalpy, Viscosity. * \param[in] solver_container - Container vector with all the solutions. @@ -389,9 +388,8 @@ class CEulerSolver : public CSolver { * \param[in] Output - boolean to determine whether to print output. * \return - The number of non-physical points. */ - unsigned long SetPrimitive_Variables(CSolver **solver_container, - CConfig *config, - bool Output); + virtual unsigned long SetPrimitive_Variables(CSolver **solver_container, + CConfig *config, bool Output); protected: diff --git a/SU2_CFD/include/solvers/CNSSolver.hpp b/SU2_CFD/include/solvers/CNSSolver.hpp index 0566e637c363..aa0bcfe59b6d 100644 --- a/SU2_CFD/include/solvers/CNSSolver.hpp +++ b/SU2_CFD/include/solvers/CNSSolver.hpp @@ -68,7 +68,6 @@ class CNSSolver final : public CEulerSolver { */ void SetRoe_Dissipation(CGeometry *geometry, CConfig *config) override; -private: /*! * \brief Compute the velocity^2, SoundSpeed, Pressure, Enthalpy, Viscosity. * \param[in] solver_container - Container vector with all the solutions. @@ -77,8 +76,7 @@ class CNSSolver final : public CEulerSolver { * \return - The number of non-physical points. */ unsigned long SetPrimitive_Variables(CSolver **solver_container, - CConfig *config, - bool Output); + CConfig *config, bool Output) override; protected: diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index c45bdc5ba616..b35a3faa5f60 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2614,7 +2614,7 @@ void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_con if(!ReducerStrategy && !Output) { LinSysRes.SetValZero(); - if (implicit && !config->GetDiscrete_Adjoint()) Jacobian.SetValZero(); + if (implicit && !disc_adjoint) Jacobian.SetValZero(); else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes } @@ -2646,8 +2646,7 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container case LEAST_SQUARES: case WEIGHTED_LEAST_SQUARES: SetPrimitive_Gradient_LS(geometry, config, true); break; - default: - break; + default: break; } /*--- Limiter computation ---*/ @@ -3713,18 +3712,21 @@ void CEulerSolver::SetUndivided_Laplacian_And_Centered_Dissipation_Sensor(CGeome } if (isPeriodic) { - /*--- Correct the sensor values across any periodic boundaries. ---*/ + /*--- Correct the Laplacian and sensor values across any periodic boundaries. ---*/ SU2_OMP_MASTER { for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); } } SU2_OMP_BARRIER - /*--- Set pressure switch for each point ---*/ + /*--- Set final pressure switch for each point ---*/ SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) @@ -3733,13 +3735,6 @@ void CEulerSolver::SetUndivided_Laplacian_And_Centered_Dissipation_Sensor(CGeome SU2_OMP_MASTER { - /*--- Correct the Laplacian values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_LAPLACIAN); - } - /*--- MPI parallelization ---*/ InitiateComms(geometry, config, UNDIVIDED_LAPLACIAN); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index d2031ca454ec..3b22cc0c292a 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -190,8 +190,7 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C case LEAST_SQUARES: case WEIGHTED_LEAST_SQUARES: SetPrimitive_Gradient_LS(geometry, config, true); break; - default: - break; + default: break; } } From 45e3cbc8a15d8541c50928c5837b07f848b9d6c6 Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Thu, 5 Mar 2020 15:32:57 +0000 Subject: [PATCH 24/28] fix some indentation in CDriver --- SU2_CFD/src/drivers/CDriver.cpp | 116 ++++++++++++++++---------------- 1 file changed, 57 insertions(+), 59 deletions(-) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 85205cffa674..738fbac3ae94 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2586,15 +2586,15 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe /*--- On the target side ---*/ nMarkerTarget = config[targetZone]->GetnMarker_All(); - for (iMarkerTarget = 0; iMarkerTarget < nMarkerTarget; iMarkerTarget++) { + for (iMarkerTarget = 0; iMarkerTarget < nMarkerTarget; iMarkerTarget++) { - /*--- If the tag GetMarker_All_ZoneInterface(iMarker) equals the index we are looping at ---*/ - if ( config[targetZone]->GetMarker_All_ZoneInterface(iMarkerTarget) == iMarkerInt ) { - /*--- We have identified the identifier for the interface marker ---*/ - markTarget = iMarkerTarget; + /*--- If the tag GetMarker_All_ZoneInterface(iMarker) equals the index we are looping at ---*/ + if ( config[targetZone]->GetMarker_All_ZoneInterface(iMarkerTarget) == iMarkerInt ) { + /*--- We have identified the identifier for the interface marker ---*/ + markTarget = iMarkerTarget; - break; - } + break; + } } #ifdef HAVE_MPI @@ -2602,36 +2602,34 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe Donor_check = -1; Target_check = -1; - /*--- We gather a vector in MASTER_NODE that determines if the boundary is not on the processor because - * of the partition or because the zone does not include it ---*/ + /*--- We gather a vector in MASTER_NODE that determines if the boundary is not on the processor because + * of the partition or because the zone does not include it ---*/ - SU2_MPI::Gather(&markDonor , 1, MPI_INT, Buffer_Recv_mark, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); + SU2_MPI::Gather(&markDonor , 1, MPI_INT, Buffer_Recv_mark, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); if (rank == MASTER_NODE) { for (iRank = 0; iRank < nProcessor; iRank++) { if( Buffer_Recv_mark[iRank] != -1 ) { - Donor_check = Buffer_Recv_mark[iRank]; - - break; - } + Donor_check = Buffer_Recv_mark[iRank]; + break; } } + } - SU2_MPI::Bcast(&Donor_check , 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); + SU2_MPI::Bcast(&Donor_check , 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Gather(&markTarget, 1, MPI_INT, Buffer_Recv_mark, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); + SU2_MPI::Gather(&markTarget, 1, MPI_INT, Buffer_Recv_mark, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); if (rank == MASTER_NODE){ for (iRank = 0; iRank < nProcessor; iRank++){ if( Buffer_Recv_mark[iRank] != -1 ){ - Target_check = Buffer_Recv_mark[iRank]; - - break; - } + Target_check = Buffer_Recv_mark[iRank]; + break; } } + } - SU2_MPI::Bcast(&Target_check, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); + SU2_MPI::Bcast(&Target_check, 1, MPI_INT, MASTER_NODE, MPI_COMM_WORLD); #else Donor_check = markDonor; @@ -2704,39 +2702,39 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe /*--- Match Zones ---*/ if (rank == MASTER_NODE) cout << "Setting coupling "; - bool conservative_interp = config[donorZone]->GetConservativeInterpolation(); + bool conservative_interp = config[donorZone]->GetConservativeInterpolation(); - /*--- Conditions for conservative interpolation are not met, we cannot fallback on the consistent approach - because CFlowTractionInterface relies on the information in config to be correct. ---*/ - if ( conservative_interp && targetZone == 0 && structural_target ) - SU2_MPI::Error("Conservative interpolation assumes the structural model mesh is evaluated second, " - "somehow this has not happened.",CURRENT_FUNCTION); + /*--- Conditions for conservative interpolation are not met, we cannot fallback on the consistent approach + because CFlowTractionInterface relies on the information in config to be correct. ---*/ + if ( conservative_interp && targetZone == 0 && structural_target ) + SU2_MPI::Error("Conservative interpolation assumes the structural model mesh is evaluated second, " + "somehow this has not happened.",CURRENT_FUNCTION); - switch (config[donorZone]->GetKindInterpolation()) { + switch (config[donorZone]->GetKindInterpolation()) { - case NEAREST_NEIGHBOR: - if ( conservative_interp && targetZone > 0 && structural_target ) { - interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); - if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " - "from opposite mesh." << endl; - } - else { + case NEAREST_NEIGHBOR: + if ( conservative_interp && targetZone > 0 && structural_target ) { + interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); + if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " + "from opposite mesh." << endl; + } + else { interpolation[donorZone][targetZone] = new CNearestNeighbor(geometry, config, donorZone, targetZone); if (rank == MASTER_NODE) cout << "using a nearest-neighbor approach." << endl; - } - break; + } + break; - case ISOPARAMETRIC: - if ( conservative_interp && targetZone > 0 && structural_target ) { - interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); - if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " - "from opposite mesh." << endl; - } - else { + case ISOPARAMETRIC: + if ( conservative_interp && targetZone > 0 && structural_target ) { + interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); + if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " + "from opposite mesh." << endl; + } + else { interpolation[donorZone][targetZone] = new CIsoparametric(geometry, config, donorZone, targetZone); if (rank == MASTER_NODE) cout << "using an isoparametric approach." << endl; - } - break; + } + break; case WEIGHTED_AVERAGE: interpolation[donorZone][targetZone] = new CSlidingMesh(geometry, config, donorZone, targetZone); @@ -2744,19 +2742,19 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe break; - case RADIAL_BASIS_FUNCTION: - if ( conservative_interp && targetZone > 0 && structural_target ) { - interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); - if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " - "from opposite mesh." << endl; - } - else { - interpolation[donorZone][targetZone] = new CRadialBasisFunction(geometry, config, - donorZone, targetZone); - if (rank == MASTER_NODE) cout << "using a radial basis function approach." << endl; - } - break; - } + case RADIAL_BASIS_FUNCTION: + if ( conservative_interp && targetZone > 0 && structural_target ) { + interpolation[donorZone][targetZone] = new CMirror(geometry, config, donorZone, targetZone); + if (rank == MASTER_NODE) cout << "using a mirror approach: matching coefficients " + "from opposite mesh." << endl; + } + else { + interpolation[donorZone][targetZone] = new CRadialBasisFunction(geometry, config, + donorZone, targetZone); + if (rank == MASTER_NODE) cout << "using a radial basis function approach." << endl; + } + break; + } /*--- Initialize the appropriate transfer strategy ---*/ if (rank == MASTER_NODE) cout << "Transferring "; From 7fb78f83f7d5c7847b6877dd34a3227d1c8005eb Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 11 Mar 2020 10:40:33 +0000 Subject: [PATCH 25/28] add hybrid options to config_template --- config_template.cfg | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/config_template.cfg b/config_template.cfg index ec58a8a5c79d..352d37bed08b 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1333,6 +1333,26 @@ UQ_URLX= 0.1 % Perturbation magnitude (float [0,1], default= 1.0) UQ_DELTA_B= 1.0 % +% --------------------- HYBRID PARALLEL (MPI+OpenMP) OPTIONS ---------------------% +% +% An advanced performance parameter for FVM solvers, a large-ish value should be best +% when relatively few threads per MPI rank are in use (~4). However, maximum parallelism +% is obtained with EDGE_COLORING_GROUP_SIZE=1, consider using this value only if SU2 +% warns about low coloring efficiency during preprocessing (performance is usually worse). +% Setting the option to 0 disables coloring and a different strategy is used instead, +% that strategy is automatically used when the coloring efficiency is less than 0.875. +% The optimum value/strategy is case-dependent. +EDGE_COLORING_GROUP_SIZE= 512 +% +% Independent "threads per MPI rank" setting for LU-SGS and ILU preconditioners. +% For problems where time is spend mostly in the solution of linear systems (e.g. elasticity, +% very high CFL central schemes), AND, if the memory bandwidth of the machine is saturated +% (4 or more cores per memory channel) better performance (via a reduction in linear iterations) +% may be possible by using a smaller value than that defined by the system or in the call to +% SU2_CFD (via the -t/--threads option). +% The default (0) means "same number of threads as for all else". +LINEAR_SOLVER_PREC_THREADS= 0 +% % ------------------------- SCREEN/HISTORY VOLUME OUTPUT --------------------------% % % Screen output fields (use 'SU2_CFD -d ' to view list of available fields) From 2aaf89f3c1af251ba6a36955ced462afce0cf75a Mon Sep 17 00:00:00 2001 From: pcarruscag Date: Wed, 11 Mar 2020 14:01:15 +0000 Subject: [PATCH 26/28] convert moving mesh part of turb solver's dual time residual to point loop --- SU2_CFD/src/solvers/CTurbSolver.cpp | 106 +++++++++++----------------- 1 file changed, 41 insertions(+), 65 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index d814efdb7348..536cabc6dba7 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -607,14 +607,18 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con const bool second_order = (config->GetTime_Marching() == DT_STEPPING_2ND); const bool incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE); + /*--- Flow solution, needed to get density. ---*/ + + CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + /*--- Store the physical time step ---*/ const su2double TimeStep = config->GetDelta_UnstTimeND(); /*--- Local variables ---*/ - unsigned short iVar, iMarker, iDim; - unsigned long iPoint, jPoint, iVertex; + unsigned short iVar, iMarker, iDim, iNeigh; + unsigned long iPoint, jPoint, iVertex, iEdge; const su2double *U_time_nM1 = nullptr, *U_time_n = nullptr, *U_time_nP1 = nullptr; su2double Volume_nM1, Volume_nP1; @@ -656,14 +660,14 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con density could also be temperature dependent, but as it is not a part of the solution vector it's neither stored for previous time steps nor updated with the solution at the end of each iteration. */ - Density_nM1 = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); - Density_nP1 = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + Density_nM1 = flowNodes->GetDensity(iPoint); + Density_n = flowNodes->GetDensity(iPoint); + Density_nP1 = flowNodes->GetDensity(iPoint); } else{ - Density_nM1 = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n1(iPoint)[0]; - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n(iPoint,0); - Density_nP1 = solver_container[FLOW_SOL]->GetNodes()->GetSolution(iPoint,0); + Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; + Density_n = flowNodes->GetSolution_time_n(iPoint,0); + Density_nP1 = flowNodes->GetSolution(iPoint,0); } for (iVar = 0; iVar < nVar; iVar++) { @@ -703,71 +707,43 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con we will loop over the edges and boundaries to compute the GCL component of the dual time source term that depends on grid velocities. ---*/ - /// TODO: Make this a point loop. - SU2_OMP_MASTER - for(auto iEdge = 0ul; iEdge < geometry->GetnEdge(); ++iEdge) { - - /*--- Get indices for nodes i & j plus the face normal ---*/ - - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Normal = geometry->edge[iEdge]->GetNormal(); - - /*--- Grid velocities stored at nodes i & j ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { GridVel_i = geometry->node[iPoint]->GetGridVel(); - GridVel_j = geometry->node[jPoint]->GetGridVel(); - - /*--- Compute the GCL term by averaging the grid velocities at the - edge mid-point and dotting with the face normal. ---*/ - - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - - /*--- Compute the GCL component of the source term for node i ---*/ - - U_time_n = nodes->GetSolution_time_n(iPoint); - - /*--- Multiply by density at node i for the SST model ---*/ + U_time_n = nodes->GetSolution_time_n(iPoint); + Density_n = 1.0; if (sst_model) { if (incompressible) - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity(iPoint); // Temporary fix else - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n(iPoint,0); - - for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(iPoint,iVar) += Density_n*U_time_n[iVar]*Residual_GCL; - } - else { - for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; + Density_n = flowNodes->GetSolution_time_n(iPoint,0); } - /*--- Compute the GCL component of the source term for node j ---*/ + for (iNeigh = 0; iNeigh < geometry->node[iPoint]->GetnNeighbor(); iNeigh++) { - U_time_n = nodes->GetSolution_time_n(jPoint); + iEdge = geometry->node[iPoint]->GetEdge(iNeigh); + Normal = geometry->edge[iEdge]->GetNormal(); - /*--- Multiply by density at node j for the SST model ---*/ + jPoint = geometry->node[iPoint]->GetPoint(iNeigh); + GridVel_j = geometry->node[jPoint]->GetGridVel(); - if (sst_model) { - if (incompressible) - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint); // Temporary fix - else - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n(jPoint)[0]; + /*--- Determine whether to consider the normal outward or inward. ---*/ + su2double dir = (geometry->edge[iEdge]->GetNode(0) == iPoint)? 0.5 : -0.5; + + Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL += dir*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + + Residual_GCL *= Density_n; for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(jPoint,iVar) -= Density_n*U_time_n[iVar]*Residual_GCL; - } - else { - for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(jPoint,iVar) -= U_time_n[iVar]*Residual_GCL; + LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; } } - SU2_OMP_BARRIER - /*--- Loop over the boundary edges ---*/ + /*--- Loop over the boundary edges ---*/ for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && @@ -800,9 +776,9 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con if (sst_model) { if (incompressible) - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity(iPoint); // Temporary fix else - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n(iPoint,0); + Density_n = flowNodes->GetSolution_time_n(iPoint,0); for (iVar = 0; iVar < nVar; iVar++) LinSysRes(iPoint,iVar) += Density_n*U_time_n[iVar]*Residual_GCL; @@ -850,14 +826,14 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con density could also be temperature dependent, but as it is not a part of the solution vector it's neither stored for previous time steps nor updated with the solution at the end of each iteration. */ - Density_nM1 = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); - Density_nP1 = solver_container[FLOW_SOL]->GetNodes()->GetDensity(iPoint); + Density_nM1 = flowNodes->GetDensity(iPoint); + Density_n = flowNodes->GetDensity(iPoint); + Density_nP1 = flowNodes->GetDensity(iPoint); } else { - Density_nM1 = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n1(iPoint)[0]; - Density_n = solver_container[FLOW_SOL]->GetNodes()->GetSolution_time_n(iPoint,0); - Density_nP1 = solver_container[FLOW_SOL]->GetNodes()->GetSolution(iPoint,0); + Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; + Density_n = flowNodes->GetSolution_time_n(iPoint,0); + Density_nP1 = flowNodes->GetSolution(iPoint,0); } for (iVar = 0; iVar < nVar; iVar++) { From fea625ede52741840d621334b07405a405a8d8c2 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 24 Mar 2020 08:24:20 +0000 Subject: [PATCH 27/28] revise logic for flow variable reconstruction with MUSCL turbulence --- SU2_CFD/src/solvers/CTurbSolver.cpp | 104 +++++++++++++++------------- 1 file changed, 54 insertions(+), 50 deletions(-) diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 536cabc6dba7..661a7fa2b856 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -96,6 +96,13 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe const bool muscl = config->GetMUSCL_Turb(); const bool limiter = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER); + /*--- Only reconstruct flow variables if MUSCL is on for flow (requires upwind) and turbulence. ---*/ + const bool musclFlow = config->GetMUSCL_Flow() && muscl && + (config->GetKind_ConvNumScheme_Flow() == SPACE_UPWIND); + /*--- Only consider flow limiters for cell-based limiters, edge-based would need to be recomputed. ---*/ + const bool limiterFlow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && + (config->GetKind_SlopeLimit_Flow() != VAN_ALBADA_EDGE); + CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); /*--- Pick one numerics object per thread. ---*/ @@ -125,14 +132,14 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe /*--- Primitive variables w/o reconstruction ---*/ - auto V_i = flowNodes->GetPrimitive(iPoint); - auto V_j = flowNodes->GetPrimitive(jPoint); + const auto V_i = flowNodes->GetPrimitive(iPoint); + const auto V_j = flowNodes->GetPrimitive(jPoint); numerics->SetPrimitive(V_i, V_j); /*--- Turbulent variables w/o reconstruction ---*/ - auto Turb_i = nodes->GetSolution(iPoint); - auto Turb_j = nodes->GetSolution(jPoint); + const auto Turb_i = nodes->GetSolution(iPoint); + const auto Turb_j = nodes->GetSolution(jPoint); numerics->SetTurbVar(Turb_i, Turb_j); /*--- Grid Movement ---*/ @@ -141,75 +148,72 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[jPoint]->GetGridVel()); - if (muscl) { - su2double *Limiter_i = nullptr, *Limiter_j = nullptr; - - auto Coord_i = geometry->node[iPoint]->GetCoord(); - auto Coord_j = geometry->node[jPoint]->GetCoord(); + if (muscl || musclFlow) { + const su2double *Limiter_i = nullptr, *Limiter_j = nullptr; - /*--- Reconstruct flow variables. ---*/ + const auto Coord_i = geometry->node[iPoint]->GetCoord(); + const auto Coord_j = geometry->node[jPoint]->GetCoord(); su2double Vector_ij[MAXNDIM] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) { Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]); } - /*--- Mean flow primitive variables using gradient reconstruction and limiters ---*/ - - auto Gradient_i = flowNodes->GetGradient_Reconstruction(iPoint); - auto Gradient_j = flowNodes->GetGradient_Reconstruction(jPoint); + if (musclFlow) { + /*--- Reconstruct mean flow primitive variables. ---*/ - if (limiter) { - Limiter_i = flowNodes->GetLimiter_Primitive(iPoint); - Limiter_j = flowNodes->GetLimiter_Primitive(jPoint); - } + auto Gradient_i = flowNodes->GetGradient_Reconstruction(iPoint); + auto Gradient_j = flowNodes->GetGradient_Reconstruction(jPoint); - for (iVar = 0; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) { - su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim]; - Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim]; - } - if (limiter) { - flowPrimVar_i[iVar] = V_i[iVar] + Limiter_i[iVar]*Project_Grad_i; - flowPrimVar_j[iVar] = V_j[iVar] + Limiter_j[iVar]*Project_Grad_j; + if (limiterFlow) { + Limiter_i = flowNodes->GetLimiter_Primitive(iPoint); + Limiter_j = flowNodes->GetLimiter_Primitive(jPoint); } - else { + + for (iVar = 0; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) { + su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim]; + Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim]; + } + if (limiterFlow) { + Project_Grad_i *= Limiter_i[iVar]; + Project_Grad_j *= Limiter_j[iVar]; + } flowPrimVar_i[iVar] = V_i[iVar] + Project_Grad_i; flowPrimVar_j[iVar] = V_j[iVar] + Project_Grad_j; } - } - - numerics->SetPrimitive(flowPrimVar_i, flowPrimVar_j); - /*--- Reconstruct turbulence variables. ---*/ + numerics->SetPrimitive(flowPrimVar_i, flowPrimVar_j); + } - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); + if (muscl) { + /*--- Reconstruct turbulence variables. ---*/ - if (limiter) { - Limiter_i = nodes->GetLimiter(iPoint); - Limiter_j = nodes->GetLimiter(jPoint); - } + auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint); + auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint); - for (iVar = 0; iVar < nVar; iVar++) { - su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim]; - Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim]; - } if (limiter) { - solution_i[iVar] = Turb_i[iVar] + Limiter_i[iVar]*Project_Grad_i; - solution_j[iVar] = Turb_j[iVar] + Limiter_j[iVar]*Project_Grad_j; + Limiter_i = nodes->GetLimiter(iPoint); + Limiter_j = nodes->GetLimiter(jPoint); } - else { + + for (iVar = 0; iVar < nVar; iVar++) { + su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim]; + Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim]; + } + if (limiter) { + Project_Grad_i *= Limiter_i[iVar]; + Project_Grad_j *= Limiter_j[iVar]; + } solution_i[iVar] = Turb_i[iVar] + Project_Grad_i; solution_j[iVar] = Turb_j[iVar] + Project_Grad_j; } - } - - numerics->SetTurbVar(solution_i, solution_j); + numerics->SetTurbVar(solution_i, solution_j); + } } /*--- Update convective residual value ---*/ From b6361f1973be838de3b866953042683b72508fab Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 25 Mar 2020 19:39:55 +0000 Subject: [PATCH 28/28] update some testcases --- TestCases/parallel_regression.py | 4 ++-- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 4 ++-- TestCases/serial_regression_AD.py | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 82b58fbf9076..4475490dd1c6 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -792,7 +792,7 @@ def main(): unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" unst_inc_turb_naca0015_sa.test_iter = 1 - unst_inc_turb_naca0015_sa.test_vals = [-2.990702, -6.861616, 1.475736, 0.419771] #last 4 columns + unst_inc_turb_naca0015_sa.test_vals = [-2.990703, -6.865923, 1.475736, 0.419770] #last 4 columns unst_inc_turb_naca0015_sa.su2_exec = "parallel_computation.py -f" unst_inc_turb_naca0015_sa.timeout = 1600 unst_inc_turb_naca0015_sa.tol = 0.00001 @@ -976,7 +976,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.590195, -1.955395] #last 4 columns + bars_SST_2D.test_vals = [13.000000, -0.598555, -1.879082] #last 4 columns bars_SST_2D.su2_exec = "SU2_CFD" bars_SST_2D.timeout = 1600 bars_SST_2D.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 69ceb1be21aa..a1c17b61e2c2 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -84,7 +84,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.230578, 0.678810, 0.181780, -0.000018] #last 4 columns + discadj_rans_naca0012_sa.test_vals = [-2.230573, 0.696562, 0.181780, -0.000018] #last 4 columns discadj_rans_naca0012_sa.su2_exec = "parallel_computation.py -f" discadj_rans_naca0012_sa.timeout = 1600 discadj_rans_naca0012_sa.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 2b2ee1dd9ae8..f57084dc5757 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -926,7 +926,7 @@ def main(): unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" unst_inc_turb_naca0015_sa.test_iter = 1 - unst_inc_turb_naca0015_sa.test_vals = [ -2.994996, -6.865786, 1.434864, 0.416627] #last 4 columns + unst_inc_turb_naca0015_sa.test_vals = [-2.994996, -6.869781, 1.434864, 0.416626] #last 4 columns unst_inc_turb_naca0015_sa.su2_exec = "SU2_CFD" unst_inc_turb_naca0015_sa.timeout = 1600 unst_inc_turb_naca0015_sa.tol = 0.00001 @@ -1134,7 +1134,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.590195, -1.955395] #last 3 columns + bars_SST_2D.test_vals = [13.000000, -0.598555, -1.879082] #last 3 columns bars_SST_2D.su2_exec = "SU2_CFD" bars_SST_2D.timeout = 1600 bars_SST_2D.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index e90b5d0e973a..70efb28b1f1d 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -84,7 +84,7 @@ def main(): discadj_rans_naca0012_sa.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" discadj_rans_naca0012_sa.test_iter = 10 - discadj_rans_naca0012_sa.test_vals = [-2.230561, 0.678834, 0.180740, -0.000018] #last 4 columns + discadj_rans_naca0012_sa.test_vals = [-2.230556, 0.696586, 0.180740, -0.000018] #last 4 columns discadj_rans_naca0012_sa.su2_exec = "SU2_CFD_AD" discadj_rans_naca0012_sa.timeout = 1600 discadj_rans_naca0012_sa.tol = 0.00001