physics.h

Go to the documentation of this file.

 
#ifndef PHYSICS_H
#define PHYSICS_H
 
#include <string>
#include <unordered_map>
 
#include <Kokkos_Core.hpp>
#include <toml.hpp>
 
#include "common.h"
#include "log.h"
 
enum class PhysicsType {
    EULER,
    NAVIER_STOKES
};
 
static const std::unordered_map<std::string, PhysicsType> PHYSICS_TYPES = {
    {"euler", PhysicsType::EULER},
    {"navier_stokes", PhysicsType::NAVIER_STOKES}
};
 
static const std::unordered_map<PhysicsType, std::string> PHYSICS_NAMES = {
    {PhysicsType::EULER, "euler"},
    {PhysicsType::NAVIER_STOKES, "navier_stokes"}
};
 
enum class ViscosityModel {
    NONE,
    CONSTANT,
    SUTHERLAND,
    POWER_LAW
};
 
struct Euler {
    rtype gamma = 1.4_r;
    rtype R = 287.0;
    rtype cp = 1004.5;
    rtype cv = 717.5;
    ViscosityModel viscosity_model = ViscosityModel::NONE;
    rtype mu_ref = 0.0;      // Viscosity (constant) or at T_mu_ref (Sutherland, power law)
    rtype T_mu_ref = 273.15_r;
    rtype S_mu = 110.4_r;     // Sutherland temperature
    rtype n_mu = 0.75_r;      // Power-law exponent
    rtype Pr = 0.72_r;
 
    static Euler from_reference(rtype gamma, rtype p_ref, rtype T_ref, rtype rho_ref);
 
    static Euler from_input(const toml::value & input);
 
    logging::Items summary() const;
 
    KOKKOS_INLINE_FUNCTION
    PhysicsType get_type() const {
        return viscosity_model == ViscosityModel::NONE ? PhysicsType::EULER : PhysicsType::NAVIER_STOKES;
    }
 
    KOKKOS_INLINE_FUNCTION
    bool is_viscous() const { return viscosity_model != ViscosityModel::NONE; }
 
    KOKKOS_INLINE_FUNCTION
    rtype viscosity(const rtype T) const {
        switch (viscosity_model) {
            case ViscosityModel::CONSTANT:
                return mu_ref;
            case ViscosityModel::SUTHERLAND:
                return mu_ref * Kokkos::pow(T / T_mu_ref, 1.5_r) * (T_mu_ref + S_mu) / (T + S_mu);
            case ViscosityModel::POWER_LAW:
                return mu_ref * Kokkos::pow(T / T_mu_ref, n_mu);
            default:
                return 0.0;
        }
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype conductivity(const rtype mu) const { return mu * cp / Pr; }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_gamma() const { return gamma; }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_energy_from_temperature(const rtype T) const { return cv * T; }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_temperature_from_energy(const rtype e) const { return e / cv; }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_density_from_pressure_temperature(const rtype p, const rtype T) const {
        return p / (R * T);
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_temperature_from_density_pressure(const rtype rho, const rtype p) const {
        return p / (rho * R);
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_pressure_from_density_temperature(const rtype rho, const rtype T) const {
        return rho * R * T;
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_pressure_from_density_energy(const rtype rho, const rtype e) const {
        return (gamma - 1.0_r) * rho * e;
    }
 
    KOKKOS_INLINE_FUNCTION
    rtype get_sound_speed_from_pressure_density(const rtype p, const rtype rho) const {
        return Kokkos::sqrt(gamma * p / rho);
    }
 
    KOKKOS_INLINE_FUNCTION
    void compute_primitives_from_conservatives(rtype * primitives,
                                               const rtype * conservatives) const {
        const rtype rho = conservatives[0];
        rtype u[N_DIM];
        FOR_I_DIM u[i] = conservatives[1 + i] / rho;
        const rtype e = conservatives[N_DIM + 1] / rho - 0.5_r * dot<N_DIM>(u, u);
        const rtype p = get_pressure_from_density_energy(rho, e);
        FOR_I_DIM primitives[i] = u[i];
        primitives[N_DIM] = p;
        primitives[N_DIM + 1] = get_temperature_from_energy(e);
        primitives[N_DIM + 2] = e + p / rho;
    }
 
    KOKKOS_INLINE_FUNCTION
    void compute_W_from_conservatives(rtype * W, const rtype * conservatives) const {
        const rtype rho = conservatives[0];
        W[0] = rho;
        FOR_I_DIM W[1 + i] = conservatives[1 + i] / rho;
        W[N_DIM + 1] = (gamma - 1.0_r) * (conservatives[N_DIM + 1] - 0.5_r * rho * dot<N_DIM>(W + 1, W + 1));
    }
 
    KOKKOS_INLINE_FUNCTION
    void compute_conservatives_from_W(rtype * conservatives, const rtype * W) const {
        conservatives[0] = W[0];
        FOR_I_DIM conservatives[1 + i] = W[0] * W[1 + i];
        conservatives[N_DIM + 1] = W[N_DIM + 1] / (gamma - 1.0_r) + 0.5_r * W[0] * dot<N_DIM>(W + 1, W + 1);
    }
};
 
#endif // PHYSICS_H