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 : #include <type_traits>
11 :
12 : #include "DataStructures/DataBox/Prefixes.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "DataStructures/Variables.hpp"
16 : #include "Evolution/BoundaryConditions/Type.hpp"
17 : #include "Evolution/Systems/ScalarWave/BoundaryConditions/BoundaryCondition.hpp"
18 : #include "Evolution/Systems/ScalarWave/Tags.hpp"
19 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
20 : #include "Options/Options.hpp"
21 : #include "Options/String.hpp"
22 : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
23 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
24 : #include "Utilities/Gsl.hpp"
25 : #include "Utilities/Serialization/CharmPupable.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 :
35 : namespace ScalarWave::BoundaryConditions {
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+ 2 \lambda_0n^j \mathcal{C}_{ji}.
102 : * \f}
103 : *
104 : * Below we assume the normal vector \f$n^i\f$ is the radial unit normal vector.
105 : * That is, we assume the outer boundary is spherical. A Sommerfeld
106 : * \cite Sommerfeld1949 radiation condition is given by
107 : *
108 : * \f{align*}{
109 : * \partial_t\Psi=n^i\Phi_i
110 : * \f}
111 : *
112 : * Or, assuming that \f$\partial_tn^i=0\f$ (or is very small),
113 : *
114 : * \f{align*}{
115 : * \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i
116 : * \f}
117 : *
118 : * The Bayliss-Turkel \cite BaylissTurkel boundary conditions are given by:
119 : *
120 : * \f{align*}{
121 : * \prod_{l=1}^m\left(\partial_t + \partial_r + \frac{2l-1}{r}\right)\Psi=0
122 : * \f}
123 : *
124 : * The first-order form is
125 : *
126 : * \f{align*}{
127 : * \partial_t\Pi^{\mathrm{BC}}=n^i\partial_t\Phi_i + \frac{1}{r}\partial_t\Psi,
128 : * \f}
129 : *
130 : * assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
131 : *
132 : * The second-order boundary condition is given by,
133 : *
134 : * \f{align*}{
135 : * \partial_t\Pi^{\mathrm{BC}}
136 : * &=\left(\partial_t\partial_r + \partial_r\partial_t +
137 : * \partial_r^2+\frac{4}{r}\partial_t
138 : * +\frac{4}{r}\partial_r + \frac{2}{r^2}\right)\Psi \notag \\
139 : * &=n^i(\partial_t\Phi_i-\partial_i\Pi) + n^i n^j\partial_i\Phi_j -
140 : * \frac{4}{r}\Pi+\frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi \notag \\
141 : * &=n^i(\gamma_2(\partial_i\Psi - \Phi_i)-2\partial_i\Pi)
142 : * + n^i n^j\partial_i\Phi_j -
143 : * \frac{4}{r}\Pi+\frac{4}{r}n^i\Phi_i + \frac{2}{r^2}\Psi,
144 : * \f}
145 : *
146 : * where \f$\partial_t\f$ is the derivative with respect to the inertial frame
147 : * time, and we are assuming \f$\partial_t n^i=0\f$ and \f$\partial_t r=0\f$.
148 : *
149 : * The moving mesh can be accounted for by using
150 : *
151 : * \f{align*}{
152 : * \partial_t r = \frac{1}{r}x^i\delta_{ij}\partial_t x^j
153 : * \f}
154 : *
155 : * \note It is not clear if \f$\partial_t\Phi_i\f$ should be replaced by
156 : * \f$-\partial_i\Pi\f$, which is the evolution equation but without the
157 : * constraint.
158 : *
159 : * \note On a moving mesh the characteristic speeds change according to
160 : * \f$\lambda\to\lambda-v^i_gn_i\f$ where \f$v^i_g\f$ is the mesh velocity.
161 : *
162 : * \note For the scalar wave system \f$\lambda_0 = \lambda_\psi\f$
163 : *
164 : * \warning The boundary conditions are implemented assuming the outer boundary
165 : * is spherical. It might be possible to generalize the condition to
166 : * non-spherical boundaries by using \f$x^i/r\f$ instead of \f$n^i\f$, but this
167 : * hasn't been tested.
168 : *
169 : * \warning The received time derivatives are in the logical frame, not the
170 : * inertial frame and would need to first be corrected by subtracting the mesh
171 : * velocity. Similarly, the time derivative correction is in the logical frame.
172 : * We instead substitute the evolution equations directly in since the scalar
173 : * wave system is quite simple.
174 : */
175 : template <size_t Dim>
176 1 : class ConstraintPreservingSphericalRadiation final
177 : : public BoundaryCondition<Dim> {
178 : public:
179 0 : struct TypeOptionTag {
180 0 : using type = detail::ConstraintPreservingSphericalRadiationType;
181 0 : static std::string name() { return "Type"; }
182 0 : static constexpr Options::String help{
183 : "Whether to impose Sommerfeld, first-order Bayliss-Turkel, or "
184 : "second-order Bayliss-Turkel spherical radiation boundary conditions."};
185 : };
186 :
187 0 : using options = tmpl::list<TypeOptionTag>;
188 0 : static constexpr Options::String help{
189 : "Constraint-preserving spherical radiation boundary conditions setting "
190 : "the time derivatives of Psi, Phi, and Pi to avoid incoming constraint "
191 : "violations, and imposing radiation boundary conditions."};
192 :
193 0 : ConstraintPreservingSphericalRadiation(
194 : detail::ConstraintPreservingSphericalRadiationType type);
195 :
196 0 : ConstraintPreservingSphericalRadiation() = default;
197 : /// \cond
198 : ConstraintPreservingSphericalRadiation(
199 : ConstraintPreservingSphericalRadiation&&) = default;
200 : ConstraintPreservingSphericalRadiation& operator=(
201 : ConstraintPreservingSphericalRadiation&&) = default;
202 : ConstraintPreservingSphericalRadiation(
203 : const ConstraintPreservingSphericalRadiation&) = default;
204 : ConstraintPreservingSphericalRadiation& operator=(
205 : const ConstraintPreservingSphericalRadiation&) = default;
206 : /// \endcond
207 0 : ~ConstraintPreservingSphericalRadiation() override = default;
208 :
209 0 : explicit ConstraintPreservingSphericalRadiation(CkMigrateMessage* msg);
210 :
211 0 : WRAPPED_PUPable_decl_base_template(
212 : domain::BoundaryConditions::BoundaryCondition,
213 : ConstraintPreservingSphericalRadiation);
214 :
215 0 : auto get_clone() const -> std::unique_ptr<
216 : domain::BoundaryConditions::BoundaryCondition> override;
217 :
218 0 : static constexpr evolution::BoundaryConditions::Type bc_type =
219 : evolution::BoundaryConditions::Type::TimeDerivative;
220 :
221 0 : void pup(PUP::er& p) override;
222 :
223 0 : using dg_interior_evolved_variables_tags =
224 : tmpl::list<ScalarWave::Tags::Psi, ScalarWave::Tags::Pi,
225 : ScalarWave::Tags::Phi<Dim>>;
226 0 : using dg_interior_temporary_tags =
227 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
228 : Tags::ConstraintGamma2>;
229 : // using dg_interior_dt_vars_tags = tmpl::list<>;
230 0 : using dg_interior_deriv_vars_tags = tmpl::list<
231 : ::Tags::deriv<ScalarWave::Tags::Psi, tmpl::size_t<Dim>, Frame::Inertial>,
232 : ::Tags::deriv<ScalarWave::Tags::Pi, tmpl::size_t<Dim>, Frame::Inertial>,
233 : ::Tags::deriv<ScalarWave::Tags::Phi<Dim>, tmpl::size_t<Dim>,
234 : Frame::Inertial>>;
235 0 : using dg_gridless_tags = tmpl::list<>;
236 :
237 0 : std::optional<std::string> dg_time_derivative(
238 : gsl::not_null<Scalar<DataVector>*> dt_psi_correction,
239 : gsl::not_null<Scalar<DataVector>*> dt_pi_correction,
240 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
241 : dt_phi_correction,
242 :
243 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
244 : face_mesh_velocity,
245 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
246 :
247 : const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
248 : const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
249 :
250 : const tnsr::I<DataVector, Dim, Frame::Inertial>& coords,
251 : const Scalar<DataVector>& gamma2,
252 :
253 : const tnsr::i<DataVector, Dim, Frame::Inertial>& d_psi,
254 : const tnsr::i<DataVector, Dim, Frame::Inertial>& d_pi,
255 : const tnsr::ij<DataVector, Dim, Frame::Inertial>& d_phi) const;
256 :
257 : private:
258 0 : detail::ConstraintPreservingSphericalRadiationType type_{
259 : detail::ConstraintPreservingSphericalRadiationType::Sommerfeld};
260 : };
261 : } // namespace ScalarWave::BoundaryConditions
262 :
263 : template <>
264 : struct Options::create_from_yaml<
265 : ScalarWave::BoundaryConditions::detail::
266 : ConstraintPreservingSphericalRadiationType> {
267 : template <typename Metavariables>
268 : static typename ScalarWave::BoundaryConditions::detail::
269 : ConstraintPreservingSphericalRadiationType
270 : create(const Options::Option& options) {
271 : return ScalarWave::BoundaryConditions::detail::
272 : convert_constraint_preserving_spherical_radiation_type_from_yaml(
273 : options);
274 : }
275 : };
|
