UpwindPenaltyCorrection.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
9 #include "Domain/FaceNormal.hpp"
10 #include "Evolution/Systems/CurvedScalarWave/Characteristics.hpp"
12 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
13 #include "Options/Options.hpp"
14 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
15 #include "Utilities/Gsl.hpp"
16 #include "Utilities/ProtocolHelpers.hpp"
17 #include "Utilities/TMPL.hpp"
18 
19 /// \cond
20 class DataVector;
21 template <typename X, typename Symm, typename IndexList>
22 class Tensor;
23 /// \endcond
24 
25 namespace CurvedScalarWave {
26 /*!
27  * \ingroup NumericalFluxesGroup
28  * \brief Computes the upwind flux for scalar-waves in curved spacetime.
29  *
30  * \details The upwind flux is given by \cite Teukolsky2015ega :
31  *
32  * \f[
33  * G = S\left(\Lambda^+ S^{-1} U^{\rm int}
34  * + \Lambda^- S^{-1} U^{\rm ext}\right)
35  * = S \left(\Lambda^+ \hat{U}^{\rm int}
36  * + \Lambda^- \hat{U}^{\rm ext}\right),
37  * \f]
38  *
39  * where
40  *
41  * - \f$G\f$ is the numerical upwind flux dotted with the interface normal;
42  * - \f$U\f$ is a vector of all evolved variables;
43  * - \f$S\f$ is a matrix whose columns are the eigenvectors of the
44  * characteristic matrix for the evolution system. It maps the
45  * evolved variables to characteristic variables \f$\hat{U}\f$, s.t.
46  * \f$\hat{U} := S^{-1}\cdot U\f$; and
47  * - \f$\Lambda^\pm\f$ are diagonal matrices containing
48  * positive / negative eigenvalues of the same matrix as its elements.
49  *
50  * The superscripts \f${\rm int}\f$ and \f${\rm ext}\f$ on \f$U\f$ indicate
51  * that the corresponding set of variables at the element interface have been
52  * taken from the _interior_ or _exterior_ of the element. Exterior of the
53  * element is naturally the interior of its neighboring element. Therefore,
54  * \f$\hat{U}^{\rm int}:= S^{-1}U^{\rm int}\f$ are the characteristic variables
55  * at the element interface computed using evolved variables from the interior
56  * of the element, and \f$\hat{U}^{\rm ext}= S^{-1}U^{\rm ext}\f$ are the
57  * same computed from evolved variables taken from the element exterior, i.e.
58  * the neighboring element.
59  * The sign of characteristic speeds indicates the direction of propagation of
60  * the corresponding characteristic field with respect to the interface normal
61  * that the field has been computed along (with negative speeds indicating
62  * incoming characteristics, and with positive speeds indicating outgoing
63  * characteristics). Therefore, \f$\Lambda^+\f$ contains characterstic speeds
64  * for variables that are outgoing from the element at the interface, and
65  * \f$\Lambda^-\f$ contains characteristic speeds for variables that are
66  * incoming to the element at the interface. An ambiguity naturally arises
67  * as to which set of evolved variables (''interior'' or ''exterior'') to use
68  * when computing these speeds. We compute both and use their average, i.e.
69  * \f$\lambda^{\rm avg} = (\lambda^{\rm int} + \lambda^{\rm ext})/2\f$, to
70  * populate \f$\Lambda^\pm\f$.
71  *
72  *
73  * This function computes the upwind flux as follows:
74  * -# Computes internal and external characteristic variables and speeds using
75  * evolved variables from
76  * both the element interior and exterior.
77  * -# Computes the average characteristic speeds and constructs
78  * \f$\Lambda^{\pm}\f$.
79  * -# Computes the upwind flux as a weigted sum of external and internal
80  * characteristic variables
81  *
82  * \f[
83  * G = S\left(w_{\rm ext} \hat{U}^{\rm ext}
84  * + w_{\rm int} \hat{U}^{\rm int}\right),
85  * \f]
86  *
87  * with weights \f$w_{\rm ext} = \Theta(-\Lambda)\cdot\Lambda\f$, and
88  * \f$w_{\rm int} = \Theta(\Lambda)\cdot\Lambda\f$, where \f$\Theta\f$ is the
89  * step function centered at zero, \f$\Lambda = \Lambda^+ + \Lambda^-\f$, and
90  * the dot operator \f$(\cdot)\f$ indicates an element-wise product.
91  *
92  * \warning With the averaging of characteristic speeds, this flux does not
93  * satisfy the generalized Rankine-Hugoniot conditions
94  * \f${G}^{\rm int}(\hat{U}^{\rm int}, \hat{U}^{\rm ext})
95  * = - {G}^{\rm ext}(\hat{U}^{\rm ext}, \hat{U}^{\rm int})\f$,
96  * which enforces that the flux leaving the element is equal to the flux
97  * entering the neighboring element, and vice versa. This condition is
98  * important for a well-balanced scheme, and so please use this flux with
99  * caution.
100  */
101 template <size_t Dim>
102 struct UpwindPenaltyCorrection : tt::ConformsTo<dg::protocols::NumericalFlux> {
103  public:
104  using options = tmpl::list<>;
105  static constexpr Options::String help = {
106  "Computes the curved scalar-wave upwind flux."};
107 
108  // clang-tidy: non-const reference
109  void pup(PUP::er& /*p*/) noexcept {} // NOLINT
110 
111  using variables_tags = tmpl::list<Pi, Phi<Dim>, Psi>;
112 
113  // This is the data needed to compute the numerical flux.
114  // `dg::SendBoundaryFluxes` calls `package_data` to store these tags in a
115  // Variables. Local and remote values of this data are then combined in the
116  // `()` operator.
117  using package_field_tags = tmpl::list<
124  using package_extra_tags = tmpl::list<>;
125 
126  // These tags on the interface of the element are passed to
127  // `package_data` to provide the data needed to compute the numerical fluxes.
128  using argument_tags = tmpl::list<
135 
136  // pseudo-interface: used internally by Algorithm infrastructure, not
137  // user-level code
138  // Following the not-null pointer to packaged_data, this function expects as
139  // arguments the databox types of the `argument_tags`.
140  void package_data(
141  gsl::not_null<Scalar<DataVector>*> packaged_pi,
142  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> packaged_phi,
143  gsl::not_null<Scalar<DataVector>*> packaged_psi,
144  gsl::not_null<Scalar<DataVector>*> packaged_lapse,
145  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*> packaged_shift,
146  gsl::not_null<tnsr::II<DataVector, Dim, Frame::Inertial>*>
147  packaged_inverse_spatial_metric,
148  gsl::not_null<Scalar<DataVector>*> packaged_gamma1,
149  gsl::not_null<Scalar<DataVector>*> packaged_gamma2,
150  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
151  packaged_interface_unit_normal,
152  const Scalar<DataVector>& pi,
153  const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
154  const Scalar<DataVector>& psi, const Scalar<DataVector>& lapse,
155  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
156  const tnsr::II<DataVector, Dim, Frame::Inertial>& inverse_spatial_metric,
157  const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
158  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
159  const noexcept;
160 
161  // pseudo-interface: used internally by Algorithm infrastructure, not
162  // user-level code
163  // The arguments are first the system::variables_tag::tags_list wrapped in
164  // Tags::NormalDotNumericalFlux as not-null pointers to write the results
165  // into, then the package_tags on the interior side of the mortar followed by
166  // the package_tags on the exterior side.
167  void operator()(
168  gsl::not_null<Scalar<DataVector>*> pi_normal_dot_numerical_flux,
169  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
170  phi_normal_dot_numerical_flux,
171  gsl::not_null<Scalar<DataVector>*> psi_normal_dot_numerical_flux,
172  const Scalar<DataVector>& pi_int,
173  const tnsr::i<DataVector, Dim, Frame::Inertial>& phi_int,
174  const Scalar<DataVector>& psi_int, const Scalar<DataVector>& lapse_int,
175  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift_int,
176  const tnsr::II<DataVector, Dim, Frame::Inertial>&
177  inverse_spatial_metric_int,
178  const Scalar<DataVector>& gamma1_int,
179  const Scalar<DataVector>& gamma2_int,
180  const tnsr::i<DataVector, Dim, Frame::Inertial>&
181  interface_unit_normal_int,
182  const Scalar<DataVector>& pi_ext,
183  const tnsr::i<DataVector, Dim, Frame::Inertial>& phi_ext,
184  const Scalar<DataVector>& psi_ext, const Scalar<DataVector>& lapse_ext,
185  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift_ext,
186  const tnsr::II<DataVector, Dim, Frame::Inertial>&
187  inverse_spatial_metric_ext,
188  const Scalar<DataVector>& gamma1_ext,
189  const Scalar<DataVector>& gamma2_ext,
190  const tnsr::i<DataVector, Dim, Frame::Inertial>&
191  interface_unit_normal_ext) const noexcept;
192 };
193 } // namespace CurvedScalarWave
domain::Tags::UnnormalizedFaceNormal
Definition: FaceNormal.hpp:112
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 .
CurvedScalarWave::Tags::ConstraintGamma1
Definition: Tags.hpp:37
Options.hpp
FaceNormal.hpp
Tags.hpp
CurvedScalarWave::Tags::ConstraintGamma2
Definition: Tags.hpp:41
CurvedScalarWave::Phi
Definition: Tags.hpp:32
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.
Tags::Normalized
Definition: Magnitude.hpp:137
cstddef
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
CurvedScalarWave::Pi
Definition: Tags.hpp:27
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...
CurvedScalarWave::Psi
Definition: Tags.hpp:23
gr::Tags::Shift
Definition: Tags.hpp:48
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
Gsl.hpp
TypeAliases.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
CurvedScalarWave::UpwindPenaltyCorrection
Computes the upwind flux for scalar-waves in curved spacetime.
Definition: UpwindPenaltyCorrection.hpp:102
gr::Tags::Lapse
Definition: Tags.hpp:52
TMPL.hpp
gr::Tags::InverseSpatialMetric
Inverse of the spatial metric.
Definition: Tags.hpp:33
CurvedScalarWave
Items related to evolving a scalar wave on a curved background.
Definition: Characteristics.hpp:18
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13