limiter.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2020 - 2024 by the ryujin authors
//
 
#pragma once
 
#include <compile_time_options.h>
 
#include "hyperbolic_system.h"
 
#include <multicomponent_vector.h>
#include <newton.h>
#include <simd.h>
 
namespace ryujin
{
  namespace EulerAEOS
  {
    template <typename ScalarNumber = double>
    class LimiterParameters : public dealii::ParameterAcceptor
    {
    public:
      LimiterParameters(const std::string &subsection = "/Limiter")
          : ParameterAcceptor(subsection)
      {
        iterations_ = 2;
        add_parameter(
            "iterations", iterations_, "Number of limiter iterations");
 
        if constexpr (std::is_same_v<ScalarNumber, double>)
          newton_tolerance_ = 1.e-10;
        else
          newton_tolerance_ = 1.e-4;
        add_parameter("newton tolerance",
                      newton_tolerance_,
                      "Tolerance for the quadratic newton stopping criterion");
 
        newton_max_iterations_ = 2;
        add_parameter("newton max iterations",
                      newton_max_iterations_,
                      "Maximal number of quadratic newton iterations performed "
                      "during limiting");
 
        relaxation_factor_ = ScalarNumber(1.);
        add_parameter("relaxation factor",
                      relaxation_factor_,
                      "Factor for scaling the relaxation window with r_i = "
                      "factor * (m_i/|Omega|)^(1.5/d).");
      }
 
      ACCESSOR_READ_ONLY(iterations);
      ACCESSOR_READ_ONLY(newton_tolerance);
      ACCESSOR_READ_ONLY(newton_max_iterations);
      ACCESSOR_READ_ONLY(relaxation_factor);
 
    private:
      unsigned int iterations_;
      ScalarNumber newton_tolerance_;
      unsigned int newton_max_iterations_;
      ScalarNumber relaxation_factor_;
    };
 
 
    template <int dim, typename Number = double>
    class Limiter
    {
    public:
 
      using View = HyperbolicSystemView<dim, Number>;
 
      using ScalarNumber = typename View::ScalarNumber;
 
      static constexpr auto problem_dimension = View::problem_dimension;
 
      using state_type = typename View::state_type;
 
      using flux_contribution_type = typename View::flux_contribution_type;
 
      using precomputed_type = typename View::precomputed_type;
 
      using PrecomputedVector = typename View::PrecomputedVector;
 
      using Parameters = LimiterParameters<ScalarNumber>;
 
 
      //
 
      static constexpr unsigned int n_bounds = 4;
 
      using Bounds = std::array<Number, n_bounds>;
 
      Limiter(const HyperbolicSystem &hyperbolic_system,
              const Parameters &parameters,
              const PrecomputedVector &precomputed_values)
          : hyperbolic_system(hyperbolic_system)
          , parameters(parameters)
          , precomputed_values(precomputed_values)
      {
      }
 
      Bounds projection_bounds_from_state(const unsigned int i,
                                          const state_type &U_i) const;
 
      Bounds combine_bounds(const Bounds &bounds_left,
                            const Bounds &bounds_right) const;
 
      Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const;
 
 
 
      void reset(const unsigned int i,
                 const state_type &U_i,
                 const flux_contribution_type &flux_i);
 
      void accumulate(const unsigned int *js,
                      const state_type &U_j,
                      const flux_contribution_type &flux_j,
                      const dealii::Tensor<1, dim, Number> &scaled_c_ij,
                      const state_type &affine_shift);
 
      Bounds bounds(const Number hd_i) const;
 
      //*}
 
      std::tuple<Number, bool> limit(const Bounds &bounds,
                                     const state_type &U,
                                     const state_type &P,
                                     const Number t_min = Number(0.),
                                     const Number t_max = Number(1.)) const;
 
    private:
 
 
      const HyperbolicSystem &hyperbolic_system;
      const Parameters &parameters;
      const PrecomputedVector &precomputed_values;
 
      state_type U_i;
      flux_contribution_type flux_i;
 
      Bounds bounds_;
 
