hyperbolic_system.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2023 - 2024 by the ryujin authors
//
 
#pragma once
 
#include "flux_library.h"
 
#include <convenience_macros.h>
#include <deal.II/base/config.h>
#include <discretization.h>
#include <multicomponent_vector.h>
#include <openmp.h>
#include <patterns_conversion.h>
#include <simd.h>
#include <state_vector.h>
 
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/tensor.h>
 
#include <array>
 
namespace ryujin
{
  namespace ScalarConservation
  {
    template <int dim, typename Number>
    class HyperbolicSystemView;
 
    class HyperbolicSystem final : public dealii::ParameterAcceptor
    {
    public:
      static inline std::string problem_name = "Scalar conservation equation";
 
      HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
 
      template <int dim, typename Number>
      auto view() const
      {
        return HyperbolicSystemView<dim, Number>{*this};
      }
 
    private:
      std::string flux_;
 
      FluxLibrary::flux_list_type flux_list_;
      using Flux = FluxLibrary::Flux;
      std::shared_ptr<Flux> selected_flux_;
 
      template <int dim, typename Number>
      friend class HyperbolicSystemView;
    }; /* HyperbolicSystem */
 
 
    template <int dim, typename Number>
    class HyperbolicSystemView
    {
    public:
      HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
          : hyperbolic_system_(hyperbolic_system)
      {
      }
 
      template <int dim2, typename Number2>
      auto view() const
      {
        return HyperbolicSystemView<dim2, Number2>{hyperbolic_system_};
      }
 
      using ScalarNumber = typename get_value_type<Number>::type;
 
 
      DEAL_II_ALWAYS_INLINE inline const std::string &flux() const
      {
        return hyperbolic_system_.flux_;
      }
 
      DEAL_II_ALWAYS_INLINE inline ScalarNumber
      derivative_approximation_delta() const
      {
        const auto &flux = hyperbolic_system_.selected_flux_;
        return ScalarNumber(flux->derivative_approximation_delta());
      }
 
 
 
      DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
      flux_function(const Number &u) const;
 
      DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
      flux_gradient_function(const Number &u) const;
 
 
 
    private:
      const HyperbolicSystem &hyperbolic_system_;
 
    public:
 
 
      static constexpr unsigned int problem_dimension = 1;
 
      using state_type = dealii::Tensor<1, problem_dimension, Number>;
 
      using flux_type =
          dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
 
      using flux_contribution_type = flux_type;
 
      static inline const auto component_names =
          std::array<std::string, problem_dimension>{"u"};
 
      static inline const auto primitive_component_names =
          std::array<std::string, problem_dimension>{"u"};
 
      static constexpr unsigned int n_precomputed_values = 2 * dim;
 
      using precomputed_type = std::array<Number, n_precomputed_values>;
 
      static inline const auto precomputed_names =
          []() -> std::array<std::string, n_precomputed_values> {
        if constexpr (dim == 1)
          return {"f", "df"};
        else if constexpr (dim == 2)
          return {"f_1", "f_2", "df_1", "df_2"};
        else if constexpr (dim == 3)
          return {"f_1", "f_2", "f_3", "df_1", "df_2", "df_3"};
        __builtin_trap();
      }();
 
      static constexpr unsigned int n_initial_precomputed_values = 0;
 
      using initial_precomputed_type =
          std::array<Number, n_initial_precomputed_values>;
 
      static inline const auto initial_precomputed_names =
          std::array<std::string, n_initial_precomputed_values>{};
 
      using StateVector = Vectors::
          StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
 
      using HyperbolicVector =
          Vectors::MultiComponentVector<ScalarNumber, problem_dimension>;
 
      using PrecomputedVector =
          Vectors::MultiComponentVector<ScalarNumber, n_precomputed_values>;
 
      using InitialPrecomputedVector =
          Vectors::MultiComponentVector<ScalarNumber,
                                        n_initial_precomputed_values>;
 
 
 
 
      static constexpr unsigned int n_precomputation_cycles = 1;
 
