SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/NewtonianEuler/BoundaryCorrections - Hllc.hpp Hit Total Coverage
Commit: 2ae2b99409ac582030d56a4560a92a3e066a7e54 Lines: 1 26 3.8 %
Date: 2022-01-15 08:40:38
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <memory>
       7             : #include <optional>
       8             : 
       9             : #include "DataStructures/DataBox/Prefixes.hpp"
      10             : #include "DataStructures/Tensor/TypeAliases.hpp"
      11             : #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp"
      12             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      13             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
      14             : #include "Options/Options.hpp"
      15             : #include "Parallel/CharmPupable.hpp"
      16             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      17             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      18             : #include "Utilities/Gsl.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : 
      21             : /// \cond
      22             : class DataVector;
      23             : namespace gsl {
      24             : template <typename T>
      25             : class not_null;
      26             : }  // namespace gsl
      27             : namespace PUP {
      28             : class er;
      29             : }  // namespace PUP
      30             : /// \endcond
      31             : 
      32             : namespace NewtonianEuler::BoundaryCorrections {
      33             : /*!
      34             :  * \brief The HLLC (Harten-Lax-van Leer-Contact) Riemann solver for the
      35             :  * NewtonianEuler system
      36             :  *
      37             :  * Let \f$U\f$ be the evolved variable, \f$F^i\f$ the flux, \f$v^i\f$ the
      38             :  * spatial velocity, and \f$n_i\f$ be the outward directed unit normal to the
      39             :  * interface. Denoting \f$F := n_i F^i\f$ and \f$v:=n_iv^i\f$, the HLLC boundary
      40             :  * correction is \cite Toro1994
      41             :  *
      42             :  * \f{align*}
      43             :  * G_\text{HLLC} = \left\{\begin{array}{lcl}
      44             :  * F_\text{int} & \text{if} & 0 \leq \lambda_\text{min} \\
      45             :  * F_\text{int} + \lambda_\text{min}(U_\text{*int} - U_\text{int}) & \text{if} &
      46             :  * \lambda_\text{min} \leq 0 \leq \lambda_\text{*} \\
      47             :  * -F_\text{ext} + \lambda_\text{max}(U_\text{*ext} - U_\text{ext}) & \text{if}
      48             :  * & \lambda_\text{*} \leq 0 \leq \lambda_\text{max} \\
      49             :  * -F_\text{ext} & \text{if} &  \lambda_\text{max} \leq 0
      50             :  * \end{array}\right\}
      51             :  * \f}
      52             :  *
      53             :  * where "int" and "ext" stand for interior and exterior.
      54             :  *
      55             :  * Intermediate ('star') states are given by
      56             :  *
      57             :  * \f{align*}
      58             :  * U_\text{*int} = \left(\frac{\lambda_\text{min} - v_\text{int}}
      59             :  * {\lambda_\text{min} - \lambda_*}\right)
      60             :  * \left[\begin{array}{c}
      61             :  * \displaystyle \rho_\text{int} \\
      62             :  * \displaystyle \rho_\text{int}[v_\text{int}^x + (\lambda_* - v_\text{int})
      63             :  * n_x^\text{int}] \\
      64             :  * \displaystyle \rho_\text{int}[v_\text{int}^y + (\lambda_* - v_\text{int})
      65             :  * n_y^\text{int}] \\
      66             :  * \displaystyle \rho_\text{int}[v_\text{int}^z + (\lambda_* - v_\text{int})
      67             :  * n_z^\text{int}] \\
      68             :  * \displaystyle E_\text{int} + p_\text{int} \frac{\lambda_* - v_\text{int}}
      69             :  * {\lambda_\text{min} - v_\text{int}} + \rho_\text{int}\lambda_*(\lambda_* -
      70             :  * v_\text{int}) \end{array}\right]
      71             :  * \f}
      72             :  *
      73             :  * and
      74             :  *
      75             :  * \f{align*}
      76             :  * U_\text{*ext} = \left(\frac{\lambda_\text{max} + v_\text{ext}}
      77             :  * {\lambda_\text{max} - \lambda_*}\right)
      78             :  * \left[\begin{array}{c}
      79             :  * \displaystyle \rho_\text{ext} \\
      80             :  * \displaystyle \rho_\text{ext}[-v_\text{ext}^x - (\lambda_* + v_\text{ext})
      81             :  * n_x^\text{ext}] \\
      82             :  * \displaystyle \rho_\text{ext}[-v_\text{ext}^y - (\lambda_* + v_\text{ext})
      83             :  * n_y^\text{ext}] \\
      84             :  * \displaystyle \rho_\text{ext}[-v_\text{ext}^z - (\lambda_* + v_\text{ext})
      85             :  * n_z^\text{ext}] \\
      86             :  * \displaystyle E_\text{ext} + p_\text{ext} \frac{\lambda_* + v_\text{ext}}
      87             :  * {\lambda_\text{max} + v_\text{ext}} + \rho_\text{ext}\lambda_*(\lambda_* +
      88             :  * v_\text{ext}) \end{array}\right].
      89             :  * \f}
      90             :  *
      91             :  * The contact wave speed \f$\lambda_*\f$ is \cite Toro2009
      92             :  *
      93             :  * \f{align*}
      94             :  * \lambda_* = \frac
      95             :  * { p_\text{ext} - p_\text{int} +
      96             :  * \rho_\text{int}v_\text{int}(\lambda_\text{min}-v_\text{int})
      97             :  * + \rho_\text{ext}v_\text{ext}(\lambda_\text{max}+v_\text{ext})}
      98             :  * { \rho_\text{int}(\lambda_\text{min}-v_\text{int}) -
      99             :  * \rho_\text{ext}(\lambda_\text{max} + v_\text{ext})}.
     100             :  * \f}
     101             :  *
     102             :  * \f$\lambda_\text{min}\f$ and \f$\lambda_\text{max}\f$ are estimated by
     103             :  * \cite Davis1988
     104             :  *
     105             :  * \f{align*}
     106             :  * \lambda_\text{min} &=
     107             :  * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}\right) \\
     108             :  * \lambda_\text{max} &=
     109             :  * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}\right)
     110             :  * \f}
     111             :  *
     112             :  * where \f$\lambda^{+}\f$ (\f$\lambda^{-}\f$) is the largest characteristic
     113             :  * speed in the outgoing (ingoing) direction for each domain.
     114             :  *
     115             :  * Note the minus signs in front of \f$\lambda^{\pm}_\text{ext}\f$, which is
     116             :  * because an outgoing speed w.r.t. the neighboring element is an ingoing speed
     117             :  * w.r.t. the local element, and vice versa. Similarly, the \f$F_{\text{ext}}\f$
     118             :  * term in \f$G_\text{HLLC}\f$ and the \f$v_\text{ext}\f$ term in
     119             :  * \f$U_\text{*ext}\f$ have a positive sign because the outward directed normal
     120             :  * of the neighboring element has the opposite sign, i.e.
     121             :  * \f$n_i^{\text{ext}}=-n_i^{\text{int}}\f$.
     122             :  *
     123             :  * For the NewtonianEuler system, \f$\lambda^\pm\f$ are given as
     124             :  *
     125             :  * \f{align*}
     126             :  * \lambda^\pm = v \pm c_s
     127             :  * \f}
     128             :  *
     129             :  * where \f$c_s\f$ is the sound speed.
     130             :  *
     131             :  * \note
     132             :  * - In the strong form the `dg_boundary_terms` function returns \f$G -
     133             :  *   F_\text{int}\f$
     134             :  * - In the implementation, we use
     135             :  *
     136             :  *   \f{align*}
     137             :  *   G_\text{HLLC} = \left\{\begin{array}{lcl}
     138             :  *   F_\text{int} + \lambda_\text{min}(U_\text{*int} - U_\text{int}) & \text{if}
     139             :  *   & 0 \leq \lambda_\text{*} \\
     140             :  *   -F_\text{ext} + \lambda_\text{max}(U_\text{*ext} - U_\text{ext}) &
     141             :  *   \text{if} & \lambda_\text{*} \leq 0 \\
     142             :  *   \end{array}\right\},
     143             :  *   \f}
     144             :  *
     145             :  *   with
     146             :  *
     147             :  *   \f{align*}
     148             :  *   \lambda_\text{min} &=
     149             :  *   \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}, 0\right) \\
     150             :  *   \lambda_\text{max} &=
     151             :  *   \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}, 0\right).
     152             :  *   \f}
     153             :  *
     154             :  *   Provided that \f$\lambda_*\f$ falls in the correct range i.e
     155             :  *
     156             :  *   \f{align*}
     157             :  *   \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}\right)
     158             :  *   < \lambda_* <
     159             :  *   \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}\right),
     160             :  *   \f}
     161             :  *
     162             :  *   this prescription recovers the original HLLC boundary correction for all
     163             :  *   four cases. For either \f$\lambda_\text{min} = 0\f$ or
     164             :  *   \f$\lambda_\text{max} = 0\f$ (i.e. all characteristics move in the same
     165             :  *   direction), boundary correction reduces to pure upwinding.
     166             :  * - Some references use \f$S\f$ instead of \f$\lambda\f$ for the
     167             :  *   signal/characteristic speeds.
     168             :  */
     169             : template <size_t Dim>
     170           1 : class Hllc final : public BoundaryCorrection<Dim> {
     171             :  private:
     172           0 :   struct InterfaceUnitNormal : db::SimpleTag {
     173           0 :     using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
     174             :   };
     175           0 :   struct NormalDotVelocity : db::SimpleTag {
     176           0 :     using type = Scalar<DataVector>;
     177             :   };
     178             : 
     179             :  public:
     180           0 :   struct LargestOutgoingCharSpeed : db::SimpleTag {
     181           0 :     using type = Scalar<DataVector>;
     182             :   };
     183           0 :   struct LargestIngoingCharSpeed : db::SimpleTag {
     184           0 :     using type = Scalar<DataVector>;
     185             :   };
     186             : 
     187           0 :   using options = tmpl::list<>;
     188           0 :   static constexpr Options::String help = {
     189             :       "Computes the HLLC boundary correction term for the "
     190             :       "Newtonian Euler/hydrodynamics system."};
     191             : 
     192           0 :   Hllc() = default;
     193           0 :   Hllc(const Hllc&) = default;
     194           0 :   Hllc& operator=(const Hllc&) = default;
     195           0 :   Hllc(Hllc&&) = default;
     196           0 :   Hllc& operator=(Hllc&&) = default;
     197           0 :   ~Hllc() override = default;
     198             : 
     199             :   /// \cond
     200             :   explicit Hllc(CkMigrateMessage* msg);
     201             :   using PUP::able::register_constructor;
     202             :   WRAPPED_PUPable_decl_template(Hllc);  // NOLINT
     203             :   /// \endcond
     204           0 :   void pup(PUP::er& p) override;  // NOLINT
     205             : 
     206           0 :   std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const override;
     207             : 
     208           0 :   using dg_package_field_tags = tmpl::list<
     209             :       Tags::MassDensityCons, Tags::MomentumDensity<Dim>, Tags::EnergyDensity,
     210             :       Tags::Pressure<DataVector>, ::Tags::NormalDotFlux<Tags::MassDensityCons>,
     211             :       ::Tags::NormalDotFlux<Tags::MomentumDensity<Dim>>,
     212             :       ::Tags::NormalDotFlux<Tags::EnergyDensity>, InterfaceUnitNormal,
     213             :       NormalDotVelocity, LargestOutgoingCharSpeed, LargestIngoingCharSpeed>;
     214           0 :   using dg_package_data_temporary_tags = tmpl::list<>;
     215           0 :   using dg_package_data_primitive_tags =
     216             :       tmpl::list<NewtonianEuler::Tags::Velocity<DataVector, Dim>,
     217             :                  NewtonianEuler::Tags::SpecificInternalEnergy<DataVector>>;
     218           0 :   using dg_package_data_volume_tags =
     219             :       tmpl::list<hydro::Tags::EquationOfStateBase>;
     220             : 
     221             :   template <size_t ThermodynamicDim>
     222           0 :   double dg_package_data(
     223             :       gsl::not_null<Scalar<DataVector>*> packaged_mass_density,
     224             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     225             :           packaged_momentum_density,
     226             :       gsl::not_null<Scalar<DataVector>*> packaged_energy_density,
     227             :       gsl::not_null<Scalar<DataVector>*> packaged_pressure,
     228             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_mass_density,
     229             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     230             :           packaged_normal_dot_flux_momentum_density,
     231             :       gsl::not_null<Scalar<DataVector>*>
     232             :           packaged_normal_dot_flux_energy_density,
     233             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     234             :           packaged_interface_unit_normal,
     235             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_velocity,
     236             :       gsl::not_null<Scalar<DataVector>*> packaged_largest_outgoing_char_speed,
     237             :       gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_char_speed,
     238             : 
     239             :       const Scalar<DataVector>& mass_density,
     240             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density,
     241             :       const Scalar<DataVector>& energy_density,
     242             : 
     243             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_mass_density,
     244             :       const tnsr::IJ<DataVector, Dim, Frame::Inertial>& flux_momentum_density,
     245             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_energy_density,
     246             : 
     247             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& velocity,
     248             :       const Scalar<DataVector>& specific_internal_energy,
     249             : 
     250             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
     251             :       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
     252             :       /*mesh_velocity*/,
     253             :       const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
     254             :       const EquationsOfState::EquationOfState<false, ThermodynamicDim>&
     255             :           equation_of_state) const;
     256             : 
     257           0 :   void dg_boundary_terms(
     258             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_mass_density,
     259             :       gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
     260             :           boundary_correction_momentum_density,
     261             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_energy_density,
     262             :       const Scalar<DataVector>& mass_density_int,
     263             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_int,
     264             :       const Scalar<DataVector>& energy_density_int,
     265             :       const Scalar<DataVector>& pressure_int,
     266             :       const Scalar<DataVector>& normal_dot_flux_mass_density_int,
     267             :       const tnsr::I<DataVector, Dim, Frame::Inertial>&
     268             :           normal_dot_flux_momentum_density_int,
     269             :       const Scalar<DataVector>& normal_dot_flux_energy_density_int,
     270             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     271             :           interface_unit_normal_int,
     272             :       const Scalar<DataVector>& normal_dot_velocity_int,
     273             :       const Scalar<DataVector>& largest_outgoing_char_speed_int,
     274             :       const Scalar<DataVector>& largest_ingoing_char_speed_int,
     275             :       const Scalar<DataVector>& mass_density_ext,
     276             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& momentum_density_ext,
     277             :       const Scalar<DataVector>& energy_density_ext,
     278             :       const Scalar<DataVector>& pressure_ext,
     279             :       const Scalar<DataVector>& normal_dot_flux_mass_density_ext,
     280             :       const tnsr::I<DataVector, Dim, Frame::Inertial>&
     281             :           normal_dot_flux_momentum_density_ext,
     282             :       const Scalar<DataVector>& normal_dot_flux_energy_density_ext,
     283             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     284             :           interface_unit_normal_ext,
     285             :       const Scalar<DataVector>& normal_dot_velocity_ext,
     286             :       const Scalar<DataVector>& largest_outgoing_char_speed_ext,
     287             :       const Scalar<DataVector>& largest_ingoing_char_speed_ext,
     288             :       dg::Formulation dg_formulation) const;
     289             : };
     290             : }  // namespace NewtonianEuler::BoundaryCorrections

Generated by: LCOV version 1.14