UpwindPenalty.hpp
1 // 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 
12 #include "Evolution/Systems/GeneralizedHarmonic/BoundaryCorrections/BoundaryCorrection.hpp"
13 #include "Evolution/Systems/GeneralizedHarmonic/ConstraintDamping/Tags.hpp"
14 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
15 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
16 #include "Options/Options.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 
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 class UpwindPenalty final : public BoundaryCorrection<Dim> {
179  private:
180  struct NormalTimesVPlus : db::SimpleTag {
181  using type = tnsr::iaa<DataVector, Dim, Frame::Inertial>;
182  };
183  struct NormalTimesVMinus : db::SimpleTag {
184  using type = tnsr::iaa<DataVector, Dim, Frame::Inertial>;
185  };
186  struct Gamma2VSpacetimeMetric : db::SimpleTag {
187  using type = tnsr::aa<DataVector, Dim, Frame::Inertial>;
188  };
189  struct CharSpeedsTensor : db::SimpleTag {
190  using type = tnsr::a<DataVector, 3, Frame::Inertial>;
191  };
192 
193  public:
194  using options = tmpl::list<>;
195  static constexpr Options::String help = {
196  "Computes the UpwindPenalty boundary correction term for the generalized "
197  "harmonic system."};
198 
199  UpwindPenalty() = default;
200  UpwindPenalty(const UpwindPenalty&) = default;
201  UpwindPenalty& operator=(const UpwindPenalty&) = default;
202  UpwindPenalty(UpwindPenalty&&) = default;
203  UpwindPenalty& operator=(UpwindPenalty&&) = default;
204  ~UpwindPenalty() override = default;
205 
206  /// \cond
207  explicit UpwindPenalty(CkMigrateMessage* msg) noexcept;
208  using PUP::able::register_constructor;
210  /// \endcond
211  void pup(PUP::er& p) override; // NOLINT
212 
213  std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const noexcept override;
214 
215  using dg_package_field_tags =
216  tmpl::list<Tags::VSpacetimeMetric<Dim, Frame::Inertial>,
219  Tags::VMinus<Dim, Frame::Inertial>, NormalTimesVPlus,
220  NormalTimesVMinus, Gamma2VSpacetimeMetric, CharSpeedsTensor>;
221  using dg_package_data_temporary_tags = tmpl::list<
226  using dg_package_data_primitive_tags = tmpl::list<>;
227  using dg_package_data_volume_tags = tmpl::list<>;
228 
229  double dg_package_data(
230  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
231  packaged_char_speed_v_spacetime_metric,
232  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
233  packaged_char_speed_v_zero,
234  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
235  packaged_char_speed_v_plus,
236  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
237  packaged_char_speed_v_minus,
238  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
239  packaged_char_speed_n_times_v_plus,
240  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
241  packaged_char_speed_n_times_v_minus,
242  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
243  packaged_char_speed_gamma2_v_spacetime_metric,
244  gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*>
245  packaged_char_speeds,
246 
247  const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
248  const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
249  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
250 
251  const Scalar<DataVector>& constraint_gamma1,
252  const Scalar<DataVector>& constraint_gamma2,
253  const Scalar<DataVector>& lapse,
254  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
255 
256  const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
257  const tnsr::I<DataVector, Dim, Frame::Inertial>& normal_vector,
258  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
259  /*mesh_velocity*/,
260  const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity)
261  const noexcept;
262 
263  void dg_boundary_terms(
264  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
265  boundary_correction_spacetime_metric,
266  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
267  boundary_correction_pi,
268  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
269  boundary_correction_phi,
270 
271  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
272  char_speed_v_spacetime_metric_int,
273  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
274  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_plus_int,
275  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_minus_int,
276  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
277  char_speed_normal_times_v_plus_int,
278  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
279  char_speed_normal_times_v_minus_int,
280  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
281  char_speed_constraint_gamma2_v_spacetime_metric_int,
282  const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_int,
283 
284  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
285  char_speed_v_spacetime_metric_ext,
286  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
287  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_plus_ext,
288  const tnsr::aa<DataVector, Dim, Frame::Inertial>& char_speed_v_minus_ext,
289  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
290  char_speed_normal_times_v_plus_ext,
291  const tnsr::iaa<DataVector, Dim, Frame::Inertial>&
292  char_speed_normal_times_v_minus_ext,
293  const tnsr::aa<DataVector, Dim, Frame::Inertial>&
294  char_speed_constraint_gamma2_v_spacetime_metric_ext,
295  const tnsr::a<DataVector, 3, Frame::Inertial>& char_speeds_ext,
296  dg::Formulation /*dg_formulation*/) const noexcept;
297 };
298 } // namespace GeneralizedHarmonic::BoundaryCorrections
GeneralizedHarmonic::pi
void pi(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > * > pi, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the conjugate momentum of the spacetime metric .
CharmPupable.hpp
GeneralizedHarmonic::BoundaryCorrections::BoundaryCorrection
The base class used to make boundary corrections factory createable so they can be specified in the i...
Definition: BoundaryCorrection.hpp:25
Options.hpp
GeneralizedHarmonic::BoundaryCorrections
Boundary corrections/numerical fluxes.
Definition: BoundaryCorrection.hpp:14
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
gr::lapse
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
GeneralizedHarmonic::Tags::VMinus
Definition: Tags.hpp:120
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
cstddef
GeneralizedHarmonic::Tags::VZero
Definition: Tags.hpp:112
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
GeneralizedHarmonic::phi
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > * > phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein's equations...
memory
gr::Tags::Shift
Definition: Tags.hpp:48
GeneralizedHarmonic::ConstraintDamping::Tags::ConstraintGamma2
Constraint dammping parameter for the generalized harmonic system (cf. ).
Definition: Tags.hpp:70
gr::shift
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute shift from spacetime metric and inverse spatial metric.
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
gr::spacetime_metric
void spacetime_metric(gsl::not_null< tnsr::aa< DataType, Dim, Frame > * > spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, Dim, Frame > &shift, const tnsr::ii< DataType, Dim, Frame > &spatial_metric) noexcept
Computes the spacetime metric from the spatial metric, lapse, and shift.
GeneralizedHarmonic::Tags::VPlus
Definition: Tags.hpp:116
optional
GeneralizedHarmonic::BoundaryCorrections::UpwindPenalty
Computes the generalized harmonic upwind multipenalty boundary correction.
Definition: UpwindPenalty.hpp:178
GeneralizedHarmonic::ConstraintDamping::Tags::ConstraintGamma1
Constraint dammping parameter for the generalized harmonic system (cf. ).
Definition: Tags.hpp:62
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
std::unique_ptr
Prefixes.hpp
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
gr::Tags::Lapse
Definition: Tags.hpp:52
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13