diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index a57b43b4ab09..abb11731a7fa 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1097,6 +1097,17 @@ class CConfig { su2double Const_DES; /*!< \brief Detached Eddy Simulation Constant. */ WINDOW_FUNCTION Kind_WindowFct; /*!< \brief Type of window (weight) function for objective functional. */ unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */ + bool StochasticBackscatter; /*!< \brief Option to include Stochastic Backscatter Model. */ + su2double SBS_Cdelta; /*!< \brief Stochastic Backscatter Model lengthscale coefficient. */ + unsigned short SBS_maxIterSmooth; /*!< \brief Maximum number of smoothing iterations for the SBS model. */ + su2double SBS_Ctau; /*!< \brief Stochastic Backscatter Model timescale coefficient. */ + su2double SBS_Cmag; /*!< \brief Stochastic Backscatter Model intensity coefficient. */ + bool stochSourceNu; /*!< \brief Option for including stochastic source term in turbulence model equation (Stochastic Backscatter Model). */ + bool stochSourceDiagnostics; /*!< \brief Option for writing diagnostics related to stochastic source terms in Langevin equations (Stochastic Backscatter Model). */ + su2double stochSourceRelax; /*!< \brief Relaxation factor for stochastic source term generation (Stochastic Backscatter Model). */ + bool enforceLES; /*!< \brief Option to enforce LES mode in hybrid RANS-LES simulations. */ + su2double LES_FilterWidth; /*!< \brief LES filter width for hybrid RANS-LES simulations. */ + bool LD2_Scheme; /*!< \brief Use the LD2 scheme (incompressible flows, combined with JST discretization). */ unsigned short Kind_RoeLowDiss; /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */ unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */ @@ -9563,6 +9574,48 @@ class CConfig { */ unsigned short GetKind_HybridRANSLES(void) const { return Kind_HybridRANSLES; } + /*! + * \brief Get if the Stochastic Backscatter Model must be activated. + * \return TRUE if the Stochastic Backscatter Model is activated. + */ + bool GetStochastic_Backscatter(void) const { return StochasticBackscatter; } + + /*! + * \brief Get if the LES mode must be enforced. + * \return TRUE if LES is enforced. + */ + bool GetEnforceLES(void) const { return enforceLES; } + + /*! + * \brief Get if the stochastic source term must be included in the turbulence model equation. + * \return TRUE if the stochastic source term is included in the turbulence model equation. + */ + bool GetStochSourceNu(void) const { return stochSourceNu; } + + /*! + * \brief Get if the diagnostics of the stochastic source terms in Langevin equations must be computed. + * \return TRUE if the diagnostics of the stochastic source terms in Langevin equations are computed. + */ + bool GetStochSourceDiagnostics(void) const { return stochSourceDiagnostics; } + + /*! + * \brief Get the relaxation factor for stochastic source term generation. + * \return Relaxation factor for stochastic source term generation. + */ + su2double GetStochSourceRelax(void) const { return stochSourceRelax; } + + /*! + * \brief Get the LES Filter Width. + * \return Value of LES Filter Width. + */ + su2double GetLES_FilterWidth(void) const { return LES_FilterWidth; } + + /*! + * \brief Get if the LD2 scheme must be employed (incompressible flows, combined with JST discretization). + * \return TRUE if LD2 scheme is enabled. + */ + bool GetLD2_Scheme(void) const { return LD2_Scheme; } + /*! * \brief Get the Kind of Roe Low Dissipation Scheme for Unsteady flows. * \return Value of Low dissipation approach. @@ -9575,6 +9628,30 @@ class CConfig { */ su2double GetConst_DES(void) const { return Const_DES; } + /*! + * \brief Get the SBS lengthscale coefficient. + * \return Value of SBS lengthscale coefficient. + */ + su2double GetSBS_Cdelta(void) const { return SBS_Cdelta; } + + /*! + * \brief Get the SBS timescale coefficient. + * \return Value of SBS timescale coefficient. + */ + su2double GetSBS_maxIterSmooth(void) const { return SBS_maxIterSmooth; } + + /*! + * \brief Get the SBS timescale coefficient. + * \return Value of SBS timescale coefficient. + */ + su2double GetSBS_Ctau(void) const { return SBS_Ctau; } + + /*! + * \brief Get the SBS intensity coefficient. + * \return Value of SBS intensity coefficient. + */ + su2double GetSBS_Cmag(void) const { return SBS_Cmag; } + /*! * \brief Get the type of tape that will be checked in a tape debug run. */ diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2b7a0e8c3961..abbf4a1cc9fd 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2686,6 +2686,7 @@ enum class MPI_QUANTITIES { MAX_LENGTH , /*!< \brief Maximum length communication. */ GRID_VELOCITY , /*!< \brief Grid velocity communication. */ SOLUTION_EDDY , /*!< \brief Turbulent solution plus eddy viscosity communication. */ + STOCH_SOURCE_LANG , /*!< \brief Stochastic source term for Langevin equations communication. */ SOLUTION_MATRIX , /*!< \brief Matrix solution communication. */ SOLUTION_MATRIXTRANS , /*!< \brief Matrix transposed solution communication. */ NEIGHBORS , /*!< \brief Neighbor point count communication (for JST). */ diff --git a/Common/include/toolboxes/random_toolbox.hpp b/Common/include/toolboxes/random_toolbox.hpp new file mode 100644 index 000000000000..aa50baf1e7f6 --- /dev/null +++ b/Common/include/toolboxes/random_toolbox.hpp @@ -0,0 +1,154 @@ +/*! + * \file random_toolbox.hpp + * \brief Collection of utility functions for random number generation. + * \version 8.3.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include + +namespace RandomToolbox { +/// \addtogroup RandomToolbox +/// @{ + +/*! + * \brief Combine two 64-bit integers into a single hash value. + * \param[in] v1 Current hash value. + * \param[in] v2 Value to mix in. + * \return Combined hash value. + */ +inline unsigned long HashCombine(unsigned long v1, unsigned long v2) { + const unsigned long prime = 1099511628211ULL; + v1 ^= v2; + v1 *= prime; + return v1; +} + +/*! + * \brief Convert a double to a 64-bit integer suitable for hashing. + * \param[in] x Double to integer. + * \return Hash value of the double (not portable). + */ +inline unsigned long ToUInt64(double x) { return std::hash{}(x); } + +/*! + * \brief Build a deterministic seed from physical time. + * \param[in] x First integer value. + * \param[in] y Second integer value. + * \return 64-bit seed value. + */ +inline unsigned long GetSeed(unsigned long x, unsigned long y) { return HashCombine(x, y); } + +/*! + * \brief Generate a standard normally-distributed random number. + * \param[in] gen Pseudo-random number generator. + * \param[in] mean Mean of the normal distribution (default 0). + * \param[in] stddev Standard deviation of the normal distribution (default 1). + * \return Normally-distributed random number. + */ +inline double GetRandomNormal(std::mt19937 gen, double mean = 0.0, double stddev = 1.0) { + std::normal_distribution rnd(mean, stddev); + return rnd(gen); +} + +/*! + * \brief Generate a uniformly-distributed random number. + * \param[in] gen Pseudo-random number generator. + * \param[in] xmin Lower boundary of the interval (default 0). + * \param[in] xmax Upper bounary of the interval (default 1). + * \return Uniformly-distributed random number. + */ +inline double GetRandomUniform(std::mt19937 gen, double xmin = 0.0, double xmax = 1.0) { + std::uniform_real_distribution rnd(xmin, xmax); + return rnd(gen); +} + +/*! + * \brief Compute modified bessel function of first kind (order 0). + * \param[in] x Argument of Bessel funtion. + * \return Value of Bessel function. + */ +inline double GetBesselZero(double x) { + double abx = fabs(x); + + if (abx < 3.75) { + double t = x / 3.75; + double p = + 1.0 + + t * t * + (3.5156229 + + t * t * (3.0899424 + t * t * (1.2067492 + t * t * (0.2659732 + t * t * (0.0360768 + t * t * 0.0045813))))); + return log(p); + } else { + double t = 3.75 / abx; + double poly = + 0.39894228 + + t * (0.01328592 + + t * (0.00225319 + + t * (-0.00157565 + + t * (0.00916281 + t * (-0.02057706 + t * (0.02635537 + t * (-0.01647633 + t * 0.00392377))))))); + return abx - 0.5 * log(abx) + log(poly); + } +} + +/*! + * \brief Compute integral involving product of three modified Bessel functions. + * \param[in] beta_x Argument in x-direction. + * \param[in] beta_y Argument in y-direction. + * \param[in] beta_z Argument in z-direction. + * \return Value of the integral. + */ +inline double GetBesselIntegral(double beta_x, double beta_y, double beta_z) { + const double A = 1.0 + 2.0 * (beta_x + beta_y + beta_z); + const double Bx = 2.0 * beta_x; + const double By = 2.0 * beta_y; + const double Bz = 2.0 * beta_z; + + const int N = 4000; + const double t_max = 40.0; + const double dt = t_max / N; + + double sum = 0.0; + + for (int i = 1; i < N; i++) { + double t = i * dt; + + double e = exp(-A * t); + + double lx = GetBesselZero(Bx * t); + double ly = GetBesselZero(By * t); + double lz = GetBesselZero(Bz * t); + + double lin = log(t) - A * t + lx + ly + lz; + + double integrand = exp(lin); + sum += integrand; + } + + return sum * dt; +} + +/// @} +} // namespace RandomToolbox diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a1dc6b06c572..c6a481098310 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2944,9 +2944,42 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: DES Constant */ addDoubleOption("DES_CONST", Const_DES, 0.65); + /* DESCRIPTION: SBS lengthscale coefficient */ + addDoubleOption("SBS_LENGTHSCALE_COEFF", SBS_Cdelta, 0.1); + + /* DESCRIPTION: Maximum number of smoothing iterations for SBS model. */ + addUnsignedShortOption("SBS_MAX_ITER_SMOOTH", SBS_maxIterSmooth, 100); + + /* DESCRIPTION: SBS timescale coefficient */ + addDoubleOption("SBS_TIMESCALE_COEFF", SBS_Ctau, 0.05); + + /* DESCRIPTION: SBS intensity coefficient */ + addDoubleOption("SBS_INTENSITY_COEFF", SBS_Cmag, 1.0); + /* DESCRIPTION: Specify Hybrid RANS/LES model */ addEnumOption("HYBRID_RANSLES", Kind_HybridRANSLES, HybridRANSLES_Map, NO_HYBRIDRANSLES); + /* DESCRIPTION: Specify if the Stochastic Backscatter Model must be activated */ + addBoolOption("STOCHASTIC_BACKSCATTER", StochasticBackscatter, false); + + /* DESCRIPTION: Specify if the LES mode must be enforced */ + addBoolOption("ENFORCE_LES", enforceLES, false); + + /* DESCRIPTION: Specify if the stochastic source term must be included in the turbulence model equation */ + addBoolOption("STOCH_SOURCE_NU", stochSourceNu, true); + + /* DESCRIPTION: Enable diagnostics of the stochastic source term in Langevin equations. */ + addBoolOption("STOCH_SOURCE_DIAGNOSTICS", stochSourceDiagnostics, false); + + /* DESCRIPTION: Relaxation factor for the stochastic source term (Stochastic Backscatter Model). */ + addDoubleOption("SBS_RELAXATION_FACTOR", stochSourceRelax, 0.0); + + /* DESCRIPTION: Filter width for LES (if negative, it is computed based on the local cell size) */ + addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0); + + /* DESCRIPTION: Specify if the LD2 scheme must be employed (incompressible flows, combined with JST discretization). */ + addBoolOption("LD2_OPTION", LD2_Scheme, false); + /* DESCRIPTION: Roe with low dissipation for unsteady flows */ addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS); @@ -6499,6 +6532,44 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { case SA_ZDES: cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break; case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break; } + if (Kind_HybridRANSLES != NO_HYBRIDRANSLES) { + if (LES_FilterWidth > 0.0) cout << "User-specified LES filter width: " << LES_FilterWidth << endl; + cout << "Stochastic Backscatter: "; + if (StochasticBackscatter) { + cout << "ON" << endl; + cout << "Backscatter intensity coefficient: " << SBS_Cmag << endl; + if (SBS_Cmag < 0.0) + SU2_MPI::Error("Backscatter intensity coefficient must be non-negative.", CURRENT_FUNCTION); + cout << "Backscatter lengthscale coefficient: " << SBS_Cdelta << endl; + if (SBS_Cdelta < 0.0) + SU2_MPI::Error("Backscatter lengthscale coefficient must be non-negative.", CURRENT_FUNCTION); + cout << "Backscatter timescale coefficient: " << SBS_Ctau << endl; + if (SBS_Ctau < 0.0) + SU2_MPI::Error("Backscatter timescale coefficient must be non-negative.", CURRENT_FUNCTION); + if (SBS_maxIterSmooth > 0) + cout << "Maximum number of iterations for implicit smoothing: " << SBS_maxIterSmooth << endl; + else + cout << "No smoothing applied to source terms in Langevin equations." << endl; + if (stochSourceNu) + cout << "Stochastic source term included in turbulence model equation." << endl; + else + cout << "Stochastic source term NOT included in turbulence model equation." << endl; + if (stochSourceRelax > 0.0) + cout << "Relaxation factor for stochastic source term: " << stochSourceRelax << endl; + else + cout << "No relaxation factor for stochastic source term." << endl; + } else { + cout << "OFF" << endl; + } + } + if (StochasticBackscatter && Kind_HybridRANSLES == NO_HYBRIDRANSLES) + SU2_MPI::Error("Stochastic Backscatter can only be activated with Hybrid RANS/LES.", CURRENT_FUNCTION); + if (enforceLES) { + if (Kind_HybridRANSLES == NO_HYBRIDRANSLES) + SU2_MPI::Error("ENFORCE_LES can only be activated with Hybrid RANS/LES.", CURRENT_FUNCTION); + else + cout << "LES enforced in the whole computational domain." << endl; + } break; case MAIN_SOLVER::NEMO_EULER: if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE) cout << "Compressible two-temperature thermochemical non-equilibrium Euler equations." << endl; @@ -7053,12 +7124,26 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "First order integration." << endl; } else { - cout << "Jameson-Schmidt-Turkel scheme (2nd order in space) for the flow inviscid terms.\n"; + if (LD2_Scheme) { + cout << "Low-Dissipation Low-Dispersion (LD2) scheme for the flow inviscid terms." << endl; + if (!(Kind_Solver==MAIN_SOLVER::INC_EULER || Kind_Solver==MAIN_SOLVER::INC_NAVIER_STOKES || Kind_Solver==MAIN_SOLVER::INC_RANS)) + SU2_MPI::Error("LD2 option available for incompressible flow simulations only.", CURRENT_FUNCTION); + if (Kind_FluidModel != CONSTANT_DENSITY) + SU2_MPI::Error("LD2 option available for constant density flow simulations only.", CURRENT_FUNCTION); + if (Energy_Equation) + cout << "WARNING: LD2 option not compatible with the energy equation. JST discretization in energy equation employed." << endl; + } else { + cout << "Jameson-Schmidt-Turkel scheme (2nd order in space) for the flow inviscid terms.\n"; + } cout << "JST viscous coefficients (2nd & 4th): " << Kappa_2nd_Flow << ", " << Kappa_4th_Flow << ".\n"; cout << "The method includes a grid stretching correction (p = 0.3)."<< endl; } } + if ((Kind_ConvNumScheme_Flow != SPACE_CENTERED || (Kind_ConvNumScheme_Flow == SPACE_CENTERED && Kind_Centered_Flow == CENTERED::LAX)) && LD2_Scheme) { + SU2_MPI::Error("LD2 option available for JST scheme only.", CURRENT_FUNCTION); + } + if (Kind_ConvNumScheme_Flow == SPACE_UPWIND) { if (Kind_Upwind_Flow == UPWIND::ROE) cout << "Roe (with entropy fix = "<< EntropyFix_Coeff <<") solver for the flow inviscid terms."<< endl; if (Kind_Upwind_Flow == UPWIND::TURKEL) cout << "Roe-Turkel solver for the flow inviscid terms."<< endl; diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index d5da1b16b411..cacc6b356585 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include "../../../Common/include/CConfig.hpp" #include "../../../Common/include/linear_algebra/blas_structure.hpp" @@ -183,6 +184,14 @@ class CNumerics { roughness_j = 0.0; /*!< \brief Roughness of the wall nearest to point j. */ su2double MeanPerturbedRSM[3][3]; /*!< \brief Perturbed Reynolds stress tensor */ + su2double stochReynStress[3][3]; /*!< \brief Stochastic contribution to Reynolds stress tensor for Backscatter Model. */ + su2double stochSource[3] = {0.0}; /*!< \brief Source term for Langevin equations in Stochastic Backscatter Model. */ + su2double + stochVar_i[3], /*!< \brief Stochastic variables at point i for Stochastic Backscatter Model. */ + stochVar_j[3]; /*!< \brief Stochastic variables at point j for Stochastic Backscatter Model. */ + su2double + lesMode_i, /*!< \brief LES sensor at point i for hybrid RANS-LES methods. */ + lesMode_j; /*!< \brief LES sensor at point j for hybrid RANS-LES methods. */ 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 */ @@ -639,6 +648,38 @@ 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 + * eddy simulations." Flow, Turbulence and Combustion 99.1 (2017): 119-150. + * \param[in] nDim - Dimension of the flow problem, 2 or 3. + * \param[in] density - Density. + * \param[in] eddyVis - Eddy viscosity. + * \param[in] velGrad - Velocity gradient matrix. + * \param[out] stochReynStress - Stochastic tensor (to be added to the Reynolds stress tensor). + */ + template + NEVERINLINE static void ComputeStochReynStress(size_t nDim, Scalar density, Scalar eddyVis, + Scalar turbKE, Vector rndVec, Scalar lesSensor, + Mat& stochReynStress, Scalar Cmag) { + + /* --- Calculate stochastic tensor --- */ + + stochReynStress[1][0] = - std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[2]; + stochReynStress[2][0] = std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[1]; + stochReynStress[2][1] = - std::nearbyint(lesSensor) * Cmag * density * turbKE * rndVec[0]; + for (size_t iDim = 0; iDim < nDim; iDim++) { + for (size_t jDim = 0; jDim <= iDim; jDim++) { + if (iDim==jDim) { + stochReynStress[iDim][jDim] = 0.0; + } else { + stochReynStress[jDim][iDim] = - stochReynStress[iDim][jDim]; + } + } + } + + } /*! * \brief Project average gradient onto normal (with or w/o correction) for viscous fluxes of scalar quantities. @@ -837,6 +878,27 @@ class CNumerics { turb_ke_j = val_turb_ke_j; } + /*! + * \brief Set the stochastic variables from Langevin equations (Stochastic Backscatter Model). + * \param[in] val_stochvar_i - Value of the stochastic variable at point i. + * \param[in] val_stochvar_j - Value of the stochastic variable at point j. + * \param[in] iDim - Index of Langevin equation. + */ + inline void SetStochVar(su2double val_stochvar_i, su2double val_stochvar_j, unsigned short iDim) { + stochVar_i[iDim] = val_stochvar_i; + stochVar_j[iDim] = val_stochvar_j; + } + + /*! + * \brief Set the LES sensor for hybrid RANS-LES methods. + * \param[in] val_lesMode_i - Value of the LES sensor at point i. + * \param[in] val_lesMode_j - Value of the LES sensor at point j. + */ + inline void SetLES_Mode(su2double val_lesMode_i, su2double val_lesMode_j) { + lesMode_i = val_lesMode_i; + lesMode_j = val_lesMode_j; + } + /*! * \brief Set the value of the distance from the nearest wall. * \param[in] val_dist_i - Value of of the distance from point i to the nearest wall. @@ -847,6 +909,15 @@ class CNumerics { dist_j = val_dist_j; } + /*! + * \brief Set the stochastic source term for the Langevin equations (Backscatter Model). + * \param[in] val_stoch_source - Value of stochastic source term. + * \param[in] iDim - Index of Langevin equation. + */ + void SetStochSource(su2double val_stoch_source, unsigned short iDim) { + stochSource[iDim] = val_stoch_source; + } + /*! * \brief Set the value of the roughness from the nearest wall. * \param[in] val_dist_i - Value of of the roughness of the nearest wall from point i diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp index efff30b19a94..e9043164f1ab 100644 --- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp +++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp @@ -61,6 +61,7 @@ class CAvgGrad_Base : public CNumerics { Mean_Cp, /*!< \brief Mean value of the specific heat capacity at constant pressure. */ Mean_turb_ke, /*!< \brief Mean value of the turbulent kinetic energy. */ Mean_TauWall, /*!< \brief Mean wall shear stress (wall functions). */ + Mean_StochVar[3], /*!< \brief Mean stochastic variables (Stochastic Backscatter Model). */ TauWall_i, TauWall_j, /*!< \brief Wall shear stress at point i and j (wall functions). */ dist_ij_2, /*!< \brief Length of the edge and face, squared */ Edge_Vector[MAXNDIM] = {0.0}; /*!< \brief Vector from point i to point j. */ @@ -199,12 +200,14 @@ class CAvgGrad_Base : public CNumerics { * \param[in] val_turb_ke - Turbulent kinetic energy * \param[in] val_laminar_viscosity - Laminar viscosity. * \param[in] val_eddy_viscosity - Eddy viscosity. + * \param[in] config - Definition of the particular problem. */ void SetStressTensor(const su2double *val_primvar, const su2double* const *val_gradprimvar, su2double val_turb_ke, su2double val_laminar_viscosity, - su2double val_eddy_viscosity); + su2double val_eddy_viscosity, + const CConfig* config); /*! * \brief Get a component of the viscous stress tensor. diff --git a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp index 91b8429e0aaf..4395a4a2092b 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_convection.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_convection.hpp @@ -52,6 +52,7 @@ class CUpwScalar : public CNumerics { const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double a0 = 0.0; /*!< \brief The maximum of the face-normal velocity and 0. */ su2double a1 = 0.0; /*!< \brief The minimum of the face-normal velocity and 0. */ + su2double qij = 0.0; /*!< \brief The face-normal velocity (Langevin equations). */ su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ @@ -140,6 +141,13 @@ class CUpwScalar : public CNumerics { a1 = fmin(0.0, q_ij); } + if (config->GetStochastic_Backscatter()) { + qij = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + qij += 0.5 * (V_i[idx.Density()]*V_i[iDim + idx.Velocity()] + V_j[idx.Density()]*V_j[iDim + idx.Velocity()]) * Normal[iDim]; + } + } + FinishResidualCalc(config); AD::SetPreaccOut(Flux, nVar); diff --git a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp index 13ee01f7eda6..d92776f2f3b7 100644 --- a/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp +++ b/SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp @@ -64,7 +64,7 @@ class CAvgGrad_Scalar : public CNumerics { const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ su2double Proj_Mean_GradScalarVar[MAXNVAR]; /*!< \brief Mean_gradScalarVar DOT normal, corrected if required. */ su2double proj_vector_ij = 0.0; /*!< \brief (Edge_Vector DOT normal)/|Edge_Vector|^2 */ - su2double Flux[MAXNVAR]; /*!< \brief Final result, diffusive flux/residual. */ + su2double Flux[MAXNVAR] = {0.0}; /*!< \brief Final result, diffusive flux/residual. */ su2double* Jacobian_i[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node i. */ su2double* Jacobian_j[MAXNVAR]; /*!< \brief Flux Jacobian w.r.t. node j. */ su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */ diff --git a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp index 03c33d381384..15635dc9fcce 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_convection.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_convection.hpp @@ -48,6 +48,11 @@ class CUpwSca_TurbSA final : public CUpwScalar { using Base::ScalarVar_i; using Base::ScalarVar_j; using Base::bounded_scalar; + using Base::V_i; + using Base::V_j; + using Base::idx; + using Base::nVar; + using Base::qij; /*! * \brief Adds any extra variables to AD. @@ -59,6 +64,18 @@ class CUpwSca_TurbSA final : public CUpwScalar { * \param[in] config - Definition of the particular problem. */ void FinishResidualCalc(const CConfig* config) override { + bool backscatter = config->GetStochastic_Backscatter(); + if (backscatter) { + for (unsigned short iVar = 1; iVar < nVar; iVar++) { + Flux[iVar] = qij * 0.5 * (ScalarVar_i[iVar] + ScalarVar_j[iVar]); + } + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + for (unsigned short jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = (iVar == jVar) ? 0.5*qij : 0.0; + Jacobian_j[iVar][jVar] = (iVar == jVar) ? 0.5*qij : 0.0; + } + } + } Flux[0] = a0*ScalarVar_i[0] + a1*ScalarVar_j[0]; Jacobian_i[0][0] = a0; Jacobian_j[0][0] = a1; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 2195f94bd21f..f9cc15a146c1 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -74,6 +74,8 @@ class CSourceBase_TurbSA : public CNumerics { /*--- Residual and Jacobian ---*/ su2double Residual, *Jacobian_i; su2double Jacobian_Buffer; /*!< \brief Static storage for the Jacobian (which needs to be pointer for return type). */ + su2double* ResidSB = nullptr; + su2double** JacobianSB_i = nullptr; const FlowIndices idx; /*!< \brief Object to manage the access to the flow primitives. */ const SA_ParsedOptions options; /*!< \brief Struct with SA options. */ @@ -112,6 +114,86 @@ class CSourceBase_TurbSA : public CNumerics { Residual += yinv * dv_axi * Volume; } + /*! + * \brief Include source-term residuals for Langevin equations (Stochastic Backscatter Model) + */ + inline void ResidualStochEquations(su2double timeStep, const su2double ct, + su2double lengthScale, su2double DES_const, + const CSAVariables& var, TIME_MARCHING time_marching) { + + const su2double& nue = ScalarVar_i[0]; + const auto& density = V_i[idx.Density()]; + + const su2double nut_small = 1.0e-12; + + const su2double nut = max(nue * var.fv1, nut_small); + + su2double tLES = ct * pow(lengthScale/DES_const, 2) / nut; + su2double tRANS = pow(lengthScale, 2) / nut; + su2double tLR = tLES / tRANS; + su2double tTurb = tLES / (lesMode_i + (1.0-lesMode_i)*tLR); + su2double tRat = timeStep / tTurb; + + su2double corrFac = 1.0; + if (time_marching == TIME_MARCHING::DT_STEPPING_2ND) { + corrFac = sqrt(0.5*(1.0+tRat)*(4.0+tRat)/(2.0+tRat)); + } else if (time_marching == TIME_MARCHING::DT_STEPPING_1ST) { + corrFac = sqrt(1.0+0.5*tRat); + } + + su2double scaleFactor = 0.0; + if (std::nearbyint(lesMode_i) == 1) scaleFactor = 1.0/tTurb * sqrt(2.0/tRat) * density * corrFac; + + ResidSB[0] = Residual; + for (unsigned short iVar = 1; iVar < nVar; iVar++ ) { + ResidSB[iVar] = scaleFactor * stochSource[iVar-1] - 1.0/tTurb * density * ScalarVar_i[iVar]; + ResidSB[iVar] *= Volume; + } + + for (unsigned short iVar = 0; iVar < nVar; iVar++ ) { + for (unsigned short jVar = 0; jVar < nVar; jVar++ ) { + JacobianSB_i[iVar][jVar] = 0.0; + } + } + + for (unsigned short iVar = 1; iVar < nVar; iVar++ ) { + JacobianSB_i[iVar][iVar] = -1.0/tTurb * density * Volume; + } + + JacobianSB_i[0][0] = Jacobian_i[0]; + + } + + /*! + * \brief Include stochastic source term in the Spalart-Allmaras turbulence model equation (Stochastic Backscatter Model). + */ + template + inline void AddStochSource(const CSAVariables& var, const MatrixType& velGrad, const su2double Cmag) { + + su2double dist2 = dist_i * dist_i; + const su2double eps = 1.0e-12; + su2double xi3 = pow(var.Ji, 3); + + su2double factor = dist2 / (2.0 * var.fv1 * ScalarVar_i[0] + eps); + factor /= (3.0 * xi3 * var.cv1_3 / pow(xi3 + var.cv1_3, 2) + var.fv1 + eps); + + su2double tke = pow(var.fv1*ScalarVar_i[0]/dist_i, 2); + + su2double R12 = (nDim==3 ? Cmag * tke * ScalarVar_i[3] : 0.0); + su2double R13 = - Cmag * tke * ScalarVar_i[2]; + su2double R23 = Cmag * tke * ScalarVar_i[1]; + + su2double RGradU = R12 * (velGrad[0][1] - velGrad[1][0]) + + (nDim==3 ? R13 * (velGrad[0][2] - velGrad[2][0]) + + R23 * (velGrad[1][2] - velGrad[2][1]) : 0.0); + + su2double source_k = - RGradU; + su2double source_nu = factor * source_k; + + Residual += source_nu * Volume; + + } + public: /*! * \brief Constructor of the class. @@ -119,7 +201,7 @@ class CSourceBase_TurbSA : public CNumerics { * \param[in] config - Definition of the particular problem. */ CSourceBase_TurbSA(unsigned short nDim, const CConfig* config) - : CNumerics(nDim, 1, config), + : CNumerics(nDim, config->GetStochastic_Backscatter() ? 1+nDim : 1, config), idx(nDim, config->GetnSpecies()), options(config->GetSAParsedOptions()), axisymmetric(config->GetAxisymmetric()), @@ -127,8 +209,30 @@ class CSourceBase_TurbSA : public CNumerics { /*--- Setup the Jacobian pointer, we need to return su2double** but we know * the Jacobian is 1x1 so we use this trick to avoid heap allocation. ---*/ Jacobian_i = &Jacobian_Buffer; + /*--- Setup the Jacobian for Stochastic Backscatter Model. ---*/ + if (config->GetStochastic_Backscatter()) { + ResidSB = new su2double [1+nDim] (); + JacobianSB_i = new su2double* [1+nDim]; + for (unsigned short iVar = 0; iVar < 1+nDim; iVar++ ) { + JacobianSB_i[iVar] = new su2double [1+nDim] (); + } + } } + /*! + * \brief Destructor of the class. + */ + ~CSourceBase_TurbSA() { + if (JacobianSB_i) { + for (unsigned short iVar = 0; iVar < 1+nDim; iVar++) { + delete [] JacobianSB_i[iVar]; + } + delete [] JacobianSB_i; + } + if (ResidSB) { + delete [] ResidSB; + } + } /*! * \brief Residual for source term integration. @@ -140,10 +244,16 @@ class CSourceBase_TurbSA : public CNumerics { const auto& laminar_viscosity = V_i[idx.LaminarViscosity()]; AD::StartPreacc(); - AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, ScalarVar_i[0], Volume, dist_i, roughness_i); + AD::SetPreaccIn(density, laminar_viscosity, StrainMag_i, Volume, dist_i, roughness_i); + AD::SetPreaccIn(ScalarVar_i, nVar); AD::SetPreaccIn(Vorticity_i, 3); AD::SetPreaccIn(PrimVar_Grad_i + idx.Velocity(), nDim, nDim); - AD::SetPreaccIn(ScalarVar_Grad_i[0], nDim); + AD::SetPreaccIn(ScalarVar_Grad_i, nVar, nDim); + + bool backscatter = config->GetStochastic_Backscatter(); + if (backscatter) { + AD::SetPreaccIn(stochSource, 3); + } /*--- Common auxiliary variables and constants of the model. ---*/ CSAVariables var; @@ -252,12 +362,45 @@ class CSourceBase_TurbSA : public CNumerics { if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma); Jacobian_i[0] *= Volume; + + /*--- Compute residual for Langevin equations (Stochastic Backscatter Model). ---*/ + + if (backscatter) { + if (std::nearbyint(lesMode_i) == 1 && config->GetStochSourceNu()) { + su2double SBS_Cmag = config->GetSBS_Cmag(); + su2double lesSensor = max(lesMode_i, lesMode_j); + su2double SBS_RelaxFactor = config->GetStochSourceRelax(); + su2double intensityCoeff = SBS_Cmag; + if (SBS_RelaxFactor > 0.0) { + su2double FS_Vel = config->GetModVel_FreeStream(); + su2double ReynoldsLength = config->GetLength_Reynolds(); + su2double timeScale = ReynoldsLength / FS_Vel; + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + su2double timeStep = config->GetTime_Step(); + su2double currentTime = (timeIter - restartIter) * timeStep; + intensityCoeff = SBS_Cmag * (1.0 - exp(-SBS_RelaxFactor * currentTime / timeScale)); + } + AddStochSource(var, PrimVar_Grad_i + idx.Velocity(), intensityCoeff); + } + const su2double DES_const = config->GetConst_DES(); + const su2double ctau = config->GetSBS_Ctau(); + ResidualStochEquations(config->GetDelta_UnstTime(), ctau, dist_i, DES_const, + var, config->GetTime_Marching()); + } } - AD::SetPreaccOut(Residual); + if (backscatter) { + AD::SetPreaccOut(ResidSB, nVar); + } else { + AD::SetPreaccOut(Residual); + } AD::EndPreacc(); - return ResidualType<>(&Residual, &Jacobian_i, nullptr); + if (backscatter) + return ResidualType<>(ResidSB, JacobianSB_i, nullptr); + else + return ResidualType<>(&Residual, &Jacobian_i, nullptr); } }; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 8dc0e4b63931..1d7322f0e11c 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -434,9 +434,10 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const bool backscatter = config->GetStochastic_Backscatter(); CVariable* turbNodes = nullptr; - if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); + if (tkeNeeded || backscatter) turbNodes = solver_container[TURB_SOL]->GetNodes(); /*--- Points, coordinates and normal vector in edge ---*/ @@ -467,6 +468,33 @@ void CFVMFlowSolverBase::Viscous_Residual_impl(unsigned long iEdge, CGeome numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(jPoint,0)); + /*--- Stochastic variables from Langevin equations (Stochastic Backscatter Model). ---*/ + + if (backscatter) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + numerics->SetStochVar(turbNodes->GetSolution(iPoint, 1 + iDim), + turbNodes->GetSolution(jPoint, 1 + iDim), iDim); + } + su2double rho, eddy_visc_i, eddy_visc_j, DES_length_i, + DES_length_j, tke_i, tke_j, lesMode_i, lesMode_j; + rho = nodes->GetDensity(iPoint); + eddy_visc_i = turbNodes->GetmuT(iPoint) / rho; + eddy_visc_j = turbNodes->GetmuT(jPoint) / rho; + DES_length_i = turbNodes->GetDES_LengthScale(iPoint); + DES_length_j = turbNodes->GetDES_LengthScale(jPoint); + lesMode_i = turbNodes->GetLES_Mode(iPoint); + lesMode_j = turbNodes->GetLES_Mode(jPoint); + const su2double tol = 1e-12; + if (DES_length_i < tol || DES_length_j < tol) { + tke_i = tke_j = 0.0; + } else { + tke_i = pow(eddy_visc_i/DES_length_i, 2); + tke_j = pow(eddy_visc_j/DES_length_j, 2); + } + numerics->SetTurbKineticEnergy(tke_i, tke_j); + numerics->SetLES_Mode(lesMode_i, lesMode_j); + } + /*--- Wall shear stress values (wall functions) ---*/ numerics->SetTau_Wall(nodes->GetTau_Wall(iPoint), diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index 7220ffc32907..0aaf50225462 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -39,7 +39,13 @@ class CTurbSASolver final : public CTurbSolver { private: - su2double nu_tilde_Engine, nu_tilde_ActDisk; + su2double* nu_tilde_Engine = nullptr; + su2double* nu_tilde_ActDisk = nullptr; + su2double* ext_AverageNu = nullptr; + su2double* Res_Wall = nullptr; + su2double** Jacobian_i = nullptr; + su2double* nu_tilde_inturb = nullptr; + su2double* nu_tilde_WF = nullptr; /*! * \brief A virtual member. @@ -51,6 +57,27 @@ class CTurbSASolver final : public CTurbSolver { CGeometry *geometry, CConfig *config); + /*! + * \brief Update the source terms of the stochastic equations (Stochastic Backscatter Model). + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition. + */ + void SetLangevinSourceTerms(CConfig *config, CGeometry* geometry); + + /*! + * \brief Set seed for Langevin equations (Stochastic Backscatter Model). + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition. + */ + void SetLangevinGen(CConfig* config, CGeometry* geometry); + + /*! + * \brief Apply Laplacian smoothing to the source terms in Langevin equations (Stochastic Backscatter Model). + * \param[in] config - Definition of the particular problem. + * \param[in] geometry - Geometrical definition. + */ + void SmoothLangevinSourceTerms(CConfig* config, CGeometry* geometry); + /*! * \brief Compute nu tilde from the wall functions. * \param[in] geometry - Geometrical definition of the problem. @@ -85,7 +112,21 @@ class CTurbSASolver final : public CTurbSolver { /*! * \brief Destructor of the class. */ - ~CTurbSASolver() = default; + ~CTurbSASolver() { + + for(unsigned short iVar = 0; iVar < nVar; ++iVar) { + delete [] Jacobian_i[iVar]; + } + delete [] Jacobian_i; + + delete [] Res_Wall; + delete [] nu_tilde_Engine; + delete [] nu_tilde_ActDisk; + delete [] ext_AverageNu; + delete [] nu_tilde_inturb; + delete [] nu_tilde_WF; + + } /*! * \brief Restart residual and compute gradients. diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index 8753093b346b..db6415a0fda4 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -152,6 +152,7 @@ class CTurbSolver : public CScalarSolver { * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over * a nonlinear iteration for stability. * \param[in] allowableRatio - Maximum percentage update in variable per iteration. + * \param[in] backscatter - Flag for Stochastic Backscatter Model. */ void ComputeUnderRelaxationFactorHelper(su2double allowableRatio); }; diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index f9bf112f9819..79fd2190e078 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/SU2_CFD/include/variables/CIncNSVariable.hpp @@ -39,7 +39,8 @@ class CIncNSVariable final : public CIncEulerVariable { private: VectorType Tau_Wall; /*!< \brief Magnitude of the wall shear stress from a wall function. */ - VectorType DES_LengthScale; + VectorType DES_LengthScale; /*!< \brief DES Length Scale. */ + VectorType LES_Mode; /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES).*/ const bool Energy; /*!< \brief Flag for Energy equation in incompressible flows. */ public: @@ -133,4 +134,17 @@ class CIncNSVariable final : public CIncEulerVariable { */ inline su2double GetDES_LengthScale(unsigned long iPoint) const override { return DES_LengthScale(iPoint); } + /*! + * \brief Set the LES sensor. + */ + inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override { + LES_Mode(iPoint) = val_les_mode; + } + + /*! + * \brief Get the LES sensor. + * \return Value of the LES sensor. + */ + inline su2double GetLES_Mode(unsigned long iPoint) const override { return LES_Mode(iPoint); } + }; diff --git a/SU2_CFD/include/variables/CNEMONSVariable.hpp b/SU2_CFD/include/variables/CNEMONSVariable.hpp index 1d399a7325ca..bdda130e14c0 100644 --- a/SU2_CFD/include/variables/CNEMONSVariable.hpp +++ b/SU2_CFD/include/variables/CNEMONSVariable.hpp @@ -52,6 +52,7 @@ class CNEMONSVariable final : public CNEMOEulerVariable { VectorType Tau_Wall; /*!< \brief Magnitude of the wall shear stress from a wall function. */ VectorType DES_LengthScale; /*!< \brief DES Length Scale. */ + VectorType LES_Mode; /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES). */ VectorType Roe_Dissipation; /*!< \brief Roe low dissipation coefficient. */ VectorType Vortex_Tilting; /*!< \brief Value of the vortex tilting variable for DES length scale computation. */ diff --git a/SU2_CFD/include/variables/CNSVariable.hpp b/SU2_CFD/include/variables/CNSVariable.hpp index 69855f3d8a67..44cc1c5c60f6 100644 --- a/SU2_CFD/include/variables/CNSVariable.hpp +++ b/SU2_CFD/include/variables/CNSVariable.hpp @@ -41,6 +41,7 @@ class CNSVariable final : public CEulerVariable { VectorType Tau_Wall; /*!< \brief Magnitude of the wall shear stress from a wall function. */ VectorType DES_LengthScale; /*!< \brief DES Length Scale. */ + VectorType LES_Mode; /*!< \brief Sensor for local simulation mode (0=RANS, 1=LES).*/ VectorType Roe_Dissipation; /*!< \brief Roe low dissipation coefficient. */ VectorType Vortex_Tilting; /*!< \brief Value of the vortex tilting variable for DES length scale computation. */ @@ -185,6 +186,19 @@ class CNSVariable final : public CEulerVariable { DES_LengthScale(iPoint) = val_des_lengthscale; } +/*! + * \brief Set the LES sensor. + */ + inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override { + LES_Mode(iPoint) = val_les_mode; + } + + /*! + * \brief Get the LES sensor. + * \return Value of the LES sensor. + */ + inline su2double GetLES_Mode(unsigned long iPoint) const override { return LES_Mode(iPoint); } + /*! * \brief Set the new solution for Roe Dissipation. * \param[in] val_delta - A scalar measure of the grid size diff --git a/SU2_CFD/include/variables/CTurbSAVariable.hpp b/SU2_CFD/include/variables/CTurbSAVariable.hpp index adf3cb5cf17b..4b625b901012 100644 --- a/SU2_CFD/include/variables/CTurbSAVariable.hpp +++ b/SU2_CFD/include/variables/CTurbSAVariable.hpp @@ -28,6 +28,7 @@ #pragma once #include "CTurbVariable.hpp" +#include /*! * \class CTurbSAVariable @@ -40,7 +41,12 @@ class CTurbSAVariable final : public CTurbVariable { private: VectorType DES_LengthScale; + VectorType LES_Mode; + MatrixType stochSource; + MatrixType stochSource_old; VectorType Vortex_Tilting; + MatrixTypeGen stochGen; + VectorType besselIntegral; public: /*! @@ -73,6 +79,49 @@ class CTurbSAVariable final : public CTurbVariable { */ inline void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) override { DES_LengthScale(iPoint) = val_des_lengthscale; } + /*! + * \brief Get the source terms for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \return Value of the source term for the stochastic equations. + */ + inline su2double GetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim) const override { return stochSource(iPoint, iDim); } + + /*! + * \brief Set the source terms for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochSource - Value of the source term for the stochastic equations. + */ + inline void SetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim, su2double val_stochSource) override { stochSource(iPoint, iDim) = val_stochSource; } + + /*! + * \brief Get the old value of the source terms for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \return Old value of the source terms for the stochastic equations. + */ + inline su2double GetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim) const override { return stochSource_old(iPoint, iDim); } + + /*! + * \brief Set the old value of source terms for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochSource_old - Old value of the source term for the stochastic equations. + */ + inline void SetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim, su2double val_stochSource_old) override { stochSource_old(iPoint, iDim) = val_stochSource_old; } + + /*! + * \brief Set the LES sensor. + */ + inline void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) override { LES_Mode(iPoint) = val_les_mode; } + + /*! + * \brief Get the LES sensor. + * \return Value of the LES sensor. + */ + inline su2double GetLES_Mode(unsigned long iPoint) const override { return LES_Mode(iPoint); } + /*! * \brief Set the vortex tilting measure for computation of the EDDES length scale * \param[in] iPoint - Point index. @@ -87,4 +136,33 @@ class CTurbSAVariable final : public CTurbVariable { */ inline su2double GetVortex_Tilting(unsigned long iPoint) const override { return Vortex_Tilting(iPoint); } + /*! + * \brief Get the seed for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \return Value of the seed for the stochastic equations. + */ + inline std::mt19937 GetLangevinGen(unsigned long iPoint, unsigned short iDim) const override { return stochGen(iPoint, iDim); } + + /*! + * \brief Set the seed for the stochastic equations. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochGen - Value of the seed for the stochastic equations. + */ + inline void SetLangevinGen(unsigned long iPoint, unsigned short iDim, std::mt19937 val_stochGen) override { stochGen(iPoint, iDim) = val_stochGen; } + + /*! + * \brief Set the integral of the product of three Bessel functions appearing in Laplacian smoothing. + * \param[in] iPoint - Point index. + * \param[in] val_integral - Value of the integral. + */ + inline void SetBesselIntegral(unsigned long iPoint, su2double val_integral) override { besselIntegral(iPoint) = val_integral; } + + /*! + * \brief Get the the integral of the product of three Bessel functions appearing in Laplacian smoothing. + * \return Value of the integral. + */ + inline su2double GetBesselIntegral(unsigned long iPoint) const override { return besselIntegral(iPoint); } + }; diff --git a/SU2_CFD/include/variables/CTurbVariable.hpp b/SU2_CFD/include/variables/CTurbVariable.hpp index c0c33efd8c90..2fbec9abd244 100644 --- a/SU2_CFD/include/variables/CTurbVariable.hpp +++ b/SU2_CFD/include/variables/CTurbVariable.hpp @@ -40,7 +40,7 @@ class CTurbVariable : public CScalarVariable { VectorType muT; /*!< \brief Eddy viscosity. */ public: - static constexpr size_t MAXNVAR = 2; + static constexpr size_t MAXNVAR = 4; VectorType turb_index; VectorType intermittency; /*!< \brief Value of the intermittency for the trans. model. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index dc905d79fcb4..5a432dedcbb7 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include "../../../Common/include/CConfig.hpp" #include "../../../Common/include/containers/container_decorators.hpp" @@ -51,6 +52,7 @@ class CVariable { protected: using VectorType = C2DContainer; using MatrixType = C2DContainer; + using MatrixTypeGen = C2DContainer; MatrixType Solution; /*!< \brief Solution of the problem. */ MatrixType Solution_Old; /*!< \brief Old solution of the problem R-K. */ @@ -397,6 +399,76 @@ class CVariable { */ inline virtual void SetDES_LengthScale(unsigned long iPoint, su2double val_des_lengthscale) {} + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + */ + inline virtual su2double GetLES_Mode(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] val_les_mode - Value of the LES sensor. + */ + inline virtual void SetLES_Mode(unsigned long iPoint, su2double val_les_mode) {} + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + */ + inline virtual su2double GetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochSource - Source term in Langevin equations. + */ + inline virtual void SetLangevinSourceTerms(unsigned long iPoint, unsigned short iDim, su2double val_stochSource) {} + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + */ + inline virtual su2double GetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochSource_old - Old value of source term in Langevin equations. + */ + inline virtual void SetLangevinSourceTermsOld(unsigned long iPoint, unsigned short iDim, su2double val_stochSource_old) {} + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + */ + inline virtual std::mt19937 GetLangevinGen(unsigned long iPoint, unsigned short iDim) const {std::mt19937 gen(123); return gen; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] val_integral - Value of the integral. + */ + inline virtual void SetBesselIntegral(unsigned long iPoint, su2double val_integral) {} + + /*! + * \brief A virtual member. + */ + inline virtual su2double GetBesselIntegral(unsigned long iPoint) const { return 0.0; } + + /*! + * \brief A virtual member. + * \param[in] iPoint - Point index. + * \param[in] iDim - Dimension index. + * \param[in] val_stochGen - Pseudo-random number generator for Langevin equations. + */ + inline virtual void SetLangevinGen(unsigned long iPoint, unsigned short iDim, std::mt19937 val_stochGen) {} + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index 3d1b85502659..ec47f3f580fd 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -311,6 +311,8 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi su2double U_i[5] = {0.0}, U_j[5] = {0.0}; su2double ProjGridVel = 0.0; + bool LD2_Scheme = config->GetLD2_Scheme(); + const su2double alpha_LD2 = 0.36; /*--- Primitive variables at point i and j ---*/ @@ -326,6 +328,31 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi for (iDim = 0; iDim < nDim; iDim++) { Velocity_i[iDim] = V_i[iDim+1]; Velocity_j[iDim] = V_j[iDim+1]; + } + + if (LD2_Scheme) { + su2double d_ij[3] = {0.0}; + for (iDim = 0; iDim < nDim; iDim++) + d_ij[iDim] = Coord_j[iDim]-Coord_i[iDim]; + su2double velGrad_i[3][3] = {{0.0}}; + su2double velGrad_j[3][3] = {{0.0}}; + for (unsigned short jDim = 0; jDim < nDim; jDim++) { + for (iDim = 0; iDim < nDim; iDim++) { + velGrad_i[iDim][jDim] = PrimVar_Grad_i[iDim+1][jDim]; + velGrad_j[iDim][jDim] = PrimVar_Grad_j[iDim+1][jDim]; + } + } + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] += alpha_LD2 * (velGrad_i[iDim][0] * d_ij[0] + + velGrad_i[iDim][1] * d_ij[1] + + velGrad_i[iDim][2] * d_ij[2]); + Velocity_j[iDim] -= alpha_LD2 * (velGrad_j[iDim][0] * d_ij[0] + + velGrad_j[iDim][1] * d_ij[1] + + velGrad_j[iDim][2] * d_ij[2]); + } + } + + for (iDim = 0; iDim < nDim; iDim++) { MeanVelocity[iDim] = 0.5*(Velocity_i[iDim]+Velocity_j[iDim]); sq_vel_i += 0.5*Velocity_i[iDim]*Velocity_i[iDim]; sq_vel_j += 0.5*Velocity_j[iDim]*Velocity_j[iDim]; diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp index 4c67cc13305c..902b8a6f0f06 100644 --- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp +++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp @@ -120,7 +120,8 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, const su2double* const *val_gradprimvar, const su2double val_turb_ke, const su2double val_laminar_viscosity, - const su2double val_eddy_viscosity) { + const su2double val_eddy_viscosity, + const CConfig* config) { const su2double Density = val_primvar[nDim+2]; @@ -140,6 +141,15 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar, // turb_ke is not considered in the stress tensor, see #797 ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, su2double(0.0)); } + + /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/ + + if (config->GetStochastic_Backscatter()) { + for (unsigned short iDim = 0 ; iDim < nDim; iDim++) + for (unsigned short jDim = 0 ; jDim < nDim; jDim++) + tau[iDim][jDim] += stochReynStress[iDim][jDim]; + } + } void CAvgGrad_Base::SetHeatFluxVector(const su2double* const* val_gradprimvar, const su2double val_eddy_viscosity, @@ -445,10 +455,33 @@ CNumerics::ResidualType<> CAvgGrad_Flow::ComputeResidual(const CConfig* config) Mean_turb_ke, MeanPerturbedRSM); } + /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/ + + if (config->GetStochastic_Backscatter()) { + for (iVar = 0; iVar < nDim; iVar++) + Mean_StochVar[iVar] = 0.5*(stochVar_i[iVar] + stochVar_j[iVar]); + su2double SBS_Cmag = config->GetSBS_Cmag(); + su2double lesSensor = max(lesMode_i, lesMode_j); + su2double SBS_RelaxFactor = config->GetStochSourceRelax(); + su2double intensityCoeff = SBS_Cmag; + if (SBS_RelaxFactor > 0.0) { + su2double FS_Vel = config->GetModVel_FreeStream(); + su2double ReynoldsLength = config->GetLength_Reynolds(); + su2double timeScale = ReynoldsLength / FS_Vel; + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + su2double timeStep = config->GetTime_Step(); + su2double currentTime = (timeIter - restartIter) * timeStep; + intensityCoeff = SBS_Cmag * (1.0 - exp(-SBS_RelaxFactor * currentTime / timeScale)); + } + ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_turb_ke, + Mean_StochVar, lesSensor, stochReynStress, intensityCoeff); + } + /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + 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); @@ -619,9 +652,33 @@ CNumerics::ResidualType<> CAvgGradInc_Flow::ComputeResidual(const CConfig* confi Mean_turb_ke, MeanPerturbedRSM); } + /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/ + + if (config->GetStochastic_Backscatter()) { + for (iVar = 0; iVar < nDim; iVar++) + Mean_StochVar[iVar] = 0.5*(stochVar_i[iVar] + stochVar_j[iVar]); + su2double SBS_Cmag = config->GetSBS_Cmag(); + su2double lesSensor = max(lesMode_i, lesMode_j); + su2double SBS_RelaxFactor = config->GetStochSourceRelax(); + su2double intensityCoeff = SBS_Cmag; + if (SBS_RelaxFactor > 0.0) { + su2double FS_Vel = config->GetModVel_FreeStream(); + su2double ReynoldsLength = config->GetLength_Reynolds(); + su2double timeScale = ReynoldsLength / FS_Vel; + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + su2double timeStep = config->GetTime_Step(); + su2double currentTime = (timeIter - restartIter) * timeStep; + intensityCoeff = SBS_Cmag * (1.0 - exp(-SBS_RelaxFactor * currentTime / timeScale)); + } + ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_turb_ke, + Mean_StochVar, lesSensor, stochReynStress, intensityCoeff); + } + /*--- Get projected flux tensor (viscous residual) ---*/ + SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + 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); @@ -944,10 +1001,33 @@ CNumerics::ResidualType<> CGeneralAvgGrad_Flow::ComputeResidual(const CConfig* c Mean_turb_ke, MeanPerturbedRSM); } + /* --- If the Stochastic Backscatter Model is active, add random contribution to stress tensor ---*/ + + if (config->GetStochastic_Backscatter()) { + for (iVar = 0; iVar < nDim; iVar++) + Mean_StochVar[iVar] = 0.5*(stochVar_i[iVar] + stochVar_j[iVar]); + su2double SBS_Cmag = config->GetSBS_Cmag(); + su2double lesSensor = max(lesMode_i, lesMode_j); + su2double SBS_RelaxFactor = config->GetStochSourceRelax(); + su2double intensityCoeff = SBS_Cmag; + if (SBS_RelaxFactor > 0.0) { + su2double FS_Vel = config->GetModVel_FreeStream(); + su2double ReynoldsLength = config->GetLength_Reynolds(); + su2double timeScale = ReynoldsLength / FS_Vel; + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + su2double timeStep = config->GetTime_Step(); + su2double currentTime = (timeIter - restartIter) * timeStep; + intensityCoeff = SBS_Cmag * (1.0 - exp(-SBS_RelaxFactor * currentTime / timeScale)); + } + ComputeStochReynStress(nDim, Mean_PrimVar[nDim+2], Mean_Eddy_Viscosity, Mean_turb_ke, + Mean_StochVar, lesSensor, stochReynStress, intensityCoeff); + } + /*--- Get projected flux tensor (viscous residual) ---*/ SetStressTensor(Mean_PrimVar, Mean_GradPrimVar, Mean_turb_ke, - Mean_Laminar_Viscosity, Mean_Eddy_Viscosity); + 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/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 2476498e84d7..b96cdc6fbf32 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -234,6 +234,8 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("GRID_VELOCITY-Z", "Grid_Velocity_z", "GRID_VELOCITY", "z-component of the grid velocity vector"); } + AddVolumeOutput("VELOCITY_DIVERGENCE", "Velocity_Divergence", "DERIVED", "Divergence of the velocity field"); + // Primitive variables AddVolumeOutput("PRESSURE", "Pressure", "PRIMITIVE", "Pressure"); AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature"); @@ -325,6 +327,13 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv SetVolumeOutputValue("ENERGY", iPoint, Node_Flow->GetSolution(iPoint, 3)); } + const auto VelocityGradient = Node_Flow->GetVelocityGradient(iPoint); + su2double divVel = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + divVel += VelocityGradient[iDim][iDim]; + } + SetVolumeOutputValue("VELOCITY_DIVERGENCE", iPoint, divVel); + if (gridMovement){ SetVolumeOutputValue("GRID_VELOCITY-X", iPoint, Node_Geo->GetGridVel(iPoint)[0]); SetVolumeOutputValue("GRID_VELOCITY-Y", iPoint, Node_Geo->GetGridVel(iPoint)[1]); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 9d408892cf6e..f6e28ce131ef 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -322,6 +322,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("GRID_VELOCITY-Z", "Grid_Velocity_z", "GRID_VELOCITY", "z-component of the grid velocity vector"); } + AddVolumeOutput("VELOCITY_DIVERGENCE", "Velocity_Divergence", "DERIVED", "Divergence of the velocity field"); + // Primitive variables AddVolumeOutput("PRESSURE_COEFF", "Pressure_Coefficient", "PRIMITIVE", "Pressure coefficient"); AddVolumeOutput("DENSITY", "Density", "PRIMITIVE", "Density"); @@ -432,6 +434,13 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("GRID_VELOCITY-Z", iPoint, Node_Geo->GetGridVel(iPoint)[2]); } + const auto VelocityGradient = Node_Flow->GetVelocityGradient(iPoint); + su2double divVel = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + divVel += VelocityGradient[iDim][iDim]; + } + SetVolumeOutputValue("VELOCITY_DIVERGENCE", iPoint, divVel); + const su2double factor = solver[FLOW_SOL]->GetReferenceDynamicPressure(); SetVolumeOutputValue("PRESSURE_COEFF", iPoint, (Node_Flow->GetPressure(iPoint) - solver[FLOW_SOL]->GetPressure_Inf())/factor); SetVolumeOutputValue("DENSITY", iPoint, Node_Flow->GetDensity(iPoint)); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 696a4c595475..44338b645a58 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -964,6 +964,14 @@ void CFlowOutput::AddHistoryOutputFields_ScalarRMS_RES(const CConfig* config) { case TURB_FAMILY::SA: /// DESCRIPTION: Root-mean square residual of nu tilde (SA model). AddHistoryOutput("RMS_NU_TILDE", "rms[nu]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of nu tilde (SA model).", HistoryFieldType::RESIDUAL); + if (config->GetStochastic_Backscatter()) { + /// DESCRIPTION: Root-mean square residual of stochastic vector x-component (Stochastic Backscatter Model). + AddHistoryOutput("RMS_STOCH_VAR-X", "rms[stoch_x]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector x-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Root-mean square residual of stochastic vector y-component (Stochastic Backscatter Model). + AddHistoryOutput("RMS_STOCH_VAR-Y", "rms[stoch_y]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector y-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Root-mean square residual of stochastic vector z-component (Stochastic Backscatter Model). + if (nDim==3) AddHistoryOutput("RMS_STOCH_VAR-Z", "rms[stoch_z]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of stochastic vector z-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + } break; case TURB_FAMILY::KW: @@ -1019,6 +1027,14 @@ void CFlowOutput::AddHistoryOutputFields_ScalarMAX_RES(const CConfig* config) { case TURB_FAMILY::SA: /// DESCRIPTION: Maximum residual of nu tilde (SA model). AddHistoryOutput("MAX_NU_TILDE", "max[nu]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of nu tilde (SA model).", HistoryFieldType::RESIDUAL); + if (config->GetStochastic_Backscatter()) { + /// DESCRIPTION: Maximum residual of stochastic vector x-component (Stochastic Backscatter Model). + AddHistoryOutput("MAX_STOCH_VAR-X", "max[stoch_x]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector x-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of stochastic vector y-component (Stochastic Backscatter Model). + AddHistoryOutput("MAX_STOCH_VAR-Y", "max[stoch_y]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector y-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + /// DESCRIPTION: Maximum residual of stochastic vector z-component (Stochastic Backscatter Model). + if (nDim==3) AddHistoryOutput("MAX_STOCH_VAR-Z", "max[stoch_z]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of stochastic vector z-component (Stochastic Backscatter Model).", HistoryFieldType::RESIDUAL); + } break; case TURB_FAMILY::KW: @@ -1160,8 +1176,21 @@ void CFlowOutput::LoadHistoryDataScalar(const CConfig* config, const CSolver* co case TURB_FAMILY::SA: SetHistoryOutputValue("RMS_NU_TILDE", log10(solver[TURB_SOL]->GetRes_RMS(0))); SetHistoryOutputValue("MAX_NU_TILDE", log10(solver[TURB_SOL]->GetRes_Max(0))); + if (config->GetStochastic_Backscatter()) { + SetHistoryOutputValue("RMS_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_RMS(1))); + SetHistoryOutputValue("RMS_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_RMS(2))); + if (nDim==3) SetHistoryOutputValue("RMS_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_RMS(3))); + SetHistoryOutputValue("MAX_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_Max(1))); + SetHistoryOutputValue("MAX_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_Max(2))); + if (nDim==3) SetHistoryOutputValue("MAX_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_Max(3))); + } if (multiZone) { SetHistoryOutputValue("BGS_NU_TILDE", log10(solver[TURB_SOL]->GetRes_BGS(0))); + if (config->GetStochastic_Backscatter()) { + SetHistoryOutputValue("BGS_STOCH_VAR-X", log10(solver[TURB_SOL]->GetRes_BGS(1))); + SetHistoryOutputValue("BGS_STOCH_VAR-Y", log10(solver[TURB_SOL]->GetRes_BGS(2))); + if (nDim==3) SetHistoryOutputValue("BGS_STOCH_VAR-Z", log10(solver[TURB_SOL]->GetRes_BGS(3))); + } } break; @@ -1483,6 +1512,16 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) { if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { AddVolumeOutput("DES_LENGTHSCALE", "DES_LengthScale", "DDES", "DES length scale value"); AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "DDES", "Wall distance value"); + AddVolumeOutput("LES_SENSOR","LES_Sensor","DDES","LES sensor value"); + if (config->GetStochastic_Backscatter()) { + AddVolumeOutput("STOCHVAR_X", "StochVar_x", "BACKSCATTER", "x-component of the stochastic vector potential"); + AddVolumeOutput("STOCHVAR_Y", "StochVar_y", "BACKSCATTER", "y-component of the stochastic vector potential"); + if (nDim==3) AddVolumeOutput("STOCHVAR_Z", "StochVar_z", "BACKSCATTER", "z-component of the stochastic vector potential"); + AddVolumeOutput("STOCHSOURCE_X", "StochSource_x", "BACKSCATTER", "x-component of the stochastic source vector"); + AddVolumeOutput("STOCHSOURCE_Y", "StochSource_y", "BACKSCATTER", "y-component of the stochastic source vector"); + if (nDim==3) AddVolumeOutput("STOCHSOURCE_Z", "StochSource_z", "BACKSCATTER", "z-component of the stochastic source vector"); + AddVolumeOutput("ENERGY_BACKSCATTER_RATIO", "Energy_Backscatter_Ratio", "BACKSCATTER", "Energy backscatter from unresolved to resolved scales (divided by the turbulent dissipation of resolved kinetic energy)"); + } } if (config->GetViscous()) { @@ -1582,6 +1621,84 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { SetVolumeOutputValue("DES_LENGTHSCALE", iPoint, Node_Flow->GetDES_LengthScale(iPoint)); SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint)); + SetVolumeOutputValue("LES_SENSOR", iPoint, Node_Flow->GetLES_Mode(iPoint)); + if (config->GetStochastic_Backscatter()) { + SetVolumeOutputValue("STOCHVAR_X", iPoint, Node_Turb->GetSolution(iPoint, 1)); + SetVolumeOutputValue("STOCHVAR_Y", iPoint, Node_Turb->GetSolution(iPoint, 2)); + if (nDim==3) SetVolumeOutputValue("STOCHVAR_Z", iPoint, Node_Turb->GetSolution(iPoint, 3)); + SetVolumeOutputValue("STOCHSOURCE_X", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 0)); + SetVolumeOutputValue("STOCHSOURCE_Y", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 1)); + if (nDim==3) SetVolumeOutputValue("STOCHSOURCE_Z", iPoint, Node_Turb->GetLangevinSourceTerms(iPoint, 2)); + const su2double rho = Node_Flow->GetDensity(iPoint); + const su2double nu_t = Node_Flow->GetEddyViscosity(iPoint) / rho; + const su2double DES_lengthscale = max(Node_Flow->GetDES_LengthScale(iPoint), 1.0E-10); + const su2double lesSensor = Node_Flow->GetLES_Mode(iPoint); + const su2double mag = config->GetSBS_Cmag(); + const su2double tke_estim = pow(nu_t/DES_lengthscale, 2.0); + const su2double csi_x = Node_Turb->GetSolution(iPoint, 1); + const su2double csi_y = Node_Turb->GetSolution(iPoint, 2); + const su2double csi_z = (nDim==3) ? Node_Turb->GetSolution(iPoint, 3) : 0.0; + const su2double R_xy = mag * tke_estim * csi_z * std::nearbyint(lesSensor); + const su2double R_xz = - mag * tke_estim * csi_y * std::nearbyint(lesSensor); + const su2double R_yz = mag * tke_estim * csi_x * std::nearbyint(lesSensor); + const auto vel_grad = Node_Flow->GetVelocityGradient(iPoint); + const su2double vel_div = vel_grad(0,0) + vel_grad(1,1) + (nDim ==3 ? vel_grad(2,2) : 0.0); + 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)); + su2double tau_zz=0.0, tau_xz=0.0, tau_yz=0.0; + if (nDim == 3){ + tau_zz = nu_t * (2*vel_grad(2,2) - (2.0/3.0)*vel_div); + tau_xz = nu_t * (vel_grad(0,2) + vel_grad(2,0)); + tau_yz = nu_t * (vel_grad(1,2) + vel_grad(2,1)); + } + const su2double energy_res_to_mod = tau_xx * vel_grad(0,0) + tau_yy * vel_grad(1,1) + + (nDim==3 ? tau_zz * vel_grad(2,2) : 0.0) + + tau_xy * (vel_grad(0,1) + vel_grad(1,0)) + + (nDim==3 ? tau_xz * (vel_grad(0,2) + vel_grad(2,0)) + + tau_yz * (vel_grad(1,2) + vel_grad(2,1)) : 0.0); + const su2double energy_backscatter = R_xy * (vel_grad(0,1) - vel_grad(1,0)) + + (nDim==3 ? R_xz * (vel_grad(0,2) - vel_grad(2,0)) + + R_yz * (vel_grad(1,2) - vel_grad(2,1)) : 0.0); + const su2double energy_backscatter_ratio = energy_backscatter / (energy_res_to_mod + 1e-12); + SetVolumeOutputValue("ENERGY_BACKSCATTER_RATIO", iPoint, energy_backscatter_ratio); + } + } + + if (config->GetTime_Domain()) { + const su2double rho = Node_Flow->GetDensity(iPoint); + const su2double nu_t = Node_Flow->GetEddyViscosity(iPoint) / rho; + const auto vel_grad = Node_Flow->GetVelocityGradient(iPoint); + const su2double vel_div = vel_grad(0,0) + vel_grad(1,1) + (nDim ==3 ? vel_grad(2,2) : 0.0); + 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); + 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); + } + if (config->GetKind_HybridRANSLES()!=NO_HYBRIDRANSLES && config->GetStochastic_Backscatter()) { + const su2double DES_lengthscale = max(Node_Flow->GetDES_LengthScale(iPoint), 1.0E-10); + const su2double lesSensor = Node_Flow->GetLES_Mode(iPoint); + const su2double mag = config->GetSBS_Cmag(); + const su2double tke_estim = pow(nu_t/DES_lengthscale, 2.0); + const su2double csi_x = Node_Turb->GetSolution(iPoint, 1); + const su2double csi_y = Node_Turb->GetSolution(iPoint, 2); + const su2double csi_z = (nDim==3) ? Node_Turb->GetSolution(iPoint, 3) : 0.0; + const su2double R_xy = mag * tke_estim * csi_z * std::nearbyint(lesSensor); + const su2double R_xz = - mag * tke_estim * csi_y * std::nearbyint(lesSensor); + const su2double R_yz = mag * tke_estim * csi_x * std::nearbyint(lesSensor); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_12", iPoint, -R_xy); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_13", iPoint, -R_xz); + SetAvgVolumeOutputValue("STOCHASTIC_REYNOLDS_STRESS_23", iPoint, -R_yz); + } } switch (config->GetKind_Species_Model()) { @@ -1661,6 +1778,13 @@ void CFlowOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver SetVolumeOutputValue("SKIN_FRICTION-Z", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 2)); SetVolumeOutputValue("HEAT_FLUX", iPoint, solver[heat_sol]->GetHeatFlux(iMarker, iVertex)); SetVolumeOutputValue("Y_PLUS", iPoint, solver[FLOW_SOL]->GetYPlus(iMarker, iVertex)); + + if (config->GetTime_Domain()) { + SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-X", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 0)); + SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-Y", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 1)); + if (nDim == 3) + SetAvgVolumeOutputValue("MEAN_SKIN_FRICTION-Z", iPoint, solver[FLOW_SOL]->GetCSkinFriction(iMarker, iVertex, 2)); + } } void CFlowOutput::AddAerodynamicCoefficients(const CConfig* config) { @@ -3968,6 +4092,11 @@ void CFlowOutput::SetTimeAveragedFields() { AddVolumeOutput("MEAN_VELOCITY-Z", "MeanVelocity_z", "TIME_AVERAGE", "Mean velocity z-component"); AddVolumeOutput("MEAN_PRESSURE", "MeanPressure", "TIME_AVERAGE", "Mean pressure"); + AddVolumeOutput("MEAN_SKIN_FRICTION-X", "MeanSkinFriction_x", "TIME_AVERAGE", "Mean skin friction x-component"); + AddVolumeOutput("MEAN_SKIN_FRICTION-Y", "MeanSkinFriction_y", "TIME_AVERAGE", "Mean skin friction y-component"); + if (nDim==3) + AddVolumeOutput("MEAN_SKIN_FRICTION-Z", "MeanSkinFriction_z", "TIME_AVERAGE", "Mean skin friction z-component"); + AddVolumeOutput("RMS_U", "RMS[u]", "TIME_AVERAGE", "RMS u"); AddVolumeOutput("RMS_V", "RMS[v]", "TIME_AVERAGE", "RMS v"); AddVolumeOutput("RMS_UV", "RMS[uv]", "TIME_AVERAGE", "RMS uv"); @@ -3984,6 +4113,19 @@ void CFlowOutput::SetTimeAveragedFields() { AddVolumeOutput("UWPRIME", "w'u'", "TIME_AVERAGE", "Mean Reynolds-stress component w'u'"); AddVolumeOutput("VWPRIME", "w'v'", "TIME_AVERAGE", "Mean Reynolds-stress component w'v'"); } + + 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"); + 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("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"); } void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *Node_Flow){ @@ -3994,7 +4136,6 @@ void CFlowOutput::LoadTimeAveragedData(unsigned long iPoint, const CVariable *No SetAvgVolumeOutputValue("MEAN_VELOCITY-Z", iPoint, Node_Flow->GetVelocity(iPoint,2)); SetAvgVolumeOutputValue("MEAN_PRESSURE", iPoint, Node_Flow->GetPressure(iPoint)); - SetAvgVolumeOutputValue("RMS_U", iPoint, pow(Node_Flow->GetVelocity(iPoint,0),2)); SetAvgVolumeOutputValue("RMS_V", iPoint, pow(Node_Flow->GetVelocity(iPoint,1),2)); SetAvgVolumeOutputValue("RMS_UV", iPoint, Node_Flow->GetVelocity(iPoint,0) * Node_Flow->GetVelocity(iPoint,1)); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 709aa7b689bb..b35471b72efa 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1134,6 +1134,7 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool jst_scheme = ((config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0)); const bool bounded_scalar = config->GetBounded_Scalar(); + const bool LD2_Scheme = config->GetLD2_Scheme(); /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of * variables, otherwise switch to the faster adjoint evaluation mode. ---*/ @@ -1171,6 +1172,16 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co numerics->SetSensor(nodes->GetSensor(iPoint), nodes->GetSensor(jPoint)); } + if (LD2_Scheme) { + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(jPoint)); + if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint) + && !geometry->nodes->GetPeriodicBoundary(jPoint))) { + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); + } else { + numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); + } + } + /*--- Grid movement ---*/ if (dynamic_grid) { diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 28747c7e1d71..b4d5859a997c 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -360,7 +360,7 @@ void CIncNSSolver::Compute_Enthalpy_Diffusion(unsigned long iEdge, CGeometry* ge unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) { unsigned long iPoint, nonPhysicalPoints = 0; - su2double eddy_visc = 0.0, turb_ke = 0.0, DES_LengthScale = 0.0; + su2double eddy_visc = 0.0, turb_ke = 0.0, DES_LengthScale = 0.0, LES_Mode = 0.0; const su2double* scalar = nullptr; const TURB_MODEL turb_model = config->GetKind_Turb_Model(); const SPECIES_MODEL species_model = config->GetKind_Species_Model(); @@ -380,6 +380,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES){ DES_LengthScale = solver_container[TURB_SOL]->GetNodes()->GetDES_LengthScale(iPoint); + LES_Mode = solver_container[TURB_SOL]->GetNodes()->GetLES_Mode(iPoint); } } @@ -400,6 +401,10 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c nodes->SetDES_LengthScale(iPoint,DES_LengthScale); + /*--- Set the LES sensor ---*/ + + nodes->SetLES_Mode(iPoint, LES_Mode); + } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index ab90660ef22e..33684c4936e8 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -150,6 +150,8 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, cons if (config->GetKind_HybridRANSLES() != NO_HYBRIDRANSLES) { su2double DES_LengthScale = solver_container[TURB_SOL]->GetNodes()->GetDES_LengthScale(iPoint); nodes->SetDES_LengthScale(iPoint, DES_LengthScale); + su2double LES_Mode = solver_container[TURB_SOL]->GetNodes()->GetLES_Mode(iPoint); + nodes->SetLES_Mode(iPoint, LES_Mode); } } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 99be4b04cc40..51baaef76acf 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1356,6 +1356,10 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar+1; MPI_TYPE = COMM_TYPE_DOUBLE; break; + case MPI_QUANTITIES::STOCH_SOURCE_LANG: + COUNT_PER_POINT = nDim; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; case MPI_QUANTITIES::SOLUTION_FEA: if (config->GetTime_Domain()) COUNT_PER_POINT = nVar*3; @@ -1483,6 +1487,10 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); bufDSend[buf_offset+nVar] = base_nodes->GetmuT(iPoint); break; + case MPI_QUANTITIES::STOCH_SOURCE_LANG: + for (iDim = 0; iDim < nDim; iDim++) + bufDSend[buf_offset+iDim] = base_nodes->GetLangevinSourceTerms(iPoint, iDim); + break; case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetUndivided_Laplacian(iPoint, iVar); @@ -1631,6 +1639,10 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); base_nodes->SetmuT(iPoint,bufDRecv[buf_offset+nVar]); break; + case MPI_QUANTITIES::STOCH_SOURCE_LANG: + for (iDim = 0; iDim < nDim; iDim++) + base_nodes->SetLangevinSourceTerms(iPoint, iDim, bufDRecv[buf_offset+iDim]); + break; case MPI_QUANTITIES::UNDIVIDED_LAPLACIAN: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]); diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 7649165d0be3..6b050fb38418 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -30,6 +30,7 @@ #include "../../include/variables/CFlowVariable.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../../Common/include/toolboxes/random_toolbox.hpp" CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned short iMesh, CFluidModel* FluidModel) @@ -54,6 +55,11 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor nDim = geometry->GetnDim(); + /*--- Add Langevin equations if the Stochastic Backscatter Model is used ---*/ + + bool backscatter = config->GetStochastic_Backscatter(); + if (backscatter) { nVar += nDim; nVarGrad = nPrimVar = nVar; } + /*--- Single grid simulation ---*/ if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) { @@ -109,17 +115,24 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor } Solution_Inf[0] = nu_tilde_Inf; + if (backscatter) { + for (unsigned short iVar = 1; iVar < nVar; iVar++) { + Solution_Inf[iVar] = 0.0; + } + } /*--- Factor_nu_Engine ---*/ - Factor_nu_Engine = config->GetNuFactor_Engine(); - nu_tilde_Engine = Factor_nu_Engine*Viscosity_Inf/Density_Inf; + Factor_nu_Engine = config->GetNuFactor_Engine(); + nu_tilde_Engine = new su2double [nVar] (); + nu_tilde_Engine[0] = Factor_nu_Engine*Viscosity_Inf/Density_Inf; if (config->GetSAParsedOptions().bc) { - nu_tilde_Engine = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf; + nu_tilde_Engine[0] = 0.005*Factor_nu_Engine*Viscosity_Inf/Density_Inf; } /*--- Factor_nu_ActDisk ---*/ - Factor_nu_ActDisk = config->GetNuFactor_Engine(); - nu_tilde_ActDisk = Factor_nu_ActDisk*Viscosity_Inf/Density_Inf; + Factor_nu_ActDisk = config->GetNuFactor_Engine(); + nu_tilde_ActDisk = new su2double [nVar] (); + nu_tilde_ActDisk[0] = Factor_nu_ActDisk*Viscosity_Inf/Density_Inf; /*--- Eddy viscosity at infinity ---*/ su2double Ji, Ji_3, fv1, cv1_3 = 7.1*7.1*7.1; @@ -129,6 +142,16 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor fv1 = Ji_3/(Ji_3+cv1_3); muT_Inf = Density_Inf*fv1*nu_tilde_Inf; + /*--- Allocate memory for boundary conditions ---*/ + ext_AverageNu = new su2double [nVar] (); + Res_Wall = new su2double [nVar] (); + Jacobian_i = new su2double* [nVar]; + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar] = new su2double [nVar] (); + } + nu_tilde_inturb = new su2double [nVar] (); + nu_tilde_WF = new su2double [nVar] (); + /*--- Initialize the solution to the far-field state everywhere. ---*/ nodes = new CTurbSAVariable(nu_tilde_Inf, muT_Inf, nPoint, nDim, nVar, config); @@ -155,8 +178,16 @@ CTurbSASolver::CTurbSASolver(CGeometry *geometry, CConfig *config, unsigned shor * due to arbitrary number of turbulence variables ---*/ Inlet_TurbVars.resize(nMarker); - for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) + for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) { Inlet_TurbVars[iMarker].resize(nVertex[iMarker],nVar) = nu_tilde_Inf; + if (backscatter) { + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { + for (unsigned short iVar = 1; iVar < nVar; iVar++) { + Inlet_TurbVars[iMarker](iVertex,iVar) = 0.0; + } + } + } + } /*--- Store the initial CFL number for all grid points. ---*/ @@ -203,6 +234,35 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe SetDES_LengthScale(solver_container, geometry, config); + /*--- Compute source terms for Langevin equations ---*/ + + bool backscatter = config->GetStochastic_Backscatter(); + unsigned long innerIter = config->GetInnerIter(); + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + + if (backscatter && innerIter==0) { + SetLangevinGen(config, geometry); + SetLangevinSourceTerms(config, geometry); + const unsigned short maxIter = config->GetSBS_maxIterSmooth(); + bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)); + if (maxIter>0) SmoothLangevinSourceTerms(config, geometry); + if (timeIter == restartIter) { + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned short iVar = 1; iVar < nVar; iVar++) { + const su2double randomSource = nodes->GetLangevinSourceTerms(iPoint, iVar-1); + nodes->SetSolution(iPoint, iVar, randomSource); + nodes->SetSolution_Old(iPoint, iVar, randomSource); + if (dual_time) { + nodes->Set_Solution_time_n(iPoint, iVar, randomSource); + nodes->Set_Solution_time_n1(iPoint, iVar, randomSource); + } + } + } + } + } + } } @@ -378,6 +438,14 @@ void CTurbSASolver::Source_Residual(CGeometry *geometry, CSolver **solver_contai numerics->SetDistance(nodes->GetDES_LengthScale(iPoint), 0.0); + /*--- Compute source terms in Langevin equations (Stochastic Basckscatter Model) ---*/ + + if (config->GetStochastic_Backscatter()) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + numerics->SetStochSource(nodes->GetLangevinSourceTerms(iPoint, iDim), iDim); + numerics->SetLES_Mode(nodes->GetLES_Mode(iPoint), 0.0); + } + } /*--- Effective Intermittency ---*/ @@ -401,7 +469,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); @@ -503,19 +571,20 @@ void CTurbSASolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_conta su2double coeff = (nu_total/sigma); su2double RoughWallBC = nodes->GetSolution(iPoint,0)/(0.03*Roughness_Height); - su2double Res_Wall;// = new su2double [nVar]; - Res_Wall = coeff*RoughWallBC*Area; - LinSysRes.SubtractBlock(iPoint, &Res_Wall); + Res_Wall[0] = coeff*RoughWallBC*Area; + LinSysRes.SubtractBlock(iPoint, Res_Wall); - su2double Jacobian_i = (laminar_viscosity /density *Area)/(0.03*Roughness_Height*sigma); - Jacobian_i += 2.0*RoughWallBC*Area/sigma; - if (implicit) Jacobian.AddVal2Diag(iPoint, -Jacobian_i); + Jacobian_i[0][0] = (laminar_viscosity /density *Area)/(0.03*Roughness_Height*sigma); + Jacobian_i[0][0] += 2.0*RoughWallBC*Area/sigma; + if (implicit) Jacobian.AddVal2Diag(iPoint, 0, -Jacobian_i[0][0]); } } } END_SU2_OMP_FOR } + + void CTurbSASolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { @@ -558,7 +627,7 @@ void CTurbSASolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, CN conv_numerics->SetPrimitive(V_domain, V_inlet); /*--- Non-dimensionalize Inlet_TurbVars if Inlet-Files are used. ---*/ - su2double Inlet_Vars[MAXNVAR]; + su2double Inlet_Vars[MAXNVAR] = {0.0}; Inlet_Vars[0] = Inlet_TurbVars[val_marker][iVertex][0]; if (config->GetInlet_Profile_From_File()) { Inlet_Vars[0] *= config->GetDensity_Ref() / config->GetViscosity_Ref(); @@ -863,7 +932,7 @@ void CTurbSASolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_cont /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde_Engine); + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_Engine); /*--- Set various other quantities in the conv_numerics class ---*/ @@ -1016,7 +1085,7 @@ void CTurbSASolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, } else { /*--- Outflow analysis ---*/ - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde_ActDisk); + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_ActDisk); } /*--- Grid Movement ---*/ @@ -1077,6 +1146,8 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c su2double extAverageNu = solver_container[FLOW_SOL]->GetExtAverageNu(val_marker, iSpan); + ext_AverageNu[0] = extAverageNu; + /*--- Loop over all the vertices on this boundary marker ---*/ SU2_OMP_FOR_STAT(OMP_MIN_SIZE) @@ -1111,7 +1182,7 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &extAverageNu); + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), ext_AverageNu); /*--- Set various other quantities in the conv_numerics class ---*/ @@ -1141,7 +1212,7 @@ void CTurbSASolver::BC_Inlet_MixingPlane(CGeometry *geometry, CSolver **solver_c /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &extAverageNu); + visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), ext_AverageNu); visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); @@ -1216,7 +1287,9 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain /*--- Set the turbulent variable states (prescribed for an inflow) ---*/ - conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde); + nu_tilde_inturb[0] = nu_tilde; + + conv_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_inturb); if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -1246,7 +1319,7 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain /*--- Turbulent variables w/o reconstruction, and its gradients ---*/ - visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde); + visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), nu_tilde_inturb); visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint), nodes->GetGradient(iPoint)); @@ -1337,7 +1410,9 @@ void CTurbSASolver::SetTurbVars_WF(CGeometry *geometry, CSolver **solver_contain if (counter > max_iter) break; } - nodes->SetSolution_Old(iPoint_Neighbor, &nu_til); + nu_tilde_WF[0] = nu_til; + + nodes->SetSolution_Old(iPoint_Neighbor, nu_tilde_WF); LinSysRes.SetBlock_Zero(iPoint_Neighbor); /*--- includes 1 in the diagonal ---*/ @@ -1396,7 +1471,9 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC su2double psi_2 = (1.0 - (cb1/(cw1*k2*fw_star))*(ft2 + (1.0 - ft2)*fv2))/(fv1 * max(1.0e-10,1.0-ft2)); psi_2 = min(100.0,psi_2); - su2double lengthScale = 0.0; + su2double lengthScale = 0.0, lesSensor = 0.0; + + const su2double LES_FilterWidth = config->GetLES_FilterWidth(); switch(kindHybridRANSLES){ case SA_DES: { @@ -1405,9 +1482,18 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC 1997 ---*/ - const su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); + su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); + if (LES_FilterWidth > 0.0){ + maxDelta = LES_FilterWidth; + } const su2double distDES = constDES * maxDelta; lengthScale = min(distDES,wallDistance); + lesSensor = (wallDistance<=distDES) ? 0.0 : 1.0; + + if (config->GetEnforceLES()) { + lengthScale = distDES; + lesSensor = 1.0; + } break; } @@ -1417,13 +1503,22 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC Theoretical and Computational Fluid Dynamics - 2006 ---*/ - const su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); + su2double maxDelta = geometry->nodes->GetMaxLength(iPoint); + if (LES_FilterWidth > 0.0){ + maxDelta = LES_FilterWidth; + } const su2double r_d = (kinematicViscosityTurb+kinematicViscosity)/(uijuij*k2*pow(wallDistance, 2.0)); const su2double f_d = 1.0-tanh(pow(8.0*r_d,3.0)); const su2double distDES = constDES * maxDelta; lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES)); + lesSensor = (wallDistance<=distDES) ? 0.0 : f_d; + + if (config->GetEnforceLES()) { + lengthScale = distDES; + lesSensor = 1.0; + } break; } @@ -1462,8 +1557,17 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC maxDelta = deltaDDES; } + if (LES_FilterWidth > 0.0){ + maxDelta = LES_FilterWidth; + } const su2double distDES = constDES * maxDelta; lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES)); + lesSensor = (wallDistance<=distDES) ? 0.0 : f_d; + + if (config->GetEnforceLES()) { + lengthScale = distDES; + lesSensor = 1.0; + } break; } @@ -1514,19 +1618,261 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC maxDelta = deltaDDES; } + if (LES_FilterWidth > 0.0){ + maxDelta = LES_FilterWidth; + } const su2double distDES = constDES * maxDelta; lengthScale = wallDistance-f_d*max(0.0,(wallDistance-distDES)); + lesSensor = (wallDistance<=distDES) ? 0.0 : f_d; + + if (config->GetEnforceLES()) { + lengthScale = distDES; + lesSensor = 1.0; + } break; } } nodes->SetDES_LengthScale(iPoint, lengthScale); + nodes->SetLES_Mode(iPoint, lesSensor); } END_SU2_OMP_FOR } +void CTurbSASolver::SetLangevinSourceTerms(CConfig *config, CGeometry* geometry) { + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++){ + for (auto iDim = 0u; iDim < nDim; iDim++){ + auto gen = nodes->GetLangevinGen(iPoint, iDim); + su2double lesSensor = nodes->GetLES_Mode(iPoint); + su2double rnd = RandomToolbox::GetRandomNormal(gen) * std::nearbyint(lesSensor); + nodes->SetLangevinSourceTermsOld(iPoint, iDim, rnd); + nodes->SetLangevinSourceTerms(iPoint, iDim, rnd); + } + } + END_SU2_OMP_FOR + +} + +void CTurbSASolver::SetLangevinGen(CConfig* config, CGeometry* geometry) { + + unsigned long timeStep = config->GetTimeIter(); + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++){ + const auto iGlobalPoint = geometry->nodes->GetGlobalIndex(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++){ + unsigned long seed = RandomToolbox::GetSeed(iGlobalPoint+1,iDim+1); + std::mt19937 gen(seed + timeStep); + nodes->SetLangevinGen(iPoint, iDim, gen); + } + } + END_SU2_OMP_FOR + +} + +void CTurbSASolver::SmoothLangevinSourceTerms(CConfig* config, CGeometry* geometry) { + + const su2double LES_FilterWidth = config->GetLES_FilterWidth(); + const su2double constDES = config->GetConst_DES(); + const su2double cDelta = config->GetSBS_Cdelta(); + const unsigned short maxIter = config->GetSBS_maxIterSmooth(); + const su2double tol = -5.0; + const su2double sourceLim = 3.0; + const su2double omega = 0.8; + unsigned long timeIter = config->GetTimeIter(); + unsigned long restartIter = config->GetRestart_Iter(); + + /*--- 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); + CompleteComms(geometry, config, MPI_QUANTITIES::STOCH_SOURCE_LANG); + + su2double localResNorm = 0.0; + unsigned long local_nPointLES = 0; + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + const su2double lesSensor = nodes->GetLES_Mode(iPoint); + if (std::nearbyint(lesSensor) == 0.0) continue; + local_nPointLES += 1; + const su2double DES_LengthScale = nodes->GetDES_LengthScale(iPoint); + su2double maxDelta = DES_LengthScale / constDES; + 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 source_i = nodes->GetLangevinSourceTerms(iPoint, iDim); + su2double source_i_old = nodes->GetLangevinSourceTermsOld(iPoint, iDim); + auto coord_i = geometry->nodes->GetCoord(iPoint); + + /*--- Assemble system matrix. ---*/ + + su2double diag = 1.0; + su2double sum = 0.0; + for (unsigned short iNode = 0; iNode < geometry->nodes->GetnPoint(iPoint); iNode++) { + auto jPoint = geometry->nodes->GetPoint(iPoint, iNode); + auto coord_j = geometry->nodes->GetCoord(jPoint); + auto iEdge = geometry->nodes->GetEdge(iPoint, iNode); + auto* normal = geometry->edges->GetNormal(iEdge); + su2double area = GeometryToolbox::Norm(nDim, normal); + su2double dx_ij = GeometryToolbox::Distance(nDim, coord_i, coord_j); + su2double source_j = nodes->GetLangevinSourceTerms(jPoint, iDim); + su2double a_ij = area/volume_iPoint * b2/dx_ij; + diag += a_ij; + sum += a_ij * source_j; + } + + /*--- Update the solution. ---*/ + + su2double source_tmp = (source_i_old + sum) / diag; + localResNorm += pow(omega * (source_tmp - source_i), 2); + source_i = (1.0-omega)*source_i + omega*source_tmp; + nodes->SetLangevinSourceTerms(iPoint, iDim, source_i); + + } + END_SU2_OMP_FOR + + /*--- 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_INT, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&localResNorm, &globalResNorm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + globalResNorm = sqrt(globalResNorm / global_nPointLES); + + if (rank == MASTER_NODE) { + if (iter == 0) { + cout << endl + << "Residual of Laplacian smoothing along dimension " << iDim+1 << "." << endl + << "---------------------------------" << endl + << " Iter RMS Residual" << endl + << "---------------------------------" << endl; + } + if (iter%10 == 0) { + cout << " " + << std::setw(5) << iter + << " " + << std::setw(12) << std::fixed << std::setprecision(6) << log10(globalResNorm) + << endl; + } + } + + if (log10(globalResNorm) < tol || iter == maxIter-1) { + + if (rank == MASTER_NODE) { + 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; + 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; + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + const su2double lesSensor = nodes->GetLES_Mode(iPoint); + if (std::nearbyint(lesSensor) == 0.0) continue; + su2double source = nodes->GetLangevinSourceTerms(iPoint, iDim); + su2double source_notSmoothed = nodes->GetLangevinSourceTermsOld(iPoint, iDim); + mean_check_old += source; + var_check_old += source * source; + mean_check_notSmoothed += source_notSmoothed; + var_check_notSmoothed += source_notSmoothed * source_notSmoothed; + su2double integral = 0.0; + if (timeIter==restartIter) { + const su2double DES_LengthScale = nodes->GetDES_LengthScale(iPoint); + su2double maxDelta = DES_LengthScale / constDES; + if (LES_FilterWidth > 0.0) maxDelta = LES_FilterWidth; + su2double b = sqrt(cDelta) * maxDelta; + su2double b2 = b * b; + auto coord_i = geometry->nodes->GetCoord(iPoint); + su2double max_dx = 0.0; + su2double max_dy = 0.0; + su2double max_dz = 0.0; + unsigned short nNeigh = geometry->nodes->GetnPoint(iPoint); + for (unsigned short iNode = 0; iNode < nNeigh; iNode++) { + auto jPoint = geometry->nodes->GetPoint(iPoint, iNode); + auto coord_j = geometry->nodes->GetCoord(jPoint); + su2double dx = fabs(coord_i[0]-coord_j[0]); + su2double dy = fabs(coord_i[1]-coord_j[1]); + su2double dz = 0.0; + if (nDim == 3) dz = fabs(coord_i[2]-coord_j[2]); + if (dx > max_dx) max_dx = dx; + if (dy > max_dy) max_dy = dy; + if (dz > max_dz) max_dz = dz; + } + su2double dx2 = max_dx * max_dx; + su2double dy2 = max_dy * max_dy; + su2double dz2 = max_dz * max_dz; + su2double bx = b2 / dx2; + su2double by = b2 / dy2; + su2double bz = 0.0; + if (nDim == 3) bz = b2 / dz2; + integral = RandomToolbox::GetBesselIntegral(bx, by, bz); + nodes->SetBesselIntegral(iPoint, integral); + } else { + integral = nodes->GetBesselIntegral(iPoint); + } + su2double scaleFactor = 1.0 / sqrt(max(integral, 1e-10)); + source *= scaleFactor; + source = min(max(source, -sourceLim), sourceLim); + mean_check_new += source; + var_check_new += source * source; + nodes->SetLangevinSourceTerms(iPoint, iDim, source); + } + if (config->GetStochSourceDiagnostics()) { + su2double mean_check_old_G = 0.0; + su2double mean_check_new_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_new, &mean_check_new_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; + mean_check_new_G /= global_nPointLES; + 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 before scaling: " << mean_check_old_G-mean_check_notSmoothed_G <<". After scaling: " << mean_check_new_G-mean_check_notSmoothed_G << "." << endl; + cout << "Variance of stochastic source term before scaling: " << var_check_old_G/var_check_notSmoothed_G <<". After scaling: " << var_check_new_G/var_check_notSmoothed_G << "." << endl; + cout << endl; + } + } + break; + + } + } + } + +} + void CTurbSASolver::SetInletAtVertex(const su2double *val_inlet, unsigned short iMarker, unsigned long iVertex) { @@ -1561,7 +1907,8 @@ void CTurbSASolver::ComputeUnderRelaxationFactor(const CConfig *config) { SA_NEG model is more robust due to allowing for negative nu_tilde, so the under-relaxation is not applied to that variant. */ - if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG) return; + if (config->GetSAParsedOptions().version == SA_OPTIONS::NEG || + config->GetStochastic_Backscatter()) return; /* Loop over the solution update given by relaxing the linear system for this nonlinear iteration. */ diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index e31109fcd662..9efffc6066ff 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -37,6 +37,7 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su StrainMag.resize(nPoint); Tau_Wall.resize(nPoint) = su2double(-1.0); DES_LengthScale.resize(nPoint) = su2double(0.0); + LES_Mode.resize(nPoint) = su2double(0.0); Max_Lambda_Visc.resize(nPoint); /*--- Allocate memory for the AuxVar and its gradient. See e.g. CIncEulerSolver::Source_Residual: * Axisymmetric: total-viscosity * y-vel / y-coord diff --git a/SU2_CFD/src/variables/CNEMONSVariable.cpp b/SU2_CFD/src/variables/CNEMONSVariable.cpp index e673ad5603b7..b379ce2752df 100644 --- a/SU2_CFD/src/variables/CNEMONSVariable.cpp +++ b/SU2_CFD/src/variables/CNEMONSVariable.cpp @@ -70,6 +70,7 @@ CNEMONSVariable::CNEMONSVariable(su2double val_pressure, StrainMag.resize(nPoint) = su2double(0.0); Tau_Wall.resize(nPoint) = su2double(-1.0); DES_LengthScale.resize(nPoint) = su2double(0.0); + LES_Mode.resize(nPoint) = su2double(0.0); Roe_Dissipation.resize(nPoint) = su2double(0.0); Vortex_Tilting.resize(nPoint) = su2double(0.0); Max_Lambda_Visc.resize(nPoint) = su2double(0.0); diff --git a/SU2_CFD/src/variables/CNSVariable.cpp b/SU2_CFD/src/variables/CNSVariable.cpp index 6c360efc995e..753f7b4feebf 100644 --- a/SU2_CFD/src/variables/CNSVariable.cpp +++ b/SU2_CFD/src/variables/CNSVariable.cpp @@ -39,6 +39,7 @@ CNSVariable::CNSVariable(su2double density, const su2double *velocity, su2double StrainMag.resize(nPoint) = su2double(0.0); Tau_Wall.resize(nPoint) = su2double(-1.0); DES_LengthScale.resize(nPoint) = su2double(0.0); + LES_Mode.resize(nPoint) = su2double(0.0); Roe_Dissipation.resize(nPoint) = su2double(0.0); Vortex_Tilting.resize(nPoint) = su2double(0.0); Max_Lambda_Visc.resize(nPoint) = su2double(0.0); diff --git a/SU2_CFD/src/variables/CTurbSAVariable.cpp b/SU2_CFD/src/variables/CTurbSAVariable.cpp index 75119c89ed5f..384cc70e51cc 100644 --- a/SU2_CFD/src/variables/CTurbSAVariable.cpp +++ b/SU2_CFD/src/variables/CTurbSAVariable.cpp @@ -28,12 +28,22 @@ #include "../../include/variables/CTurbSAVariable.hpp" - CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) : CTurbVariable(npoint, ndim, nvar, config) { - Solution_Old = Solution = val_nu_tilde; + /*--- Initialize solution (check if the Stochastic Backscatter Model is active) ---*/ + bool backscatter = config->GetStochastic_Backscatter(); + if (!backscatter) { + Solution_Old = Solution = val_nu_tilde; + } else { + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + Solution_Old(iPoint, 0) = Solution(iPoint, 0) = val_nu_tilde; + for (unsigned short iVar = 1; iVar < nVar; iVar++) { + Solution_Old(iPoint, iVar) = Solution(iPoint, iVar) = 0.0; + } + } + } muT.resize(nPoint) = val_muT; @@ -47,7 +57,13 @@ CTurbSAVariable::CTurbSAVariable(su2double val_nu_tilde, su2double val_muT, unsi } DES_LengthScale.resize(nPoint) = su2double(0.0); + LES_Mode.resize(nPoint) = su2double(0.0); Vortex_Tilting.resize(nPoint); + stochSource.resize(nPoint, nDim) = su2double(0.0); + stochSource_old.resize(nPoint, nDim) = su2double(0.0); + stochGen.resize(nPoint, nDim); + besselIntegral.resize(nPoint); + } void CTurbSAVariable::SetVortex_Tilting(unsigned long iPoint, CMatrixView PrimGrad_Flow, diff --git a/TestCases/backscatter/DIHT/backscatter_DIHT.cfg b/TestCases/backscatter/DIHT/backscatter_DIHT.cfg new file mode 100644 index 000000000000..6fd2aa724855 --- /dev/null +++ b/TestCases/backscatter/DIHT/backscatter_DIHT.cfg @@ -0,0 +1,107 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Decaying Isotropic Homogeneous Turbulence (DIHT) % +% Author: Angelo Passariello % +% Institution: University of Naples Federico II % +% Date: 2025.10.01 % +% File Version 8.3.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% + +SOLVER= RANS +KIND_TURB_MODEL= SA +HYBRID_RANSLES= SA_DES +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ------------- STOCHASTIC BACKSCATTER MODEL PARAMETERS -----------------------% + +STOCHASTIC_BACKSCATTER= YES +SBS_CTAU= 0.001 + +% ------------- DECAYING ISOTROPIC HOMOGENEOUS TURBULENCE ---------------------% + +DIHT= YES +DIHT_DOMAIN_LENGTH= (0.5588, 0.5588, 0.5588) +DIHT_NPOINT= (64, 64, 64) +DIHT_NUM_MODES= 800 + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% + +MACH_NUMBER= 0.001 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 10129 +REYNOLDS_LENGTH= 0.5588 +FREESTREAM_NU_FACTOR= 1.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% + +REF_LENGTH= 1.0 +REF_AREA= 1.0 + +% ------------------------- UNSTEADY SIMULATION -------------------------------% + +TIME_DOMAIN= NO +%TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER +%TIME_STEP= 0.001 +%MAX_TIME= 20 +%UNST_CFL_NUMBER= 0.0 +INNER_ITER= 0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% + +MARKER_PERIODIC= (xmin, xmax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588, 0.0, 0.0, ymin, ymax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588, 0.0, zmin, zmax, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5588) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% + +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +TIME_ITER= 1 +LINEAR_SOLVER_PREC= ILU +LOW_MACH_PREC= YES +%LOW_MACH_CORR= YES + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% + +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +MUSCL_TURB= NO +VENKAT_LIMITER_COEFF= 0.05 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= SLAU2 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +TIME_DISCRE_TURB= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -15 +CONV_STARTITER= 0 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= cube_64.su2 +MESH_FORMAT= SU2 +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow +VOLUME_FILENAME= flow +VOLUME_OUTPUT= DDES, BACKSCATTER +SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_MOMENTUM-Z, RMS_NU_TILDE, RMS_STOCH_VAR_X, RMS_STOCH_VAR_Y, RMS_STOCH_VAR_Z) +HISTORY_OUTPUT= (TIME_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y, RMS_MOMENTUM-Z, RMS_NU_TILDE) +HISTORY_WRT_FREQ_INNER= 1000000 +OUTPUT_WRT_FREQ= 1000000 +OUTPUT_FILES= (RESTART, PARAVIEW) diff --git a/TestCases/ddes/flatplate/ddes_flatplate.cfg b/TestCases/ddes/flatplate/ddes_flatplate.cfg index 5c1d0804b3b4..03b780bcb067 100644 --- a/TestCases/ddes/flatplate/ddes_flatplate.cfg +++ b/TestCases/ddes/flatplate/ddes_flatplate.cfg @@ -13,9 +13,11 @@ % SOLVER= RANS KIND_TURB_MODEL= SA -HYBRID_RANSLES= SA_EDDES +HYBRID_RANSLES= SA_DES +STOCHASTIC_BACKSCATTER= YES MATH_PROBLEM= DIRECT RESTART_SOL= NO +RESTART_ITER= 100 % ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% % @@ -41,7 +43,7 @@ TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER % % U_inf = 69.4448 - dt*=0.02 - dt=0.000288 TIME_STEP= 0.000288 -MAX_TIME= 20.0 +MAX_TIME= 5 UNST_CFL_NUMBER= 0.0 INNER_ITER= 20 @@ -61,7 +63,7 @@ CFL_NUMBER= 10.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -TIME_ITER= 10 +TIME_ITER= 500 % ----------------------- SLOPE LIMITER DEFINITION ----------------------------% % @@ -103,5 +105,5 @@ VOLUME_ADJ_FILENAME= adjoint GRAD_OBJFUNC_FILENAME= of_grad SURFACE_FILENAME= surface_flow SURFACE_ADJ_FILENAME= surface_adjoint -OUTPUT_WRT_FREQ= 1000 +OUTPUT_WRT_FREQ= 1000000 SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_NU_TILDE, LIFT, DRAG, TOTAL_HEATFLUX) diff --git a/config_template.cfg b/config_template.cfg index 5bb7f5ead4f1..4d206f2e9b9b 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -180,7 +180,37 @@ HYBRID_RANSLES= SA_DDES % % DES Constant (0.65) DES_CONST= 0.65 - +% +% User-specified LES filter width (if negative, it is computed based on the local cell size, default -1.0) +LES_FILTER_WIDTH= -1.0 +% +% Stochastic Backscatter Model (NO, YES) +STOCHASTIC_BACKSCATTER= NO +% +% Backscatter lengthscale coefficient (0.1) +SBS_LENGTHSCALE_COEFF= 0.1 +% +% Maximum number of smoothing iterations for SBS model (100) +SBS_MAX_ITER_SMOOTH= 100 +% +% Backscatter timescale coefficient (0.05) +SBS_TIMESCALE_COEFF= 0.05 +% +% Backscatter intensity coefficient (1.0) +SBS_INTENSITY_COEFF= 1.0 +% +% Relaxation factor for the stochastic source term in the main balance equations (0.0) +SBS_RELAXATION_FACTOR= 0.0 +% +% Include stochastic source in turbulence model equation (NO, YES) +STOCH_SOURCE_NU= YES +% +% Include diagnostics of the stochastic source term in Langevin equations (NO, YES) +STOCH_SOURCE_DIAGNOSTICS= NO +% +% Enforce LES mode in Hybrid RANS-LES Simulations (NO, YES) +ENFORCE_LES= NO +% % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% % % Mach number (non-dimensional, based on the free-stream values) @@ -1627,6 +1657,9 @@ SMOOTH_GEOMETRY= 0 % SW, MSW, FDS, SLAU, SLAU2, L2ROE, LMROE) CONV_NUM_METHOD_FLOW= ROE % +% Option to employ the Low-Dissipation Low-Dispersion (LD2) scheme for incompressible flows (NO, YES) +LD2_OPTION= NO +% % Roe Low Dissipation function for Hybrid RANS/LES simulations (FD, NTS, NTS_DUCROS) ROE_LOW_DISSIPATION= FD %