From fd86820508966ad201c1786a43b73489aa55ab2f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 24 Apr 2026 18:24:02 +0100 Subject: [PATCH 1/2] address comments 1 --- AUTHORS.md | 1 + SU2_CFD/include/numerics/CNumerics.hpp | 15 ++++---- .../include/numerics_simd/CNumericsSIMD.cpp | 4 ++- SU2_CFD/src/output/CFlowOutput.cpp | 36 +++++++++---------- SU2_CFD/src/solvers/CTurbSASolver.cpp | 18 +++++----- 5 files changed, 38 insertions(+), 36 deletions(-) diff --git a/AUTHORS.md b/AUTHORS.md index f6580fadb892..554d925a8e97 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -52,6 +52,7 @@ Amit Sachdeva Ana Lourenco Andrew Burkett Andrew Wendorff +Angelo Passariello Aniket C. Aranake Antonio Rubino Arne Bachmann diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e910c6c5a487..db2c5d8f35af 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -192,9 +192,8 @@ class CNumerics { su2double lesMode_i = 0.0, /*!< \brief LES sensor at point i for hybrid RANS-LES methods. */ lesMode_j = 0.0; /*!< \brief LES sensor at point j for hybrid RANS-LES methods. */ - unsigned short - sbsInBox_i = 0, /*!< \brief Sensor to assess if point i lies inside the box where the Stochastic Backscatter Model is active. */ - sbsInBox_j = 0; /*!< \brief Sensor to assess if point j lies inside the box where the Stochastic Backscatter Model is active. */ + int8_t + sbsInBox_i = 0; /*!< \brief Sensor to assess if point i lies inside the box where the Stochastic Backscatter Model is active. */ SST_ParsedOptions sstParsedOptions; /*!< \brief additional options for the SST turbulence model */ unsigned short Eig_Val_Comp; /*!< \brief Component towards which perturbation is perfromed */ su2double uq_delta_b; /*!< \brief Magnitude of perturbation */ @@ -651,7 +650,7 @@ class CNumerics { } } } - + /*! * \brief Compute a random contribution to the Reynolds stress tensor (Stochastic Backscatter Model). * \details See: Kok, Johan C. "A stochastic backscatter model for grey-area mitigation in detached @@ -664,12 +663,12 @@ class CNumerics { * \param[out] stochReynStress - Stochastic tensor (to be added to the Reynolds stress tensor). */ template - inline void ComputeStochReynStress(su2double density, su2double tke, const Vec& rndVec, + inline void ComputeStochReynStress(su2double density, su2double tke, const Vec& rndVec, su2double Cmag, Mat& stochReynStress) { /* --- Calculate stochastic tensor --- */ - su2double stochLim = 3.0; + su2double stochLim = 3.0; stochReynStress[0][0] = 0.0; stochReynStress[1][1] = 0.0; @@ -918,11 +917,9 @@ class CNumerics { /*! * \brief Set the sensor to locate the box where the Stochastic Backscatter Model is active. * \param[in] val_sbsInBox_i - 1 if point i lies inside the box where the model is active. - * \param[in] val_sbsInBox_j - 1 if point j lies inside the box where the model is active. */ - inline void SetSbsInBoxSensor(unsigned short val_sbsInBox_i, unsigned short val_sbsInBox_j) { + inline void SetSbsInBoxSensor(int8_t val_sbsInBox_i) { sbsInBox_i = val_sbsInBox_i; - sbsInBox_j = val_sbsInBox_j; } /*! diff --git a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp index 68afef986934..e05435d2d269 100644 --- a/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp +++ b/SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp @@ -74,7 +74,6 @@ CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CV obj = new CLaxScheme(config, iMesh, turbVars); break; case CENTERED::JST: - case CENTERED::LD2: // Just to silence compiler warnings (LD2 implemented in the incompressible solver only). obj = new CJSTScheme(config, iMesh, turbVars); break; case CENTERED::JST_KE: @@ -83,6 +82,9 @@ CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CV case CENTERED::JST_MAT: obj = new CJSTmatScheme(config, iMesh, turbVars); break; + case CENTERED::LD2: + /*--- LD2 implemented only in the incompressible solver. ---*/ + break; } return obj; } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 7b6ae71b4873..495d956e0ab5 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -4136,19 +4136,19 @@ void CFlowOutput::SetTimeAveragedFields(const CConfig *config) { } if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) { - AddVolumeOutput("MODELED_REYNOLDS_STRESS_11", "ModeledReynoldsStress_11", "TIME_AVERAGE", "Modeled Reynolds stress xx-component"); - AddVolumeOutput("MODELED_REYNOLDS_STRESS_22", "ModeledReynoldsStress_22", "TIME_AVERAGE", "Modeled Reynolds stress yy-component"); - AddVolumeOutput("MODELED_REYNOLDS_STRESS_12", "ModeledReynoldsStress_12", "TIME_AVERAGE", "Modeled Reynolds stress xy-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_XX", "ModeledReynoldsStress_XX", "TIME_AVERAGE", "Modeled Reynolds stress xx-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_YY", "ModeledReynoldsStress_YY", "TIME_AVERAGE", "Modeled Reynolds stress yy-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_XY", "ModeledReynoldsStress_XY", "TIME_AVERAGE", "Modeled Reynolds stress xy-component"); if (nDim == 3){ - AddVolumeOutput("MODELED_REYNOLDS_STRESS_33", "ModeledReynoldsStress_33", "TIME_AVERAGE", "Modeled Reynolds stress zz-component"); - AddVolumeOutput("MODELED_REYNOLDS_STRESS_13", "ModeledReynoldsStress_13", "TIME_AVERAGE", "Modeled Reynolds stress xz-component"); - AddVolumeOutput("MODELED_REYNOLDS_STRESS_23", "ModeledReynoldsStress_23", "TIME_AVERAGE", "Modeled Reynolds stress yz-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_ZZ", "ModeledReynoldsStress_ZZ", "TIME_AVERAGE", "Modeled Reynolds stress zz-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_XZ", "ModeledReynoldsStress_XZ", "TIME_AVERAGE", "Modeled Reynolds stress xz-component"); + AddVolumeOutput("MODELED_REYNOLDS_STRESS_YZ", "ModeledReynoldsStress_YZ", "TIME_AVERAGE", "Modeled Reynolds stress yz-component"); } if (config->GetSBSParam().StochasticBackscatter) { - AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_12", "StochasticReynoldsStress_12", "TIME_AVERAGE", "Stochastic Reynolds stress xy-component"); - AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_13", "StochasticReynoldsStress_13", "TIME_AVERAGE", "Stochastic Reynolds stress xz-component"); - AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_23", "StochasticReynoldsStress_23", "TIME_AVERAGE", "Stochastic Reynolds stress yz-component"); + AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_XY", "StochasticReynoldsStress_XY", "TIME_AVERAGE", "Stochastic Reynolds stress xy-component"); + AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_XZ", "StochasticReynoldsStress_XZ", "TIME_AVERAGE", "Stochastic Reynolds stress xz-component"); + AddVolumeOutput("STOCHASTIC_REYNOLDS_STRESS_YZ", "StochasticReynoldsStress_YZ", "TIME_AVERAGE", "Stochastic Reynolds stress yz-component"); } } } @@ -4201,16 +4201,16 @@ void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *No const su2double tau_xx = nu_t * (2*vel_grad(0,0) - (2.0/3.0)*vel_div); const su2double tau_yy = nu_t * (2*vel_grad(1,1) - (2.0/3.0)*vel_div); const su2double tau_xy = nu_t * (vel_grad(0,1) + vel_grad(1,0)); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_11", iPoint, -tau_xx); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_22", iPoint, -tau_yy); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_12", iPoint, -tau_xy); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_XX", iPoint, -tau_xx); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_YY", iPoint, -tau_yy); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_XY", iPoint, -tau_xy); if (nDim == 3){ const su2double tau_zz = nu_t * (2*vel_grad(2,2) - (2.0/3.0)*vel_div); const su2double tau_xz = nu_t * (vel_grad(0,2) + vel_grad(2,0)); const su2double tau_yz = nu_t * (vel_grad(1,2) + vel_grad(2,1)); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_33", iPoint, -tau_zz); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_13", iPoint, -tau_xz); - SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_23", iPoint, -tau_yz); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_ZZ", iPoint, -tau_zz); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_XZ", iPoint, -tau_xz); + SetAvgVolumeOutputValue("MODELED_REYNOLDS_STRESS_YZ", iPoint, -tau_yz); } if (config->GetSBSParam().StochasticBackscatter) { @@ -4226,9 +4226,9 @@ void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *No const su2double R_xy = - mag * tke_estim * csi_z; const su2double R_xz = + mag * tke_estim * csi_y; const su2double R_yz = - mag * tke_estim * csi_x; - SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_12", iPoint, -R_xy); - SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_13", iPoint, -R_xz); - SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_23", iPoint, -R_yz); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_XY", iPoint, -R_xy); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_XY", iPoint, -R_xz); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_YZ", iPoint, -R_yz); } } } diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index dbad3f42e69d..d67b1e117c44 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -427,7 +427,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai for (unsigned short iDim = 0; iDim < nDim; iDim++) numerics->SetStochSource(nodes->GetLangevinSourceTerms(iPoint, iDim), iDim); numerics->SetLES_Mode(nodes->GetLES_Mode(iPoint), 0.0); - numerics->SetSbsInBoxSensor(nodes->GetSbsInBox(iPoint), 0); + numerics->SetSbsInBoxSensor(nodes->GetSbsInBox(iPoint)); } } @@ -453,7 +453,7 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai if (transition_BC || config->GetKind_Trans_Model() != TURB_TRANS_MODEL::NONE) { nodes->SetIntermittency(iPoint,numerics->GetIntermittencyEff()); } - + /*--- Subtract residual and the Jacobian ---*/ LinSysRes.SubtractBlock(iPoint, residual); @@ -1636,7 +1636,8 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC } void CTurbSASolver::SetBackscatterInBox(CConfig *config, CGeometry *geometry) { - + SU2_ZONE_SCOPED + auto sbsBoxBounds = config->GetSBSParam().StochBackscatterBoxBounds; SU2_OMP_FOR_STAT(omp_chunk_size) @@ -1653,6 +1654,7 @@ void CTurbSASolver::SetBackscatterInBox(CConfig *config, CGeometry *geometry) { } void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) { + SU2_ZONE_SCOPED const su2double threshold = config->GetSBSParam().stochFdThreshold; const su2double dummySource = 1e3; @@ -1710,9 +1712,9 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet /*--- Start SOR algorithm for the Laplacian smoothing. ---*/ for (unsigned short iDim = 0; iDim < nDim; iDim++) { - + for (unsigned short iter = 0; iter < maxIter; iter++) { - + /*--- MPI communication. ---*/ InitiateComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG); @@ -1778,7 +1780,7 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet << "---------------------------------" << endl; } if (iter%10 == 0) { - cout << " " + cout << " " << std::setw(5) << iter << " " << std::setw(12) << std::fixed << std::setprecision(6) << log10(globalResNorm) @@ -1789,14 +1791,14 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet if (log10(globalResNorm) < tol || iter == maxIter-1) { if (rank == MASTER_NODE) { - cout << " " + cout << " " << std::setw(5) << iter << " " << std::setw(12) << std::fixed << ::setprecision(6) << log10(globalResNorm) << endl; cout << "---------------------------------" << endl; } - + /*--- Scale source terms for variance preservation. ---*/ su2double var_check_old = 0.0; From 4c4bc203c411746254ebf7620981efb341ba0d10 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 25 Apr 2026 10:40:15 +0100 Subject: [PATCH 2/2] fix OMP issues --- SU2_CFD/include/variables/CVariable.hpp | 1 - SU2_CFD/src/numerics/flow/flow_diffusion.cpp | 6 +- SU2_CFD/src/solvers/CTurbSASolver.cpp | 167 ++++++++++--------- 3 files changed, 90 insertions(+), 84 deletions(-) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index f3a6bec8f66d..82cf99645d11 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -34,7 +34,6 @@ #include #include #include -#include #include "../../../Common/include/CConfig.hpp" #include "../../../Common/include/containers/container_decorators.hpp" diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 9d687b4acb06..f2d1d9910ac9 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -477,7 +477,7 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, config); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -655,7 +655,7 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, config); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); @@ -985,7 +985,7 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity,config); + Mean_Laminar_Viscosity, Mean_Eddy_Viscosity, config); if (config->GetSAParsedOptions().qcr2000) AddQCR(nDim, &Mean_GradPrimVar[1], tau); if (Mean_TauWall > 0) AddTauWall(UnitNormal, Mean_TauWall); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index d67b1e117c44..952a245072a7 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -1679,8 +1679,8 @@ void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) } END_SU2_OMP_FOR - SU2_OMP_FOR_DYN(omp_chunk_size) for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) for (unsigned long iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { unsigned long iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) { @@ -1690,8 +1690,8 @@ void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) } } } + END_SU2_OMP_FOR } - END_SU2_OMP_FOR InitiateComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG); CompleteComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG); @@ -1699,6 +1699,11 @@ void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) } void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geometry) { + SU2_ZONE_SCOPED + + static su2double globalResNorm; + static unsigned long global_nPointLES; + static std::array globalChecks; const su2double LES_FilterWidth = config->GetLES_FilterWidth(); const su2double cDelta = config->GetSBSParam().SBS_Cdelta; @@ -1723,16 +1728,15 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet su2double localResNorm = 0.0; unsigned long local_nPointLES = 0; - SU2_OMP_FOR_DYN(omp_chunk_size) + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { su2double source_i_old = nodes->GetLangevinSourceTermsOld(iPoint, iDim); if (source_i_old > 3.0*sourceLim) continue; local_nPointLES += 1; su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); if (LES_FilterWidth > 0.0) maxDelta = LES_FilterWidth; - su2double b = sqrt(cDelta) * maxDelta; - su2double b2 = b * b; - su2double volume_iPoint = geometry->nodes->GetVolume(iPoint); + su2double b2 = cDelta * pow(maxDelta, 2); + su2double volume_iPoint = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); su2double source_i = nodes->GetLangevinSourceTerms(iPoint, iDim); auto coord_i = geometry->nodes->GetCoord(iPoint); @@ -1761,23 +1765,24 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet nodes->SetLangevinSourceTerms(iPoint, iDim, source_i); } - END_SU2_OMP_FOR + END_SU2_OMP_SAFE_GLOBAL_ACCESS /*--- Stop integration if residual drops below tolerance. ---*/ - su2double globalResNorm = 0.0; - unsigned long global_nPointLES = 0; - SU2_MPI::Allreduce(&local_nPointLES, &global_nPointLES, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&localResNorm, &globalResNorm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - globalResNorm = (global_nPointLES==0) ? su2double(0.0) : sqrt(globalResNorm / global_nPointLES); + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + SU2_MPI::Allreduce(&local_nPointLES, &global_nPointLES, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&localResNorm, &globalResNorm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + globalResNorm = (global_nPointLES==0) ? su2double(0.0) : sqrt(globalResNorm / global_nPointLES); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + SU2_OMP_MASTER if (rank == MASTER_NODE) { if (iter == 0) { - cout << endl - << "Residual of Laplacian smoothing along dimension " << iDim+1 << "." << endl - << "---------------------------------" << endl - << " Iter RMS Residual" << endl - << "---------------------------------" << endl; + cout << "\nResidual of Laplacian smoothing along dimension " << iDim+1 + << "\n---------------------------------" + << "\n Iter RMS Residual" + << "\n---------------------------------" << endl; } if (iter%10 == 0) { cout << " " @@ -1787,9 +1792,11 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet << endl; } } + END_SU2_OMP_MASTER if (log10(globalResNorm) < tol || iter == maxIter-1) { + SU2_OMP_MASTER if (rank == MASTER_NODE) { cout << " " << std::setw(5) << iter @@ -1798,55 +1805,49 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet << endl; cout << "---------------------------------" << endl; } + END_SU2_OMP_MASTER /*--- Scale source terms for variance preservation. ---*/ - su2double var_check_old = 0.0; - su2double mean_check_old = 0.0; - su2double var_check_new = 0.0; - su2double mean_check_new = 0.0; - su2double var_check_notSmoothed = 0.0; - su2double mean_check_notSmoothed = 0.0; + su2double mean_check_old = 0.0, var_check_old = 0.0; + su2double mean_check_new = 0.0, var_check_new = 0.0; + su2double mean_check_notSmoothed = 0.0, var_check_notSmoothed = 0.0; + + SU2_OMP_FOR_(schedule(static, omp_chunk_size) SU2_NOWAIT) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { su2double source_notSmoothed = nodes->GetLangevinSourceTermsOld(iPoint, iDim); if (source_notSmoothed > 3.0*sourceLim) continue; su2double source = nodes->GetLangevinSourceTerms(iPoint, iDim); mean_check_old += source; - var_check_old += source * source; + var_check_old += pow(source, 2); mean_check_notSmoothed += source_notSmoothed; - var_check_notSmoothed += source_notSmoothed * source_notSmoothed; + var_check_notSmoothed += pow(source_notSmoothed, 2); su2double integral = 0.0; if (timeIter==restartIter) { su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); if (LES_FilterWidth > 0.0) maxDelta = LES_FilterWidth; - su2double b = sqrt(cDelta) * maxDelta; - su2double b2 = b * b; + su2double b2 = cDelta * pow(maxDelta, 2); auto coord_i = geometry->nodes->GetCoord(iPoint); - su2double maxDelta_vec[3] = {0.0}; - su2double maxDelta_tmp = 0.0; - unsigned short nNeigh = geometry->nodes->GetnPoint(iPoint); - for (unsigned short iNode = 0; iNode < nNeigh; iNode++) { - auto jPoint = geometry->nodes->GetPoint(iPoint, iNode); + su2double maxDelta_vec[MAXNDIM] = {0.0}; + su2double maxDelta2_tmp = 0.0; + for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) { auto coord_j = geometry->nodes->GetCoord(jPoint); - su2double dist_ij = GeometryToolbox::Distance(nDim, coord_j, coord_i); - if (dist_ij > maxDelta_tmp) { - maxDelta_tmp = dist_ij; + su2double dist2_ij = GeometryToolbox::SquaredDistance(nDim, coord_j, coord_i); + if (dist2_ij > maxDelta2_tmp) { + maxDelta2_tmp = dist2_ij; GeometryToolbox::Distance(nDim, coord_j, coord_i, maxDelta_vec); } } su2double max_dist_ij_normal = 0.0; - for (unsigned short iNode = 0; iNode < nNeigh; iNode++) { - auto jPoint = geometry->nodes->GetPoint(iPoint, iNode); + for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) { auto coord_j = geometry->nodes->GetCoord(jPoint); - su2double dist_ij2 = GeometryToolbox::SquaredDistance(nDim, coord_j, coord_i); - su2double dist_ij_vec[3] = {0.0}; + su2double dist_ij_vec[MAXNDIM] = {0.0}; GeometryToolbox::Distance(nDim, coord_j, coord_i, dist_ij_vec); + su2double dist2_ij = GeometryToolbox::SquaredNorm(nDim, dist_ij_vec); su2double dist_ij_parallel = GeometryToolbox::DotProduct(nDim, dist_ij_vec, maxDelta_vec); dist_ij_parallel /= maxDelta; - su2double dist_ij_normal = sqrt(max(dist_ij2 - dist_ij_parallel*dist_ij_parallel, 0.0)); - if (dist_ij_normal > max_dist_ij_normal) { - max_dist_ij_normal = dist_ij_normal; - } + su2double dist_ij_normal = sqrt(max(dist2_ij - pow(dist_ij_parallel, 2), 0.0)); + max_dist_ij_normal = max(max_dist_ij_normal, dist_ij_normal); } su2double dI = maxDelta; su2double dJ = max_dist_ij_normal; @@ -1866,53 +1867,59 @@ void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geomet source *= scaleFactor; if (source < -sourceLim || source > sourceLim) source = 0.0; mean_check_new += source; - var_check_new += source * source; + var_check_new += pow(source, 2); nodes->SetLangevinSourceTerms(iPoint, iDim, source); } - su2double mean_check_new_G = 0.0; - SU2_MPI::Allreduce(&mean_check_new, &mean_check_new_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - mean_check_new_G /= global_nPointLES; + END_SU2_OMP_FOR + + SU2_OMP_SAFE_GLOBAL_ACCESS(globalChecks = {0, 0, 0, 0, 0, 0};) + + atomicAdd(mean_check_old, globalChecks[0]); + atomicAdd(var_check_old, globalChecks[1]); + atomicAdd(mean_check_notSmoothed, globalChecks[2]); + atomicAdd(var_check_notSmoothed, globalChecks[3]); + atomicAdd(mean_check_new, globalChecks[4]); + atomicAdd(var_check_new, globalChecks[5]); + + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + auto tmp = globalChecks; + SU2_MPI::Allreduce(tmp.data(), globalChecks.data(), tmp.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS + + const auto invDenom = 1.0 / max(global_nPointLES, 1ul); + mean_check_new = globalChecks[4] * invDenom; + + SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { su2double source_notSmoothed = nodes->GetLangevinSourceTermsOld(iPoint, iDim); su2double source = nodes->GetLangevinSourceTerms(iPoint, iDim); if (source_notSmoothed > 3.0*sourceLim) continue; - source -= mean_check_new_G; - nodes->SetLangevinSourceTerms(iPoint, iDim, source); + nodes->SetLangevinSourceTerms(iPoint, iDim, source - mean_check_new); } - if (config->GetSBSParam().stochSourceDiagnostics) { - su2double mean_check_old_G = 0.0; - su2double mean_check_notSmoothed_G = 0.0; - su2double var_check_old_G = 0.0; - su2double var_check_new_G = 0.0; - su2double var_check_notSmoothed_G = 0.0; - SU2_MPI::Allreduce(&mean_check_old, &mean_check_old_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&mean_check_notSmoothed, &mean_check_notSmoothed_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&var_check_old, &var_check_old_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&var_check_new, &var_check_new_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - SU2_MPI::Allreduce(&var_check_notSmoothed, &var_check_notSmoothed_G, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - mean_check_old_G /= global_nPointLES; - var_check_old_G /= global_nPointLES; - var_check_old_G -= mean_check_old_G * mean_check_old_G; - var_check_new_G /= global_nPointLES; - var_check_new_G -= mean_check_new_G * mean_check_new_G; - mean_check_notSmoothed_G /= global_nPointLES; - var_check_notSmoothed_G /= global_nPointLES; - var_check_notSmoothed_G -= mean_check_notSmoothed_G * mean_check_notSmoothed_G; - if (rank == MASTER_NODE) { - cout << "Mean of stochastic source term in Langevin equations: " << endl; - cout << " Uncorrelated --> " << mean_check_notSmoothed_G << endl; - cout << " Smoothed before scaling --> " << mean_check_old_G << endl; - cout << " Smoothed after scaling --> " << mean_check_new_G << " (subtracted from stochastic field to guarantee zero mean)" << endl; - cout << "Variance of stochastic source term in Langevin equations: " << endl; - cout << " Uncorrelated --> " << var_check_notSmoothed_G << endl; - cout << " Smoothed before scaling --> " << var_check_old_G << endl; - cout << " Smoothed after scaling --> " << var_check_new_G << endl; - cout << endl; - } + END_SU2_OMP_FOR + + SU2_OMP_MASTER + if (rank == MASTER_NODE && config->GetSBSParam().stochSourceDiagnostics) { + mean_check_old = globalChecks[0] * invDenom; + var_check_old = globalChecks[1] * invDenom - pow(mean_check_old, 2); + mean_check_notSmoothed = globalChecks[2] * invDenom; + var_check_notSmoothed = globalChecks[3] * invDenom - pow(mean_check_notSmoothed, 2); + var_check_new = globalChecks[5] * invDenom - pow(mean_check_new, 2); + + cout << "Mean of stochastic source term in Langevin equations:"; + cout << "\n Uncorrelated --> " << mean_check_notSmoothed; + cout << "\n Smoothed before scaling --> " << mean_check_old; + cout << "\n Smoothed after scaling --> " << mean_check_new << " (subtracted from stochastic field to guarantee zero mean)"; + cout << "\nVariance of stochastic source term in Langevin equations:"; + cout << "\n Uncorrelated --> " << var_check_notSmoothed; + cout << "\n Smoothed before scaling --> " << var_check_old; + cout << "\n Smoothed after scaling --> " << var_check_new << '\n' << endl; } + END_SU2_OMP_MASTER + /*--- Converged or maximum number of iterations reached. ---*/ break; - } } }