ConstraintPreservingSphericalRadiation.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 #include <pup.h>
9 #include <string>
10 #include <type_traits>
11 
13 #include "DataStructures/DataVector.hpp"
16 #include "Evolution/BoundaryConditions/Type.hpp"
17 #include "Evolution/Systems/ScalarWave/BoundaryConditions/BoundaryCondition.hpp"
20 #include "Options/Options.hpp"
22 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
23 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
24 #include "Time/Tags.hpp"
25 #include "Utilities/Gsl.hpp"
26 #include "Utilities/TMPL.hpp"
27 
28 /// \cond
29 namespace domain::Tags {
30 template <size_t Dim, typename Frame>
31 struct Coordinates;
32 } // namespace domain::Tags
33 /// \endcond
34 
36 namespace detail {
37 /// The type of spherical radiation boundary condition to impose
38 enum class ConstraintPreservingSphericalRadiationType {
39  /// Impose \f$(\partial_t + \partial_r)\Psi=0\f$
40  Sommerfeld,
41  /// Impose \f$(\partial_t + \partial_r + r^{-1})\Psi=0\f$
42  FirstOrderBaylissTurkel,
43  /// Imposes a second-order Bayliss-Turkel boundary condition
44  SecondOrderBaylissTurkel
45 };
46 
47 ConstraintPreservingSphericalRadiationType
48 convert_constraint_preserving_spherical_radiation_type_from_yaml(
49  const Options::Option& options);
50 } // namespace detail
51 
52 /*!
53  * \brief Constraint-preserving spherical radiation boundary condition that
54  * seeks to avoid ingoing constraint violations and radiation.
55  *
56  * The constraint-preserving part of the boundary condition imposes the
57  * following condition on the time derivatives of the characteristic fields:
58  *
59  * \f{align*}{
60  * d_tw^\Psi&\to d_tw^{\Psi}+\lambda_{\Psi}n^i\mathcal{C}_i, \\
61  * d_tw^0_i&\to d_tw^{0}_i+\lambda_{0}n^jP^k_i\mathcal{C}_{ik},
62  * \f}
63  *
64  * where
65  *
66  * \f{align*}{
67  * P^k{}_i=\delta^k_i-n^kn_i
68  * \f}
69  *
70  * projects a tensor onto the spatial surface to which \f$n_i\f$ is normal, and
71  * \f$d_t w\f$ is the evolved to characteristic field transformation applied to
72  * the time derivatives of the evolved fields. That is,
73  *
74  * \f{align*}{
75  * d_t w^\Psi&=\partial_t \Psi, \\
76  * d_t w_i^0&=(\delta^k_i-n^k n_i)\partial_t \Phi_k, \\
77  * d_t w^{\pm}&=\partial_t\Pi\pm n^k\partial_t\Phi_k - \gamma_2\partial_t\Psi.
78  * \f}
79  *
80  * The constraints are defined as:
81  *
82  * \f{align*}{
83  * \mathcal{C}_i&=\partial_i\Psi - \Phi_i=0, \\
84  * \mathcal{C}_{ij}&=\partial_{[i}\Phi_{j]}=0
85  * \f}
86  *
87  * Radiation boundary conditions impose a condition on \f$\Pi\f$ or its time
88  * derivative. We denote the boundary condition value of the time derivative of
89  * \f$\Pi\f$ by \f$\partial_t\Pi^{\mathrm{BC}}\f$. With this, we can impose
90  * boundary conditions on the time derivatives of the evolved variables as
91  * follows:
92  *
93  * \f{align*}{
94  * \partial_{t} \Psi&\to\partial_{t}\Psi +
95  * \lambda_\Psi n^i \mathcal{C}_i, \\
96  * \partial_{t}\Pi&\to\partial_{t}\Pi-\left(\partial_t\Pi
97  * - \partial_t\Pi^{\mathrm{BC}}\right)
98  * +\gamma_2\lambda_\Psi n^i \mathcal{C}_i
99  * =\partial_t\Pi^{\mathrm{BC}}
100  * +\gamma_2\lambda_\Psi n^i \mathcal{C}_i, \\
101  * \partial_{t}\Phi_i&\to\partial_{t}\Phi_i+
102  * \lambda_0n^jP^k{}_i\mathcal{C}_{jk}
103  * = \partial_{t}\Phi_i+ \lambda_0n^j \mathcal{C}_{ji}.
104  * \f}
105  *
106  * Below we assume the normal vector \f$n^i\f$ is the radial unit normal vector.
107  * That is, we assume the outer boundary is spherical. A Sommerfeld
108  * \cite Sommerfeld1949 radiation condition is given by
109  *
110  * \f{align*}{
111  * \partial_t\Psi=n^i\Phi_i
112  * \f}
113  *
114  * Or, assuming that \f$\partial_tn^i=0\f$ (or is very small),
115  *
116  * \f{align*}{
117  * \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i
118  * \f}
119  *
120  * The Bayliss-Turkel \cite BaylissTurkel boundary conditions are given by:
121  *
122  * \f{align*}{
123  * \prod_{l=1}^m\left(\partial_t + \partial_r + \frac{2l-1}{r}\right)\Psi=0
124  * \f}
125  *
126  * The first-order form is
127  *
128  * \f{align*}{
129  * \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i + \frac{1}{r}\partial_t\Psi,
130  * \f}
131  *
132  * assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
133  *
134  * The second-order boundary condition is given by,
135  *
136  * \f{align*}{
137  * \partial_t\Pi^{\mathrm{BC}}
138  * &=\left(\partial_t\partial_r + \partial_r\partial_t +
139  * \partial_r^2+\frac{4}{r}\partial_t
140  * +\frac{4}{r}\partial_r + \frac{2}{r^2}\right)\Psi \notag \\
141  * &=n^i(\partial_t\Phi_i-\partial_i\Pi) + n^i n^j\partial_i\Phi_j -
142  * \frac{4}{r}\Pi+\frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi,
143  * \f}
144  *
145  * assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
146  *
147  * The moving mesh can be accounted for by using
148  *
149  * \f{align*}{
150  * \partial_t r = \frac{1}{r}x^i\delta_{ij}\partial_t x^j
151  * \f}
152  *
153  * \note It is not clear if \f$\partial_t\Phi_i\f$ should be replaced by
154  * \f$-\partial_i\Pi\f$, which is the evolution equation but without the
155  * constraint.
156  *
157  * \note On a moving mesh the characteristic speeds change according to
158  * \f$\lambda\to\lambda-v^i_gn_i\f$ where \f$v^i_g\f$ is the mesh velocity.
159  *
160  * \note For the scalar wave system \f$\lambda_0 = \lambda_\psi\f$
161  *
162  * \warning The boundary conditions are implemented assuming the outer boundary
163  * is spherical. It might be possible to generalize the condition to
164  * non-spherical boundaries by using \f$x^i/r\f$ instead of \f$n^i\f$, but this
165  * hasn't been tested.
166  */
167 template <size_t Dim>
169  : public BoundaryCondition<Dim> {
170  public:
171  struct TypeOptionTag {
172  using type = detail::ConstraintPreservingSphericalRadiationType;
173  static std::string name() noexcept { return "Type"; }
174  static constexpr Options::String help{
175  "Whether to impose Sommerfeld, first-order Bayliss-Turkel, or "
176  "second-order Bayliss-Turkel spherical radiation boundary conditions."};
177  };
178 
179  using options = tmpl::list<TypeOptionTag>;
180  static constexpr Options::String help{
181  "Constraint-preserving spherical radiation boundary conditions setting "
182  "the time derivatives of Psi, Phi, and Pi to avoid incoming constraint "
183  "violations, and imposing radiation boundary conditions."};
184 
186  detail::ConstraintPreservingSphericalRadiationType type) noexcept;
187 
189  /// \cond
191  ConstraintPreservingSphericalRadiation&&) noexcept = default;
193  ConstraintPreservingSphericalRadiation&&) noexcept = default;
195  const ConstraintPreservingSphericalRadiation&) = default;
197  const ConstraintPreservingSphericalRadiation&) = default;
198  /// \endcond
199  ~ConstraintPreservingSphericalRadiation() override = default;
200 
202  CkMigrateMessage* msg) noexcept;
203 
204  WRAPPED_PUPable_decl_base_template(
205  domain::BoundaryConditions::BoundaryCondition,
207 
208  auto get_clone() const noexcept -> std::unique_ptr<
209  domain::BoundaryConditions::BoundaryCondition> override;
210 
211  static constexpr evolution::BoundaryConditions::Type bc_type =
212  evolution::BoundaryConditions::Type::TimeDerivative;
213 
214  void pup(PUP::er& p) override;
215 
216  using dg_interior_evolved_variables_tags =
217  tmpl::list<ScalarWave::Pi, ScalarWave::Phi<Dim>, ScalarWave::Psi>;
218  using dg_interior_temporary_tags =
219  tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
220  Tags::ConstraintGamma2>;
221  using dg_interior_dt_vars_tags =
222  tmpl::list<::Tags::dt<ScalarWave::Pi>, ::Tags::dt<ScalarWave::Phi<Dim>>,
223  ::Tags::dt<ScalarWave::Psi>>;
224  using dg_interior_deriv_vars_tags = tmpl::list<
225  ::Tags::deriv<ScalarWave::Pi, tmpl::size_t<Dim>, Frame::Inertial>,
226  ::Tags::deriv<ScalarWave::Psi, tmpl::size_t<Dim>, Frame::Inertial>,
227  ::Tags::deriv<ScalarWave::Phi<Dim>, tmpl::size_t<Dim>, Frame::Inertial>>;
228  using dg_gridless_tags = tmpl::list<>;
229 
230  std::optional<std::string> dg_time_derivative(
231  gsl::not_null<Scalar<DataVector>*> dt_pi_correction,
232  gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
233  dt_phi_correction,
234  gsl::not_null<Scalar<DataVector>*> dt_psi_correction,
235  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
236  face_mesh_velocity,
237  const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
238  const Scalar<DataVector>& pi,
239  const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
240  const Scalar<DataVector>& psi,
241  const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
242  const Scalar<DataVector>& gamma2, const Scalar<DataVector>& dt_pi,
243  const tnsr::i<DataVector, Dim, Frame::Inertial>& dt_phi,
244  const Scalar<DataVector>& dt_psi,
245  const tnsr::i<DataVector, Dim, Frame::Inertial>& d_pi,
246  const tnsr::i<DataVector, Dim, Frame::Inertial>& d_psi,
247  const tnsr::ij<DataVector, Dim, Frame::Inertial>& d_phi) const noexcept;
248 
249  private:
250  detail::ConstraintPreservingSphericalRadiationType type_{
251  detail::ConstraintPreservingSphericalRadiationType::Sommerfeld};
252 };
253 } // namespace ScalarWave::BoundaryConditions
254 
255 template <>
257  ScalarWave::BoundaryConditions::detail::
258  ConstraintPreservingSphericalRadiationType> {
259  template <typename Metavariables>
260  static typename ScalarWave::BoundaryConditions::detail::
261  ConstraintPreservingSphericalRadiationType
262  create(const Options::Option& options) {
263  return ScalarWave::BoundaryConditions::detail::
264  convert_constraint_preserving_spherical_radiation_type_from_yaml(
265  options);
266  }
267 };
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 .
std::string
CharmPupable.hpp
ScalarWave::TimeDerivative
Compute the time derivatives for scalar wave system.
Definition: TimeDerivative.hpp:26
evolution
Functionality for evolving hyperbolic partial differential equations.
Definition: RunEventsAndDenseTriggers.hpp:30
Options.hpp
ScalarWave::Phi
Definition: Tags.hpp:31
Tags.hpp
ScalarWave::BoundaryConditions::ConstraintPreservingSphericalRadiation::TypeOptionTag
Definition: ConstraintPreservingSphericalRadiation.hpp:171
Options::Option
Definition: Options.hpp:108
ScalarWave::BoundaryConditions::BoundaryCondition
The base class off of which all boundary conditions must inherit.
Definition: BoundaryCondition.hpp:27
Options::create_from_yaml
Definition: MinmodType.hpp:11
ScalarWave::BoundaryConditions
Boundary conditions for the scalar wave system.
Definition: BoundaryCondition.hpp:24
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
Variables.hpp
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
Tensor.hpp
ScalarWave
Items related to evolving the scalar wave equation.
Definition: BoundaryCondition.hpp:24
ScalarWave::Psi
Definition: Tags.hpp:20
optional
PartialDerivatives.hpp
Prefixes.hpp
ScalarWave::BoundaryConditions::ConstraintPreservingSphericalRadiation
Constraint-preserving spherical radiation boundary condition that seeks to avoid ingoing constraint v...
Definition: ConstraintPreservingSphericalRadiation.hpp:168
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
type_traits
ScalarWave::Pi
Definition: Tags.hpp:25
TMPL.hpp
domain::Tags
Tags for the domain.
Definition: FaceNormal.hpp:107
string