diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 13c266e4c38e..2f142f4fc4ee 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4332,6 +4332,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i } } + if(Kind_TimeIntScheme_Turb != EULER_IMPLICIT && + Kind_TimeIntScheme_Turb != EULER_EXPLICIT){ + SU2_MPI::Error("Only TIME_DISCRE_TURB = EULER_IMPLICIT, EULER_EXPLICIT have been implemented.", CURRENT_FUNCTION); + } + if (nIntCoeffs == 0) { nIntCoeffs = 2; Int_Coeffs = new su2double[2]; Int_Coeffs[0] = 0.25; Int_Coeffs[1] = 0.5; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 183aea428d3c..d28cbfc564db 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -775,13 +775,8 @@ class CFVMFlowSolverBase : public CSolver { const su2double RK_FuncCoeff[] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0}; const su2double RK_TimeCoeff[] = {0.5, 0.5, 1.0, 1.0}; - /*--- Set shared residual variables to 0 and declare - * local ones for current thread to work on. ---*/ - - SetResToZero(); - + /*--- Local residual variables for current thread ---*/ su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; unsigned long idxMax[MAXNVAR] = {0}; /*--- Update the solution and residuals ---*/ @@ -832,23 +827,13 @@ class CFVMFlowSolverBase : public CSolver { } /*--- Update residual information for current thread. ---*/ - resRMS[iVar] += Res*Res; - if (fabs(Res) > resMax[iVar]) { - resMax[iVar] = fabs(Res); - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + ResidualReductions_PerThread(iPoint, iVar, Res, resRMS, resMax, idxMax); } } END_SU2_OMP_FOR /*--- Reduce residual information over all threads in this rank. ---*/ - SU2_OMP_CRITICAL - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Residual_RMS[iVar] += resRMS[iVar]; - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - END_SU2_OMP_CRITICAL - SU2_OMP_BARRIER + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); + } /*--- MPI solution ---*/ @@ -857,9 +842,6 @@ class CFVMFlowSolverBase : public CSolver { CompleteComms(geometry, config, SOLUTION); if (!adjoint) { - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); - /*--- For verification cases, compute the global error metrics. ---*/ ComputeVerificationError(geometry, config); } @@ -893,12 +875,8 @@ class CFVMFlowSolverBase : public CSolver { const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - /*--- Set shared residual variables to 0 and declare local ones for current thread to work on. ---*/ - - SetResToZero(); - + /*--- Local residual variables for current thread ---*/ su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; unsigned long idxMax[MAXNVAR] = {0}; /*--- Add pseudotime term to Jacobian. ---*/ @@ -948,26 +926,15 @@ class CFVMFlowSolverBase : public CSolver { LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]); LinSysSol[total_index] = 0.0; - su2double Res = fabs(LinSysRes[total_index]); - resRMS[iVar] += Res*Res; - if (Res > resMax[iVar]) { - resMax[iVar] = Res; - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint, iVar, LinSysRes[total_index], resRMS, resMax, idxMax); } } END_SU2_OMP_FOR - SU2_OMP_CRITICAL - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Residual_RMS[iVar] += resRMS[iVar]; - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - END_SU2_OMP_CRITICAL - SU2_OMP_BARRIER - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); + } /*! diff --git a/SU2_CFD/include/solvers/CScalarSolver.hpp b/SU2_CFD/include/solvers/CScalarSolver.hpp index 4f7cc2492da9..a95f3b3de8b2 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.hpp +++ b/SU2_CFD/include/solvers/CScalarSolver.hpp @@ -259,6 +259,14 @@ class CScalarSolver : public CSolver { */ void CompleteImplicitIteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) final; + /*! + * \brief Update the solution using the explicit Euler scheme. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + void ExplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, CConfig* config) final; + /*! * \brief Update the solution using an implicit solver. * \param[in] geometry - Geometrical definition of the problem. @@ -287,7 +295,8 @@ class CScalarSolver : public CSolver { * \param[in] val_iter - Current external iteration number. * \param[in] val_update_geo - Flag for updating coords and grid velocity. */ - virtual void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, bool val_update_geo) override = 0; + void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, + bool val_update_geo) override = 0; /*! * \brief Scalar solvers support OpenMP+MPI. diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 8cc4e044ee91..5d5986426df1 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -79,7 +79,8 @@ CScalarSolver::~CScalarSolver() { template void CScalarSolver::Upwind_Residual(CGeometry* geometry, CSolver** solver_container, - CNumerics** numerics_container, CConfig* config, unsigned short iMesh) { + CNumerics** numerics_container, CConfig* config, + unsigned short iMesh) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool muscl = config->GetMUSCL_Turb(); const bool limiter = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER); @@ -256,7 +257,7 @@ void CScalarSolver::SumEdgeFluxes(CGeometry* geometry) { template void CScalarSolver::BC_Periodic(CGeometry* geometry, CSolver** solver_container, CNumerics* numerics, - CConfig* config) { + CConfig* config) { /*--- Complete residuals for periodic boundary conditions. We loop over the periodic BCs in matching pairs so that, in the event that there are adjacent periodic markers, the repeated points will have their residuals @@ -271,7 +272,7 @@ void CScalarSolver::BC_Periodic(CGeometry* geometry, CSolver** sol template void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, CSolver** solver_container, - CConfig* config) { + CConfig* config) { const auto flowNodes = solver_container[FLOW_SOL]->GetNodes(); /*--- Set shared residual variables to 0 and declare @@ -280,7 +281,6 @@ void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, SetResToZero(); su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; unsigned long idxMax[MAXNVAR] = {0}; /*--- Build implicit system ---*/ @@ -308,31 +308,19 @@ void CScalarSolver::PrepareImplicitIteration(CGeometry* geometry, LinSysRes[total_index] = -LinSysRes[total_index]; LinSysSol[total_index] = 0.0; - su2double Res = fabs(LinSysRes[total_index]); - resRMS[iVar] += Res * Res; - if (Res > resMax[iVar]) { - resMax[iVar] = Res; - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint, iVar, LinSysRes[total_index], resRMS, resMax, idxMax); } } END_SU2_OMP_FOR - SU2_OMP_CRITICAL - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - Residual_RMS[iVar] += resRMS[iVar]; - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - END_SU2_OMP_CRITICAL - SU2_OMP_BARRIER - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); } template void CScalarSolver::CompleteImplicitIteration(CGeometry* geometry, CSolver** solver_container, - CConfig* config) { + CConfig* config) { const bool compressible = (config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE); const auto flowNodes = solver_container[FLOW_SOL]->GetNodes(); @@ -380,7 +368,7 @@ void CScalarSolver::CompleteImplicitIteration(CGeometry* geometry, template void CScalarSolver::ImplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, - CConfig* config) { + CConfig* config) { PrepareImplicitIteration(geometry, solver_container, config); /*--- Solve or smooth the linear system. ---*/ @@ -404,10 +392,41 @@ void CScalarSolver::ImplicitEuler_Iteration(CGeometry* geometry, C CompleteImplicitIteration(geometry, solver_container, config); } +template +void CScalarSolver::ExplicitEuler_Iteration(CGeometry* geometry, CSolver** solver_container, + CConfig* config) { + const auto flowNodes = solver_container[FLOW_SOL]->GetNodes(); + + /*--- Local residual variables for current thread ---*/ + su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; + unsigned long idxMax[MAXNVAR] = {0}; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + const su2double dt = nodes->GetLocalCFL(iPoint) / flowNodes->GetLocalCFL(iPoint) * flowNodes->GetDelta_Time(iPoint); + nodes->SetDelta_Time(iPoint, dt); + const su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint, iVar, LinSysRes(iPoint, iVar), resRMS, resMax, idxMax); + /*--- Explicit Euler step: ---*/ + LinSysSol(iPoint, iVar) = -dt / Vol * LinSysRes(iPoint, iVar); + } + } + END_SU2_OMP_FOR + + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS, resMax, idxMax); + + /*--- Use LinSysSol for solution update. ---*/ + CompleteImplicitIteration(geometry, solver_container, config); +} + template void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSolver** solver_container, CConfig* config, - unsigned short iRKStep, unsigned short iMesh, - unsigned short RunTime_EqSystem) { + unsigned short iRKStep, unsigned short iMesh, + unsigned short RunTime_EqSystem) { const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST); const bool second_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index b239474da7ca..729c15d6027f 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -435,19 +435,6 @@ class CSolver { unsigned short iMesh, unsigned short RunTime_EqSystem) { } - /*! - * \brief Set the RMS and MAX residual to zero. - */ - inline void SetResToZero() { - SU2_OMP_MASTER { - for (auto& r : Residual_RMS) r = 0; - for (auto& r : Residual_Max) r = 0; - for (auto& p : Point_Max) p = 0; - } - END_SU2_OMP_MASTER - SU2_OMP_BARRIER - } - /*! * \brief Get the maximal residual, this is useful for the convergence history. * \param[in] val_var - Index of the variable. @@ -455,25 +442,6 @@ class CSolver { */ inline su2double GetRes_RMS(unsigned short val_var) const { return Residual_RMS[val_var]; } - /*! - * \brief Adds the maximal residual, this is useful for the convergence history (overload). - * \param[in] val_var - Index of the variable. - * \param[in] val_residual - Value of the residual to store in the position val_var. - * \param[in] val_point - Value of the point index for the max residual. - * \param[in] val_coord - Location (x, y, z) of the max residual point. - */ - inline void AddRes_Max(unsigned short val_var, - su2double val_residual, - unsigned long val_point, - const su2double* val_coord) { - if (val_residual > Residual_Max[val_var]) { - Residual_Max[val_var] = val_residual; - Point_Max[val_var] = val_point; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Point_Max_Coord[val_var][iDim] = val_coord[iDim]; - } - } - /*! * \brief Get the maximal residual, this is useful for the convergence history. * \param[in] val_var - Index of the variable. @@ -488,25 +456,6 @@ class CSolver { */ inline su2double GetRes_BGS(unsigned short val_var) const { return Residual_BGS[val_var]; } - /*! - * \brief Adds the maximal residual for BGS subiterations. - * \param[in] val_var - Index of the variable. - * \param[in] val_residual - Value of the residual to store in the position val_var. - * \param[in] val_point - Value of the point index for the max residual. - * \param[in] val_coord - Location (x, y, z) of the max residual point. - */ - inline void AddRes_Max_BGS(unsigned short val_var, - su2double val_residual, - unsigned long val_point, - const su2double* val_coord) { - if (val_residual > Residual_Max_BGS[val_var]) { - Residual_Max_BGS[val_var] = val_residual; - Point_Max_BGS[val_var] = val_point; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Point_Max_Coord_BGS[val_var][iDim] = val_coord[iDim]; - } - } - /*! * \brief Get the maximal residual for BGS subiterations. * \param[in] val_var - Index of the variable. @@ -4326,4 +4275,94 @@ class CSolver { void SetVerificationSolution(unsigned short nDim, unsigned short nVar, CConfig *config); + + /*! + * \brief "Add" residual at (iPoint,iVar) to residual variables local to the thread. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \param[in] res - Residual at (iPoint,iVar), e.g. LinSysRes(iPoint,iVar) + * \param[in,out] resRMS - increases by pow(Residual, 2) + * \param[in,out] resMax - increases to max(resMax, Residual) + * \param[in,out] idxMax - changes when resMax increases + */ + static inline void ResidualReductions_PerThread(unsigned long iPoint, unsigned short iVar, su2double res, su2double* resRMS, su2double* resMax, + unsigned long* idxMax) { + res = fabs(res); + resRMS[iVar] += res * res; + if (res > resMax[iVar]) { + resMax[iVar] = res; + idxMax[iVar] = iPoint; + } + } + + /*! + * \brief "Add" local residual variables of all threads to compute global residual variables. + */ + inline void ResidualReductions_FromAllThreads(const CGeometry* geometry, const CConfig* config, const su2double* resRMS, const su2double* resMax, + const unsigned long* idxMax){ + SetResToZero(); + + SU2_OMP_CRITICAL + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Residual_RMS[iVar] += resRMS[iVar]; + AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), geometry->nodes->GetCoord(idxMax[iVar])); + } + END_SU2_OMP_CRITICAL + SU2_OMP_BARRIER + + /*--- Compute the root mean square residual ---*/ + SetResidual_RMS(geometry, config); + } + + /*! + * \brief Set the RMS and MAX residual to zero. + */ + inline void SetResToZero() { + SU2_OMP_MASTER { + for (auto& r : Residual_RMS) r = 0; + for (auto& r : Residual_Max) r = 0; + for (auto& p : Point_Max) p = 0; + } + END_SU2_OMP_MASTER + SU2_OMP_BARRIER + } + + /*! + * \brief Adds the maximal residual, this is useful for the convergence history. + * \param[in] val_var - Index of the variable. + * \param[in] val_residual - Value of the residual to store in the position val_var. + * \param[in] val_point - Value of the point index for the max residual. + * \param[in] val_coord - Location (x, y, z) of the max residual point. + */ + inline void AddRes_Max(unsigned short val_var, + su2double val_residual, + unsigned long val_point, + const su2double* val_coord) { + if (val_residual > Residual_Max[val_var]) { + Residual_Max[val_var] = val_residual; + Point_Max[val_var] = val_point; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Point_Max_Coord[val_var][iDim] = val_coord[iDim]; + } + } + + /*! + * \brief Adds the maximal residual for BGS subiterations. + * \param[in] val_var - Index of the variable. + * \param[in] val_residual - Value of the residual to store in the position val_var. + * \param[in] val_point - Value of the point index for the max residual. + * \param[in] val_coord - Location (x, y, z) of the max residual point. + */ + inline void AddRes_Max_BGS(unsigned short val_var, + su2double val_residual, + unsigned long val_point, + const su2double* val_coord) { + if (val_residual > Residual_Max_BGS[val_var]) { + Residual_Max_BGS[val_var] = val_residual; + Point_Max_BGS[val_var] = val_point; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Point_Max_Coord_BGS[val_var][iDim] = val_coord[iDim]; + } + } + }; diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 363b2be39b42..f0915208673b 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -305,9 +305,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi const su2double relax = (config->GetInnerIter()==0) ? 1.0 : config->GetRelaxation_Factor_Adjoint(); /*--- Thread-local residual variables. ---*/ - su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; unsigned long idxMax[MAXNVAR] = {0}; /*--- Set the old solution and compute residuals. ---*/ @@ -329,13 +327,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi nodes->AddSolution(iPoint, iVar, relax*residual); if (iPoint < nPointDomain) { - /*--- Update residual information for current thread. ---*/ - resRMS[iVar] += residual*residual; - if (fabs(residual) > resMax[iVar]) { - resMax[iVar] = fabs(residual); - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint,iVar,residual,resRMS,resMax,idxMax); } } } @@ -344,18 +337,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/ if (CrossTerm) return; - SetResToZero(); - - /*--- Reduce residual information over all threads in this rank. ---*/ - SU2_OMP_CRITICAL - for (auto iVar = 0u; iVar < nVar; iVar++) { - Residual_RMS[iVar] += resRMS[iVar]; - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - END_SU2_OMP_CRITICAL - SU2_OMP_BARRIER - - SetResidual_RMS(geometry, config); + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS,resMax,idxMax); SU2_OMP_MASTER { SetIterLinSolver(direct_solver->System.GetIterations()); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 6aa942874a56..28f2857b4974 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1847,41 +1847,23 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CConfig *config, CNumerics const auto ResidualAux = computeLinearResidual(Jacobian, LinSysSol, LinSysRes); - /*--- Set maximum residual to zero. ---*/ - - SetResToZero(); - SU2_OMP_PARALLEL { /*--- Compute the residual. ---*/ - su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; unsigned long idxMax[MAXNVAR] = {0}; SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { for (auto iVar = 0ul; iVar < nVar; iVar++) { - su2double Res = fabs(ResidualAux(iPoint, iVar)); - resRMS[iVar] += Res*Res; - if (Res > resMax[iVar]) { - resMax[iVar] = Res; - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + /*--- "Add" residual at (iPoint,iVar) to local residual variables. ---*/ + ResidualReductions_PerThread(iPoint, iVar, ResidualAux(iPoint, iVar), resRMS, resMax, idxMax); } } END_SU2_OMP_FOR - SU2_OMP_CRITICAL - for (auto iVar = 0ul; iVar < nVar; iVar++) { - Residual_RMS[iVar] += resRMS[iVar]; - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - END_SU2_OMP_CRITICAL - SU2_OMP_BARRIER - /*--- Compute the root mean square residual. ---*/ - SetResidual_RMS(geometry, config); + /*--- "Add" residuals from all threads to global residual variables. ---*/ + ResidualReductions_FromAllThreads(geometry, config, resRMS,resMax,idxMax); } END_SU2_OMP_PARALLEL diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index 3f5aed4d0227..2d1bedf06e46 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -171,7 +171,6 @@ CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iM const bool euler_implicit = (config->GetKind_TimeIntScheme_Heat() == EULER_IMPLICIT); SetImplicitPeriodic(euler_implicit); - /*--- MPI solution ---*/ InitiateComms(geometry, config, SOLUTION); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index f4287ea0aa97..2b0aa4d4d68f 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -164,11 +164,6 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar) = nu_tilde_Inf; - /*--- The turbulence models are always solved implicitly, so set the - implicit flag in case we have periodic BCs. ---*/ - - SetImplicitPeriodic(true); - /*--- Store the initial CFL number for all grid points. ---*/ const su2double CFL = config->GetCFL(MGLevel)*config->GetCFLRedCoeff_Turb(); diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 0d526c7cdfb0..0cf75b3d0290 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -167,11 +167,6 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh } } - /*--- The turbulence models are always solved implicitly, so set the - implicit flag in case we have periodic BCs. ---*/ - - SetImplicitPeriodic(true); - /*--- Store the initial CFL number for all grid points. ---*/ const su2double CFL = config->GetCFL(MGLevel)*config->GetCFLRedCoeff_Turb(); diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp index 1f4f689c9d53..230585a8abbb 100644 --- a/SU2_CFD/src/solvers/CTurbSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSolver.cpp @@ -33,7 +33,11 @@ /*--- Explicit instantiation of the parent class of CTurbSolver. ---*/ template class CScalarSolver; -CTurbSolver::CTurbSolver(CGeometry* geometry, CConfig *config, bool conservative) : CScalarSolver(geometry, config, conservative) { } +CTurbSolver::CTurbSolver(CGeometry* geometry, CConfig *config, bool conservative) + : CScalarSolver(geometry, config, conservative) { + /*--- Store if an implicit scheme is used, for use during periodic boundary conditions. ---*/ + SetImplicitPeriodic(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT); +} CTurbSolver::~CTurbSolver() { for (auto& mat : SlidingState) { diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index dd491243c7ca..5bef7dee6a63 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -210,6 +210,14 @@ def main(): turb_naca0012_sst_fixedvalues.test_vals = [-5.192502, -9.575898, -1.568269, 1.022571, 0.040527, -2.384329] test_list.append(turb_naca0012_sst_fixedvalues) + # NACA0012 (SST, explicit Euler for flow and turbulence equations) + turb_naca0012_sst_expliciteuler = TestCase('turb_naca0012_sst_expliciteuler') + turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" + turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" + turb_naca0012_sst_expliciteuler.test_iter = 10 + turb_naca0012_sst_expliciteuler.test_vals = [-3.532228, -3.157766, 3.364025, 1.124824, 0.501717, -float("inf")] + test_list.append(turb_naca0012_sst_expliciteuler) + # PROPELLER propeller = TestCase('propeller') propeller.cfg_dir = "rans/propeller" diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ccdb0f6a8308..14902a094d67 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -369,6 +369,17 @@ def main(): turb_naca0012_sst_fixedvalues.tol = 0.00001 test_list.append(turb_naca0012_sst_fixedvalues) + # NACA0012 (SST, explicit Euler for flow and turbulence equations) + turb_naca0012_sst_expliciteuler = TestCase('turb_naca0012_sst_expliciteuler') + turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" + turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" + turb_naca0012_sst_expliciteuler.test_iter = 10 + turb_naca0012_sst_expliciteuler.test_vals = [-3.532228, -3.157766, 3.364025, 1.124824, 0.501717, -float("inf")] + turb_naca0012_sst_expliciteuler.su2_exec = "parallel_computation.py -f" + turb_naca0012_sst_expliciteuler.timeout = 3200 + turb_naca0012_sst_expliciteuler.tol = 0.00001 + test_list.append(turb_naca0012_sst_expliciteuler) + # PROPELLER propeller = TestCase('propeller') propeller.cfg_dir = "rans/propeller" diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_expliciteuler.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_expliciteuler.cfg new file mode 100644 index 000000000000..8cef50bf7b24 --- /dev/null +++ b/TestCases/rans/naca0012/turb_NACA0012_sst_expliciteuler.cfg @@ -0,0 +1,79 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Testing explicit solution of flow and turbulence equations % +% Author: Max Aehle % +% Institution: TU Kaiserslautern % +% Date: Nov 17th, 2021 % +% File Version 7.2.1 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +SOLVER= RANS +KIND_TURB_MODEL= SST +MATH_PROBLEM= DIRECT +RESTART_SOL= NO +READ_BINARY_RESTART= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +MACH_NUMBER= 0.15 +AOA= 10.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 6.0E6 +REYNOLDS_LENGTH= 1.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +REF_LENGTH= 1.0 +REF_AREA= 1.0 +REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +MARKER_HEATFLUX= ( airfoil, 0.0 ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING= ( airfoil ) +MARKER_MONITORING= ( airfoil ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +NUM_METHOD_GRAD_RECON= LEAST_SQUARES +CFL_NUMBER= 0.1 +MAX_DELTA_TIME= 1E10 +ITER= 99999 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +CONV_NUM_METHOD_FLOW= ROE +USE_VECTORIZATION= YES +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_EXPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= NO +SLOPE_LIMITER_TURB= NONE +TIME_DISCRE_TURB= EULER_EXPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +CONV_FIELD= RMS_DENSITY +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +MESH_FILENAME= n0012_113-33.su2 +MESH_FORMAT= SU2 +% +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG, LINSOL_RESIDUAL) +% +TABULAR_FORMAT= CSV +CONV_FILENAME= history +% +OUTPUT_FILES= (RESTART_ASCII, PARAVIEW, SURFACE_PARAVIEW) +OUTPUT_WRT_FREQ= 10000 +RESTART_FILENAME= restart_flow.dat +SOLUTION_FILENAME= solution_flow_sst_expliciteuler.dat +VOLUME_FILENAME= flow +SURFACE_FILENAME= surface_flow diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_fixedvalues.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_fixedvalues.cfg index 3dc75033aacc..7e1de50b4362 100644 --- a/TestCases/rans/naca0012/turb_NACA0012_sst_fixedvalues.cfg +++ b/TestCases/rans/naca0012/turb_NACA0012_sst_fixedvalues.cfg @@ -10,199 +10,81 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% -% -% Physical governing equations (EULER, NAVIER_STOKES, -% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY, -% POISSON_EQUATION) SOLVER= RANS -% -% Specify turbulent model (NONE, SA, SA_NEG, SST) KIND_TURB_MODEL= SST -% -% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT) MATH_PROBLEM= DIRECT -% -% Restart solution (NO, YES) RESTART_SOL= NO -% -% Read binary restart files (YES, NO) READ_BINARY_RESTART= NO -% % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% -% -% Mach number (non-dimensional, based on the free-stream values) MACH_NUMBER= 0.15 -% -% Angle of attack (degrees, only for compressible flows) AOA= 10.0 -% -% Free-stream temperature (288.15 K by default) FREESTREAM_TEMPERATURE= 300.0 -% -% Reynolds number (non-dimensional, based on the free-stream values) REYNOLDS_NUMBER= 6.0E6 -% -% Reynolds length (1 m by default) REYNOLDS_LENGTH= 1.0 -% -% Free-stream turbulence intensity FREESTREAM_TURBULENCEINTENSITY= 0.0008165 -% -% Free-stream ratio between turbulent and laminar viscosity FREESTREAM_TURB2LAMVISCRATIO= 1.2 -% -% Fix turbulence quantities to far-field values at some upstream half-space TURB_FIXED_VALUES= YES -% -% Shift of the fixed values half-space in unit far-field velocity vectors TURB_FIXED_VALUES_DOMAIN= -1.0 % ---------------------- REFERENCE VALUE DEFINITION ---------------------------% -% -% Reference origin for moment computation REF_ORIGIN_MOMENT_X = 0.25 REF_ORIGIN_MOMENT_Y = 0.00 REF_ORIGIN_MOMENT_Z = 0.00 -% -% Reference length for pitching, rolling, and yawing non-dimensional moment REF_LENGTH= 1.0 -% -% Reference area for force coefficients (0 implies automatic calculation) REF_AREA= 1.0 -% -% Flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE, -% FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE) REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE % -------------------- BOUNDARY CONDITION DEFINITION --------------------------% -% -% Navier-Stokes wall boundary marker(s) (NONE = no marker) MARKER_HEATFLUX= ( airfoil, 0.0 ) -% -% Farfield boundary marker(s) (NONE = no marker) MARKER_FAR= ( farfield ) -% -% Marker(s) of the surface to be plotted or designed MARKER_PLOTTING= ( airfoil ) -% -% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated MARKER_MONITORING= ( airfoil ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% -% -% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES NUM_METHOD_GRAD_RECON= LEAST_SQUARES -% -% Courant-Friedrichs-Lewy condition of the finest grid CFL_NUMBER= 1000.0 -% -% Max Delta time MAX_DELTA_TIME= 1E10 -% -% Number of total iterations ITER= 99999 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% -% -% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) LINEAR_SOLVER= FGMRES -% -% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI) LINEAR_SOLVER_PREC= ILU -% -% Minimum error of the linear solver for implicit formulations LINEAR_SOLVER_ERROR= 1E-10 -% -% Max number of iterations of the linear solver for the implicit formulation LINEAR_SOLVER_ITER= 20 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% -% -% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, -% TURKEL_PREC, MSW) CONV_NUM_METHOD_FLOW= ROE USE_VECTORIZATION= YES -% -% Spatial numerical order integration (1ST_ORDER, 2ND_ORDER, 2ND_ORDER_LIMITER) MUSCL_FLOW= YES -% -% Slope limiter (VENKATAKRISHNAN, MINMOD) SLOPE_LIMITER_FLOW= NONE -% -% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) TIME_DISCRE_FLOW= EULER_IMPLICIT % -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% -% -% Convective numerical method (SCALAR_UPWIND) CONV_NUM_METHOD_TURB= SCALAR_UPWIND -% -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. -% Required for 2nd order upwind schemes (NO, YES) MUSCL_TURB= NO -% -% Slope limiter (VENKATAKRISHNAN, MINMOD) SLOPE_LIMITER_TURB= NONE -% -% Time discretization (EULER_IMPLICIT) TIME_DISCRE_TURB= EULER_IMPLICIT -% -% Reduction factor of the CFL coefficient in the turbulence problem CFL_REDUCTION_TURB= 1.0 % --------------------------- CONVERGENCE PARAMETERS --------------------------% -% -% Convergence field CONV_FIELD= RMS_DENSITY -% -% Min value of the residual (log10 of the residual) CONV_RESIDUAL_MINVAL= -12 -% -% Start convergence criteria at iteration number CONV_STARTITER= 10 -% -% Number of elements to apply the criteria CONV_CAUCHY_ELEMS= 100 -% -% Epsilon to control the series convergence CONV_CAUCHY_EPS= 1E-6 -% % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% -% -% Mesh input file MESH_FILENAME= n0012_113-33.su2 -% -% Mesh input file format (SU2, CGNS, NETCDF_ASCII) MESH_FORMAT= SU2 -% -% Mesh output file MESH_OUT_FILENAME= mesh_out.su2 -% -% Restart flow input file SOLUTION_FILENAME= solution_flow_sst_fixedvalues.dat -% -% Output file format (PARAVIEW, TECPLOT, STL) TABULAR_FORMAT= CSV -% -% Output file convergence history (w/o extension) CONV_FILENAME= history -% -% Output file restart flow RESTART_FILENAME= restart_flow.dat -% -% Output file flow (w/o extension) variables VOLUME_FILENAME= flow -% -% Output file surface flow coefficient (w/o extension) SURFACE_FILENAME= surface_flow -% -% Writing solution file frequency OUTPUT_WRT_FREQ= 10000 -% -% -% Screen output fields SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG, LINSOL_RESIDUAL) OUTPUT_FILES= (RESTART_ASCII, PARAVIEW, SURFACE_PARAVIEW) diff --git a/config_template.cfg b/config_template.cfg index 7bf02ba3cbb4..358e9773fd5f 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1209,7 +1209,7 @@ KIND_MATRIX_COLORING= GREEDY_COLORING % Convective numerical method (SCALAR_UPWIND) CONV_NUM_METHOD_TURB= SCALAR_UPWIND % -% Time discretization (EULER_IMPLICIT) +% Time discretization (EULER_IMPLICIT, EULER_EXPLICIT) TIME_DISCRE_TURB= EULER_IMPLICIT % % Reduction factor of the CFL coefficient in the turbulence problem