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/ScalarWave/BoundaryCorrections/BoundaryCorrection.hpp"
      12             : #include "Evolution/Systems/ScalarWave/Tags.hpp"
      13             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "Utilities/Serialization/CharmPupable.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : namespace gsl {
      21             : template <typename T>
      22             : class not_null;
      23             : }  // namespace gsl
      24             : namespace PUP {
      25             : class er;
      26             : }  // namespace PUP
      27             : /// \endcond
      28             : 
      29             : namespace ScalarWave::BoundaryCorrections {
      30             : /*!
      31             :  * \brief Computes the scalar wave upwind multipenalty boundary
      32             :  * correction.
      33             :  *
      34             :  * This implements the upwind multipenalty boundary correction term
      35             :  * \f$D_\beta\f$. The general form is given by:
      36             :  *
      37             :  * \f{align*}{
      38             :  *   D_\beta =
      39             :  *   T_{\beta\hat{\beta}}^{\mathrm{ext}}
      40             :  *   \Lambda^{\mathrm{ext},-}_{\hat{\beta}\hat{\alpha}}
      41             :  *   v^{\mathrm{ext}}_{\hat{\alpha}}
      42             :  *   -T_{\beta\hat{\beta}}^{\mathrm{int}}
      43             :  *   \Lambda^{\mathrm{int},-}_{\hat{\beta}\hat{\alpha}}
      44             :  *   v^{\mathrm{int}}_{\hat{\alpha}}.
      45             :  * \f}
      46             :  *
      47             :  * We denote the evolved fields by \f$u_{\alpha}\f$, the characteristic fields
      48             :  * by \f$v_{\hat{\alpha}}\f$, and implicitly sum over reapeated indices.
      49             :  * \f$T_{\beta\hat{\beta}}\f$ transforms characteristic fields to evolved
      50             :  * fields, while \f$\Lambda_{\hat{\beta}\hat{\alpha}}^-\f$ is a diagonal matrix
      51             :  * with only the negative characteristic speeds, and has zeros on the diagonal
      52             :  * for all other entries. The int and ext superscripts denote quantities on the
      53             :  * internal and external side of the mortar. Note that Eq. (6.3) of
      54             :  * \cite Teukolsky2015ega is not exactly what's implemented since that boundary
      55             :  * term does not consistently treat both sides of the interface on the same
      56             :  * footing.
      57             :  *
      58             :  * For the scalar wave system the correction is:
      59             :  *
      60             :  * \f{align}{
      61             :  *   D_{\Psi} &= \lambda_{v^{\Psi}}^{\mathrm{ext},-}
      62             :  *                v^{\mathrm{ext},\Psi}
      63             :  *                - \lambda_{v^{\Psi}}^{\mathrm{int},-}
      64             :  *                v^{\mathrm{int},\Psi}, \\
      65             :  *   D_{\Pi} &= \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{ext},-}
      66             :  *             v^{\mathrm{ext},+} +
      67             :  *             \lambda_{v^-}^{\mathrm{ext},-}
      68             :  *             v^{\mathrm{ext},-}\right)
      69             :  *             + \lambda_{v^\Psi}^{\mathrm{ext},-}\gamma_2
      70             :  *             v^{\mathrm{ext},\Psi}
      71             :  *             \notag \\
      72             :  *           &-\frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
      73             :  *             v^{\mathrm{int},+} +
      74             :  *             \lambda_{v^-}^{\mathrm{int},-}
      75             :  *             v^{\mathrm{int},-}\right)
      76             :  *             - \lambda_{v^\Psi}^{\mathrm{int},-}\gamma_2
      77             :  *             v^{\mathrm{int},\Psi} , \\
      78             :  *   D_{\Phi_{i}}
      79             :  *              &= \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{ext},-}
      80             :  *                v^{\mathrm{ext},+}
      81             :  *                - \lambda_{v^-}^{\mathrm{ext},-}
      82             :  *                v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
      83             :  *                + \lambda_{v^0}^{\mathrm{ext},-}
      84             :  *                v^{\mathrm{ext},0}_{i}
      85             :  *                \notag \\
      86             :  *              &-
      87             :  *                \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
      88             :  *                v^{\mathrm{int},+}
      89             :  *                - \lambda_{v^-}^{\mathrm{int},-}
      90             :  *                v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
      91             :  *                - \lambda_{v^0}^{\mathrm{int},-}
      92             :  *                v^{\mathrm{int},0}_{i},
      93             :  * \f}
      94             :  *
      95             :  * with characteristic fields
      96             :  *
      97             :  * \f{align}{
      98             :  *   v^{\Psi} &= \Psi, \\
      99             :  *   v^{0}_{i} &= (\delta^k_i-n^k n_i)\Phi_{k}, \\
     100             :  *   v^{\pm} &= \Pi\pm n^i\Phi_{i} -\gamma_2 \Psi,
     101             :  * \f}
     102             :  *
     103             :  * and characteristic speeds
     104             :  *
     105             :  * \f{align}{
     106             :  *   \lambda_{v^\Psi} =& -v^i_g n_i, \\
     107             :  *   \lambda_{v^0} =& -v^i_g n_i, \\
     108             :  *   \lambda_{v^\pm} =& \pm 1 - v^i_g n_i,
     109             :  * \f}
     110             :  *
     111             :  * where \f$v_g^i\f$ is the mesh velocity and \f$n_i\f$ is the outward directed
     112             :  * unit normal covector to the interface. We have also defined
     113             :  *
     114             :  * \f{align}{
     115             :  *   \lambda^{\pm}_{\hat{\alpha}} =
     116             :  *    \left\{
     117             :  *        \begin{array}{ll}
     118             :  *          \lambda_{\hat{\alpha}} &
     119             :  *            \mathrm{if}\;\pm\lambda_{\hat{\alpha}}> 0 \\
     120             :  *          0 & \mathrm{otherwise}
     121             :  *        \end{array}\right.
     122             :  * \f}
     123             :  *
     124             :  * In the implementation we store the speeds in a rank-3 tensor with the zeroth
     125             :  * component being \f$\lambda_{v^\Psi}\f$, the first being \f$\lambda_{v^+}\f$
     126             :  * and the second being \f$\lambda_{v^-}\f$.
     127             :  *
     128             :  * Note that we have assumed \f$n_i^{\mathrm{ext}}\f$ points in the same
     129             :  * direction as \f$n_i^{\mathrm{int}}\f$, but in the code they point in opposite
     130             :  * directions. If \f$n_i^{\mathrm{ext}}\f$ points in the opposite direction the
     131             :  * external speeds have their sign flipped and the \f$\pm\f$ fields and their
     132             :  * speeds reverse roles (i.e. the \f$v^{\mathrm{ext},+}\f$ field is now flowing
     133             :  * into the element, while \f$v^{\mathrm{ext},-}\f$ flows out). In our
     134             :  * implementation this reversal actually cancels out, and we have the following
     135             :  * equations:
     136             :  *
     137             :  * \f{align}{
     138             :  *   D_{\Psi} &= -\lambda_{v^{\Psi}}^{\mathrm{ext},+}
     139             :  *                v^{\mathrm{ext},\Psi}
     140             :  *                - \lambda_{v^{\Psi}}^{\mathrm{int},-}
     141             :  *                v^{\mathrm{int},\Psi}, \\
     142             :  *   D_{\Pi}
     143             :  *              &= \frac{1}{2}\left(-\lambda_{v^+}^{\mathrm{ext},+}
     144             :  *                v^{\mathrm{ext},+} -
     145             :  *                \lambda_{v^-}^{\mathrm{ext},+}
     146             :  *                v^{\mathrm{ext},-}\right)
     147             :  *                - \lambda_{v^\Psi}^{\mathrm{ext},+}\gamma_2
     148             :  *                v^{\mathrm{ext},\Psi}
     149             :  *                \notag \\
     150             :  *              &-\frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
     151             :  *                v^{\mathrm{int},+} +
     152             :  *                \lambda_{v^-}^{\mathrm{int},-}
     153             :  *                v^{\mathrm{int},-}\right)
     154             :  *                - \lambda_{v^\Psi}^{\mathrm{int},-}\gamma_2
     155             :  *                v^{\mathrm{int},\Psi} , \\
     156             :  *   D_{\Phi_{i}}
     157             :  *              &= \frac{1}{2}\left(-\lambda_{v^+}^{\mathrm{ext},+}
     158             :  *                v^{\mathrm{ext},+}
     159             :  *                + \lambda_{v^-}^{\mathrm{ext},+}
     160             :  *                v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
     161             :  *                - \lambda_{v^0}^{\mathrm{ext},+}
     162             :  *                v^{\mathrm{ext},0}_{i}
     163             :  *                \notag \\
     164             :  *              &-
     165             :  *                \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
     166             :  *                v^{\mathrm{int},+}
     167             :  *                - \lambda_{v^-}^{\mathrm{int},-}
     168             :  *                v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
     169             :  *                - \lambda_{v^0}^{\mathrm{int},-}
     170             :  *                v^{\mathrm{int},0}_{i},
     171             :  * \f}
     172             :  */
     173             : template <size_t Dim>
     174           1 : class UpwindPenalty final : public BoundaryCorrection<Dim> {
     175             :  private:
     176           0 :   struct NormalTimesVPlus : db::SimpleTag {
     177           0 :     using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
     178             :   };
     179           0 :   struct NormalTimesVMinus : db::SimpleTag {
     180           0 :     using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
     181             :   };
     182           0 :   struct Gamma2VPsi : db::SimpleTag {
     183           0 :     using type = Scalar<DataVector>;
     184             :   };
     185           0 :   struct CharSpeedsTensor : db::SimpleTag {
     186           0 :     using type = tnsr::i<DataVector, 3, Frame::Inertial>;
     187             :   };
     188             : 
     189             :  public:
     190           0 :   using options = tmpl::list<>;
     191           0 :   static constexpr Options::String help = {
     192             :       "Computes the UpwindPenalty boundary correction term for the scalar wave "
     193             :       "system."};
     194             : 
     195           0 :   UpwindPenalty() = default;
     196           0 :   UpwindPenalty(const UpwindPenalty&) = default;
     197           0 :   UpwindPenalty& operator=(const UpwindPenalty&) = default;
     198           0 :   UpwindPenalty(UpwindPenalty&&) = default;
     199           0 :   UpwindPenalty& operator=(UpwindPenalty&&) = default;
     200           0 :   ~UpwindPenalty() override = default;
     201             : 
     202             :   /// \cond
     203             :   explicit UpwindPenalty(CkMigrateMessage* msg);
     204             :   using PUP::able::register_constructor;
     205             :   WRAPPED_PUPable_decl_template(UpwindPenalty);  // NOLINT
     206             :   /// \endcond
     207           0 :   void pup(PUP::er& p) override;  // NOLINT
     208             : 
     209           0 :   std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const override;
     210             : 
     211           0 :   using dg_package_field_tags =
     212             :       tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
     213             :                  NormalTimesVPlus, NormalTimesVMinus, Gamma2VPsi,
     214             :                  CharSpeedsTensor>;
     215           0 :   using dg_package_data_temporary_tags = tmpl::list<Tags::ConstraintGamma2>;
     216           0 :   using dg_package_data_volume_tags = tmpl::list<>;
     217           0 :   using dg_boundary_terms_volume_tags = tmpl::list<>;
     218             : 
     219           0 :   double dg_package_data(
     220             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_psi,
     221             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     222             :           packaged_char_speed_v_zero,
     223             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_plus,
     224             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_minus,
     225             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     226             :           packaged_char_speed_n_times_v_plus,
     227             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     228             :           packaged_char_speed_n_times_v_minus,
     229             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_gamma2_v_psi,
     230             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
     231             :           packaged_char_speeds,
     232             : 
     233             :       const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
     234             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
     235             : 
     236             :       const Scalar<DataVector>& constraint_gamma2,
     237             : 
     238             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
     239             :       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
     240             :       /*mesh_velocity*/,
     241             :       const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
     242             : 
     243           0 :   void dg_boundary_terms(
     244             :       gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
     245             :       gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
     246             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     247             :           phi_boundary_correction,
     248             : 
     249             :       const Scalar<DataVector>& char_speed_v_psi_int,
     250             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
     251             :       const Scalar<DataVector>& char_speed_v_plus_int,
     252             :       const Scalar<DataVector>& char_speed_v_minus_int,
     253             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     254             :           char_speed_normal_times_v_plus_int,
     255             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     256             :           char_speed_normal_times_v_minus_int,
     257             :       const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_int,
     258             :       const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_int,
     259             : 
     260             :       const Scalar<DataVector>& char_speed_v_psi_ext,
     261             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
     262             :       const Scalar<DataVector>& char_speed_v_plus_ext,
     263             :       const Scalar<DataVector>& char_speed_v_minus_ext,
     264             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     265             :           char_speed_minus_normal_times_v_plus_ext,
     266             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     267             :           char_speed_minus_normal_times_v_minus_ext,
     268             :       const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_ext,
     269             :       const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_ext,
     270             :       dg::Formulation /*dg_formulation*/) const;
     271             : };
     272             : }  // namespace ScalarWave::BoundaryCorrections