      template <typename DISPATCH, typename SPARSITY>
      void precomputation_loop(unsigned int cycle,
                               const DISPATCH &dispatch_check,
                               const SPARSITY &sparsity_simd,
                               StateVector &state_vector,
                               unsigned int left,
                               unsigned int right,
                               const bool skip_constrained_dofs = true) const;
 
 
 
      static Number state(const state_type &U);
 
      Number square_entropy(const Number &u) const;
 
      Number square_entropy_derivative(const Number &u) const;
 
      Number kruzkov_entropy(const Number &k, const Number &u) const;
 
      Number kruzkov_entropy_derivative(const Number &k, const Number &u) const;
 
      bool is_admissible(const state_type & /*U*/) const
      {
        return true;
      }
 
 
 
      template <typename Lambda>
      state_type
      apply_boundary_conditions(const dealii::types::boundary_id id,
                                const state_type &U,
                                const dealii::Tensor<1, dim, Number> &normal,
                                const Lambda &get_dirichlet_data) const;
 
 
 
      dealii::Tensor<1, dim, Number>
      construct_flux_tensor(const precomputed_type &precomputed_state) const;
 
      dealii::Tensor<1, dim, Number> construct_flux_gradient_tensor(
          const precomputed_type &precomputed_state) const;
 
      flux_contribution_type
      flux_contribution(const PrecomputedVector &pv,
                        const InitialPrecomputedVector & /*piv*/,
                        const unsigned int i,
                        const state_type & /*U_i*/) const;
 
      flux_contribution_type
      flux_contribution(const PrecomputedVector &pv,
                        const InitialPrecomputedVector & /*piv*/,
                        const unsigned int *js,
                        const state_type & /*U_j*/) const;
 
      state_type
      flux_divergence(const flux_contribution_type &flux_i,
                      const flux_contribution_type &flux_j,
                      const dealii::Tensor<1, dim, Number> &c_ij) const;
 
      static constexpr bool have_high_order_flux = false;
 
      state_type high_order_flux_divergence(
          const flux_contribution_type &,
          const flux_contribution_type &,
          const dealii::Tensor<1, dim, Number> &c_ij) const = delete;
 
 
 
      static constexpr bool have_source_terms = false;
 
      state_type nodal_source(const PrecomputedVector &pv,
                              const unsigned int i,
                              const state_type &U_i,
                              const ScalarNumber tau) const = delete;
 
      state_type nodal_source(const PrecomputedVector &pv,
                              const unsigned int *js,
                              const state_type &U_j,
                              const ScalarNumber tau) const = delete;
 
 
 
      template <typename ST>
      state_type expand_state(const ST &state) const
      {
        return state;
      }
 
      state_type from_primitive_state(const state_type &primitive_state) const
      {
        return primitive_state;
      }
 
      state_type to_primitive_state(const state_type &state) const
      {
        return state;
      }
 
