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/BoundaryCorrection.hpp"
13 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
14 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/ConstraintDampingTags.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 evolution::BoundaryCorrection {
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> 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::Tags::ConstraintGamma1, ::gh::Tags::ConstraintGamma2,
222 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, Dim>>;
223 0 : using dg_package_data_primitive_tags = tmpl::list<>;
224 0 : using dg_package_data_volume_tags = tmpl::list<>;
225 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
226 :
227 0 : double dg_package_data(
228 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
229 : packaged_char_speed_v_spacetime_metric,
230 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
231 : packaged_char_speed_v_zero,
232 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
233 : packaged_char_speed_v_plus,
234 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
235 : packaged_char_speed_v_minus,
236 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
237 : packaged_char_speed_n_times_v_plus,
238 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
239 : packaged_char_speed_n_times_v_minus,
240 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
241 : packaged_char_speed_gamma2_v_spacetime_metric,
242 : gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
243 : packaged_char_speeds,
244 :
245 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
246 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
247 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
248 :
249 : const Scalar<DataVector>& constraint_gamma1,
250 : const Scalar<DataVector>& constraint_gamma2,
251 : const Scalar<DataVector>& lapse,
252 : const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
253 :
254 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
255 : const tnsr::I<DataVector, Dim, Frame::Inertial>& normal_vector,
256 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
257 : /*mesh_velocity*/,
258 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
259 :
260 0 : void dg_boundary_terms(
261 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
262 : boundary_correction_spacetime_metric,
263 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
264 : boundary_correction_pi,
265 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
266 : boundary_correction_phi,
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_normal_times_v_plus_ext,
288 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
289 : char_speed_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 : dg::Formulation /*dg_formulation*/) const;
294 : };
295 :
296 : template <size_t Dim>
297 0 : bool operator==(const UpwindPenalty<Dim>& lhs, const UpwindPenalty<Dim>& rhs);
298 : template <size_t Dim>
299 0 : bool operator!=(const UpwindPenalty<Dim>& lhs, const UpwindPenalty<Dim>& rhs);
300 : } // namespace gh::BoundaryCorrections
|
