Line data Source code
1 0 : // 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 :
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Evolution/BoundaryConditions/Type.hpp"
15 : #include "Evolution/Systems/CurvedScalarWave/BoundaryConditions/BoundaryCondition.hpp"
16 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
18 : #include "Options/String.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
20 : #include "Utilities/Gsl.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : namespace domain::Tags {
26 : template <size_t Dim, typename Frame>
27 : struct Coordinates;
28 : } // namespace domain::Tags
29 : /// \endcond
30 :
31 : namespace CurvedScalarWave::BoundaryConditions {
32 :
33 : /*!
34 : * \brief Implements constraint-preserving boundary conditions with a second
35 : * order Bayliss-Turkel radiation boundary condition.
36 : *
37 : * \details The Bayliss-Turkel boundary conditions are technically only valid in
38 : * flat space and should therefore only be used at boundaries where the
39 : * background spacetime is approximately Minkwoski such as (sufficiently far
40 : * out) outer boundaries for asymptotically flat spacetimes. Small reflections
41 : * are still likely to occur.
42 : *
43 : * The constraint-preserving part of the boundary conditions are set on the time
44 : * derivatives of all evolved fields. The physical Bayliss-Turkel boundary
45 : * conditions are additionally set onto the time derivative of \f$\Pi\f$.
46 : *
47 : * The constraints are defined as follows:
48 : *
49 : * \f{align*}{
50 : * \mathcal{C}_i&=\partial_i\Psi - \Phi_i=0, \\
51 : * \mathcal{C}_{ij}&=\partial_{[i}\Phi_{j]}=0
52 : * \f}
53 : * Inspection of the constraint evolution system (Eqs. 29-30 in
54 : * \cite Holst2004wt) shows that the constraints themselves are characteristic
55 : * fields. We can derive constraint boundary conditions the same way
56 : * \cite Kidder2004rw does for the Einstein equations:
57 : *
58 : * We express the constraints in terms of (evolution) characeristic fields and
59 : * demand that the normal component of the constraint has to be zero when
60 : * flowing into the boundary i.e. there are no constraints flowing into our
61 : * numerical domain:
62 : *
63 : * \f{align*}{
64 : * 0 &= n^i \mathcal{C}_i &= n^i \partial_i w^\Psi - \frac{1}{2}(w^{+} - w^-) +
65 : * n^i
66 : * w_i^0 \\
67 : * (n^i \partial_i w^\Psi)_{BC} &= \frac{1}{2}(w^{+} - w^-) - n^i w_i^0
68 : * \f}
69 : *
70 : * and
71 : *
72 : * \f{align*}{
73 : * 0 &= 2 n^i \mathcal{C}_{ij} = n^i \partial_i w^0_j + \frac{1}{2}n^i n_j
74 : * (\partial_i w^+ - \partial_i w^-) - \frac{1}{2}(\partial_j w^+ - \partial_j
75 : * w^-) - n^i \partial_j w^0_i \\
76 : * (n^i \partial_i w^0_j)_{BC} &= - \frac{1}{2}n^i
77 : * n_j (\partial_i w^+ - \partial_i w^-) + \frac{1}{2}(\partial_j w^+ +
78 : * \partial_j w^-) + n^i \partial_j w^0_i \f}
79 : *
80 : * This condition is applied to the time derivative using the Bjorhus condition
81 : * \cite Bjorhus1995 :
82 : * \f{align*}{
83 : * \partial_t u^\alpha + A^{i \alpha}_\beta \partial_i u^\beta &= F^\alpha \\
84 : * e^{\hat{\alpha}}_\alpha (\partial_t u^\alpha + A^{i \alpha}_\beta
85 : * \partial_i u^\beta) &= e^{\hat{\alpha}}_\alpha F^\alpha \\
86 : * d_t u^{\hat{\alpha}} + e^{\hat{\alpha}}_\alpha A^{i
87 : * \alpha}_\beta(P^k_i + n^k n_i) \partial_k u^\beta &= e^{\hat{\alpha}}_\alpha
88 : * F^\alpha \\
89 : * d_t u^{\hat{\alpha}} + \lambda_{(\hat{\alpha})} n^k \partial_k
90 : * u^{\hat{\alpha}} + e^{\hat{\alpha}}_\alpha A^{i \alpha}_\beta P^k_i
91 : * \partial_k u^\beta &= e^{\hat{\alpha}}_\alpha F^\alpha
92 : * \f}
93 : *
94 : * Defining the volume time derivative of the characteristic fields as:
95 : * \f{equation*}{
96 : * D_t u^{\hat{\alpha}} \equiv e^{\hat{\alpha}}_\alpha (- A^{i \alpha}_\beta
97 : * \partial_i u^\beta + F^\alpha)
98 : * \f}
99 : *
100 : * The boundary conditions are now formulated as follows:
101 : *
102 : * \f{equation*}{
103 : * d_t u^{\hat{\alpha}} = D_t u^{\hat{\alpha}} + \lambda_{(\hat{\alpha})}
104 : * (n^i\partial_i u^{\hat{\alpha}} - (n^i\partial_i u^{\hat{\alpha}})_{BC})
105 : * \f}
106 : *
107 : * Using the condition that there are no incoming constraint fields, this gives:
108 : *
109 : * \f{align*}{
110 : * d_t Z^1 &= D_t w^\Psi + \lambda_\Psi n^i \mathcal{C}_i \\
111 : * d_t Z^2_i &= D_t w^0_i + 2 \lambda_0 n^i \mathcal{C}_{ij}
112 : * \f}
113 : *
114 : * The Bayliss-Turkel boundary conditions are given by:
115 : *
116 : * \f{align*}{
117 : * \prod_{l=1}^m\left(\partial_t + \partial_r + \frac{2l-1}{r}\right)\Psi=0,
118 : * \f}
119 : *
120 : * which we expand here to second order (\f$m=2\f$) to derive conditions for
121 : * \f$\partial_t\Pi^{\mathrm{BC}}\f$:
122 : *
123 : * \f{align*}{
124 : * \partial_t\Pi^{\mathrm{BC}}
125 : * &=\left(\partial_t\partial_r + \partial_r\partial_t +
126 : * \partial_r^2+\frac{4}{r}\partial_t
127 : * +\frac{4}{r}\partial_r + \frac{2}{r^2}\right)\Psi \notag \\
128 : * &=\left((2n^i + \beta^i) \partial_t \Phi_i + n^i n^j\partial_i\Phi_j +
129 : * \frac{4}{r}\partial_t\Psi + \frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi
130 : * \right) / \alpha.
131 : * \f}
132 : *
133 : *
134 : * This derivation makes the following assumptions:
135 : *
136 : * - The lapse, shift, normal vector and radius are time-independent,
137 : * \f$\partial_t \alpha = \partial_t \beta^i = \partial_t n^i = \partial_t r =
138 : * 0\f$. If necessary, these time derivatives can be accounted for in the future
139 : * by inserting the appropriate terms in a straightforward manner.
140 : *
141 : * - The outer boundary is spherical. It might be possible to generalize this
142 : * condition but we have not tried this.
143 : *
144 : * * The boundary conditions to the time derivative of the evolved variables
145 : * are then given by:
146 : *
147 : * The full boundary conditions, as applied to the time derivative of each
148 : * evolved field are then given by:
149 : * \f{align*}{ \partial_{t}
150 : * \Psi&\to\partial_{t}\Psi +
151 : * \lambda_\Psi n^i \mathcal{C}_i, \\
152 : * \partial_{t}\Pi&\to\partial_{t}\Pi-\left(\partial_t\Pi
153 : * - \partial_t\Pi^{\mathrm{BC}}\right)
154 : * +\gamma_2\lambda_\Psi n^i \mathcal{C}_i
155 : * =\partial_t\Pi^{\mathrm{BC}}
156 : * +\gamma_2\lambda_\Psi n^i \mathcal{C}_i, \\
157 : * \partial_{t}\Phi_i&\to\partial_{t}\Phi_i+ 2 \lambda_0 n^j \mathcal{C}_{ji}.
158 : * \f}
159 : *
160 : */
161 : template <size_t Dim>
162 1 : class ConstraintPreservingSphericalRadiation final
163 : : public BoundaryCondition<Dim> {
164 : public:
165 0 : using options = tmpl::list<>;
166 0 : static constexpr Options::String help{
167 : "Constraint-preserving boundary conditions with a second order "
168 : "Bayliss-Turkel radiation boundary condition."};
169 0 : ConstraintPreservingSphericalRadiation() = default;
170 0 : ConstraintPreservingSphericalRadiation(
171 : ConstraintPreservingSphericalRadiation&&) = default;
172 0 : ConstraintPreservingSphericalRadiation& operator=(
173 : ConstraintPreservingSphericalRadiation&&) = default;
174 0 : ConstraintPreservingSphericalRadiation(
175 : const ConstraintPreservingSphericalRadiation&) = default;
176 0 : ConstraintPreservingSphericalRadiation& operator=(
177 : const ConstraintPreservingSphericalRadiation&) = default;
178 0 : ~ConstraintPreservingSphericalRadiation() override = default;
179 :
180 0 : explicit ConstraintPreservingSphericalRadiation(CkMigrateMessage* msg);
181 :
182 0 : WRAPPED_PUPable_decl_base_template(
183 : domain::BoundaryConditions::BoundaryCondition,
184 : ConstraintPreservingSphericalRadiation);
185 :
186 0 : auto get_clone() const -> std::unique_ptr<
187 : domain::BoundaryConditions::BoundaryCondition> override;
188 :
189 0 : static constexpr evolution::BoundaryConditions::Type bc_type =
190 : evolution::BoundaryConditions::Type::TimeDerivative;
191 :
192 0 : void pup(PUP::er& p) override;
193 :
194 0 : using dg_interior_evolved_variables_tags =
195 : tmpl::list<Tags::Psi, Tags::Phi<Dim>>;
196 0 : using dg_interior_temporary_tags =
197 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
198 : Tags::ConstraintGamma1, Tags::ConstraintGamma2,
199 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, Dim>>;
200 0 : using dg_interior_dt_vars_tags =
201 : tmpl::list<::Tags::dt<Tags::Psi>, ::Tags::dt<Tags::Pi>,
202 : ::Tags::dt<Tags::Phi<Dim>>>;
203 0 : using dg_interior_deriv_vars_tags = tmpl::list<
204 : ::Tags::deriv<Tags::Psi, tmpl::size_t<Dim>, Frame::Inertial>,
205 : ::Tags::deriv<Tags::Pi, tmpl::size_t<Dim>, Frame::Inertial>,
206 : ::Tags::deriv<Tags::Phi<Dim>, tmpl::size_t<Dim>, Frame::Inertial>>;
207 0 : using dg_gridless_tags = tmpl::list<>;
208 :
209 0 : std::optional<std::string> dg_time_derivative(
210 : gsl::not_null<Scalar<DataVector>*> dt_psi_correction,
211 : gsl::not_null<Scalar<DataVector>*> dt_pi_correction,
212 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
213 : dt_phi_correction,
214 : const std::optional<tnsr::I<DataVector, Dim>>& face_mesh_velocity,
215 : const tnsr::i<DataVector, Dim>& normal_covector,
216 : const tnsr::I<DataVector, Dim>& normal_vector,
217 : const Scalar<DataVector>& psi, const tnsr::i<DataVector, Dim>& phi,
218 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
219 : const Scalar<DataVector>& gamma1, const Scalar<DataVector>& gamma2,
220 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
221 : const Scalar<DataVector>& logical_dt_psi,
222 : const Scalar<DataVector>& logical_dt_pi,
223 : const tnsr::i<DataVector, Dim>& logical_dt_phi,
224 : const tnsr::i<DataVector, Dim>& d_psi,
225 : const tnsr::i<DataVector, Dim>& d_pi,
226 : const tnsr::ij<DataVector, Dim>& d_phi) const;
227 : };
228 : } // namespace CurvedScalarWave::BoundaryConditions
|
