characteristic.h

Go to the documentation of this file.

 
#ifndef CHARACTERISTIC_H
#define CHARACTERISTIC_H
 
#include <Kokkos_Core.hpp>
 
#include "boundary.h"
#include "common.h"
#include "gradient.h"
 
struct CharacteristicTransverseFunctor {
    BoundaryData boundaries;
    Kokkos::View<rtype *[N_DIM]> face_coords;
    Kokkos::View<rtype *[N_DIM]> face_normals;
    Kokkos::View<rtype *> face_area;
    Kokkos::View<int32_t *[2]> cells_of_face;
    Kokkos::View<rtype *[N_CONSERVATIVE]> W_cells;
 
    KOKKOS_INLINE_FUNCTION
    void face_state(const uint32_t f, rtype * W) const {
        FOR_I_CONSERVATIVE W[i] = W_cells(cells_of_face(f, 0), i);
    }
 
    KOKKOS_INLINE_FUNCTION
    void operator()(const uint32_t k) const {
        constexpr uint8_t E = N_DIM + 1;
        const uint32_t f = boundaries.char_faces(k);
        rtype n[N_DIM], n_vec[N_DIM];
        FOR_I_DIM n_vec[i] = face_normals(f, i);
        unit<N_DIM>(n_vec, n);
        rtype W_f[N_CONSERVATIVE];
        face_state(f, W_f);
 
        // Neighbors in the plane, and a point along the normal that sets the normal derivatives to zero
        LSQSystem lsq;
        const rtype ell = (N_DIM == 2) ? face_area(f) : Kokkos::sqrt(face_area(f));
        rtype dx_n[N_DIM];
        FOR_I_DIM dx_n[i] = ell * n[i];
        lsq.add(dx_n, W_f, W_f);
        const uint32_t n_neighbors = boundaries.char_offsets(k + 1) - boundaries.char_offsets(k);
        if (n_neighbors + 1 < N_DIM || boundaries.char_edge(k)) {
            for (uint8_t m = 0; m < 3; m++) boundaries.char_transverse(k, m) = 0.0_r;
            return;
        }
        for (uint32_t j = boundaries.char_offsets(k); j < boundaries.char_offsets(k + 1); j++) {
            const uint32_t g = boundaries.char_faces(boundaries.char_neighbors(j));
            rtype dx[N_DIM], W_g[N_CONSERVATIVE];
            FOR_I_DIM dx[i] = boundaries.char_neighbor_dx(j, i);
            const rtype dx_n_g = dot<N_DIM>(dx, n);
            FOR_I_DIM dx[i] -= dx_n_g * n[i];
            face_state(g, W_g);
            lsq.add(dx, W_f, W_g);
        }
        rtype g[N_CONSERVATIVE][N_DIM];
        lsq.solve(g);
 
        const rtype u_n = dot<N_DIM>(W_f + 1, n);
        rtype u_t[N_DIM];
        FOR_I_DIM u_t[i] = W_f[1 + i] - u_n * n[i];
        rtype div_tangential = 0.0_r, u_t_grad_u_n = 0.0_r;
        FOR_I_DIM {
            div_tangential += g[1 + i][i];
            for (uint8_t j = 0; j < N_DIM; j++) {
                div_tangential -= n[i] * g[1 + i][j] * n[j];
                u_t_grad_u_n += n[i] * g[1 + i][j] * u_t[j];
            }
        }
        boundaries.char_transverse(k, 0) = dot<N_DIM>(u_t, g[E]);
        boundaries.char_transverse(k, 1) = W_f[0] * div_tangential;
        boundaries.char_transverse(k, 2) = W_f[0] * u_t_grad_u_n;
    }
};
 
struct CharacteristicStateFunctor {
    BoundaryData boundaries;
    Kokkos::View<rtype *[N_DIM]> face_normals;
    Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution;
    Kokkos::View<rtype **[2][2]> face_thermo;  // Gas mixtures: [gamma, e0]; empty for a single gas
    Kokkos::View<rtype *> quad_weights;        // 2D
    Kokkos::View<rtype **> face_weights;       // 3D: (face, q)
    rtype dt;
    bool initialize;
 