      Number rho_relaxation_numerator;
      Number rho_relaxation_denominator;
      Number s_interp_max;
 
    };
 
 
    /*
     * -------------------------------------------------------------------------
     * Inline definitions
     * -------------------------------------------------------------------------
     */
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    Limiter<dim, Number>::projection_bounds_from_state(
        const unsigned int i, const state_type &U_i) const -> Bounds
    {
      const auto view = hyperbolic_system.view<dim, Number>();
      const auto rho_i = view.density(U_i);
      const auto &[p_i, gamma_min_i, s_i, eta_i] =
          precomputed_values.template get_tensor<Number, precomputed_type>(i);
 
      return {/*rho_min*/ rho_i,
              /*rho_max*/ rho_i,
              /*s_min*/ s_i,
              /*gamma_min*/ gamma_min_i};
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto Limiter<dim, Number>::combine_bounds(
        const Bounds &bounds_left, const Bounds &bounds_right) const -> Bounds
    {
      const auto &[rho_min_l, rho_max_l, s_min_l, gamma_min_l] = bounds_left;
      const auto &[rho_min_r, rho_max_r, s_min_r, gamma_min_r] = bounds_right;
 
      return {std::min(rho_min_l, rho_min_r),
              std::max(rho_max_l, rho_max_r),
              std::min(s_min_l, s_min_r),
              std::min(gamma_min_l, gamma_min_r)};
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    Limiter<dim, Number>::fully_relax_bounds(const Bounds &bounds,
                                             const Number &hd) const -> Bounds
    {
      const auto view = hyperbolic_system.view<dim, Number>();
 
      const auto &[rho_min, rho_max, s_min, gamma_min] = bounds;
 
      auto relaxed_bounds = bounds;
      auto &[rho_min_relaxed, rho_max_relaxed, s_min_relaxed, g_m_relaxed] =
          relaxed_bounds;
 
      /* Use r = factor * (m_i / |Omega|) ^ (1.5 / d): */
 
      Number r = std::sqrt(hd);                              // in 3D: ^ 3/6
      if constexpr (dim == 2)                                //
        r = dealii::Utilities::fixed_power<3>(std::sqrt(r)); // in 2D: ^ 3/4
      else if constexpr (dim == 1)                           //
        r = dealii::Utilities::fixed_power<3>(r);            // in 1D: ^ 3/2
      r *= parameters.relaxation_factor();
 
      constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
      rho_min_relaxed *= std::max(Number(1.) - r, Number(eps));
      rho_max_relaxed *= (Number(1.) + r);
      s_min_relaxed *= std::max(Number(1.) - r, Number(eps));
 
      /*
       * If we have a maximum compressibility constant, b, the maximum
       * bound for rho changes. See @cite ryujin-2023-4 for how to define
       * rho_max.
       */
 
      const auto numerator = (gamma_min + Number(1.)) * rho_max;
      const auto interpolation_b = view.eos_interpolation_b();
      const auto denominator =
          gamma_min - Number(1.) + ScalarNumber(2.) * interpolation_b * rho_max;
      const auto rho_compressibility_bound = numerator / denominator;
 
      rho_max_relaxed = std::min(rho_compressibility_bound, rho_max_relaxed);
 
      return relaxed_bounds;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline void
    Limiter<dim, Number>::reset(const unsigned int i,
                                const state_type &new_U_i,
                                const flux_contribution_type &new_flux_i)
    {
      U_i = new_U_i;
      flux_i = new_flux_i;
 
      /* Bounds: */
 
      auto &[rho_min, rho_max, s_min, gamma_min] = bounds_;
 
      rho_min = Number(std::numeric_limits<ScalarNumber>::max());
      rho_max = Number(0.);
      s_min = Number(std::numeric_limits<ScalarNumber>::max());
 
      const auto &[p_i, gamma_min_i, s_i, eta_i] =
          precomputed_values.template get_tensor<Number, precomputed_type>(i);
 
      gamma_min = gamma_min_i;
 
      /* Relaxation: */
 
      rho_relaxation_numerator = Number(0.);
      rho_relaxation_denominator = Number(0.);
      s_interp_max = Number(0.);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline void Limiter<dim, Number>::accumulate(
        const unsigned int *js,
        const state_type &U_j,
        const flux_contribution_type &flux_j,
        const dealii::Tensor<1, dim, Number> &scaled_c_ij,
        const state_type &affine_shift)
    {
      // TODO: Currently we only apply the affine_shift to U_ij_bar (which
      // then enters all bounds), but we do not modify s_interp and
      // rho_relaxation. When actually adding a source term to the Euler
      // equations verify that this does the right thing.
      Assert(std::max(affine_shift.norm(), Number(0.)) == Number(0.),
             dealii::ExcNotImplemented());
 
      const auto view = hyperbolic_system.view<dim, Number>();
 
      /* Bounds: */
      auto &[rho_min, rho_max, s_min, gamma_min] = bounds_;
 
      const auto rho_i = view.density(U_i);
      const auto rho_j = view.density(U_j);
 
      /* bar state shifted by an affine shift: */
      const auto U_ij_bar =
          ScalarNumber(0.5) * (U_i + U_j) -
          ScalarNumber(0.5) * contract(add(flux_j, -flux_i), scaled_c_ij) +
          affine_shift;
 
      const auto rho_ij_bar = view.density(U_ij_bar);
 
      /* Density bounds: */
 
      rho_min = std::min(rho_min, rho_ij_bar);
      rho_max = std::max(rho_max, rho_ij_bar);
 
      /* Density relaxation: */
 
      /* Use a uniform weight. */
      const auto beta_ij = Number(1.);
      rho_relaxation_numerator += beta_ij * (rho_i + rho_j);
      rho_relaxation_denominator += std::abs(beta_ij);
 
      /* Surrogate entropy bounds and relaxation: */
 
      if (view.compute_strict_bounds()) {
        /*
         * Compute strict bounds precisely as outlined in @cite ryujin-2023-4
         *
         * This means, we compute
         *  - the surrogate entropy at dof j with the gamma_min of index i,
         *  - the surrogate entropy of the bar state U_ij_bar
         *  - an interpolated surrogate entropy at (U_i + U_j) / 2 for
         *    bounds relaxation:
         */
 
        const auto s_j = view.surrogate_specific_entropy(U_j, gamma_min);
 
        const auto s_ij_bar =
            view.surrogate_specific_entropy(U_ij_bar, gamma_min);
 
        const Number s_interp = view.surrogate_specific_entropy(
            (U_i + U_j) * ScalarNumber(.5), gamma_min);
 
        s_min = std::min(s_min, s_j);
        s_min = std::min(s_min, s_ij_bar);
        s_interp_max = std::max(s_interp_max, s_interp);
 
      } else {
        /*
         * Compute a cheaper bound solely relying on the diagonal s_j
         * (computed with gamma_min_j) and the surrogate entropy s_ij_bar
         * of the bar state. We use the s_ij_bar for computing the bounds
         * relaxation as well.
         */
        const auto [p_j, gamma_min_j, s_j, eta_j] =
            precomputed_values.template get_tensor<Number, precomputed_type>(
                js);
 
        const auto s_ij_bar =
            view.surrogate_specific_entropy(U_ij_bar, gamma_min);
 
        s_min = std::min(s_min, s_j);
        s_min = std::min(s_min, s_ij_bar);
        s_interp_max = std::max(s_interp_max, s_ij_bar);
      }
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    Limiter<dim, Number>::bounds(const Number hd_i) const -> Bounds
    {
      const auto &[rho_min, rho_max, s_min, gamma_min] = bounds_;
 
      auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i);
      auto &[rho_min_relaxed, rho_max_relaxed, s_min_relaxed, g_m_relaxed] =
          relaxed_bounds;
 
      /* Apply a stricter window: */
 
      constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
 
      const auto rho_relaxation =
          ScalarNumber(2. * parameters.relaxation_factor()) *
          std::abs(rho_relaxation_numerator) /
          (std::abs(rho_relaxation_denominator) + Number(eps));
 
      const auto entropy_relaxation =
          parameters.relaxation_factor() * (s_interp_max - s_min);
 
      rho_min_relaxed = std::max(rho_min_relaxed, rho_min - rho_relaxation);
      rho_max_relaxed = std::min(rho_max_relaxed, rho_max + rho_relaxation);
      s_min_relaxed = std::max(s_min_relaxed, s_min - entropy_relaxation);
 
      return relaxed_bounds;
    }
  } // namespace EulerAEOS
} // namespace ryujin