SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/ScalarWave - UpwindPenaltyCorrection.hpp Hit Total Coverage
Commit: b1342d46f40e2d46bbd11d0cef68fd973031a24b Lines: 1 20 5.0 %
Date: 2020-09-24 20:24:42
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 <cstddef>
       7             : #include <string>
       8             : 
       9             : #include "DataStructures/DataBox/Tag.hpp"
      10             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      11             : #include "DataStructures/Tensor/TypeAliases.hpp"
      12             : #include "Domain/FaceNormal.hpp"
      13             : #include "Evolution/Systems/ScalarWave/Tags.hpp"
      14             : #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
      15             : #include "Options/Options.hpp"
      16             : #include "Utilities/ForceInline.hpp"
      17             : #include "Utilities/ProtocolHelpers.hpp"
      18             : #include "Utilities/TMPL.hpp"
      19             : 
      20             : /// \cond
      21             : template <typename>
      22             : class Variables;
      23             : 
      24             : class DataVector;
      25             : 
      26             : namespace gsl {
      27             : template <typename T>
      28             : class not_null;
      29             : }  // namespace gsl
      30             : 
      31             : namespace Tags {
      32             : template <typename>
      33             : struct NormalDotFlux;
      34             : template <typename>
      35             : struct Normalized;
      36             : }  // namespace Tags
      37             : 
      38             : namespace PUP {
      39             : class er;
      40             : }  // namespace PUP
      41             : /// \endcond
      42             : 
      43             : // IWYU pragma: no_forward_declare Tensor
      44             : 
      45             : namespace ScalarWave {
      46             : /*!
      47             :  * \brief Computes the scalar wave upwind multipenalty boundary
      48             :  * correction.
      49             :  *
      50             :  * This implements the upwind multipenalty boundary correction term
      51             :  * \f$D_\alpha\f$. The general form is given by:
      52             :  *
      53             :  * \f{align*}{
      54             :  *   \label{eq:pnpm upwind boundary term characteristics}
      55             :  *   D_\beta =
      56             :  *   T_{\beta\hat{\beta}}^{\mathrm{ext}}
      57             :  *   \Lambda^{\mathrm{ext},-}_{\hat{\beta}\hat{\alpha}}
      58             :  *   v^{\mathrm{ext}}_{\hat{\alpha}}
      59             :  *   -T_{\beta\hat{\beta}}^{\mathrm{int}}
      60             :  *   \Lambda^{\mathrm{int},-}_{\hat{\beta}\hat{\alpha}}
      61             :  *   v^{\mathrm{int}}_{\hat{\alpha}},
      62             :  * \f}
      63             :  *
      64             :  * We denote the evolved fields by \f$u_{\alpha}\f$, the characteristic fields
      65             :  * by \f$v_{\hat{\alpha}}\f$, and implicitly sum over reapeated indices.
      66             :  * \f$T_{\alpha\hat{\alpha}}\f$ transforms characteristic fields to evolved
      67             :  * fields, while \f$\Lambda_{\hat{\alpha}\hat{\beta}}^-\f$ is a diagonal matrix
      68             :  * with only the negative characteristic speeds. The int and ext superscripts
      69             :  * denote quantities on the internal and external side of the mortar. Note that
      70             :  * Eq. (6.3) of \cite Teukolsky2015ega is not exactly what's implemented since
      71             :  * that boundary term does not consistently treat both sides of the interface on
      72             :  * the same footing.
      73             :  *
      74             :  * For the scalar wave system the correction is:
      75             :  *
      76             :  * \f{align}{
      77             :  *   D_{\Psi} &= \tilde{\lambda}_{v^{\Psi}}^{\mathrm{ext}}
      78             :  *                v^{\mathrm{ext},\Psi}
      79             :  *                - \tilde{\lambda}_{v^{\Psi}}^{\mathrm{int}}
      80             :  *                v^{\mathrm{int},\Psi}, \\
      81             :  *   D_{\Pi} &= \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{ext}}
      82             :  *             v^{\mathrm{ext},+} +
      83             :  *             \tilde{\lambda}_{v^-}^{\mathrm{ext}}
      84             :  *             v^{\mathrm{ext},-}\right)
      85             :  *             + \tilde{\lambda}_{v^\Psi}^\mathrm{ext}\gamma_2
      86             :  *             v^{\mathrm{ext},\Psi}
      87             :  *             \notag \\
      88             :  *           &-\frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
      89             :  *             v^{\mathrm{int},+} +
      90             :  *             \tilde{\lambda}_{v^-}^{\mathrm{int}}
      91             :  *             v^{\mathrm{int},-}\right)
      92             :  *             - \tilde{\lambda}_{v^\Psi}^\mathrm{int}\gamma_2
      93             :  *             v^{\mathrm{int},\Psi} , \\
      94             :  *   D_{\Phi_{i}}
      95             :  *              &= \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{ext}}
      96             :  *                v^{\mathrm{ext},+}
      97             :  *                - \tilde{\lambda}_{v^-}^{\mathrm{ext}}
      98             :  *                v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
      99             :  *                + \tilde{\lambda}_{v^0}^{\mathrm{ext}}
     100             :  *                v^{\mathrm{ext},0}_{i}
     101             :  *                \notag \\
     102             :  *              &-
     103             :  *                \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
     104             :  *                v^{\mathrm{int},+}
     105             :  *                - \tilde{\lambda}_{v^-}^{\mathrm{int}}
     106             :  *                v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
     107             :  *                - \tilde{\lambda}_{v^0}^{\mathrm{int}}
     108             :  *                v^{\mathrm{int},0}_{i},
     109             :  * \f}
     110             :  *
     111             :  * with characteristic fields
     112             :  *
     113             :  * \f{align}{
     114             :  *   v^{\Psi} &= \Psi, \\
     115             :  *   v^{0}_{i} &= (\delta^k_i-n^k n_i)\Phi_{k}, \\
     116             :  *   v^{\pm} &= \Pi\pm n^i\Phi_{i} -\gamma_2 \Psi,
     117             :  * \f}
     118             :  *
     119             :  * and characteristic speeds
     120             :  *
     121             :  * \f{align}{
     122             :  *   \lambda_{v^\Psi} =& -v^i_g n_i, \\
     123             :  *   \lambda_{v^0} =& -v^i_g n_i, \\
     124             :  *   \lambda_{v^\pm} =& \pm 1 - v^i_g n_i,
     125             :  * \f}
     126             :  *
     127             :  * where \f$v_g^i\f$ is the mesh velocity. We have also defined
     128             :  *
     129             :  * \f{align}{
     130             :  *   \tilde{\lambda}_{\hat{\alpha}} =
     131             :  *    \left\{
     132             :  *        \begin{array}{ll}
     133             :  *          \lambda_{\hat{\alpha}} &
     134             :  *            \mathrm{if}\;\lambda_{\hat{\alpha}}\le 0.0 \\
     135             :  *          0 & \mathrm{otherwise}
     136             :  *        \end{array}\right.
     137             :  * \f}
     138             :  *
     139             :  * Note that we have assumed \f$n_i^{\mathrm{ext}}\f$ points in the same
     140             :  * direction as \f$n_i^{\mathrm{int}}\f$. If \f$n_i^{\mathrm{ext}}\f$ points in
     141             :  * the opposite direction the external speeds have their sign flipped and the
     142             :  * \f$\pm\f$ fields and their speeds reverse roles. Specifically, in the code we
     143             :  * have:
     144             :  *
     145             :  * \f{align}{
     146             :  *   D_{\Psi} &= \bar{\lambda}_{v^{\Psi}}^{\mathrm{ext}}
     147             :  *                v^{\mathrm{ext},\Psi}
     148             :  *                - \tilde{\lambda}_{v^{\Psi}}^{\mathrm{int}}
     149             :  *                v^{\mathrm{int},g}, \\
     150             :  *   D_{\Pi}
     151             :  *              &= \frac{1}{2}\left(\bar{\lambda}_{v^+}^{\mathrm{ext}}
     152             :  *                v^{\mathrm{ext},+} +
     153             :  *                \bar{\lambda}_{v^-}^{\mathrm{ext}}
     154             :  *                v^{\mathrm{ext},-}\right)
     155             :  *                + \bar{\lambda}_{v^\Psi}^\mathrm{ext}\gamma_2
     156             :  *                v^{\mathrm{ext},\Psi}
     157             :  *                \notag \\
     158             :  *              &-\frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
     159             :  *                v^{\mathrm{int},+} +
     160             :  *                \tilde{\lambda}_{v^-}^{\mathrm{int}}
     161             :  *                v^{\mathrm{int},-}\right)
     162             :  *                - \tilde{\lambda}_{v^\Psi}^\mathrm{int}\gamma_2
     163             :  *                v^{\mathrm{int},\Psi} , \\
     164             :  *   D_{\Phi_{i}}
     165             :  *              &= \frac{1}{2}\left(\bar{\lambda}_{v^+}^{\mathrm{ext}}
     166             :  *                v^{\mathrm{ext},+}
     167             :  *                - \bar{\lambda}_{v^-}^{\mathrm{ext}}
     168             :  *                v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
     169             :  *                + \bar{\lambda}_{v^0}^{\mathrm{ext}}
     170             :  *                v^{\mathrm{ext},0}_{i}
     171             :  *                \notag \\
     172             :  *              &-
     173             :  *                \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
     174             :  *                v^{\mathrm{int},+}
     175             :  *                - \tilde{\lambda}_{v^-}^{\mathrm{int}}
     176             :  *                v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
     177             :  *                - \tilde{\lambda}_{v^0}^{\mathrm{int}}
     178             :  *                v^{\mathrm{int},0}_{i},
     179             :  * \f}
     180             :  *
     181             :  * where
     182             :  *
     183             :  * \f{align}{
     184             :  *   \bar{\lambda}_{\hat{\alpha}} =
     185             :  *    \left\{
     186             :  *        \begin{array}{ll}
     187             :  *          -\lambda_{\hat{\alpha}} &
     188             :  *            \mathrm{if}\;-\lambda_{\hat{\alpha}}\le 0.0 \\
     189             :  *          0 & \mathrm{otherwise}
     190             :  *        \end{array}\right.
     191             :  * \f}
     192             :  */
     193             : template <size_t Dim>
     194           1 : struct UpwindPenaltyCorrection : tt::ConformsTo<dg::protocols::NumericalFlux> {
     195             :  private:
     196           0 :   struct NormalTimesVPlus : db::SimpleTag {
     197           0 :     using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
     198             :   };
     199           0 :   struct NormalTimesVMinus : db::SimpleTag {
     200           0 :     using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
     201             :   };
     202           0 :   struct Gamma2VPsi : db::SimpleTag {
     203           0 :     using type = Scalar<DataVector>;
     204             :   };
     205           0 :   struct CharSpeedsTensor : db::SimpleTag {
     206           0 :     using type = tnsr::a<DataVector, 3, Frame::Inertial>;
     207             :   };
     208             : 
     209             :  public:
     210           0 :   using options = tmpl::list<>;
     211           0 :   static constexpr Options::String help = {
     212             :       "Computes the upwind penalty boundary correction for a scalar wave "
     213             :       "system. It requires no options."};
     214           0 :   static std::string name() noexcept { return "UpwindPenalty"; }
     215             : 
     216             :   // clang-tidy: non-const reference
     217           0 :   void pup(PUP::er& /*p*/) noexcept {}  // NOLINT
     218             : 
     219           0 :   using variables_tags = tmpl::list<Pi, Phi<Dim>, Psi>;
     220             : 
     221           0 :   using package_field_tags =
     222             :       tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
     223             :                  NormalTimesVPlus, NormalTimesVMinus, Gamma2VPsi,
     224             :                  CharSpeedsTensor>;
     225           0 :   using package_extra_tags = tmpl::list<>;
     226             : 
     227           0 :   using argument_tags =
     228             :       tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
     229             :                  Tags::CharacteristicSpeeds<Dim>, Tags::ConstraintGamma2,
     230             :                  ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>>;
     231             : 
     232           0 :   void package_data(
     233             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_psi,
     234             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     235             :           packaged_char_speed_v_zero,
     236             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_plus,
     237             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_minus,
     238             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     239             :           packaged_char_speed_n_times_v_plus,
     240             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     241             :           packaged_char_speed_n_times_v_minus,
     242             :       gsl::not_null<Scalar<DataVector>*> packaged_char_speed_gamma2_v_psi,
     243             :       gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
     244             :           packaged_char_speeds,
     245             : 
     246             :       const Scalar<DataVector>& v_psi,
     247             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero,
     248             :       const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
     249             :       const std::array<DataVector, 4>& char_speeds,
     250             :       const Scalar<DataVector>& constraint_gamma2,
     251             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
     252             :       const noexcept;
     253             : 
     254           0 :   void operator()(
     255             :       gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
     256             :       gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
     257             :           phi_boundary_correction,
     258             :       gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
     259             : 
     260             :       const Scalar<DataVector>& char_speed_v_psi_int,
     261             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
     262             :       const Scalar<DataVector>& char_speed_v_plus_int,
     263             :       const Scalar<DataVector>& char_speed_v_minus_int,
     264             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     265             :           char_speed_normal_times_v_plus_int,
     266             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     267             :           char_speed_normal_times_v_minus_int,
     268             :       const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_int,
     269             :       const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
     270             : 
     271             :       const Scalar<DataVector>& char_speed_v_psi_ext,
     272             :       const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
     273             :       const Scalar<DataVector>& char_speed_v_plus_ext,
     274             :       const Scalar<DataVector>& char_speed_v_minus_ext,
     275             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     276             :           char_speed_minus_normal_times_v_plus_ext,
     277             :       const tnsr::i<DataVector, Dim, Frame::Inertial>&
     278             :           char_speed_minus_normal_times_v_minus_ext,
     279             :       const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_ext,
     280             :       const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext)
     281             :       const noexcept;
     282             : };
     283             : }  // namespace ScalarWave

Generated by: LCOV version 1.14