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