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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft

AngPass wants to merge 55 commits into su2code:develop from AngPass:develop

+1,552 −75

Common/include/CConfig.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -1097,6 +1097,16 @@ 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. */
       unsigned short Kind_RoeLowDiss;      /*!< \brief Kind of Roe scheme with low dissipation for unsteady flows. */
       unsigned short nSpanWiseSections; /*!< \brief number of span-wise sections */
@@ Expand Down Expand Up / @@ -9595,6 +9605,42 @@ 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 the Kind of Roe Low Dissipation Scheme for Unsteady flows.
        * \return Value of Low dissipation approach.
@@ Expand All / @@ -9607,6 +9653,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.
+       */
+      unsigned short 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.
        */
@@ Expand Down @@

Common/include/option_structure.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -820,14 +820,16 @@ enum class CENTERED { @@
       JST,            /*!< \brief Jameson-Smith-Turkel centered numerical method. */
       LAX,            /*!< \brief Lax-Friedrich centered numerical method. */
       JST_MAT,        /*!< \brief JST with matrix dissipation. */
-      JST_KE          /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */
+      JST_KE,         /*!< \brief Kinetic Energy preserving Jameson-Smith-Turkel centered numerical method. */
+      LD2             /*!< \brief Low-Dissipation Low-Dispersion (LD2) centered scheme. */
     };
     static const MapType<std::string, CENTERED> Centered_Map = {
       MakePair("NONE", CENTERED::NONE)
       MakePair("JST", CENTERED::JST)
       MakePair("JST_KE", CENTERED::JST_KE)
       MakePair("JST_MAT", CENTERED::JST_MAT)
       MakePair("LAX-FRIEDRICH", CENTERED::LAX)
+      MakePair("LD2", CENTERED::LD2)
     };
@@ Expand Down Expand Up / @@ -2697,6 +2699,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). */
@@ Expand Down @@

Common/include/toolboxes/random_toolbox.hpp

-Original file line number
+Diff line change
@@ -0,0 +1,153 @@
+    /*!
+     * \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 <cstdint>
+    namespace RandomToolbox {
+    /// \addtogroup RandomToolbox
+    /// @{
+    /*!
+     * \brief SplitMix64 hash function for 64-bit integers.
+     * \param[in] x Input value to hash.
+     * \return Hashed 64-bit output.
+     */
+    static inline uint64_t splitmix64(uint64_t x) {
+        x += 0x9e3779b97f4a7c15ULL;                  // golden ratio offset
+        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; // first mixing step
+        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; // second mixing step
+        return x ^ (x >> 31);                        // final avalanche
+    }
+    /*!
+     * \brief Generate a deterministic 64-bit hash from three integers.
+     * \param[in] nodeIndex Global node index.
+     * \param[in] iDim Dimension index.
+     * \param[in] timeIter Current time iteration of the simulation.
+     * \return 64-bit hash value.
+     */
+    inline uint64_t GetHash(unsigned long nodeIndex, unsigned short iDim, unsigned long timeIter) {
+        uint64_t x = nodeIndex;
+        x ^= splitmix64(iDim);
+        x ^= splitmix64(timeIter);
+        return splitmix64(x);
+    }
+    /*!
+     * \brief Convert a 64-bit hash into a uniform double in (0,1].
+     * Uses the top 53 bits of the hash to fill the mantissa of a double.
+     * Ensures the result is never zero, suitable for Box-Muller transform.
+     * \param[in] x 64-bit hash.
+     * \return Uniform double in the interval (0,1].
+     */
+    inline double HashToUniform(uint64_t x) {
+        constexpr double inv53 = 1.0 / 9007199254740992.0; // 1/2^53
+        uint64_t uInt = x >> 11;                           // top 53 bits
+        return (uInt + 1) * inv53;                         // map to (0,1]
+    }
+    /*!
+     * \brief Generate a standard normal random number from a 64-bit hash.
+     * Uses two deterministic uniforms derived from the hash and its bitwise NOT
+     * as inputs to the Box-Muller transform.
+     * \param[in] x 64-bit hash.
+     * \return Standard normal random number (mean=0, stddev=1).
+     */
+    inline double HashToNormal(uint64_t x) {
+        constexpr double pi = 3.14159265358979323846;
+        double u = HashToUniform(x);    // first uniform
+        double v = HashToUniform(~x);   // second uniform (bitwise NOT)
+        double r = sqrt(-2.0 * log(u));
+        double theta = 2.0 * pi * v;
+        return r * cos(theta);          // one normal sample
+    }
+    /*!
+     * \brief Generate a deterministic standard normal number for a cell, dimension, and timestep.
+     *
+     * Combines hashing and Box-Muller in one function.
+     *
+     * \param[in] nodeIndex Global node index.
+     * \param[in] dim Dimension index.
+     * \param[in] timeIter Simulation timestep (1-based).
+     * \return Standard normal random number.
+     */
+    inline double GetNormal(unsigned long nodeIndex, unsigned long dim, unsigned long timeIter) {
+        uint64_t hash = GetHash(nodeIndex, dim, timeIter);
+        return HashToNormal(hash);
+    }
+    /*!
+     * \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 = abx / 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 the product of three modified Bessel functions.
+     *  Useful for scaling the smoothed stochastic source terms in Langevin equations.
+     * \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 = 20.0;
+      const double dt = t_max / N;
+      double sum = 0.0;
+      for (int i = 1; i <= N; i++) {
+        double t = i * dt;
+        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);
+        if (i==N) integrand *= 0.5;
+        sum += integrand;
+      }
+      return sum * dt;
+    }
+    /// @}
+    }  // namespace RandomToolbox

Common/src/CConfig.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -2960,9 +2960,39 @@ 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.02);
+      /* 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, false);
+      /* 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: Roe with low dissipation for unsteady flows */
       addEnumOption("ROE_LOW_DISSIPATION", Kind_RoeLowDiss, RoeLowDiss_Map, NO_ROELOWDISS);
@@ Expand Down Expand Up @@
               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;
@@ Expand Down Expand Up @@
               cout << "Lax viscous coefficients (1st): " << Kappa_1st_Flow << ".\n";
               cout << "First order integration." << endl;
             }
+            else if (Kind_Centered_Flow == CENTERED::LD2) {
+              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 scheme not yet implemented for the compressible flow solver.", CURRENT_FUNCTION);
+              if (Kind_FluidModel != CONSTANT_DENSITY)
+                SU2_MPI::Error("LD2 scheme available for constant density flows only.", CURRENT_FUNCTION);
+              if (Energy_Equation)
+                cout << "WARNING: Current implementation of the LD2 scheme not compatible with the energy equation. JST employed in energy equation instead." << 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;
+                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_UPWIND) {
@@ Expand Down @@

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

Diff view

Diff view

There are no files selected for viewing

Uh oh!

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

Are you sure you want to change the base?

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

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!