face_reconstruction.h

Go to the documentation of this file.

 
#ifndef FACE_RECONSTRUCTION_H
#define FACE_RECONSTRUCTION_H
 
#include <array>
#include <memory>
#include <unordered_map>
#include <vector>
 
#include <toml.hpp>
 
#include "common_typedef.h"
#include "mesh.h"
#include "quadrature.h"
#include "boundary.h"
#include "gradient.h"
#include "log.h"
#include "teno.h"
 
enum class FaceReconstructionType {
    FIRST_ORDER,
    MUSCL,
    TENO,
};
 
static const std::unordered_map<std::string, FaceReconstructionType> FACE_RECONSTRUCTION_TYPES = {
    {"FO", FaceReconstructionType::FIRST_ORDER},
    {"MUSCL", FaceReconstructionType::MUSCL},
    {"TENO", FaceReconstructionType::TENO},
};
 
static const std::unordered_map<FaceReconstructionType, std::string> FACE_RECONSTRUCTION_NAMES = {
    {FaceReconstructionType::FIRST_ORDER, "FO"},
    {FaceReconstructionType::MUSCL, "MUSCL"},
    {FaceReconstructionType::TENO, "TENO"},
};
 
class GrowingBuffer;
 
class FaceReconstruction {
    public:
        FaceReconstruction();
 
        virtual ~FaceReconstruction();
 
        virtual void init(const toml::value & input) = 0;
 
        virtual logging::Items summary() const;
 
        void set_mesh(std::shared_ptr<Mesh> mesh);
 
        void set_boundaries(const BoundaryData & boundaries);
 
        FaceReconstructionType get_type() const { return type; }
 
        virtual uint8_t n_face_quadrature_points() const = 0;
 
        virtual void calc_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                      Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) = 0;
 
        virtual std::vector<uint32_t> cells_independent_of_halo(uint32_t n_owned) const;
 
        virtual void calc_cell_face_values(const Kokkos::DefaultExecutionSpace & exec,
                                           Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                           Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution,
                                           Kokkos::View<uint32_t *> cells);
 
        virtual void finish_cell_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                             Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution);
 
        virtual bool cell_gradients(Kokkos::View<rtype *[N_CONSERVATIVE]> /*solution*/,
                                    Kokkos::View<rtype *[N_CONSERVATIVE][N_DIM]> /*gradients*/,
                                    uint32_t /*n_cells*/) {
            return false;
        }
        
        virtual void axisymmetric_source(Kokkos::View<rtype *[N_CONSERVATIVE]> /*solution*/,
                                         Kokkos::View<rtype *, Kokkos::LayoutStride> /*mu*/,
                                         Kokkos::View<rtype *> /*source*/, uint32_t /*n_cells*/) {}
 
        void init_face_quadrature_3d(uint8_t degree);
 
        void init_face_quadrature_2d_curved();
 
        Quadrature quadrature_face;
        // 3D, and 2D curved meshes: (face, q, dim) and (face, q). Weights sum to 2 per face
        // (twice the area over face_area on curved faces) and are zero on padding points
        Kokkos::View<rtype ***> face_quad_points;
        Kokkos::View<rtype **> face_quad_weights;
        Kokkos::View<rtype ***> face_quad_normals;  // curved meshes: (face, q, dim) unit normals out of cell 0
    protected:
        FaceReconstructionType type;
        std::shared_ptr<Mesh> mesh;
        BoundaryData boundaries;
    private:
};
 
class FirstOrder : public FaceReconstruction {
    public:
        FirstOrder();
 
        ~FirstOrder();
 
        void init(const toml::value & input) override;
 
        uint8_t n_face_quadrature_points() const override;
 
        void calc_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                              Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) override;
    protected:
    private:
};
 
enum class LimiterType {
    NONE,
    BARTH_JESPERSEN,
    VENKATAKRISHNAN,
};
 
static const std::unordered_map<std::string, LimiterType> LIMITER_TYPES = {
    {"none", LimiterType::NONE},
    {"barth_jespersen", LimiterType::BARTH_JESPERSEN},
    {"venkatakrishnan", LimiterType::VENKATAKRISHNAN},
};
 
static const std::unordered_map<LimiterType, std::string> LIMITER_NAMES = {
    {LimiterType::NONE, "none"},
    {LimiterType::BARTH_JESPERSEN, "barth_jespersen"},
    {LimiterType::VENKATAKRISHNAN, "venkatakrishnan"},
};
 
