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/CurvedScalarWave/BoundaryCorrections/BoundaryCorrection.hpp"
      12             : #include "Evolution/Systems/CurvedScalarWave/Characteristics.hpp"
      13             : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
      14             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
      15             : #include "Options/String.hpp"
      16             : #include "PointwiseFunctions/GeneralRelativity/Tags.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 CurvedScalarWave::BoundaryCorrections {
      32             : /*!
      33             :  * \brief Computes the upwind multipenalty boundary correction for scalar wave
      34             :  * in curved spacetime.
      35             :  *
      36             :  * \details This implements the upwind multipenalty boundary correction term.
      37             :  * The general form is given by:
      38             :  *
      39             :  * \f[
      40             :  * G = T^{\rm ext}\Lambda^{\rm ext,-} {T^{\rm ext}}^{-1} U^{\rm ext}
      41             :  *      + T^{\rm int}\Lambda^{\rm int,+} {T^{\rm int}}^{-1} U^{\rm int}
      42             :  *   = T^{\rm ext}\Lambda^{\rm ext,-} V^{\rm ext}
      43             :  *      + T^{\rm int}\Lambda^{\rm int,+} V^{\rm int},
      44             :  * \f]
      45             :  *
      46             :  * where
      47             :  *
      48             :  *  - \f$G\f$ is a vector of numerical upwind fluxes dotted with the interface
      49             :  *       normal for all evolved variables;
      50             :  *  - \f$U\f$ is a vector of all evolved variables;
      51             :  *  - \f$T\f$ is a matrix whose columns are the eigenvectors of the
      52             :  *       characteristic matrix for the evolution system. It maps the
      53             :  *       evolved variables to characteristic variables \f$V\f$, s.t.
      54             :  *       \f$V := T^{-1}\cdot U\f$; and
      55             :  *  - \f$\Lambda^\pm\f$ are diagonal matrices containing positive / negative
      56             :  *       eigenvalues of \f$T\f$ as its diagonal elements, with the rest set to
      57             :  *       \f$ 0\f$.
      58             :  *
      59             :  * The superscripts \f${\rm int}\f$ and \f${\rm ext}\f$ indicate that the
      60             :  * corresponding variable at the element interface comes from the _interior_ or
      61             :  * _exterior_ of the element. Exterior of the element is naturally the interior
      62             :  * of its neighboring element. The sign of characteristic speeds indicate the
      63             :  * direction of propagation of the corresponding characteristic field with
      64             :  * respect to the interface normal that the field has been computed along, with
      65             :  * negative speeds indicating incoming characteristics, and positive speeds
      66             :  * indicating outgoing characteristics. The expressions implemented here differ
      67             :  * from Eq.(63) of \cite Teukolsky2015ega as that boundary term does not
      68             :  * consistently treat both sides of the interface on the same footing. Unlike in
      69             :  * \cite Teukolsky2015ega, in code the interface normal vector on the interior
      70             :  * side and the one on the exterior side point in opposite directions, and the
      71             :  * characteristic speeds end up having different signs. This class therefore
      72             :  * computes:
      73             :  *
      74             :  * \f[
      75             :  * G = T^{\rm ext}w_{\rm ext} V^{\rm ext} - T^{\rm int}w_{\rm int} V^{\rm int},
      76             :  * \f]
      77             :  *
      78             :  * with weights \f$w_{\rm ext} = -\Theta(\Lambda^{\rm ext})\cdot\Lambda^{\rm
      79             :  * ext}\f$, and \f$w_{\rm int} = \Theta(-\Lambda^{\rm int})\cdot\Lambda^{\rm
      80             :  * int}\f$, where \f$\Theta\f$ is the Heaviside function centered at zero,
      81             :  * \f$\Lambda = \Lambda^+ + \Lambda^-\f$, and the dot operator \f$(\cdot)\f$
      82             :  * indicates an element-wise product.
      83             :  */
      84             : template <size_t Dim>
      85           1 : class UpwindPenalty final : public BoundaryCorrection<Dim> {
      86             :  private:
      87           0 :   struct CharSpeedsTensor : db::SimpleTag {
      88           0 :     using type = tnsr::a<DataVector, 3, Frame::Inertial>;
      89             :   };
      90             : 
      91             :  public:
      92           0 :   using options = tmpl::list<>;
      93           0 :   static constexpr Options::String help = {
      94             :       "Computes the UpwindPenalty boundary correction term for the scalar wave "
      95             :       "system in curved spacetime."};
      96             : 
      97           0 :   UpwindPenalty() = default;
      98           0 :   UpwindPenalty(const UpwindPenalty&) = default;
      99           0 :   UpwindPenalty& operator=(const UpwindPenalty&) = default;
     100           0 :   UpwindPenalty(UpwindPenalty&&) = default;
     101           0 :   UpwindPenalty& operator=(UpwindPenalty&&) = default;
     102           0 :   ~UpwindPenalty() override = default;
     103             : 
     104             :   /// \cond
     105             :   explicit UpwindPenalty(CkMigrateMessage* msg);
     106             :   using PUP::able::register_constructor;
     107             :   WRAPPED_PUPable_decl_template(UpwindPenalty);  // NOLINT
     108             :   /// \endcond
     109           0 :   void pup(PUP::er& p) override;  // NOLINT
     110             : 
     111           0 :   std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const override;
     112             : 
     113           0 :   using dg_package_field_tags =
     114             :       tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
     115             :                  Tags::ConstraintGamma2,
     116             :                  ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<
     117             :                      Dim, Frame::Inertial>>,
     118             :                  CharSpeedsTensor>;
     119           0 :   using dg_package_data_temporary_tags =
     120             :       tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, Dim>,
     121             :                  Tags::ConstraintGamma1, Tags::ConstraintGamma2>;
     122           0 :   using dg_package_data_volume_tags = tmpl::list<>;
     123             : 
     124           0 :   double dg_package_data(
     125             :       gsl::not_null<Scalar<DataVector>*> packaged_v_psi,
     126             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> packaged_v_zero,
     127             :       gsl::not_null<Scalar<DataVector>*> packaged_v_plus,
     128             :       gsl::not_null<Scalar<DataVector>*> packaged_v_minus,
     129             :       gsl::not_null<Scalar<DataVector>*> packaged_gamma2,
     130             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     131             :           packaged_interface_unit_normal,
     132             :       gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
     133             :           packaged_char_speeds,
     134             : 
     135             :       const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
     136             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
     137             : 
     138             :       const Scalar<DataVector>& lapse,
     139             :       const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
     140             :       const Scalar<DataVector>& constraint_gamma1,
     141             :       const Scalar<DataVector>& constraint_gamma2,
     142             : 
     143             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
     144             :       const tnsr::I<DataVector, Dim, Frame::Inertial>&
     145             :           interface_unit_normal_vector,
     146             :       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
     147             :       /*mesh_velocity*/,
     148             :       const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
     149             : 
     150           0 :   void dg_boundary_terms(
     151             :       gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
     152             :       gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
     153             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     154             :           phi_boundary_correction,
     155             : 
     156             :       const Scalar<DataVector>& v_psi_int,
     157             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero_int,
     158             :       const Scalar<DataVector>& v_plus_int,
     159             :       const Scalar<DataVector>& v_minus_int,
     160             :       const Scalar<DataVector>& gamma2_int,
     161             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     162             :           interface_unit_normal_int,
     163             :       const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
     164             : 
     165             :       const Scalar<DataVector>& v_psi_ext,
     166             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero_ext,
     167             :       const Scalar<DataVector>& v_plus_ext,
     168             :       const Scalar<DataVector>& v_minus_ext,
     169             :       const Scalar<DataVector>& gamma2_ext,
     170             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     171             :           interface_unit_normal_ext,
     172             :       const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext,
     173             :       dg::Formulation /*dg_formulation*/) const;
     174             : };
     175             : }  // namespace CurvedScalarWave::BoundaryCorrections