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/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp"
      12             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
      13             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : #include "Utilities/Serialization/CharmPupable.hpp"
      18             : #include "Utilities/TMPL.hpp"
      19             : 
      20             : /// \cond
      21             : class DataVector;
      22             : namespace gsl {
      23             : template <typename T>
      24             : class not_null;
      25             : }  // namespace gsl
      26             : namespace PUP {
      27             : class er;
      28             : }  // namespace PUP
      29             : /// \endcond
      30             : 
      31             : namespace grmhd::ValenciaDivClean::BoundaryCorrections {
      32             : /*!
      33             :  * \brief An HLL Riemann solver
      34             :  *
      35             :  * Let \f$U\f$ be the evolved variable, \f$F^i\f$ the flux, and \f$n_i\f$ be the
      36             :  * outward directed unit normal to the interface. Denoting \f$F := n_i F^i\f$,
      37             :  * the HLL boundary correction is \cite Harten1983
      38             :  *
      39             :  * \f{align*}
      40             :  * G_\text{HLL} = \frac{\lambda_\text{max} F_\text{int} +
      41             :  * \lambda_\text{min} F_\text{ext}}{\lambda_\text{max} - \lambda_\text{min}}
      42             :  * - \frac{\lambda_\text{min}\lambda_\text{max}}{\lambda_\text{max} -
      43             :  *   \lambda_\text{min}} \left(U_\text{int} - U_\text{ext}\right)
      44             :  * \f}
      45             :  *
      46             :  * where "int" and "ext" stand for interior and exterior.
      47             :  * \f$\lambda_\text{min}\f$ and \f$\lambda_\text{max}\f$ are defined as
      48             :  *
      49             :  * \f{align*}
      50             :  * \lambda_\text{min} &=
      51             :  * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}, 0\right) \\
      52             :  * \lambda_\text{max} &=
      53             :  * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}, 0\right)
      54             :  * \f}
      55             :  *
      56             :  * where \f$\lambda^{+}\f$ (\f$\lambda^{-}\f$) is the largest characteristic
      57             :  * speed in the outgoing (ingoing) direction. Note the minus signs in front of
      58             :  * \f$\lambda^{\pm}_\text{ext}\f$, which is because an outgoing speed w.r.t. the
      59             :  * neighboring element is an ingoing speed w.r.t. the local element, and vice
      60             :  * versa. Similarly, the \f$F_{\text{ext}}\f$ term in \f$G_\text{HLL}\f$ has a
      61             :  * positive sign because the outward directed normal of the neighboring element
      62             :  * has the opposite sign, i.e. \f$n_i^{\text{ext}}=-n_i^{\text{int}}\f$.
      63             :  *
      64             :  * The characteristic/signal speeds are given in the documentation for
      65             :  * `grmhd::ValenciaDivClean::characteristic_speeds()`. Since the fluid is
      66             :  * travelling slower than the speed of light, the speeds we are interested in
      67             :  * are
      68             :  *
      69             :  * \f{align*}{
      70             :  *   \lambda^{\pm}&=\pm\alpha-\beta^i n_i,
      71             :  * \f}
      72             :  *
      73             :  * which correspond to the divergence cleaning field.
      74             :  *
      75             :  * \note
      76             :  * - In the strong form the `dg_boundary_terms` function returns
      77             :  *   \f$G - F_\text{int}\f$
      78             :  * - For either \f$\lambda_\text{min} = 0\f$ or \f$\lambda_\text{max} = 0\f$
      79             :  *   (i.e. all characteristics move in the same direction) the HLL boundary
      80             :  *   correction reduces to pure upwinding.
      81             :  * - Some references use \f$S\f$ instead of \f$\lambda\f$ for the
      82             :  *   signal/characteristic speeds
      83             :  * - It may be possible to use the slower speeds for the magnetic field and
      84             :  *   fluid part of the system in order to make the flux less dissipative for
      85             :  *   those variables.
      86             :  */
      87           1 : class Hll final : public BoundaryCorrection {
      88             :  public:
      89           0 :   struct LargestOutgoingCharSpeed : db::SimpleTag {
      90           0 :     using type = Scalar<DataVector>;
      91             :   };
      92           0 :   struct LargestIngoingCharSpeed : db::SimpleTag {
      93           0 :     using type = Scalar<DataVector>;
      94             :   };
      95             : 
      96           0 :   using options = tmpl::list<>;
      97           0 :   static constexpr Options::String help = {
      98             :       "Computes the HLL boundary correction term for the GRMHD system."};
      99             : 
     100           0 :   Hll() = default;
     101           0 :   Hll(const Hll&) = default;
     102           0 :   Hll& operator=(const Hll&) = default;
     103           0 :   Hll(Hll&&) = default;
     104           0 :   Hll& operator=(Hll&&) = default;
     105           0 :   ~Hll() override = default;
     106             : 
     107             :   /// \cond
     108             :   explicit Hll(CkMigrateMessage* /*unused*/);
     109             :   using PUP::able::register_constructor;
     110             :   WRAPPED_PUPable_decl_template(Hll);  // NOLINT
     111             :   /// \endcond
     112           0 :   void pup(PUP::er& p) override;  // NOLINT
     113             : 
     114           0 :   std::unique_ptr<BoundaryCorrection> get_clone() const override;
     115             : 
     116           0 :   using dg_package_field_tags =
     117             :       tmpl::list<Tags::TildeD, Tags::TildeYe, Tags::TildeTau,
     118             :                  Tags::TildeS<Frame::Inertial>, Tags::TildeB<Frame::Inertial>,
     119             :                  Tags::TildePhi, ::Tags::NormalDotFlux<Tags::TildeD>,
     120             :                  ::Tags::NormalDotFlux<Tags::TildeYe>,
     121             :                  ::Tags::NormalDotFlux<Tags::TildeTau>,
     122             :                  ::Tags::NormalDotFlux<Tags::TildeS<Frame::Inertial>>,
     123             :                  ::Tags::NormalDotFlux<Tags::TildeB<Frame::Inertial>>,
     124             :                  ::Tags::NormalDotFlux<Tags::TildePhi>,
     125             :                  LargestOutgoingCharSpeed, LargestIngoingCharSpeed>;
     126           0 :   using dg_package_data_temporary_tags =
     127             :       tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>>;
     128           0 :   using dg_package_data_primitive_tags = tmpl::list<>;
     129           0 :   using dg_package_data_volume_tags = tmpl::list<>;
     130           0 :   using dg_boundary_terms_volume_tags = tmpl::list<>;
     131             : 
     132           0 :   static double dg_package_data(
     133             :       gsl::not_null<Scalar<DataVector>*> packaged_tilde_d,
     134             :       gsl::not_null<Scalar<DataVector>*> packaged_tilde_ye,
     135             :       gsl::not_null<Scalar<DataVector>*> packaged_tilde_tau,
     136             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> packaged_tilde_s,
     137             :       gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> packaged_tilde_b,
     138             :       gsl::not_null<Scalar<DataVector>*> packaged_tilde_phi,
     139             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_d,
     140             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_ye,
     141             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_tau,
     142             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
     143             :           packaged_normal_dot_flux_tilde_s,
     144             :       gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
     145             :           packaged_normal_dot_flux_tilde_b,
     146             :       gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_phi,
     147             :       gsl::not_null<Scalar<DataVector>*> packaged_largest_outgoing_char_speed,
     148             :       gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_char_speed,
     149             : 
     150             :       const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
     151             :       const Scalar<DataVector>& tilde_tau,
     152             :       const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
     153             :       const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
     154             :       const Scalar<DataVector>& tilde_phi,
     155             : 
     156             :       const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_d,
     157             :       const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_ye,
     158             :       const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_tau,
     159             :       const tnsr::Ij<DataVector, 3, Frame::Inertial>& flux_tilde_s,
     160             :       const tnsr::IJ<DataVector, 3, Frame::Inertial>& flux_tilde_b,
     161             :       const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_phi,
     162             : 
     163             :       const Scalar<DataVector>& lapse,
     164             :       const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
     165             : 
     166             :       const tnsr::i<DataVector, 3, Frame::Inertial>& normal_covector,
     167             :       const tnsr::I<DataVector, 3, Frame::Inertial>& normal_vector,
     168             :       const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
     169             :       /*mesh_velocity*/,
     170             :       const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity);
     171             : 
     172           0 :   static void dg_boundary_terms(
     173             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_d,
     174             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_ye,
     175             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_tau,
     176             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
     177             :           boundary_correction_tilde_s,
     178             :       gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
     179             :           boundary_correction_tilde_b,
     180             :       gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_phi,
     181             :       const Scalar<DataVector>& tilde_d_int,
     182             :       const Scalar<DataVector>& tilde_ye_int,
     183             :       const Scalar<DataVector>& tilde_tau_int,
     184             :       const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_int,
     185             :       const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b_int,
     186             :       const Scalar<DataVector>& tilde_phi_int,
     187             :       const Scalar<DataVector>& normal_dot_flux_tilde_d_int,
     188             :       const Scalar<DataVector>& normal_dot_flux_tilde_ye_int,
     189             :       const Scalar<DataVector>& normal_dot_flux_tilde_tau_int,
     190             :       const tnsr::i<DataVector, 3, Frame::Inertial>&
     191             :           normal_dot_flux_tilde_s_int,
     192             :       const tnsr::I<DataVector, 3, Frame::Inertial>&
     193             :           normal_dot_flux_tilde_b_int,
     194             :       const Scalar<DataVector>& normal_dot_flux_tilde_phi_int,
     195             :       const Scalar<DataVector>& largest_outgoing_char_speed_int,
     196             :       const Scalar<DataVector>& largest_ingoing_char_speed_int,
     197             :       const Scalar<DataVector>& tilde_d_ext,
     198             :       const Scalar<DataVector>& tilde_ye_ext,
     199             :       const Scalar<DataVector>& tilde_tau_ext,
     200             :       const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_ext,
     201             :       const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b_ext,
     202             :       const Scalar<DataVector>& tilde_phi_ext,
     203             :       const Scalar<DataVector>& normal_dot_flux_tilde_d_ext,
     204             :       const Scalar<DataVector>& normal_dot_flux_tilde_ye_ext,
     205             :       const Scalar<DataVector>& normal_dot_flux_tilde_tau_ext,
     206             :       const tnsr::i<DataVector, 3, Frame::Inertial>&
     207             :           normal_dot_flux_tilde_s_ext,
     208             :       const tnsr::I<DataVector, 3, Frame::Inertial>&
     209             :           normal_dot_flux_tilde_b_ext,
     210             :       const Scalar<DataVector>& normal_dot_flux_tilde_phi_ext,
     211             :       const Scalar<DataVector>& largest_outgoing_char_speed_ext,
     212             :       const Scalar<DataVector>& largest_ingoing_char_speed_ext,
     213             :       dg::Formulation dg_formulation);
     214             : };
     215             : 
     216           0 : bool operator==(const Hll& lhs, const Hll& rhs);
     217           0 : bool operator!=(const Hll& lhs, const Hll& rhs);
     218             : }  // namespace grmhd::ValenciaDivClean::BoundaryCorrections