Bjorhus.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <memory>
9 #include <optional>
10 #include <pup.h>
11 #include <string>
12 #include <type_traits>
13 
15 #include "DataStructures/DataVector.hpp"
18 #include "Evolution/BoundaryConditions/Type.hpp"
19 #include "Evolution/Systems/GeneralizedHarmonic/BoundaryConditions/BoundaryCondition.hpp"
20 #include "Evolution/Systems/GeneralizedHarmonic/ConstraintDamping/Tags.hpp"
21 #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
22 #include "Options/Options.hpp"
24 #include "PointwiseFunctions/AnalyticData/Tags.hpp"
25 #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
26 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
27 #include "Time/Tags.hpp"
28 #include "Utilities/Gsl.hpp"
29 #include "Utilities/TMPL.hpp"
30 
31 /// \cond
32 namespace domain::Tags {
33 template <size_t Dim, typename Frame>
34 struct Coordinates;
35 } // namespace domain::Tags
36 /// \endcond
37 
38 namespace GeneralizedHarmonic::BoundaryConditions::detail {
39 enum class ConstraintPreservingBjorhusType {
40  ConstraintPreserving,
41  ConstraintPreservingPhysical
42 };
43 
44 ConstraintPreservingBjorhusType
45 convert_constraint_preserving_bjorhus_type_from_yaml(
46  const Options::Option& options);
47 } // namespace GeneralizedHarmonic::BoundaryConditions::detail
48 
50 /*!
51  * \brief Sets constraint preserving boundary conditions using the Bjorhus
52  * method.
53  *
54  * \details Boundary conditions for the generalized harmonic evolution system
55  * can be divided in to three parts, constraint-preserving, physical and gauge
56  * boundary conditions.
57  *
58  * The generalized harmonic (GH) evolution system is a first-order reduction of
59  * Einstein equations brought about by the imposition of GH gauge. This
60  * introduces constraints on the free (evolved) variables in addition to the
61  * standard Hamiltonian and momentum constraints. The constraint-preserving
62  * portion of the boundary conditions is designed to prevent the influx of
63  * constraint violations from external faces of the evolution domain, by damping
64  * them away on a controlled and short time-scale. These conditions are imposed
65  * as corrections to the projections of the right-hand-sides of the GH evolution
66  * equations (i.e. using Bjorhus' method \cite Bjorhus1995), and are
67  * written down in Eq. (63) - (65) of \cite Lindblom2005qh . The gauge degrees
68  * of freedom are controlled by imposing a Sommerfeld-type condition (\f$L=0\f$
69  * member of the hierarchy derived in \cite BaylissTurkel) that allow gauge
70  * perturbations to pass through the boundary without strong reflections. These
71  * assume a spherical outer boundary, and can be written down as in Eq. (25) of
72  * \cite Rinne2007ui . Finally, the physical boundary conditions control the
73  * influx of inward propagating gravitational-wave solutions from the external
74  * boundaries. These are derived by considering the evolution system of the Weyl
75  * curvature tensor, and controlling the inward propagating characteristics of
76  * the system that are proportional to the Newman-Penrose curvature spinor
77  * components \f$\Psi_4\f$ and \f$\Psi_0\f$. Here we use Eq. (68) of
78  * \cite Lindblom2005qh to disallow any incoming waves. It is to be noted that
79  * all the above conditions are also imposed on characteristic modes with speeds
80  * exactly zero.
81  *
82  * This class provides two choices of combinations of the above corrections:
83  * - `ConstraintPreserving` : this imposes the constraint-preserving and
84  * gauge-controlling corrections;
85  * - `ConstraintPreservingPhysical` : this additionally restricts the influx of
86  * any physical gravitational waves from the outer boundary, in addition to
87  * preventing the influx of constraint violations and gauge perturbations.
88  *
89  * We refer to `Bjorhus::constraint_preserving_bjorhus_corrections_dt_v_psi()`,
90  * `Bjorhus::constraint_preserving_bjorhus_corrections_dt_v_zero()`,
91  * `Bjorhus::constraint_preserving_bjorhus_corrections_dt_v_minus()`, and
92  * `Bjorhus::constraint_preserving_physical_bjorhus_corrections_dt_v_minus()`
93  * for the further details on implementation.
94  *
95  * \note These boundary conditions assume a spherical outer boundary. Also, we
96  * do not yet have an option to inject incoming gravitational waves at the outer
97  * boundary.
98  */
99 template <size_t Dim>
101  public:
102  struct TypeOptionTag {
103  using type = detail::ConstraintPreservingBjorhusType;
104  static std::string name() noexcept { return "Type"; }
105  static constexpr Options::String help{
106  "Whether to impose ConstraintPreserving, with or without physical "
107  "terms for VMinus."};
108  };
109 
110  using options = tmpl::list<TypeOptionTag>;
111  static constexpr Options::String help{
112  "ConstraintPreservingBjorhus boundary conditions setting the value of the"
113  "time derivatives of the spacetime metric, Phi and Pi to expressions that"
114  "prevent the influx of constraint violations and reflections."};
115  static std::string name() noexcept { return "ConstraintPreservingBjorhus"; }
116 
117  ConstraintPreservingBjorhus(
118  detail::ConstraintPreservingBjorhusType type) noexcept;
119 
120  ConstraintPreservingBjorhus() = default;
121  /// \cond
122  ConstraintPreservingBjorhus(ConstraintPreservingBjorhus&&) noexcept = default;
123  ConstraintPreservingBjorhus& operator=(
124  ConstraintPreservingBjorhus&&) noexcept = default;
125  ConstraintPreservingBjorhus(const ConstraintPreservingBjorhus&) = default;
126  ConstraintPreservingBjorhus& operator=(const ConstraintPreservingBjorhus&) =
127  default;
128  /// \endcond
129  ~ConstraintPreservingBjorhus() override = default;
130 
131  explicit ConstraintPreservingBjorhus(CkMigrateMessage* msg) noexcept;
132 
133  WRAPPED_PUPable_decl_base_template(
134  domain::BoundaryConditions::BoundaryCondition,
135  ConstraintPreservingBjorhus);
136 
137  auto get_clone() const noexcept -> std::unique_ptr<
138  domain::BoundaryConditions::BoundaryCondition> override;
139 
140  static constexpr evolution::BoundaryConditions::Type bc_type =
141  evolution::BoundaryConditions::Type::TimeDerivative;
142 
143  void pup(PUP::er& p) override;
144 
145  using dg_interior_evolved_variables_tags =
146  tmpl::list<gr::Tags::SpacetimeMetric<Dim, Frame::Inertial, DataVector>,
147  Tags::Pi<Dim, Frame::Inertial>,
148  Tags::Phi<Dim, Frame::Inertial>>;
149  using dg_interior_temporary_tags = tmpl::list<
150  domain::Tags::Coordinates<Dim, Frame::Inertial>,
151  ConstraintDamping::Tags::ConstraintGamma1,
152  ConstraintDamping::Tags::ConstraintGamma2, gr::Tags::Lapse<DataVector>,
153  gr::Tags::Shift<Dim, Frame::Inertial, DataVector>,
154  gr::Tags::InverseSpacetimeMetric<Dim, Frame::Inertial, DataVector>,
155  gr::Tags::SpacetimeNormalVector<Dim, Frame::Inertial, DataVector>,
156  gr::Tags::SpacetimeNormalOneForm<Dim, Frame::Inertial, DataVector>,
157  Tags::ThreeIndexConstraint<Dim, Frame::Inertial>,
158  Tags::GaugeH<Dim, Frame::Inertial>,
159  Tags::SpacetimeDerivGaugeH<Dim, Frame::Inertial>>;
160  using dg_interior_dt_vars_tags = tmpl::list<
161  ::Tags::dt<gr::Tags::SpacetimeMetric<Dim, Frame::Inertial, DataVector>>,
162  ::Tags::dt<Tags::Pi<Dim, Frame::Inertial>>,
163  ::Tags::dt<Tags::Phi<Dim, Frame::Inertial>>>;
164  using dg_interior_deriv_vars_tags = tmpl::list<
165  ::Tags::deriv<gr::Tags::SpacetimeMetric<Dim, Frame::Inertial, DataVector>,
166  tmpl::size_t<Dim>, Frame::Inertial>,
167  ::Tags::deriv<Tags::Pi<Dim, Frame::Inertial>, tmpl::size_t<Dim>,
168  Frame::Inertial>,
169  ::Tags::deriv<Tags::Phi<Dim, Frame::Inertial>, tmpl::size_t<Dim>,
170  Frame::Inertial>>;
171  using dg_gridless_tags = tmpl::list<>;
172 
173  std::optional<std::string> dg_time_derivative(
174  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
175  dt_spacetime_metric_correction,
176  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
177  dt_pi_correction,
178  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
179  dt_phi_correction,
180 
181  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
182  face_mesh_velocity,
183  const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
184  const tnsr::I<DataVector, Dim, Frame::Inertial>& /*normal_vector*/,
185  // c.f. dg_interior_evolved_variables_tags
186  const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
187  const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
188  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
189  // c.f. dg_interior_temporary_tags
190  const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
191  const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
192  const Scalar<DataVector>& lapse,
193  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
194  const tnsr::AA<DataVector, Dim, Frame::Inertial>&
196  const tnsr::A<DataVector, Dim, Frame::Inertial>&
197  spacetime_unit_normal_vector,
198  const tnsr::a<DataVector, Dim, Frame::Inertial>&
199  spacetime_unit_normal_one_form,
200  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& three_index_constraint,
201  const tnsr::a<DataVector, Dim, Frame::Inertial>& gauge_source,
202  const tnsr::ab<DataVector, Dim, Frame::Inertial>&
203  spacetime_deriv_gauge_source,
204  // c.f. dg_interior_dt_vars_tags
205  const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_spacetime_metric,
206  const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_pi,
207  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& dt_phi,
208  // c.f. dg_interior_deriv_vars_tags
209  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_spacetime_metric,
210  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_pi,
211  const tnsr::ijaa<DataVector, Dim, Frame::Inertial>& d_phi) const noexcept;
212 
213  private:
214  void compute_intermediate_vars(
215  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
216  unit_interface_normal_vector,
217  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
219  gsl::not_null<tnsr::II<DataVector, Dim, Frame::Inertial>*>
220  inverse_spatial_metric,
221  gsl::not_null<tnsr::ii<DataVector, Dim, Frame::Inertial>*>
223  gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
224  incoming_null_one_form,
225  gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
226  outgoing_null_one_form,
227  gsl::not_null<tnsr::A<DataVector, Dim, Frame::Inertial>*>
228  incoming_null_vector,
229  gsl::not_null<tnsr::A<DataVector, Dim, Frame::Inertial>*>
230  outgoing_null_vector,
231  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*> projection_ab,
232  gsl::not_null<tnsr::Ab<DataVector, Dim, Frame::Inertial>*> projection_Ab,
233  gsl::not_null<tnsr::AA<DataVector, Dim, Frame::Inertial>*> projection_AB,
234  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
235  char_projected_rhs_dt_v_psi,
236  gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
237  char_projected_rhs_dt_v_zero,
238  gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
239  char_projected_rhs_dt_v_minus,
240  gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
241  constraint_char_zero_plus,
242  gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
243  constraint_char_zero_minus,
244  gsl::not_null<std::array<DataVector, 4>*> char_speeds,
245 
246  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
247  face_mesh_velocity,
248  const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
249  const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
250  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
251  const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
252  const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
253  const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
254  const Scalar<DataVector>& lapse,
255  const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
256  const tnsr::AA<DataVector, Dim, Frame::Inertial>&
258  const tnsr::A<DataVector, Dim, Frame::Inertial>&
259  spacetime_unit_normal_vector,
260  const tnsr::a<DataVector, Dim, Frame::Inertial>&
261  spacetime_unit_normal_one_form,
262  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& three_index_constraint,
263  const tnsr::a<DataVector, Dim, Frame::Inertial>& gauge_source,
264  const tnsr::ab<DataVector, Dim, Frame::Inertial>&
265  spacetime_deriv_gauge_source,
266  const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_pi,
267  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& dt_phi,
268  const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_spacetime_metric,
269  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_pi,
270  const tnsr::ijaa<DataVector, Dim, Frame::Inertial>& d_phi,
271  const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_spacetime_metric)
272  const noexcept;
273 
274  detail::ConstraintPreservingBjorhusType type_{
275  detail::ConstraintPreservingBjorhusType::ConstraintPreservingPhysical};
276 };
277 } // namespace GeneralizedHarmonic::BoundaryConditions
278 
279 template <>
280 struct Options::create_from_yaml<GeneralizedHarmonic::BoundaryConditions::
281  detail::ConstraintPreservingBjorhusType> {
282  template <typename Metavariables>
283  static typename GeneralizedHarmonic::BoundaryConditions::detail::
284  ConstraintPreservingBjorhusType
285  create(const Options::Option& options) {
286  return GeneralizedHarmonic::BoundaryConditions::detail::
287  convert_constraint_preserving_bjorhus_type_from_yaml(options);
288  }
289 };
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 .
GeneralizedHarmonic::three_index_constraint
tnsr::iaa< DataType, SpatialDim, Frame > three_index_constraint(const tnsr::iaa< DataType, SpatialDim, Frame > &d_spacetime_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the generalized-harmonic 3-index constraint.
std::string
CharmPupable.hpp
GeneralizedHarmonic::four_index_constraint
tnsr::iaa< DataType, SpatialDim, Frame > four_index_constraint(const tnsr::ijaa< DataType, SpatialDim, Frame > &d_phi) noexcept
Computes the generalized-harmonic 4-index constraint.
evolution
Functionality for evolving hyperbolic partial differential equations.
Definition: RunEventsAndDenseTriggers.hpp:30
Options.hpp
GeneralizedHarmonic::BoundaryConditions
Boundary conditions for the generalized harmonic system.
Definition: Bjorhus.hpp:38
GeneralizedHarmonic
Items related to evolving the first-order generalized harmonic system.
Definition: Bjorhus.hpp:38
GeneralizedHarmonic::BoundaryConditions::ConstraintPreservingBjorhus
Sets constraint preserving boundary conditions using the Bjorhus method.
Definition: Bjorhus.hpp:100
Options::Option
Definition: Options.hpp:108
GeneralizedHarmonic::BoundaryConditions::ConstraintPreservingBjorhus::TypeOptionTag
Definition: Bjorhus.hpp:102
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::extrinsic_curvature
void extrinsic_curvature(gsl::not_null< tnsr::ii< DataType, SpatialDim, Frame > * > ex_curv, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_normal_vector, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes extrinsic curvature from generalized harmonic variables and the spacetime normal vector.
cstddef
GeneralizedHarmonic::gauge_source
void gauge_source(gsl::not_null< tnsr::a< DataType, SpatialDim, Frame > * > gauge_source_h, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const Scalar< DataType > &trace_extrinsic_curvature, const tnsr::i< DataType, SpatialDim, Frame > &trace_christoffel_last_indices) noexcept
Computes generalized harmonic gauge source function.
Options::create_from_yaml
Definition: MinmodType.hpp:11
GeneralizedHarmonic::TimeDerivative
Compute the RHS of the Generalized Harmonic formulation of Einstein's equations.
Definition: TimeDerivative.hpp:87
array
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
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
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Frame
Definition: IndexType.hpp:36
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.
gr
Definition: GaugeWave.hpp:28
Tensor.hpp
GeneralizedHarmonic::BoundaryConditions::BoundaryCondition
The base class off of which all boundary conditions must inherit.
Definition: BoundaryCondition.hpp:27
optional
Prefixes.hpp
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
type_traits
gr::inverse_spacetime_metric
void inverse_spacetime_metric(gsl::not_null< tnsr::AA< DataType, SpatialDim, Frame > * > inverse_spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute inverse spacetime metric from inverse spatial metric, lapse and shift.
TMPL.hpp
domain::Tags
Tags for the domain.
Definition: FaceNormal.hpp:107
string