Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <memory>
10 : #include <optional>
11 : #include <pup.h>
12 : #include <string>
13 : #include <type_traits>
14 :
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "DataStructures/DataVector.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Evolution/BoundaryConditions/Type.hpp"
20 : #include "Evolution/Systems/GeneralizedHarmonic/BoundaryConditions/BoundaryCondition.hpp"
21 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
22 : #include "Options/Auto.hpp"
23 : #include "Options/Options.hpp"
24 : #include "Options/String.hpp"
25 : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
26 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
27 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/ConstraintDampingTags.hpp"
28 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
29 : #include "PointwiseFunctions/MathFunctions/MathFunction.hpp"
30 : #include "Utilities/Gsl.hpp"
31 : #include "Utilities/Serialization/CharmPupable.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 : /// \cond
35 : namespace domain::Tags {
36 : template <size_t Dim, typename Frame>
37 : struct Coordinates;
38 : } // namespace domain::Tags
39 : namespace Tags {
40 : struct Time;
41 : } // namespace Tags
42 : /// \endcond
43 :
44 1 : namespace gh::BoundaryConditions::detail {
45 : enum class ConstraintPreservingBjorhusType {
46 : ConstraintPreserving,
47 : ConstraintPreservingPhysical
48 : };
49 :
50 : ConstraintPreservingBjorhusType
51 : convert_constraint_preserving_bjorhus_type_from_yaml(
52 : const Options::Option& options);
53 : } // namespace gh::BoundaryConditions::detail
54 :
55 : namespace gh::BoundaryConditions {
56 : /*!
57 : * \brief Sets constraint preserving boundary conditions using the Bjorhus
58 : * method.
59 : *
60 : * \details Boundary conditions for the generalized harmonic evolution system
61 : * can be divided in to three parts, constraint-preserving, physical and gauge
62 : * boundary conditions.
63 : *
64 : * The generalized harmonic (GH) evolution system is a first-order reduction of
65 : * Einstein equations brought about by the imposition of GH gauge. This
66 : * introduces constraints on the free (evolved) variables in addition to the
67 : * standard Hamiltonian and momentum constraints. The constraint-preserving
68 : * portion of the boundary conditions is designed to prevent the influx of
69 : * constraint violations from external faces of the evolution domain, by damping
70 : * them away on a controlled and short time-scale. These conditions are imposed
71 : * as corrections to the characteristic projections of the right-hand-sides of
72 : * the GH evolution equations (i.e. using Bjorhus' method \cite Bjorhus1995),
73 : * as written down in Eq. (63) - (65) of \cite Lindblom2005qh . In addition to
74 : * these equations, the fourth projection is simply frozen in the unlikely case
75 : * its coordinate speed becomes negative, i.e. \f$d_t u^{\hat{1}+}{}_{ab}=0\f$
76 : * (in the notation of \cite Lindblom2005qh). The gauge degrees
77 : * of freedom are controlled by imposing a Sommerfeld-type condition (\f$L=0\f$
78 : * member of the hierarchy derived in \cite BaylissTurkel) that allow gauge
79 : * perturbations to pass through the boundary without strong reflections. These
80 : * assume a spherical outer boundary, and can be written down as in Eq. (25) of
81 : * \cite Rinne2007ui . Finally, the physical boundary conditions control the
82 : * influx of inward propagating gravitational-wave solutions from the external
83 : * boundaries. These are derived by considering the evolution system of the Weyl
84 : * curvature tensor, and controlling the inward propagating characteristics of
85 : * the system that are proportional to the Newman-Penrose curvature spinor
86 : * components \f$\Psi_4\f$ and \f$\Psi_0\f$. Here we use Eq. (68) of
87 : * \cite Lindblom2005qh to disallow any incoming waves. It is to be noted that
88 : * all the above conditions are also imposed on characteristic modes with speeds
89 : * exactly zero.
90 : *
91 : * An optional injected incoming-wave contribution can be specified through
92 : * `IncomingWaveProfile`. The injected strain-rate tensor is
93 : * \f[
94 : * \dot{h}_{ab} = \dot{f}(t)
95 : * \left(\hat{x}_a \hat{x}_b + \hat{y}_a \hat{y}_b - 2 \hat{z}_a
96 : * \hat{z}_b\right),
97 : * \f]
98 : * where the configured profile is interpreted directly as \f$\dot{f}(t)\f$ and
99 : * \f$\hat{x}_a\f$, \f$\hat{y}_a\f$ and \f$\hat{z}_a\f$ are the components of
100 : * the coordinate basis vectors. This formula is taken from
101 : * \cite Lindblom2005qh. It should be considered an approximate perturbation
102 : * rather than an exact gravitational wave which would have to be constructed at
103 : * null-infinity. Note that the profile of the injected wave is specified at the
104 : * outer boundary and its amplitude should be adjusted according to the usual
105 : * 1/r scaling.
106 : *
107 : * This class provides two choices of combinations of the above corrections:
108 : * - `ConstraintPreserving` : this imposes the constraint-preserving and
109 : * gauge-controlling corrections;
110 : * - `ConstraintPreservingPhysical` : this additionally restricts the influx of
111 : * any physical gravitational waves from the outer boundary, in addition to
112 : * preventing the influx of constraint violations and gauge perturbations.
113 : *
114 : * We refer to `Bjorhus::constraint_preserving_corrections_dt_v_psi()`,
115 : * `Bjorhus::constraint_preserving_corrections_dt_v_zero()`,
116 : * `Bjorhus::constraint_preserving_gauge_corrections_dt_v_minus()`, and
117 : * `Bjorhus::constraint_preserving_gauge_physical_corrections_dt_v_minus()`
118 : * for the further details on implementation.
119 : *
120 : * \note These boundary conditions assume a spherical outer boundary.
121 : */
122 : template <size_t Dim>
123 1 : class ConstraintPreservingBjorhus final : public BoundaryCondition<Dim> {
124 : public:
125 0 : struct TypeOptionTag {
126 0 : using type = detail::ConstraintPreservingBjorhusType;
127 0 : static std::string name() { return "Type"; }
128 0 : static constexpr Options::String help{
129 : "Whether to impose ConstraintPreserving, with or without physical "
130 : "terms for VMinus."};
131 : };
132 :
133 0 : struct IncomingWaveProfileOptionTag {
134 0 : using type =
135 : Options::Auto<std::unique_ptr<::MathFunction<1, Frame::Inertial>>,
136 : Options::AutoLabel::None>;
137 0 : static std::string name() { return "IncomingWaveProfile"; }
138 0 : static constexpr Options::String help{
139 : "Optional incoming-wave profile interpreted as the first time "
140 : "derivative for the injected physical wave. See the "
141 : "ConstraintPreservingBjorhus class documentation for the injected-wave "
142 : "formula. Specify `None` to disable injection. This option is only "
143 : "supported in 3D."};
144 : };
145 :
146 0 : using options = tmpl::flatten<tmpl::list<
147 : TypeOptionTag,
148 : tmpl::conditional_t<Dim == 3, tmpl::list<IncomingWaveProfileOptionTag>,
149 : tmpl::list<>>>>;
150 0 : static constexpr Options::String help{
151 : "ConstraintPreservingBjorhus boundary conditions setting the value of the"
152 : "time derivatives of the spacetime metric, Phi and Pi to expressions that"
153 : "prevent the influx of constraint violations and reflections."};
154 0 : static std::string name() { return "ConstraintPreservingBjorhus"; }
155 :
156 0 : explicit ConstraintPreservingBjorhus(
157 : detail::ConstraintPreservingBjorhusType type,
158 : std::optional<std::unique_ptr<::MathFunction<1, Frame::Inertial>>>
159 : incoming_wave_profile = std::nullopt);
160 :
161 0 : ConstraintPreservingBjorhus() = default;
162 : /// \cond
163 : ConstraintPreservingBjorhus(ConstraintPreservingBjorhus&&) = default;
164 : ConstraintPreservingBjorhus& operator=(ConstraintPreservingBjorhus&&) =
165 : default;
166 : ConstraintPreservingBjorhus(const ConstraintPreservingBjorhus&);
167 : ConstraintPreservingBjorhus& operator=(const ConstraintPreservingBjorhus&);
168 : /// \endcond
169 0 : ~ConstraintPreservingBjorhus() override = default;
170 :
171 0 : explicit ConstraintPreservingBjorhus(CkMigrateMessage* msg);
172 :
173 0 : WRAPPED_PUPable_decl_base_template(
174 : domain::BoundaryConditions::BoundaryCondition,
175 : ConstraintPreservingBjorhus);
176 :
177 0 : auto get_clone() const -> std::unique_ptr<
178 : domain::BoundaryConditions::BoundaryCondition> override;
179 :
180 0 : static constexpr evolution::BoundaryConditions::Type bc_type =
181 : evolution::BoundaryConditions::Type::TimeDerivative;
182 :
183 0 : void pup(PUP::er& p) override;
184 :
185 0 : using dg_interior_evolved_variables_tags =
186 : tmpl::list<gr::Tags::SpacetimeMetric<DataVector, Dim>,
187 : Tags::Pi<DataVector, Dim>, Tags::Phi<DataVector, Dim>>;
188 0 : using dg_interior_temporary_tags =
189 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
190 : Tags::ConstraintGamma1, Tags::ConstraintGamma2,
191 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, Dim>,
192 : gr::Tags::InverseSpacetimeMetric<DataVector, Dim>,
193 : gr::Tags::SpacetimeNormalVector<DataVector, Dim>,
194 : Tags::ThreeIndexConstraint<DataVector, Dim>,
195 : Tags::GaugeH<DataVector, Dim>,
196 : Tags::SpacetimeDerivGaugeH<DataVector, Dim>>;
197 0 : using dg_interior_dt_vars_tags =
198 : tmpl::list<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, Dim>>,
199 : ::Tags::dt<Tags::Pi<DataVector, Dim>>,
200 : ::Tags::dt<Tags::Phi<DataVector, Dim>>>;
201 0 : using dg_interior_deriv_vars_tags =
202 : tmpl::list<::Tags::deriv<gr::Tags::SpacetimeMetric<DataVector, Dim>,
203 : tmpl::size_t<Dim>, Frame::Inertial>,
204 : ::Tags::deriv<Tags::Pi<DataVector, Dim>, tmpl::size_t<Dim>,
205 : Frame::Inertial>,
206 : ::Tags::deriv<Tags::Phi<DataVector, Dim>, tmpl::size_t<Dim>,
207 : Frame::Inertial>>;
208 0 : using dg_gridless_tags = tmpl::list<::Tags::Time>;
209 :
210 0 : std::optional<std::string> dg_time_derivative(
211 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
212 : dt_spacetime_metric_correction,
213 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
214 : dt_pi_correction,
215 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
216 : dt_phi_correction,
217 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
218 : face_mesh_velocity,
219 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
220 : const tnsr::I<DataVector, Dim, Frame::Inertial>& /*normal_vector*/,
221 : // c.f. dg_interior_evolved_variables_tags
222 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
223 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
224 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
225 : // c.f. dg_interior_temporary_tags
226 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
227 : const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
228 : const Scalar<DataVector>& lapse,
229 : const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
230 : const tnsr::AA<DataVector, Dim, Frame::Inertial>&
231 : inverse_spacetime_metric,
232 : const tnsr::A<DataVector, Dim, Frame::Inertial>&
233 : spacetime_unit_normal_vector,
234 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& three_index_constraint,
235 : const tnsr::a<DataVector, Dim, Frame::Inertial>& gauge_source,
236 : const tnsr::ab<DataVector, Dim, Frame::Inertial>&
237 : spacetime_deriv_gauge_source,
238 : // c.f. dg_interior_dt_vars_tags
239 : const tnsr::aa<DataVector, Dim, Frame::Inertial>&
240 : logical_dt_spacetime_metric,
241 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& logical_dt_pi,
242 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& logical_dt_phi,
243 : // c.f. dg_interior_deriv_vars_tags
244 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_spacetime_metric,
245 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_pi,
246 : const tnsr::ijaa<DataVector, Dim, Frame::Inertial>& d_phi,
247 : double time = std::numeric_limits<double>::signaling_NaN()) const;
248 :
249 : private:
250 0 : void compute_intermediate_vars(
251 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
252 : unit_interface_normal_vector,
253 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
254 : four_index_constraint,
255 : gsl::not_null<tnsr::II<DataVector, Dim, Frame::Inertial>*>
256 : inverse_spatial_metric,
257 : gsl::not_null<tnsr::ii<DataVector, Dim, Frame::Inertial>*>
258 : extrinsic_curvature,
259 : gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
260 : incoming_null_one_form,
261 : gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
262 : outgoing_null_one_form,
263 : gsl::not_null<tnsr::A<DataVector, Dim, Frame::Inertial>*>
264 : incoming_null_vector,
265 : gsl::not_null<tnsr::A<DataVector, Dim, Frame::Inertial>*>
266 : outgoing_null_vector,
267 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*> projection_ab,
268 : gsl::not_null<tnsr::Ab<DataVector, Dim, Frame::Inertial>*> projection_Ab,
269 : gsl::not_null<tnsr::AA<DataVector, Dim, Frame::Inertial>*> projection_AB,
270 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
271 : char_projected_rhs_dt_v_psi,
272 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
273 : char_projected_rhs_dt_v_zero,
274 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
275 : char_projected_rhs_dt_v_plus,
276 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
277 : char_projected_rhs_dt_v_minus,
278 : gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
279 : constraint_char_zero_plus,
280 : gsl::not_null<tnsr::a<DataVector, Dim, Frame::Inertial>*>
281 : constraint_char_zero_minus,
282 : gsl::not_null<std::array<DataVector, 4>*> char_speeds,
283 :
284 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
285 : face_mesh_velocity,
286 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
287 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
288 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
289 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
290 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
291 : const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
292 : const Scalar<DataVector>& lapse,
293 : const tnsr::I<DataVector, Dim, Frame::Inertial>& shift,
294 : const tnsr::AA<DataVector, Dim, Frame::Inertial>&
295 : inverse_spacetime_metric,
296 : const tnsr::A<DataVector, Dim, Frame::Inertial>&
297 : spacetime_unit_normal_vector,
298 : const tnsr::a<DataVector, Dim, Frame::Inertial>&
299 : spacetime_unit_normal_one_form,
300 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& three_index_constraint,
301 : const tnsr::a<DataVector, Dim, Frame::Inertial>& gauge_source,
302 : const tnsr::ab<DataVector, Dim, Frame::Inertial>&
303 : spacetime_deriv_gauge_source,
304 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_pi,
305 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& dt_phi,
306 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& dt_spacetime_metric,
307 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_pi,
308 : const tnsr::ijaa<DataVector, Dim, Frame::Inertial>& d_phi,
309 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& d_spacetime_metric)
310 : const;
311 :
312 0 : detail::ConstraintPreservingBjorhusType type_{
313 : detail::ConstraintPreservingBjorhusType::ConstraintPreservingPhysical};
314 0 : std::unique_ptr<::MathFunction<1, Frame::Inertial>> incoming_wave_profile_{};
315 : };
316 : } // namespace gh::BoundaryConditions
317 :
318 : template <>
319 : struct Options::create_from_yaml<
320 : gh::BoundaryConditions::detail::ConstraintPreservingBjorhusType> {
321 : template <typename Metavariables>
322 : static
323 : typename gh::BoundaryConditions::detail::ConstraintPreservingBjorhusType
324 : create(const Options::Option& options) {
325 : return gh::BoundaryConditions::detail::
326 : convert_constraint_preserving_bjorhus_type_from_yaml(options);
327 : }
328 : };
|
