UpwindPenaltyCorrection.hpp
1 // 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"
12 #include "Domain/FaceNormal.hpp"
14 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
15 #include "Options/Options.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 struct UpwindPenaltyCorrection : tt::ConformsTo<dg::protocols::NumericalFlux> {
195  private:
196  struct NormalTimesVPlus : db::SimpleTag {
197  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
198  };
199  struct NormalTimesVMinus : db::SimpleTag {
200  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
201  };
202  struct Gamma2VPsi : db::SimpleTag {
203  using type = Scalar<DataVector>;
204  };
205  struct CharSpeedsTensor : db::SimpleTag {
206  using type = tnsr::a<DataVector, 3, Frame::Inertial>;
207  };
208 
209  public:
210  using options = tmpl::list<>;
211  static constexpr OptionString help = {
212  "Computes the upwind penalty boundary correction for a scalar wave "
213  "system. It requires no options."};
214  static std::string name() noexcept { return "UpwindPenalty"; }
215 
216  // clang-tidy: non-const reference
217  void pup(PUP::er& /*p*/) noexcept {} // NOLINT
218 
219  using variables_tags = tmpl::list<Pi, Phi<Dim>, Psi>;
220 
221  using package_field_tags =
222  tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
223  NormalTimesVPlus, NormalTimesVMinus, Gamma2VPsi,
224  CharSpeedsTensor>;
225  using package_extra_tags = tmpl::list<>;
226 
227  using argument_tags =
228  tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
231 
232  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  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
ScalarWave::Tags::ConstraintGamma2
Definition: Tags.hpp:37
std::string
ScalarWave::Tags::VPlus
Definition: Tags.hpp:78
Options.hpp
FaceNormal.hpp
Tags.hpp
db::SimpleTag
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
Tags::Normalized
Definition: Magnitude.hpp:145
cstddef
std::array
ScalarWave::UpwindPenaltyCorrection
Computes the scalar wave upwind multipenalty boundary correction.
Definition: UpwindPenaltyCorrection.hpp:194
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
ScalarWave::Tags::CharacteristicSpeeds
Definition: Tags.hpp:89
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
TypeAliases.hpp
ForceInline.hpp
ScalarWave
Items related to evolving the scalar wave equation.
Definition: Characteristics.cpp:16
ScalarWave::Psi
Definition: Tags.hpp:20
ScalarWave::Tags::VMinus
Definition: Tags.hpp:82
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
gsl
Implementations from the Guideline Support Library.
Definition: Gsl.hpp:80
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string