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/BoundaryCorrection.hpp"
12 : #include "Evolution/Systems/ScalarWave/Tags.hpp"
13 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
14 : #include "Options/String.hpp"
15 : #include "Utilities/Serialization/CharmPupable.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : namespace gsl {
21 : template <typename T>
22 : class not_null;
23 : } // namespace gsl
24 : namespace PUP {
25 : class er;
26 : } // namespace PUP
27 : /// \endcond
28 :
29 : namespace ScalarWave::BoundaryCorrections {
30 : /*!
31 : * \brief Computes the scalar wave upwind multipenalty boundary
32 : * correction.
33 : *
34 : * This implements the upwind multipenalty boundary correction term
35 : * \f$D_\beta\f$. The general form is given by:
36 : *
37 : * \f{align*}{
38 : * D_\beta =
39 : * T_{\beta\hat{\beta}}^{\mathrm{ext}}
40 : * \Lambda^{\mathrm{ext},-}_{\hat{\beta}\hat{\alpha}}
41 : * v^{\mathrm{ext}}_{\hat{\alpha}}
42 : * -T_{\beta\hat{\beta}}^{\mathrm{int}}
43 : * \Lambda^{\mathrm{int},-}_{\hat{\beta}\hat{\alpha}}
44 : * v^{\mathrm{int}}_{\hat{\alpha}}.
45 : * \f}
46 : *
47 : * We denote the evolved fields by \f$u_{\alpha}\f$, the characteristic fields
48 : * by \f$v_{\hat{\alpha}}\f$, and implicitly sum over reapeated indices.
49 : * \f$T_{\beta\hat{\beta}}\f$ transforms characteristic fields to evolved
50 : * fields, while \f$\Lambda_{\hat{\beta}\hat{\alpha}}^-\f$ is a diagonal matrix
51 : * with only the negative characteristic speeds, and has zeros on the diagonal
52 : * for all other entries. The int and ext superscripts denote quantities on the
53 : * internal and external side of the mortar. Note that Eq. (6.3) of
54 : * \cite Teukolsky2015ega is not exactly what's implemented since that boundary
55 : * term does not consistently treat both sides of the interface on the same
56 : * footing.
57 : *
58 : * For the scalar wave system the correction is:
59 : *
60 : * \f{align}{
61 : * D_{\Psi} &= \lambda_{v^{\Psi}}^{\mathrm{ext},-}
62 : * v^{\mathrm{ext},\Psi}
63 : * - \lambda_{v^{\Psi}}^{\mathrm{int},-}
64 : * v^{\mathrm{int},\Psi}, \\
65 : * D_{\Pi} &= \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{ext},-}
66 : * v^{\mathrm{ext},+} +
67 : * \lambda_{v^-}^{\mathrm{ext},-}
68 : * v^{\mathrm{ext},-}\right)
69 : * + \lambda_{v^\Psi}^{\mathrm{ext},-}\gamma_2
70 : * v^{\mathrm{ext},\Psi}
71 : * \notag \\
72 : * &-\frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
73 : * v^{\mathrm{int},+} +
74 : * \lambda_{v^-}^{\mathrm{int},-}
75 : * v^{\mathrm{int},-}\right)
76 : * - \lambda_{v^\Psi}^{\mathrm{int},-}\gamma_2
77 : * v^{\mathrm{int},\Psi} , \\
78 : * D_{\Phi_{i}}
79 : * &= \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{ext},-}
80 : * v^{\mathrm{ext},+}
81 : * - \lambda_{v^-}^{\mathrm{ext},-}
82 : * v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
83 : * + \lambda_{v^0}^{\mathrm{ext},-}
84 : * v^{\mathrm{ext},0}_{i}
85 : * \notag \\
86 : * &-
87 : * \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
88 : * v^{\mathrm{int},+}
89 : * - \lambda_{v^-}^{\mathrm{int},-}
90 : * v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
91 : * - \lambda_{v^0}^{\mathrm{int},-}
92 : * v^{\mathrm{int},0}_{i},
93 : * \f}
94 : *
95 : * with characteristic fields
96 : *
97 : * \f{align}{
98 : * v^{\Psi} &= \Psi, \\
99 : * v^{0}_{i} &= (\delta^k_i-n^k n_i)\Phi_{k}, \\
100 : * v^{\pm} &= \Pi\pm n^i\Phi_{i} -\gamma_2 \Psi,
101 : * \f}
102 : *
103 : * and characteristic speeds
104 : *
105 : * \f{align}{
106 : * \lambda_{v^\Psi} =& -v^i_g n_i, \\
107 : * \lambda_{v^0} =& -v^i_g n_i, \\
108 : * \lambda_{v^\pm} =& \pm 1 - v^i_g n_i,
109 : * \f}
110 : *
111 : * where \f$v_g^i\f$ is the mesh velocity and \f$n_i\f$ is the outward directed
112 : * unit normal covector to the interface. We have also defined
113 : *
114 : * \f{align}{
115 : * \lambda^{\pm}_{\hat{\alpha}} =
116 : * \left\{
117 : * \begin{array}{ll}
118 : * \lambda_{\hat{\alpha}} &
119 : * \mathrm{if}\;\pm\lambda_{\hat{\alpha}}> 0 \\
120 : * 0 & \mathrm{otherwise}
121 : * \end{array}\right.
122 : * \f}
123 : *
124 : * In the implementation we store the speeds in a rank-3 tensor with the zeroth
125 : * component being \f$\lambda_{v^\Psi}\f$, the first being \f$\lambda_{v^+}\f$
126 : * and the second being \f$\lambda_{v^-}\f$.
127 : *
128 : * Note that we have assumed \f$n_i^{\mathrm{ext}}\f$ points in the same
129 : * direction as \f$n_i^{\mathrm{int}}\f$, but in the code they point in opposite
130 : * directions. If \f$n_i^{\mathrm{ext}}\f$ points in the opposite direction the
131 : * external speeds have their sign flipped and the \f$\pm\f$ fields and their
132 : * speeds reverse roles (i.e. the \f$v^{\mathrm{ext},+}\f$ field is now flowing
133 : * into the element, while \f$v^{\mathrm{ext},-}\f$ flows out). In our
134 : * implementation this reversal actually cancels out, and we have the following
135 : * equations:
136 : *
137 : * \f{align}{
138 : * D_{\Psi} &= -\lambda_{v^{\Psi}}^{\mathrm{ext},+}
139 : * v^{\mathrm{ext},\Psi}
140 : * - \lambda_{v^{\Psi}}^{\mathrm{int},-}
141 : * v^{\mathrm{int},\Psi}, \\
142 : * D_{\Pi}
143 : * &= \frac{1}{2}\left(-\lambda_{v^+}^{\mathrm{ext},+}
144 : * v^{\mathrm{ext},+} -
145 : * \lambda_{v^-}^{\mathrm{ext},+}
146 : * v^{\mathrm{ext},-}\right)
147 : * - \lambda_{v^\Psi}^{\mathrm{ext},+}\gamma_2
148 : * v^{\mathrm{ext},\Psi}
149 : * \notag \\
150 : * &-\frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
151 : * v^{\mathrm{int},+} +
152 : * \lambda_{v^-}^{\mathrm{int},-}
153 : * v^{\mathrm{int},-}\right)
154 : * - \lambda_{v^\Psi}^{\mathrm{int},-}\gamma_2
155 : * v^{\mathrm{int},\Psi} , \\
156 : * D_{\Phi_{i}}
157 : * &= \frac{1}{2}\left(-\lambda_{v^+}^{\mathrm{ext},+}
158 : * v^{\mathrm{ext},+}
159 : * + \lambda_{v^-}^{\mathrm{ext},+}
160 : * v^{\mathrm{ext},-}\right)n_i^{\mathrm{ext}}
161 : * - \lambda_{v^0}^{\mathrm{ext},+}
162 : * v^{\mathrm{ext},0}_{i}
163 : * \notag \\
164 : * &-
165 : * \frac{1}{2}\left(\lambda_{v^+}^{\mathrm{int},-}
166 : * v^{\mathrm{int},+}
167 : * - \lambda_{v^-}^{\mathrm{int},-}
168 : * v^{\mathrm{int},-}\right)n_i^{\mathrm{int}}
169 : * - \lambda_{v^0}^{\mathrm{int},-}
170 : * v^{\mathrm{int},0}_{i},
171 : * \f}
172 : */
173 : template <size_t Dim>
174 1 : class UpwindPenalty final : public evolution::BoundaryCorrection {
175 : private:
176 0 : struct NormalTimesVPlus : db::SimpleTag {
177 0 : using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
178 : };
179 0 : struct NormalTimesVMinus : db::SimpleTag {
180 0 : using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
181 : };
182 0 : struct Gamma2VPsi : db::SimpleTag {
183 0 : using type = Scalar<DataVector>;
184 : };
185 0 : struct CharSpeedsTensor : db::SimpleTag {
186 0 : using type = tnsr::i<DataVector, 3, Frame::Inertial>;
187 : };
188 :
189 : public:
190 0 : using options = tmpl::list<>;
191 0 : static constexpr Options::String help = {
192 : "Computes the UpwindPenalty boundary correction term for the scalar wave "
193 : "system."};
194 :
195 0 : UpwindPenalty() = default;
196 0 : UpwindPenalty(const UpwindPenalty&) = default;
197 0 : UpwindPenalty& operator=(const UpwindPenalty&) = default;
198 0 : UpwindPenalty(UpwindPenalty&&) = default;
199 0 : UpwindPenalty& operator=(UpwindPenalty&&) = default;
200 0 : ~UpwindPenalty() override = default;
201 :
202 : /// \cond
203 : explicit UpwindPenalty(CkMigrateMessage* msg);
204 : using PUP::able::register_constructor;
205 : WRAPPED_PUPable_decl_template(UpwindPenalty); // NOLINT
206 : /// \endcond
207 0 : void pup(PUP::er& p) override; // NOLINT
208 :
209 0 : std::unique_ptr<BoundaryCorrection> get_clone() const override;
210 :
211 0 : using dg_package_field_tags =
212 : tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
213 : NormalTimesVPlus, NormalTimesVMinus, Gamma2VPsi,
214 : CharSpeedsTensor>;
215 0 : using dg_package_data_temporary_tags = tmpl::list<Tags::ConstraintGamma2>;
216 0 : using dg_package_data_volume_tags = tmpl::list<>;
217 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
218 :
219 0 : double dg_package_data(
220 : gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_psi,
221 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
222 : packaged_char_speed_v_zero,
223 : gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_plus,
224 : gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_minus,
225 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
226 : packaged_char_speed_n_times_v_plus,
227 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
228 : packaged_char_speed_n_times_v_minus,
229 : gsl::not_null<Scalar<DataVector>*> packaged_char_speed_gamma2_v_psi,
230 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
231 : packaged_char_speeds,
232 :
233 : const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
234 : const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
235 :
236 : const Scalar<DataVector>& constraint_gamma2,
237 :
238 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
239 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
240 : /*mesh_velocity*/,
241 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
242 :
243 0 : void dg_boundary_terms(
244 : gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
245 : gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
246 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
247 : phi_boundary_correction,
248 :
249 : const Scalar<DataVector>& char_speed_v_psi_int,
250 : const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
251 : const Scalar<DataVector>& char_speed_v_plus_int,
252 : const Scalar<DataVector>& char_speed_v_minus_int,
253 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
254 : char_speed_normal_times_v_plus_int,
255 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
256 : char_speed_normal_times_v_minus_int,
257 : const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_int,
258 : const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_int,
259 :
260 : const Scalar<DataVector>& char_speed_v_psi_ext,
261 : const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
262 : const Scalar<DataVector>& char_speed_v_plus_ext,
263 : const Scalar<DataVector>& char_speed_v_minus_ext,
264 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
265 : char_speed_minus_normal_times_v_plus_ext,
266 : const tnsr::i<DataVector, Dim, Frame::Inertial>&
267 : char_speed_minus_normal_times_v_minus_ext,
268 : const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_ext,
269 : const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_ext,
270 : dg::Formulation /*dg_formulation*/) const;
271 : };
272 : } // namespace ScalarWave::BoundaryCorrections
|