class MUSCL : public FaceReconstruction {
    public:
        MUSCL();
        ~MUSCL();
        void init(const toml::value & input) override;
        logging::Items summary() const override;
        uint8_t n_face_quadrature_points() const override;
        void calc_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                              Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) override;
        std::vector<uint32_t> cells_independent_of_halo(uint32_t n_owned) const override;
        void calc_cell_face_values(const Kokkos::DefaultExecutionSpace & exec,
                                   Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                   Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution,
                                   Kokkos::View<uint32_t *> cells) override;
        void finish_cell_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                     Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) override;
 
        LimiterType limiter = LimiterType::VENKATAKRISHNAN;
        rtype venkat_K = 5.0;
        Kokkos::View<rtype *[N_CONSERVATIVE][N_DIM]> gradients;
        Kokkos::View<rtype *[N_CONSERVATIVE]> limiters;
        LSQVertexGradientFunctor gradient;  // Its linear vertex fit serves tetrahedra
};
 
class TENO : public FaceReconstruction {
    public:
        TENO();
        ~TENO();
        void init(const toml::value & input) override;
        logging::Items summary() const override;
        uint8_t n_face_quadrature_points() const override;
        void calc_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                              Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) override;
        std::vector<uint32_t> cells_independent_of_halo(uint32_t n_owned) const override;
        void calc_cell_face_values(const Kokkos::DefaultExecutionSpace & exec,
                                   Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                   Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution,
                                   Kokkos::View<uint32_t *> cells) override;
        void finish_cell_face_values(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                     Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution) override;
        bool cell_gradients(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                            Kokkos::View<rtype *[N_CONSERVATIVE][N_DIM]> gradients, uint32_t n_cells) override;
        void axisymmetric_source(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                 Kokkos::View<rtype *, Kokkos::LayoutStride> mu, Kokkos::View<rtype *> source,
                                 uint32_t n_cells) override;
 
        uint8_t degree = 4;
        uint8_t n_dof_large = 0;
        uint16_t n_stencil_large = 0;
        uint16_t n_stencil_small = 10;
        rtype stencil_factor = 2.0;
        rtype sigma_threshold = 1.0e-3_r;
        rtype sigma_upper = 1.0e-2_r;
        rtype C_T = -1.0;  // Fixed cutoff; negative selects the adaptive cutoff
        bool characteristic = true;
        bool bound_preserving = false;
        // Stencils take mirror images across curved wall and symmetry faces too
        // (across their planes); by default they stay one-sided there
        bool curved_mirrors = false;
        // Largest degree of the central polynomials of cells whose stencils reach a
        // curved wall without mirror images: one-sided fits of higher degree grow there
        uint8_t curved_wall_degree = 3;
        rtype max_condition = 1.0e8;
        uint8_t slice_shift = teno::SLICE_SHIFT;  // log2 of the cells per slice of the packed stencils
        // Pseudo-inverses stored in single precision (opt-in: halves the
        // largest tables, and changes results at single-precision round-off)
        bool single_tables = false;
 
        // Per-cell precomputed data, for the reconstructed cells [0, n_reconstructed)
        Kokkos::View<rtype *> scale;                       // h = sqrt(V)
        Kokkos::View<rtype **> basis_mean;                 // (cell, l): mean of phi_l over the cell
        Kokkos::View<uint16_t *> stencil_large_size;       // (cell)
        teno::PackedStencils stencil_large;                // pseudo-inverse width: nk
        Kokkos::View<uint16_t **> stencil_small_size;      // (cell, face); 0 if stencil invalid
        teno::PackedStencils stencil_small;                // the faces' stencils one after another; width NK_SMALL
        Kokkos::View<rtype **> si_matrix;                  // 2D: (cell, upper_index(l, m)): symmetric
        teno::Moments moments;                             // 3D: (cell, k), from which the troubled passes form it
        teno::SmoothnessTerms si_terms;                    // 3D
        Kokkos::View<rtype *> troubled;                    // (local cell): sigma, for diagnostics
        Kokkos::View<rtype ***> troubled_coeffs;           // (queue slot, l, var): scratch for the troubled passes
        Kokkos::View<rtype ****> troubled_small_coeffs;    // (queue slot, sector, l, var): scratch for the troubled passes
        // Queued cells per round of the troubled passes, as a fraction of the
        // reconstructed cells (scratch slots); later rounds recompute the
        // central coefficients the smooth pass could not keep
        rtype troubled_capacity = N_DIM == 3 ? 0.5 : 1.0;
        Kokkos::View<uint32_t *> troubled_cells;           // queue of troubled cells
        Kokkos::View<uint32_t> n_troubled;
        // Vertex-neighbor layers each reconstructed cell's stencil search visited
        // (host); a distributed run needs this many complete layers around the cell
        std::vector<uint8_t> gather_depth;
        // Axisymmetric runs: rule on the reference triangle (q, [xi, eta, w]) of the geometric source
        Kokkos::View<rtype *[3]> cell_rule;
        // Device memory of the packed stencils when they were streamed to the device (CUDA)
        std::vector<std::shared_ptr<GrowingBuffer>> table_memory;
 
        struct Stencils {
            std::vector<uint64_t> offsets;  // (cell + 1)
            std::vector<int32_t> cells;
            std::vector<int32_t> faces;
        };
        Stencils large_stencils() const;
 
        std::vector<uint8_t> stencil_cells() const;
 
        std::array<double, 9> stencil_metric(uint32_t i) const;
 
        static uint8_t cached_halo_layers(const toml::value & input);
 
        void save_cache(uint8_t halo_layers = 0);
 
        void rebuild_cells(const std::vector<uint32_t> & cells);
 
        std::vector<uint32_t> cells_within_reach(const std::vector<uint32_t> & changed) const;
 
        // 3D setup: on the host only (tests), and the cells per device batch (0: from the device memory)
        bool setup_on_host = false;
        uint32_t setup_batch_cells = 0;
 
        // Gas mixtures: the flow block is reconstructed in primitive variables
        // W with characteristic projections on the primitive system, whose
        // sound speed comes from each cell's frozen gamma; troubled cells are
        // also found by jumps of the mixture molar mass; and the stencil
        // choice of the entropy field on each face side (SELECT_LARGE or a
        // bitmask of the kept sector stencils) is left for the species
        static constexpr uint8_t SELECT_LARGE = 0xFF;
        bool primitive = false;
        Kokkos::View<rtype *, Kokkos::LayoutStride> cell_gamma;  // (cell)
        Kokkos::View<rtype *> cell_molar_mass;                   // (cell)
        Kokkos::View<uint8_t **> selection;                      // (face, side)
 
        void set_mixture(Kokkos::View<rtype *, Kokkos::LayoutStride> gamma, Kokkos::View<rtype *> molar_mass);
 
    private:
        template <uint8_t DEG, bool PRIM, typename TABLE>
        void launch_source(Kokkos::View<rtype *[N_CONSERVATIVE]> solution, Kokkos::View<rtype *, Kokkos::LayoutStride> mu,
                           Kokkos::View<rtype *> source, uint32_t n_cells);
        template <uint8_t DEG, bool PRIM, typename TABLE>
        void launch_gradients(Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                              Kokkos::View<rtype *[N_CONSERVATIVE][N_DIM]> gradients, uint32_t n_cells);
        template <uint8_t DEG, bool PRIM, typename TABLE>
        void launch_reconstruction(const Kokkos::DefaultExecutionSpace & exec,
                                   Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                                   Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution,
                                   Kokkos::View<uint32_t *> cells, bool troubled_pass);
        void dispatch(const Kokkos::DefaultExecutionSpace & exec, Kokkos::View<rtype *[N_CONSERVATIVE]> solution,
                      Kokkos::View<rtype **[2][N_CONSERVATIVE]> face_solution, Kokkos::View<uint32_t *> cells,
                      bool troubled_pass);
 
        void read_options(const toml::value & input);
        void compute_stencils_and_matrices(const std::vector<uint32_t> * subset = nullptr);
        void compute_stencils_and_matrices_3d(const std::vector<uint32_t> * subset);
        void allocate_scratch();
        uint64_t options_key() const;
        uint64_t cache_key() const;
        bool load_cache();
        void round_to_single();
 
        std::string cache_file;  // this rank's
        bool cache_single = false;  // pseudo-inverses and smoothness-indicator matrices (3D: moments) cached in single precision
        bool cache_loaded = false;
 
        uint32_t largest_stencil = 0;        // Largest central stencil on any rank (cells)
        int64_t n_sector_unavailable = -1;   // Small sector stencils cut by boundaries, -1 if unknown
        std::string cache_status;          // Stencil cache outcome, empty without a cache file
        std::string setup_status;          // Where and how fast the tables were set up
        std::string setup_producer = "host";  // Execution space that set up the tables, or "host"
};
 
#endif // FACE_RECONSTRUCTION_H