[WIP] Stochastic Backscatter Model for Grey Area Mitigation in Hybrid RANS-LES Simulations

Common/include/CConfig.hpp

@@ -1088,6 +1088,13 @@ 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 enforceLES;                     /*!< \brief Option to enforce LES mode in hybrid RANS-LES simulations. */
+  su2double LES_FilterWidth;           /*!< \brief LES filter width for hybrid RANS-LES simulations. */
   unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
 
   unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
@@ -9518,6 +9525,24 @@ 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 the LES Filter Width.
+   * \return Value of LES Filter Width.
+   */
+  su2double GetLES_FilterWidth(void) const { return LES_FilterWidth; }
+
   /*!
    * \brief Get the Kind of Roe Low Dissipation Scheme for Unsteady flows.
    * \return Value of Low dissipation approach.
@@ -9530,6 +9555,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.
    */

Common/include/option_structure.hpp

@@ -2645,6 +2645,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). */

Common/include/toolboxes/random_toolbox.hpp

@@ -0,0 +1,149 @@
+/*!
+ * \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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <random>
+#include <functional>
+
+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<double>{}(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<double> 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<double> 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;
+        t = t * t;
+        return 1.0 + t*(3.5156229 +
+             t*(3.0899424 +
+             t*(1.2067492 +
+             t*(0.2659732 +
+             t*(0.0360768 +
+             t*0.0045813)))));
+    } else {
+        double t = 3.75/abx;
+        double ans = (exp(abx)/sqrt(abx)) *
+            (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 ans;
+    }
+}
+
+/*!
+ * \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 delta_t = t_max / N;
+    double sum = 0.0;
+    for (int i = 0; i < N; i++) {
+        double t   = i * delta_t;
+        double I0x = GetBesselZero(Bx * t);
+        double I0y = GetBesselZero(By * t);
+        double I0z = GetBesselZero(Bz * t);
+        double integrand = t * exp(-A * t)
+                             * (I0x * I0y * I0z);
+        sum += integrand;
+    }
+    return sum * delta_t;
+}
+
+/// @}
+}  // namespace RandomToolbox

Common/src/CConfig.cpp

@@ -2921,9 +2921,30 @@ 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: Filter width for LES (if negative, it is computed based on the local cell size) */
+  addDoubleOption("LES_FILTER_WIDTH", LES_FilterWidth, -1.0);
+
   /* DESCRIPTION: Roe with low dissipation for unsteady flows */
   addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS);
 
@@ -6474,6 +6495,36 @@ 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;
+          } 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;