UpwindPenalty.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <memory>
7 #include <optional>
8 
11 #include "Evolution/Systems/ScalarWave/BoundaryCorrections/BoundaryCorrection.hpp"
13 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
14 #include "Options/Options.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 
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 class UpwindPenalty final : public BoundaryCorrection<Dim> {
175  private:
176  struct NormalTimesVPlus : db::SimpleTag {
177  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
178  };
179  struct NormalTimesVMinus : db::SimpleTag {
180  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
181  };
182  struct Gamma2VPsi : db::SimpleTag {
183  using type = Scalar<DataVector>;
184  };
185  struct CharSpeedsTensor : db::SimpleTag {
186  using type = tnsr::i<DataVector, 3, Frame::Inertial>;
187  };
188 
189  public:
190  using options = tmpl::list<>;
191  static constexpr Options::String help = {
192  "Computes the UpwindPenalty boundary correction term for the scalar wave "
193  "system."};
194 
195  UpwindPenalty() = default;
196  UpwindPenalty(const UpwindPenalty&) = default;
197  UpwindPenalty& operator=(const UpwindPenalty&) = default;
198  UpwindPenalty(UpwindPenalty&&) = default;
199  UpwindPenalty& operator=(UpwindPenalty&&) = default;
200  ~UpwindPenalty() override = default;
201 
202  /// \cond
203  explicit UpwindPenalty(CkMigrateMessage* msg) noexcept;
204  using PUP::able::register_constructor;
206  /// \endcond
207  void pup(PUP::er& p) override; // NOLINT
208 
209  std::unique_ptr<BoundaryCorrection<Dim>> get_clone() const noexcept override;
210 
211  using dg_package_field_tags =
212  tmpl::list<Tags::VPsi, Tags::VZero<Dim>, Tags::VPlus, Tags::VMinus,
213  NormalTimesVPlus, NormalTimesVMinus, Gamma2VPsi,
214  CharSpeedsTensor>;
215  using dg_package_data_temporary_tags = tmpl::list<Tags::ConstraintGamma2>;
216  using dg_package_data_volume_tags = tmpl::list<>;
217 
218  double dg_package_data(
219  gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_psi,
220  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
221  packaged_char_speed_v_zero,
222  gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_plus,
223  gsl::not_null<Scalar<DataVector>*> packaged_char_speed_v_minus,
224  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
225  packaged_char_speed_n_times_v_plus,
226  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
227  packaged_char_speed_n_times_v_minus,
228  gsl::not_null<Scalar<DataVector>*> packaged_char_speed_gamma2_v_psi,
229  gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
230  packaged_char_speeds,
231 
232  const Scalar<DataVector>& pi,
233  const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
234  const Scalar<DataVector>& psi,
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)
242  const noexcept;
243 
244  void dg_boundary_terms(
245  gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
246  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
247  phi_boundary_correction,
248  gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
249 
250  const Scalar<DataVector>& char_speed_v_psi_int,
251  const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_int,
252  const Scalar<DataVector>& char_speed_v_plus_int,
253  const Scalar<DataVector>& char_speed_v_minus_int,
254  const tnsr::i<DataVector, Dim, Frame::Inertial>&
255  char_speed_normal_times_v_plus_int,
256  const tnsr::i<DataVector, Dim, Frame::Inertial>&
257  char_speed_normal_times_v_minus_int,
258  const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_int,
259  const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_int,
260 
261  const Scalar<DataVector>& char_speed_v_psi_ext,
262  const tnsr::i<DataVector, Dim, Frame::Inertial>& char_speed_v_zero_ext,
263  const Scalar<DataVector>& char_speed_v_plus_ext,
264  const Scalar<DataVector>& char_speed_v_minus_ext,
265  const tnsr::i<DataVector, Dim, Frame::Inertial>&
266  char_speed_minus_normal_times_v_plus_ext,
267  const tnsr::i<DataVector, Dim, Frame::Inertial>&
268  char_speed_minus_normal_times_v_minus_ext,
269  const Scalar<DataVector>& char_speed_constraint_gamma2_v_psi_ext,
270  const tnsr::i<DataVector, 3, Frame::Inertial>& char_speeds_ext,
271  dg::Formulation /*dg_formulation*/) const noexcept;
272 };
273 } // namespace ScalarWave::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
ScalarWave::Tags::VPlus
Definition: Tags.hpp:78
Options.hpp
ScalarWave::BoundaryCorrections::BoundaryCorrection
The base class used to make boundary corrections factory createable so they can be specified in the i...
Definition: BoundaryCorrection.hpp:24
Tags.hpp
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
ScalarWave::BoundaryCorrections
Boundary corrections/numerical fluxes.
Definition: BoundaryCorrection.hpp:13
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
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
ScalarWave::Tags::VMinus
Definition: Tags.hpp:82
optional
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: ReadSpecThirdOrderPiecewisePolynomial.hpp:11
ScalarWave::BoundaryCorrections::UpwindPenalty
Computes the scalar wave upwind multipenalty boundary correction.
Definition: UpwindPenalty.hpp:174
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13