      template <typename Lambda>
      state_type apply_galilei_transform(const state_type &state,
                                         const Lambda & /*lambda*/) const
      {
        return state;
      }
 
    }; /* HyperbolicSystemView */
 
 
    /*
     * -------------------------------------------------------------------------
     * Inline definitions
     * -------------------------------------------------------------------------
     */
 
 
    inline HyperbolicSystem::HyperbolicSystem(const std::string &subsection)
        : ParameterAcceptor(subsection)
    {
      flux_ = "burgers";
      add_parameter("flux",
                    flux_,
                    "The scalar flux. Valid names are given by any of the "
                    "subsections defined below");
 
      /*
       * And finally populate the flux list with all flux configurations
       * defined in the FluxLibrary namespace:
       */
      FluxLibrary::populate_flux_list(flux_list_, subsection);
 
      const auto populate_functions = [this]() {
        bool initialized = false;
        for (auto &it : flux_list_)
 
          /* Populate flux functions: */
          if (it->name() == flux_) {
            selected_flux_ = it;
            it->parse_parameters_call_back();
            problem_name = "Scalar conservation equation (" + it->name() +
                           ": " + it->flux_formula() + ")";
            initialized = true;
            break;
          }
 
        AssertThrow(initialized,
                    dealii::ExcMessage(
                        "Could not find a flux description with name \"" +
                        flux_ + "\""));
      };
 
      ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
      populate_functions();
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
    HyperbolicSystemView<dim, Number>::flux_function(const Number &u) const
    {
      const auto &flux = hyperbolic_system_.selected_flux_;
      dealii::Tensor<1, dim, Number> result;
 
      /* This access by calling into value() repeatedly is terrible: */
      for (unsigned int k = 0; k < dim; ++k) {
        if constexpr (std::is_same_v<ScalarNumber, Number>) {
          result[k] = flux->value(u, k);
        } else {
          for (unsigned int s = 0; s < Number::size(); ++s) {
            result[k][s] = flux->value(u[s], k);
          }
        }
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
    HyperbolicSystemView<dim, Number>::flux_gradient_function(
        const Number &u) const
    {
      const auto &flux = hyperbolic_system_.selected_flux_;
      dealii::Tensor<1, dim, Number> result;
 
      /* This access by calling into value() repeatedly is terrible: */
      for (unsigned int k = 0; k < dim; ++k) {
        if constexpr (std::is_same_v<ScalarNumber, Number>) {
          result[k] = flux->gradient(u, k);
        } else {
          for (unsigned int s = 0; s < Number::size(); ++s) {
            result[k][s] = flux->gradient(u[s], k);
          }
        }
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    template <typename DISPATCH, typename SPARSITY>
    DEAL_II_ALWAYS_INLINE inline void
    HyperbolicSystemView<dim, Number>::precomputation_loop(
        unsigned int cycle [[maybe_unused]],
        const DISPATCH &dispatch_check,
        const SPARSITY &sparsity_simd,
        StateVector &state_vector,
        unsigned int left,
        unsigned int right,
        const bool skip_constrained_dofs /*= true*/) const
    {
      Assert(cycle == 0, dealii::ExcInternalError());
 
      const auto &U = std::get<0>(state_vector);
      auto &precomputed = std::get<1>(state_vector);
 
      /* We are inside a thread parallel context */
 
      unsigned int stride_size = get_stride_size<Number>;
 
      RYUJIN_OMP_FOR
      for (unsigned int i = left; i < right; i += stride_size) {
 
        /* Skip constrained degrees of freedom: */
        const unsigned int row_length = sparsity_simd.row_length(i);
        if (skip_constrained_dofs && row_length == 1)
          continue;
 
        dispatch_check(i);
 
        const auto U_i = U.template get_tensor<Number>(i);
        const auto u_i = state(U_i);
 
        const auto f_i = flux_function(u_i);
        const auto df_i = flux_gradient_function(u_i);
 
        precomputed_type prec_i;
 
        for (unsigned int k = 0; k < n_precomputed_values / 2; ++k) {
          prec_i[k] = f_i[k];
          prec_i[dim + k] = df_i[k];
        }
 
        precomputed.template write_tensor<Number>(prec_i, i);
      }
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::state(const state_type &U)
    {
      return U[0];
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::square_entropy(const Number &u) const
    {
      return ScalarNumber(0.5) * u * u;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::square_entropy_derivative(
        const Number &u) const
    {
      return u;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::kruzkov_entropy(const Number &k,
                                                       const Number &u) const
    {
      return std::abs(k - u);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::kruzkov_entropy_derivative(
        const Number &k, const Number &u) const
    {
      constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
      // return sgn(u-k):
      return dealii::compare_and_apply_mask<gte>(u, k, Number(1.), Number(-1.));
    }
 
 
    template <int dim, typename Number>
    template <typename Lambda>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::apply_boundary_conditions(
        dealii::types::boundary_id id,
        const state_type &U,
        const dealii::Tensor<1, dim, Number> & /*normal*/,
        const Lambda &get_dirichlet_data) const -> state_type
    {
      state_type result = U;
 
      if (id == Boundary::dirichlet) {
        result = get_dirichlet_data();
 
      } else if (id == Boundary::dirichlet_momentum) {
        AssertThrow(false,
                    dealii::ExcMessage(
                        "Invalid boundary ID »Boundary::dirichlet_momentum«, "
                        "enforcing Dirichlet boundary conditions on a momentum "
                        "is not possible for scalar conservation equations."));
 
      } else if (id == Boundary::slip) {
        AssertThrow(
            false,
            dealii::ExcMessage("Invalid boundary ID »Boundary::slip«, slip "
                               "boundary conditions are unavailable for scalar "
                               "conservation equations."));
        __builtin_trap();
 
      } else if (id == Boundary::no_slip) {
        AssertThrow(
            false,
            dealii::ExcMessage("Invalid boundary ID »Boundary::no_slip«, "
                               "no-slip boundary conditions are unavailable "
                               "for scalar conservation equations."));
        __builtin_trap();
 
      } else if (id == Boundary::dynamic) {
        AssertThrow(
            false,
            dealii::ExcMessage("Invalid boundary ID »Boundary::dynamic«, "
                               "dynamic boundary conditions are unavailable "
                               "for scalar conservation equations."));
        __builtin_trap();
 
      } else {
        AssertThrow(false, dealii::ExcNotImplemented());
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
    HyperbolicSystemView<dim, Number>::construct_flux_tensor(
        const precomputed_type &precomputed) const
    {
      dealii::Tensor<1, dim, Number> result;
 
      if constexpr (dim == 1) {
        const auto &[f, df] = precomputed;
        result[0] = f;
 
      } else if constexpr (dim == 2) {
        const auto &[f_1, f_2, df_1, df_2] = precomputed;
        result[0] = f_1;
        result[1] = f_2;
 
      } else if constexpr (dim == 3) {
        const auto &[f_1, f_2, f_3, df_1, df_2, df_3] = precomputed;
        result[0] = f_1;
        result[1] = f_2;
        result[2] = f_3;
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
    HyperbolicSystemView<dim, Number>::construct_flux_gradient_tensor(
        const precomputed_type &precomputed) const
    {
      dealii::Tensor<1, dim, Number> result;
 
      if constexpr (dim == 1) {
        const auto &[f, df] = precomputed;
        result[0] = df;
 
      } else if constexpr (dim == 2) {
        const auto &[f_1, f_2, df_1, df_2] = precomputed;
        result[0] = df_1;
        result[1] = df_2;
 
      } else if constexpr (dim == 3) {
        const auto &[f_1, f_2, f_3, df_1, df_2, df_3] = precomputed;
        result[0] = df_1;
        result[1] = df_2;
        result[2] = df_3;
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::flux_contribution(
        const PrecomputedVector &pv,
        const InitialPrecomputedVector & /*piv*/,
        const unsigned int i,
        const state_type & /*U_i*/) const -> flux_contribution_type
    {
      /* The flux contribution is a rank 2 tensor, thus a little dance: */
      flux_contribution_type result;
      result[0] = construct_flux_tensor(
          pv.template get_tensor<Number, precomputed_type>(i));
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::flux_contribution(
        const PrecomputedVector &pv,
        const InitialPrecomputedVector & /*piv*/,
        const unsigned int *js,
        const state_type & /*U_j*/) const -> flux_contribution_type
    {
      /* The flux contribution is a rank 2 tensor, thus a little dance: */
      flux_contribution_type result;
      result[0] = construct_flux_tensor(
          pv.template get_tensor<Number, precomputed_type>(js));
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::flux_divergence(
        const flux_contribution_type &flux_i,
        const flux_contribution_type &flux_j,
        const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type
    {
      return -contract(add(flux_i, flux_j), c_ij);
    }
 
  } // namespace ScalarConservation
} // namespace ryujin