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 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
124 :
125 0 : double dg_package_data(
126 : gsl::not_null<Scalar<DataVector>*> packaged_v_psi,
127 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> packaged_v_zero,
128 : gsl::not_null<Scalar<DataVector>*> packaged_v_plus,
129 : gsl::not_null<Scalar<DataVector>*> packaged_v_minus,
130 : gsl::not_null<Scalar<DataVector>*> packaged_gamma2,
131 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
132 : packaged_interface_unit_normal,
133 : gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
134 : packaged_char_speeds,
135 :
136 : const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
137 : const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
138 :
139 : const Scalar<DataVector>& lapse,
140 : const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
141 : const Scalar<DataVector>& constraint_gamma1,
142 : const Scalar<DataVector>& constraint_gamma2,
143 :
144 : const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
145 : const tnsr::I<DataVector, Dim, Frame::Inertial>&
146 : interface_unit_normal_vector,
147 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
148 : /*mesh_velocity*/,
149 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
150 :
151 0 : void dg_boundary_terms(
152 : gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
153 : gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
154 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
155 : phi_boundary_correction,
156 :
157 : const Scalar<DataVector>& v_psi_int,
158 : const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero_int,
159 : const Scalar<DataVector>& v_plus_int,
160 : const Scalar<DataVector>& v_minus_int,
161 : const Scalar<DataVector>& gamma2_int,
162 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
163 : interface_unit_normal_int,
164 : const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
165 :
166 : const Scalar<DataVector>& v_psi_ext,
167 : const tnsr::i<DataVector, Dim, Frame::Inertial>& v_zero_ext,
168 : const Scalar<DataVector>& v_plus_ext,
169 : const Scalar<DataVector>& v_minus_ext,
170 : const Scalar<DataVector>& gamma2_ext,
171 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
172 : interface_unit_normal_ext,
173 : const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext,
174 : dg::Formulation /*dg_formulation*/) const;
175 : };
176 :
177 : template <size_t Dim>
178 0 : bool operator==(const UpwindPenalty<Dim>& lhs, const UpwindPenalty<Dim>& rhs);
179 : template <size_t Dim>
180 0 : bool operator!=(const UpwindPenalty<Dim>& lhs, const UpwindPenalty<Dim>& rhs);
181 : } // namespace CurvedScalarWave::BoundaryCorrections
|
