riemann_solver.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 "hyperbolic_system.h"
 
#include <simd.h>
 
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
 
namespace ryujin
{
  namespace Euler
  {
    template <typename ScalarNumber = double>
    class RiemannSolverParameters : public dealii::ParameterAcceptor
    {
    public:
      RiemannSolverParameters(const std::string &subsection = "/RiemannSolver")
          : ParameterAcceptor(subsection)
      {
        if constexpr (std::is_same<ScalarNumber, double>::value)
          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_ = 0;
        add_parameter("newton max iterations",
                      newton_max_iterations_,
                      "Maximal number of quadratic newton iterations performed "
                      "during limiting");
      }
 
      ACCESSOR_READ_ONLY(newton_tolerance);
      ACCESSOR_READ_ONLY(newton_max_iterations);
 
    private:
      ScalarNumber newton_tolerance_;
      unsigned int newton_max_iterations_;
    };
 
 
    template <int dim, typename Number = double>
    class RiemannSolver
    {
    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;
 
      static constexpr unsigned int riemann_data_size = 4;
 
      using primitive_type = std::array<Number, riemann_data_size>;
 
      using precomputed_type = typename View::precomputed_type;
 
      using PrecomputedVector = typename View::PrecomputedVector;
 
      using Parameters = RiemannSolverParameters<ScalarNumber>;
 
 
 
      RiemannSolver(const HyperbolicSystem &hyperbolic_system,
                    const Parameters &parameters,
                    const PrecomputedVector &precomputed_values)
          : hyperbolic_system(hyperbolic_system)
          , parameters(parameters)
          , precomputed_values(precomputed_values)
      {
      }
 
      Number compute(const primitive_type &riemann_data_i,
                     const primitive_type &riemann_data_j) const;
 
      Number compute(const state_type &U_i,
                     const state_type &U_j,
                     const unsigned int i,
                     const unsigned int *js,
                     const dealii::Tensor<1, dim, Number> &n_ij) const;
 
 
    protected:
 
      Number f(const primitive_type &riemann_data, const Number p_star) const;
 
 
      Number df(const primitive_type &riemann_data, const Number &p_star) const;
 
 
      Number phi(const primitive_type &riemann_data_i,
                 const primitive_type &riemann_data_j,
                 const Number p_in) const;
 
 
      Number dphi(const primitive_type &riemann_data_i,
                  const primitive_type &riemann_data_j,
                  const Number &p) const;
 
 
      Number phi_of_p_max(const primitive_type &riemann_data_i,
                          const primitive_type &riemann_data_j) const;
 
 
      Number lambda1_minus(const primitive_type &riemann_data,
                           const Number p_star) const;
 
 
      Number lambda3_plus(const primitive_type &primitive_state,
                          const Number p_star) const;
 
 
      std::array<Number, 2> compute_gap(const primitive_type &riemann_data_i,
                                        const primitive_type &riemann_data_j,
                                        const Number p_1,
                                        const Number p_2) const;
 
 
      Number compute_lambda(const primitive_type &riemann_data_i,
                            const primitive_type &riemann_data_j,
                            const Number p_star) const;
 
 
      Number p_star_two_rarefaction(const primitive_type &riemann_data_i,
                                    const primitive_type &riemann_data_j) const;
 
      Number p_star_failsafe(const primitive_type &riemann_data_i,
                             const primitive_type &riemann_data_j) const;
 
 
      primitive_type
      riemann_data_from_state(const state_type &U,
                              const dealii::Tensor<1, dim, Number> &n_ij) const;
 
    private:
      const HyperbolicSystem &hyperbolic_system;
      const Parameters &parameters;
      const PrecomputedVector &precomputed_values;
    };
  } // namespace Euler
} // namespace ryujin