wave_speed_estimator.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2020 - 2025 by the ryujin authors
//
 
#pragma once
 
#include <compile_time_options.h>
 
#include "hyperbolic_system.h"
 
#include <observer_pointer.h>
#include <simd.h>
 
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
 
namespace ryujin
{
  namespace EulerAEOS
  {
    template <int dim, typename Number = double>
    class WaveSpeedEstimatorView;
 
    template <typename ScalarNumber = double>
    class WaveSpeedEstimator : public dealii::ParameterAcceptor
    {
    public:
 
      template <int dim, typename Number = double>
      using View = WaveSpeedEstimatorView<dim, Number>;
 
 
 
      WaveSpeedEstimator(const HyperbolicSystem &hyperbolic_system,
                         const std::string &subsection = "/WaveSpeedEstimator")
          : ParameterAcceptor(subsection)
          , hyperbolic_system_(&hyperbolic_system)
      {
      }
 
 
 
      template <int dim, typename Number>
      auto view() const
      {
        return View<dim, Number>{
            hyperbolic_system_->template view<dim, Number>(), *this};
      }
 
    private:
 
 
      dealii::ObserverPointer<const HyperbolicSystem> hyperbolic_system_;
 
    };
 
 
    template <int dim, typename Number>
    class WaveSpeedEstimatorView
    {
    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 = 5;
 
      using primitive_type = typename std::array<Number, riemann_data_size>;
 
      using precomputed_type = typename View::precomputed_type;
 
      using PrecomputedVectorView = typename View::PrecomputedVectorView;
 
 
 
      WaveSpeedEstimatorView(
          const View &view,
          const WaveSpeedEstimator<ScalarNumber> &wave_speed_estimator)
          : view_(view)
          , wave_speed_estimator_(wave_speed_estimator)
      {
      }
 
      Number compute(const primitive_type &riemann_data_i,
                     const primitive_type &riemann_data_j) const;
 
      Number compute(const PrecomputedVectorView &pv,
                     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 c(const Number &gamma_Z) const;
 
      Number
      alpha(const Number &rho, const Number &gamma, const Number &a) const;
 
      Number p_star_RS_full(const primitive_type &riemann_data_i,
                            const primitive_type &riemann_data_j) const;
 
      Number p_star_SS_full(const primitive_type &riemann_data_i,
                            const primitive_type &riemann_data_j) const;
 
      /*
       * Compute only the failsafe the failsafe bound for \f$\tilde
       * p_2^\ast\f$ (5.11) in @cite ClaytonGuermondPopov-2022
       *
       * Cost: 0x pow, 3x division, 3x sqrt
       */
      Number p_star_failsafe(const primitive_type &riemann_data_i,
                             const primitive_type &riemann_data_j) const;
 
      /*
       * Compute a simultaneous upper bound on (5.7) second formula for
       * \tilde p_2^\ast (5.8) first formula for \tilde p_1^\ast (5.11)
       * formula for \tilde p_2^\ast in @cite ClaytonGuermondPopov-2022
       *
       * Cost: 3x pow, 9x division, 2x sqrt
       *
       * @todo improve documentation
       */
      Number p_star_interpolated(const primitive_type &riemann_data_i,
                                 const primitive_type &riemann_data_j) const;
 
 
#ifndef DOXYGEN
      /*
       * FIXME
       */
      Number f(const primitive_type &riemann_data, const Number p_star) const;
 
 
      /*
       * FIXME
       */
      Number phi(const primitive_type &riemann_data_i,
                 const primitive_type &riemann_data_j,
                 const Number p_in) const;
#endif
 
 
      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;
 
 
      Number compute_lambda(const primitive_type &riemann_data_i,
                            const primitive_type &riemann_data_j,
                            const Number p_star) const;
 
 
      primitive_type
      riemann_data_from_state(const state_type &U,
                              const Number &p,
                              const dealii::Tensor<1, dim, Number> &n_ij) const;
 
    private:
 
 
      const View view_;
      const WaveSpeedEstimator<ScalarNumber> &wave_speed_estimator_;
 
    };
  } // namespace EulerAEOS
} /* namespace ryujin */