les.h

Go to the documentation of this file.

 
#ifndef LES_H
#define LES_H
 
#include <string>
#include <unordered_map>
 
#include <Kokkos_Core.hpp>
#include <toml.hpp>
 
#include "common.h"
#include "log.h"
 
enum class SGSModel : uint8_t { SMAGORINSKY, WALE, VREMAN, SIGMA };
 
static const std::unordered_map<std::string, SGSModel> SGS_MODELS = {
    {"smagorinsky", SGSModel::SMAGORINSKY},
    {"wale", SGSModel::WALE},
    {"vreman", SGSModel::VREMAN},
    {"sigma", SGSModel::SIGMA},
};
 
struct LES {
    SGSModel model = SGSModel::SIGMA;
    rtype C = 1.35_r;
    rtype Pr_t = 0.9_r;
    rtype Sc_t = 0.9_r;
    bool scotti = false;
    bool dynamic = false;  // C from the Germano identity over the whole domain, once per step (see Solver)  // the eddy viscosity's width V^(1/3) f(a_1, a_2) (Scotti, Meneveau & Lilly 1993)
 
    static LES from_input(const toml::value & input);
 
    static bool in_input(const toml::value & input) { return input.contains("les"); }
 
    static double scotti_factor(double h1, double h2, double h3);
 
    static rtype default_constant(SGSModel model);
 
    logging::Items summary() const;
 
    KOKKOS_INLINE_FUNCTION
    static double operator_of(const SGSModel model, const double g[3][3]) {
        switch (model) {
            case SGSModel::SMAGORINSKY:
                return smagorinsky(g);
            case SGSModel::WALE:
                return wale(g);
            case SGSModel::VREMAN:
                return vreman(g);
            default:
                return sigma(g);
        }
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype nu_t(const double g[3][3], const rtype delta) const {
        const double c = static_cast<double>(C);
        const double coefficient = model == SGSModel::VREMAN ? c : c * c;
        const double d = static_cast<double>(delta);
        return static_cast<rtype>(coefficient * d * d * operator_of(model, g));
    }
 
    KOKKOS_INLINE_FUNCTION
    static double strain_squared(const double g[3][3]) {
        double s = 0.0;
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                const double s_ij = 0.5 * (g[i][j] + g[j][i]);
                s += s_ij * s_ij;
            }
        }
        return s;
    }
 
    KOKKOS_INLINE_FUNCTION
    static double smagorinsky(const double g[3][3]) { return Kokkos::sqrt(2.0 * strain_squared(g)); }
 
    KOKKOS_INLINE_FUNCTION
    static double wale(const double g[3][3]) {
        double g2[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                g2[i][j] = 0.0;
                for (int k = 0; k < 3; k++) g2[i][j] += g[i][k] * g[k][j];
            }
        }
        const double trace = (g2[0][0] + g2[1][1] + g2[2][2]) / 3.0;
        double sd = 0.0;
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                const double s = 0.5 * (g2[i][j] + g2[j][i]) - (i == j ? trace : 0.0);
                sd += s * s;
            }
        }
        const double ss = strain_squared(g);
        const double denominator = Kokkos::pow(ss, 2.5) + Kokkos::pow(sd, 1.25);
        return denominator > 0.0 ? Kokkos::pow(sd, 1.5) / denominator : 0.0;
    }
 
    KOKKOS_INLINE_FUNCTION
    static double vreman(const double g[3][3]) {
        double b[3][3], aa = 0.0;
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                // b_ij = sum_m a_mi a_mj = sum_m g_im g_jm
                b[i][j] = 0.0;
                for (int m = 0; m < 3; m++) b[i][j] += g[i][m] * g[j][m];
                aa += g[i][j] * g[i][j];
            }
        }
        const double B = b[0][0] * b[1][1] - b[0][1] * b[0][1] + b[0][0] * b[2][2] - b[0][2] * b[0][2] +
                         b[1][1] * b[2][2] - b[1][2] * b[1][2];
        return aa > 0.0 ? Kokkos::sqrt(Kokkos::fmax(B, 0.0) / aa) : 0.0;
    }
 
    KOKKOS_INLINE_FUNCTION
    static void singular_values(const double g[3][3], double s[3]) {
        double G[3][3];
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                G[i][j] = 0.0;
                for (int k = 0; k < 3; k++) G[i][j] += g[k][i] * g[k][j];
            }
        }
        const double I1 = G[0][0] + G[1][1] + G[2][2];
        const double I2 = G[0][0] * G[1][1] - G[0][1] * G[0][1] + G[0][0] * G[2][2] - G[0][2] * G[0][2] +
                          G[1][1] * G[2][2] - G[1][2] * G[1][2];
        const double det = g[0][0] * (g[1][1] * g[2][2] - g[1][2] * g[2][1]) -
                           g[0][1] * (g[1][0] * g[2][2] - g[1][2] * g[2][0]) +
                           g[0][2] * (g[1][0] * g[2][1] - g[1][1] * g[2][0]);
        const double I3 = det * det;
        double l1 = I1 / 3.0;
        const double a1 = I1 * I1 / 9.0 - I2 / 3.0;
        if (a1 > 0.0) {
            const double a2 = I1 * I1 * I1 / 27.0 - I1 * I2 / 6.0 + I3 / 2.0;
            const double a3 = Kokkos::acos(Kokkos::fmin(1.0, Kokkos::fmax(-1.0, a2 / Kokkos::pow(a1, 1.5)))) / 3.0;
            l1 += 2.0 * Kokkos::sqrt(a1) * Kokkos::cos(a3);
        }
        s[0] = s[1] = s[2] = 0.0;
        if (!(l1 > 0.0)) return;
        const double product = I3 / l1;
        const double sum = Kokkos::fmax((I2 - product) / l1, 0.0);
        const double l2 = 0.5 * (sum + Kokkos::sqrt(Kokkos::fmax(sum * sum - 4.0 * product, 0.0)));
        const double l3 = l2 > 0.0 ? product / l2 : 0.0;
        s[0] = Kokkos::sqrt(l1);
        s[1] = Kokkos::sqrt(Kokkos::fmin(l2, l1));
        s[2] = Kokkos::sqrt(Kokkos::fmin(l3, l2));
    }
 
    KOKKOS_INLINE_FUNCTION
    static double sigma(const double g[3][3]) {
        double s[3];
        singular_values(g, s);
        if (!(s[0] > 0.0)) return 0.0;
        return Kokkos::fmax(s[2] * (s[0] - s[1]) * (s[1] - s[2]) / (s[0] * s[0]), 0.0);
    }
};
 
#endif // LES_H