    KOKKOS_INLINE_FUNCTION
    void operator()(const uint32_t k) const {
        using Relax = BoundaryCondition::Relax;
        constexpr uint8_t E = N_DIM + 1;
        const uint32_t f = boundaries.char_faces(k);
        const BoundaryCondition & bc = boundaries.bcs(boundaries.face_bc(f));
        uint8_t n_quad;
        if constexpr (N_DIM == 2) {
            n_quad = static_cast<uint8_t>(quad_weights.extent(0));
        } else {
            n_quad = static_cast<uint8_t>(face_weights.extent(1));
        }
        rtype W[N_CONSERVATIVE] = {};
        rtype gamma = 0.0_r, w_sum = 0.0_r;
        for (uint8_t q = 0; q < n_quad; q++) {
            rtype w;
            if constexpr (N_DIM == 2) {
                w = quad_weights(q);
            } else {
                w = face_weights(f, q);
            }
            if (w == 0.0_r) continue;
            w_sum += w;
            FOR_I_CONSERVATIVE W[i] += w * face_solution(f, q, 0, i);
            gamma += w * (face_thermo.extent(0) > 0 ? face_thermo(f, q, 0, 0) : boundaries.gamma);
        }
        FOR_I_CONSERVATIVE W[i] /= w_sum;
        gamma /= w_sum;
        rtype n[N_DIM], n_vec[N_DIM];
        FOR_I_DIM n_vec[i] = face_normals(f, i);
        unit<N_DIM>(n_vec, n);
 
        const rtype c2 = gamma * W[E] / W[0];
        const rtype c = Kokkos::sqrt(c2);
        const rtype Z = W[0] * c;
        const rtype u_n = dot<N_DIM>(W + 1, n);
        const rtype w_out = W[E] + Z * u_n;
        if (initialize || u_n >= c) {
            boundaries.char_state(k, 0) = W[E];
            boundaries.char_state(k, 1) = u_n;
            if (u_n >= c) return;
        }
        const bool inlet = bc.type == BoundaryType::NSCBC_INLET;
        rtype u_target[N_DIM], du_target[N_DIM];
        FOR_I_DIM {
            u_target[i] = boundaries.char_target(k, i);
            du_target[i] = boundaries.char_target_next(k, i) - u_target[i];
        }
        if (inlet && u_n <= -c) {
            boundaries.char_state(k, 0) = bc.data[E];
            boundaries.char_state(k, 1) = dot<N_DIM>(u_target, n);
            return;
        }
        rtype w_in = boundaries.char_state(k, 0) - Z * boundaries.char_state(k, 1);
        const rtype p_b = 0.5_r * (w_out + w_in);
        const rtype u_b = 0.5_r * (w_out - w_in) / Z;
        const rtype M2 = dot<N_DIM>(W + 1, W + 1) / c2;
        const rtype K = Kokkos::fmin(bc.relax[Relax::ACOUSTIC] * c * Kokkos::fmax(1.0_r - M2, 0.0_r), 1.0_r / dt);
        const rtype relaxation = inlet ? K * Z * (u_b - dot<N_DIM>(u_target, n)) : -K * (p_b - bc.data[E]);
        // None on reversed faces of outlets, which nothing leaves through
        const rtype beta = (!inlet && u_n < 0.0_r)             ? 0.0_r
                           : (bc.relax[Relax::BETA] < 0.0_r) ? Kokkos::fmin(Kokkos::sqrt(M2), 1.0_r)
                                                             : bc.relax[Relax::BETA];
        const rtype T_in = boundaries.char_transverse(k, 0) + c2 * boundaries.char_transverse(k, 1) -
                           c * boundaries.char_transverse(k, 2);
        w_in += dt * (relaxation - beta * T_in);
        // A target that varies in time carries its own incoming wave, -Z du_n,t (Guezennec & Poinsot 2009)
        if (inlet) w_in -= Z * dot<N_DIM>(du_target, n);
        boundaries.char_state(k, 0) = 0.5_r * (w_out + w_in);
        boundaries.char_state(k, 1) = 0.5_r * (w_out - w_in) / Z;
    }
};
 
#endif // CHARACTERISTIC_H