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