UpwindPenaltyCorrection.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <array>
7 #include <cstddef>
8 #include <string>
9
10 #include "DataStructures/DataBox/Tag.hpp"
11 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
13 #include "Domain/FaceNormal.hpp"
14 #include "Evolution/Systems/GeneralizedHarmonic/ConstraintDamping/Tags.hpp"
15 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
16 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
17 #include "Options/Options.hpp"
18 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
19 #include "Utilities/Gsl.hpp"
20 #include "Utilities/ProtocolHelpers.hpp"
21 #include "Utilities/TMPL.hpp"
22
23 /// \cond
24 class DataVector;
25 namespace PUP {
26 class er;
27 } // namespace PUP
28 /// \endcond
29
30 namespace GeneralizedHarmonic {
31 /*!
32  * \brief Computes the generalized harmonic upwind multipenalty boundary
33  * correction.
34  *
35  * This implements the upwind multipenalty boundary correction term
36  * \f$D_\beta\f$. We index the evolved variables using unhatted Greek letters,
37  * and the characteristic variables using hatted Greek letters. The general form
38  * of \f$D_\beta\f$ is given by:
39  *
40  * \f{align*}{
41  * \label{eq:pnpm upwind boundary term characteristics}
42  * D_\beta =
43  * T_{\beta\hat{\beta}}^{\mathrm{ext}}
44  * \Lambda^{\mathrm{ext},-}_{\hat{\beta}\hat{\alpha}}
45  * v^{\mathrm{ext}}_{\hat{\alpha}}
46  * -T_{\beta\hat{\beta}}^{\mathrm{int}}
47  * \Lambda^{\mathrm{int},-}_{\hat{\beta}\hat{\alpha}}
48  * v^{\mathrm{int}}_{\hat{\alpha}}.
49  * \f}
50  *
51  * Note that Eq. (6.3) of \cite Teukolsky2015ega is not exactly what's
52  * implemented since the boundary term given by Eq. (6.3) does not consistently
53  * treat both sides of the interface on the same footing.
54  *
55  * For the first-order generalized harmonic system the correction is:
56  *
57  * \f{align}{
58  * \label{eq:pnpm upwind gh metric}
59  * D_{g_{ab}} &= \tilde{\lambda}_{v^{g}}^{\mathrm{ext}}
60  * v^{\mathrm{ext},g}_{ab}
61  * - \tilde{\lambda}_{v^{g}}^{\mathrm{int}}
62  * v^{\mathrm{int},g}_{ab}, \\
63  * \label{eq:pnpm upwind gh pi}
64  * D_{\Pi_{ab}}
65  * &= \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{ext}}
66  * v^{\mathrm{ext},+}_{ab} +
67  * \tilde{\lambda}_{v^-}^{\mathrm{ext}}
68  * v^{\mathrm{ext},-}_{ab}\right)
69  * + \tilde{\lambda}_{v^g}^\mathrm{ext}\gamma_2
70  * v^{\mathrm{ext},g}_{ab}
71  * \notag \\
72  * &-\frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
73  * v^{\mathrm{int},+}_{ab} +
74  * \tilde{\lambda}_{v^-}^{\mathrm{int}}
75  * v^{\mathrm{int},-}_{ab}\right)
76  * - \tilde{\lambda}_{v^g}^\mathrm{int}\gamma_2
77  * v^{\mathrm{int},g}_{ab} , \\
78  * \label{eq:pnpm upwind gh phi}
79  * D_{\Phi_{iab}}
80  * &= \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{ext}}
81  * v^{\mathrm{ext},+}_{ab}
82  * - \tilde{\lambda}_{v^-}^{\mathrm{ext}}
83  * v^{\mathrm{ext},-}_{ab}\right)n_i^{\mathrm{ext}}
84  * + \tilde{\lambda}_{v^0}^{\mathrm{ext}}
85  * v^{\mathrm{ext},0}_{iab}
86  * \notag \\
87  * &-
88  * \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
89  * v^{\mathrm{int},+}_{ab}
90  * - \tilde{\lambda}_{v^-}^{\mathrm{int}}
91  * v^{\mathrm{int},-}_{ab}\right)n_i^{\mathrm{int}}
92  * - \tilde{\lambda}_{v^0}^{\mathrm{int}}
93  * v^{\mathrm{int},0}_{iab},
94  * \f}
95  *
96  * with characteristic fields
97  *
98  * \f{align}{
99  * \label{eq:pnpm gh v g}
100  * v^{g}_{ab} &= g_{ab}, \\
101  * \label{eq:pnpm gh v 0}
102  * v^{0}_{iab} &= (\delta^k_i-n^k n_i)\Phi_{kab}, \\
103  * \label{eq:pnpm gh v pm}
104  * v^{\pm}_{ab} &= \Pi_{ab}\pm n^i\Phi_{iab} -\gamma_2 g_{ab},
105  * \f}
106  *
107  * and characteristic speeds
108  *
109  * \f{align}{
110  * \lambda_{v^g} =& -(1+\gamma_1)\beta^i n_i -v^i_g n_i, \\
111  * \lambda_{v^0} =& -\beta^i n_i -v^i_g n_i, \\
112  * \lambda_{v^\pm} =& \pm \alpha - \beta^i n_i - v^i_g n_i,
113  * \f}
114  *
115  * where \f$v_g^i\f$ is the mesh velocity. We have also defined
116  *
117  * \f{align}{
118  * \tilde{\lambda}_{\hat{\alpha}} =
119  * \left\{
120  * \begin{array}{ll}
121  * \lambda_{\hat{\alpha}} &
122  * \mathrm{if}\;\lambda_{\hat{\alpha}}\le 0.0 \\
123  * 0 & \mathrm{otherwise}
124  * \end{array}\right.
125  * \f}
126  *
127  * Note that we have assumed \f$n_i^{\mathrm{ext}}\f$ points in the same
128  * direction as \f$n_i^{\mathrm{int}}\f$, but in the code we have them point in
129  * opposite directions. If \f$n_i^{\mathrm{ext}}\f$ points in the opposite
130  * direction the external speeds have their sign flipped and the \f$\pm\f$
131  * fields and speeds reverse roles. Specifically, in the code we have:
132  *
133  * \f{align}{
134  * \label{eq:pnpm upwind gh metric code}
135  * D_{g_{ab}} &= \bar{\lambda}_{v^{g}}^{\mathrm{ext}}
136  * v^{\mathrm{ext},g}_{ab}
137  * - \tilde{\lambda}_{v^{g}}^{\mathrm{int}}
138  * v^{\mathrm{int},g}_{ab}, \\
139  * \label{eq:pnpm upwind gh pi code}
140  * D_{\Pi_{ab}}
141  * &= \frac{1}{2}\left(\bar{\lambda}_{v^+}^{\mathrm{ext}}
142  * v^{\mathrm{ext},+}_{ab} +
143  * \bar{\lambda}_{v^-}^{\mathrm{ext}}
144  * v^{\mathrm{ext},-}_{ab}\right)
145  * + \bar{\lambda}_{v^g}^\mathrm{ext}\gamma_2
146  * v^{\mathrm{ext},g}_{ab}
147  * \notag \\
148  * &-\frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
149  * v^{\mathrm{int},+}_{ab} +
150  * \tilde{\lambda}_{v^-}^{\mathrm{int}}
151  * v^{\mathrm{int},-}_{ab}\right)
152  * - \tilde{\lambda}_{v^g}^\mathrm{int}\gamma_2
153  * v^{\mathrm{int},g}_{ab} , \\
154  * \label{eq:pnpm upwind gh phi code}
155  * D_{\Phi_{iab}}
156  * &= \frac{1}{2}\left(\bar{\lambda}_{v^+}^{\mathrm{ext}}
157  * v^{\mathrm{ext},+}_{ab}
158  * - \bar{\lambda}_{v^-}^{\mathrm{ext}}
159  * v^{\mathrm{ext},-}_{ab}\right)n_i^{\mathrm{ext}}
160  * + \bar{\lambda}_{v^0}^{\mathrm{ext}}
161  * v^{\mathrm{ext},0}_{iab}
162  * \notag \\
163  * &-
164  * \frac{1}{2}\left(\tilde{\lambda}_{v^+}^{\mathrm{int}}
165  * v^{\mathrm{int},+}_{ab}
166  * - \tilde{\lambda}_{v^-}^{\mathrm{int}}
167  * v^{\mathrm{int},-}_{ab}\right)n_i^{\mathrm{int}}
168  * - \tilde{\lambda}_{v^0}^{\mathrm{int}}
169  * v^{\mathrm{int},0}_{iab},
170  * \f}
171  *
172  * where
173  *
174  * \f{align}{
175  * \bar{\lambda}_{\hat{\alpha}} =
176  * \left\{
177  * \begin{array}{ll}
178  * -\lambda_{\hat{\alpha}} &
179  * \mathrm{if}\;-\lambda_{\hat{\alpha}}\le 0.0 \\
180  * 0 & \mathrm{otherwise}
181  * \end{array}\right.
182  * \f}
183  */
184 template <size_t Dim>
185 struct UpwindPenaltyCorrection : tt::ConformsTo<dg::protocols::NumericalFlux> {
187  using type = tnsr::iaa<DataVector, Dim, Frame::Inertial>;
188  };
190  using type = tnsr::iaa<DataVector, Dim, Frame::Inertial>;
191  };
193  using type = tnsr::aa<DataVector, Dim, Frame::Inertial>;
194  };
196  // use a tensor for sending the char speeds so they all go into a Variables
197  // to reduce memory allocations.
198  using type = tnsr::a<DataVector, 3, Frame::Inertial>;
199  };
200
201  public:
202  using options = tmpl::list<>;
203  static constexpr Options::String help = {
204  "Computes the upwind penalty flux for the first-order generalized "
205  "harmonic system. It requires no options."};
206  static std::string name() noexcept { return "UpwindPenalty"; }
207
208  // clang-tidy: non-const reference
209  void pup(PUP::er& /*p*/) noexcept {} // NOLINT
210
211  using variables_tags =
212  tmpl::list<gr::Tags::SpacetimeMetric<Dim, Frame::Inertial, DataVector>,
213  Tags::Pi<Dim, Frame::Inertial>,
214  Tags::Phi<Dim, Frame::Inertial>>;
215
216  using package_field_tags =
217  tmpl::list<Tags::VSpacetimeMetric<Dim, Frame::Inertial>,
218  Tags::VZero<Dim, Frame::Inertial>,
219  Tags::VPlus<Dim, Frame::Inertial>,
220  Tags::VMinus<Dim, Frame::Inertial>, CharSpeedNormalTimesVPlus,
221  CharSpeedNormalTimesVMinus, CharSpeedGamma2VSpacetimeMetric,
222  TensorCharSpeeds>;
223  using package_extra_tags = tmpl::list<>;
224
225  using argument_tags = tmpl::list<
226  Tags::VSpacetimeMetric<Dim, Frame::Inertial>,
227  Tags::VZero<Dim, Frame::Inertial>, Tags::VPlus<Dim, Frame::Inertial>,
228  Tags::VMinus<Dim, Frame::Inertial>,
229  Tags::CharacteristicSpeeds<Dim, Frame::Inertial>,
232
233  void package_data(
234  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
235  packaged_char_speed_v_spacetime_metric,
236  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
237  packaged_char_speed_v_zero,
238  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
239  packaged_char_speed_v_plus,
240  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
241  packaged_char_speed_v_minus,
242  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
243  packaged_char_speed_n_times_v_plus,
244  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
245  packaged_char_speed_n_times_v_minus,
246  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
247  packaged_char_speed_gamma2_v_spacetime_metric,
248  gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
249  packaged_char_speeds,
250
251  const tnsr::aa<DataVector, Dim, Frame::Inertial>& v_spacetime_metric,
252  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& v_zero,
253  const tnsr::aa<DataVector, Dim, Frame::Inertial>& v_plus,
254  const tnsr::aa<DataVector, Dim, Frame::Inertial>& v_minus,
255  const std::array<DataVector, 4>& char_speeds,
256  const Scalar<DataVector>& constraint_gamma2,
257  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
258  const noexcept;
259
260  void operator()(
261  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
262  spacetime_metric_boundary_correction,
263  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
264  pi_boundary_correction,
265  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
266  phi_boundary_correction,
267
268  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
269  char_speed_v_spacetime_metric_int,
270  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
271  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_plus_int,
272  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_minus_int,
273  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
274  char_speed_normal_times_v_plus_int,
275  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
276  char_speed_normal_times_v_minus_int,
277  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
278  char_speed_constraint_gamma2_v_spacetime_metric_int,
279  const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
280
281  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
282  char_speed_v_spacetime_metric_ext,
283  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
284  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_plus_ext,
285  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_minus_ext,
286  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
287  char_speed_minus_normal_times_v_plus_ext,
288  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
289  char_speed_minus_normal_times_v_minus_ext,
290  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
291  char_speed_constraint_gamma2_v_spacetime_metric_ext,
292  const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext)
293  const noexcept;
294 };
295 } // namespace GeneralizedHarmonic
GeneralizedHarmonic::UpwindPenaltyCorrection::TensorCharSpeeds
Definition: UpwindPenaltyCorrection.hpp:195
GeneralizedHarmonic::UpwindPenaltyCorrection::CharSpeedGamma2VSpacetimeMetric
Definition: UpwindPenaltyCorrection.hpp:192
std::string
Options.hpp
FaceNormal.hpp
GeneralizedHarmonic
Items related to evolving the first-order generalized harmonic system.
Definition: Bjorhus.hpp:38
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::Normalized
Definition: Magnitude.hpp:137
cstddef
array
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
GeneralizedHarmonic::ConstraintDamping::Tags::ConstraintGamma2
Constraint dammping parameter for the generalized harmonic system (cf. ).
Definition: Tags.hpp:70
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
GeneralizedHarmonic::UpwindPenaltyCorrection::CharSpeedNormalTimesVPlus
Definition: UpwindPenaltyCorrection.hpp:186
GeneralizedHarmonic::UpwindPenaltyCorrection::CharSpeedNormalTimesVMinus
Definition: UpwindPenaltyCorrection.hpp:189
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
GeneralizedHarmonic::UpwindPenaltyCorrection
Computes the generalized harmonic upwind multipenalty boundary correction.
Definition: UpwindPenaltyCorrection.hpp:185
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
string