Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/DataBox/Tag.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Domain/Tags.hpp"
12 : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Gauges.hpp"
13 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
14 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
15 : #include "Time/Tags.hpp"
16 : #include "Utilities/ContainerHelpers.hpp"
17 : #include "Utilities/TMPL.hpp"
18 :
19 : /// \cond
20 : class DataVector;
21 : namespace gsl {
22 : template <class T>
23 : class not_null;
24 : } // namespace gsl
25 : /// \endcond
26 :
27 : namespace gh {
28 : namespace gauges {
29 : /*!
30 : * \brief Damped harmonic gauge source function and its spacetime derivative.
31 : *
32 : * \details The gauge condition has been taken from \cite Szilagyi2009qz and
33 : * \cite Deppe2018uye. We provide both a "rollon" version
34 : * (`damped_harmonic_rollon`), and a "non-rollon" version (`damped_harmonic`).
35 : * In the non-rollon version the rollon function \f$R(t)=1\f$.
36 : *
37 : * \warning Only the non-rollon version can be used with a moving mesh.
38 : *
39 : * The covariant form of the source function \f$H_a\f$ is written as:
40 : *
41 : * \f{align*}
42 : * H_a := [1 - R(t)] H_a^\mathrm{init} +
43 : * [\mu_{L1} \mathrm{log}(\sqrt{g}/N) + \mu_{L2} \mathrm{log}(1/N)] t_a
44 : * - \mu_S g_{ai} N^i / N
45 : * \f}
46 : *
47 : * where \f$N, N^k\f$ are the lapse and shift respectively, \f$t_a\f$ is the
48 : * unit normal one-form to the spatial slice, and \f$g_{ab}\f$ is
49 : * the spatial metric (obtained by projecting the spacetime metric onto the
50 : * 3-slice, i.e. \f$g_{ab} = \psi_{ab} + t_a t_b\f$). The prefactors are:
51 : *
52 : * \f{align*}
53 : * \mu_{L1} &= A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{L1}}, \\
54 : * \mu_{L2} &= A_{L2} R(t) W(x^i) \mathrm{log}(1/N)^{e_{L2}}, \\
55 : * \mu_{S} &= A_{S} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{S}},
56 : * \f}
57 : *
58 : * temporal roll-on function \f$ R(t)\f$ is:
59 : *
60 : * \f{align*}
61 : * \begin{array}{ll}
62 : * R(t) & = 0, & t< t_{0} \\
63 : * & = 1 - \exp[-((t - t_{0})/ \sigma_t)^4], & t\geq t_{0} \\
64 : * \end{array}
65 : * \f}
66 : *
67 : * and the spatial weight function is:
68 : *
69 : * \f{align*}
70 : * W(x^i) = \exp[-(r/\sigma_r)^2].
71 : * \f}
72 : *
73 : * This weight function can be written with multiple constant factors in the
74 : * exponent in literature \cite Deppe2018uye, but we absorb them all into
75 : * \f$ \sigma_r\f$ here. The coordinate \f$ r\f$ is the Euclidean radius
76 : * in Inertial coordinates.
77 : *
78 : * Note that for the last three terms in \f$H_a\f$ (with \f$ X = \{L1, L2, S\}
79 : * \f$):
80 : * - Amplitude factors \f$ A_{X} \f$ are taken as input here as `amp_coef_X`
81 : * - Exponents \f$ e_X\f$ are taken as input here as `exp_X`.
82 : * - Spatial weight function \f$W\f$ is specified completely by
83 : * \f$\sigma_r\f$, which is taken as input here as `sigma_r`.
84 : *
85 : * Also computes spacetime derivatives, i.e. \f$\partial_a H_b\f$, of the damped
86 : * harmonic source function H. Using notation from damped_harmonic_h(), we
87 : * rewrite the same as:
88 : *
89 : * \f{align*}
90 : * \partial_a H_b =& \partial_a T_1 + \partial_a T_2 + \partial_a T_3, \\
91 : * H_a =& T_1 + T_2 + T_3,
92 : * \f}
93 : *
94 : * where:
95 : *
96 : * \f{align*}
97 : * T_1 =& [1 - R(t)] H_a^\mathrm{init}, \\
98 : * T_2 =& [\mu_{L1} \mathrm{log}(\sqrt{g}/N) + \mu_{L2} \mathrm{log}(1/N)] t_a,
99 : * \\
100 : * T_3 =& - \mu_S g_{ai} N^i / N.
101 : * \f}
102 : *
103 : * Derivation:
104 : *
105 : * \f$\blacksquare\f$ For \f$ T_1 \f$, the derivatives are:
106 : * \f{align*}
107 : * \partial_a T_1 = (1 - R(t))
108 : * \partial_a H_b^\mathrm{init}
109 : * - H_b^\mathrm{init} \partial_a R.
110 : * \f}
111 : *
112 : * \f$\blacksquare\f$ Write \f$ T_2 \equiv (\mu_1 + \mu_2) t_b \f$. Then:
113 : *
114 : * \f{align*}
115 : * \partial_a T_2 =& (\partial_a \mu_1 + \partial_a \mu_2) t_b \\
116 : * +& (\mu_1 + \mu_2) \partial_a t_b,
117 : * \f}
118 : *
119 : * where
120 : *
121 : * \f{align*}
122 : * \partial_a t_b =& \left(-\partial_a N, 0, 0, 0\right) \\
123 : *
124 : * \partial_a \mu_1
125 : * =& \partial_a [A_{L1} R(t) W(x^i) \mathrm{log}(\sqrt{g}/N)^{e_{L1} +
126 : * 1}], \\
127 : * =& A_{L1} R(t) W(x^i) \partial_a [\mathrm{log}(\sqrt{g}/N)^{e_{L1} +
128 : * 1}] \\
129 : * +& A_{L1} \mathrm{log}(\sqrt{g}/N)^{e_{L1} + 1} \partial_a [R(t)
130 : * W(x^i)],\\
131 : *
132 : * \partial_a \mu_2
133 : * =& \partial_a [A_{L2} R(t) W(x^i) \mathrm{log}(1/N)^{e_{L2} + 1}], \\
134 : * =& A_{L2} R(t) W(x^i) \partial_a [\mathrm{log}(1/N)^{e_{L2} + 1}] \\
135 : * +& A_{L2} \mathrm{log}(1/N)^{e_{L2} + 1} \partial_a [R(t) W(x^i)],
136 : * \f}
137 : *
138 : * where \f$\partial_a [R W] = \left(\partial_0 R(t), \partial_i
139 : * W(x^j)\right)\f$.
140 : *
141 : * \f$\blacksquare\f$ Finally, the derivatives of \f$ T_3 \f$ are:
142 : *
143 : * \f[
144 : * \partial_a T_3 = -\partial_a(\mu_S/N) g_{bi} N^i
145 : * -(\mu_S/N) \partial_a(g_{bi}) N^i
146 : * -(\mu_S/N) g_{bi}\partial_a N^i,
147 : * \f]
148 : *
149 : * where
150 : *
151 : * \f{align*}
152 : * \partial_a(\mu_S / N) =& (1/N)\partial_a \mu_S
153 : * - \frac{\mu_S}{N^2}\partial_a N, \,\,\mathrm{and}\\
154 : * \partial_a \mu_S =& \partial_a [A_S R(t) W(x^i)
155 : * \mathrm{log}(\sqrt{g}/N)^{e_S}], \\
156 : * =& A_S R(t) W(x^i) \partial_a
157 : * [\mathrm{log}(\sqrt{g}/N)^{e_S}] \\
158 : * +& A_S \mathrm{log}(\sqrt{g} / N)^{e_S} \partial_a [R(t)
159 : * W(x^i)].
160 : * \f}
161 : */
162 : template <size_t SpatialDim, typename Frame>
163 1 : void damped_harmonic_rollon(
164 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
165 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
166 : const tnsr::a<DataVector, SpatialDim, Frame>& gauge_h_init,
167 : const tnsr::ab<DataVector, SpatialDim, Frame>& dgauge_h_init,
168 : const Scalar<DataVector>& lapse,
169 : const tnsr::I<DataVector, SpatialDim, Frame>& shift,
170 : const tnsr::a<DataVector, SpatialDim, Frame>&
171 : spacetime_unit_normal_one_form,
172 : const tnsr::A<DataVector, SpatialDim, Frame>& spacetime_unit_normal,
173 : const Scalar<DataVector>& sqrt_det_spatial_metric,
174 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
175 : const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
176 : const Scalar<DataVector>& half_pi_two_normals,
177 : const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
178 : const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
179 : const tnsr::aa<DataVector, SpatialDim, Frame>& pi,
180 : const tnsr::iaa<DataVector, SpatialDim, Frame>& phi, double time,
181 : const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
182 : double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
183 : double rollon_start_time, double rollon_width, double sigma_r);
184 :
185 : /*!
186 : * \copydoc damped_harmonic_rollon()
187 : */
188 : template <size_t SpatialDim, typename Frame>
189 1 : void damped_harmonic(
190 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*> gauge_h,
191 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame>*> d4_gauge_h,
192 : const Scalar<DataVector>& lapse,
193 : const tnsr::I<DataVector, SpatialDim, Frame>& shift,
194 : const tnsr::a<DataVector, SpatialDim, Frame>&
195 : spacetime_unit_normal_one_form,
196 : const tnsr::A<DataVector, SpatialDim, Frame>& spacetime_unit_normal,
197 : const Scalar<DataVector>& sqrt_det_spatial_metric,
198 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
199 : const tnsr::abb<DataVector, SpatialDim, Frame>& d4_spacetime_metric,
200 : const Scalar<DataVector>& half_pi_two_normals,
201 : const tnsr::i<DataVector, SpatialDim, Frame>& half_phi_two_normals,
202 : const tnsr::aa<DataVector, SpatialDim, Frame>& spacetime_metric,
203 : const tnsr::aa<DataVector, SpatialDim, Frame>& pi,
204 : const tnsr::iaa<DataVector, SpatialDim, Frame>& phi,
205 : const tnsr::I<DataVector, SpatialDim, Frame>& coords, double amp_coef_L1,
206 : double amp_coef_L2, double amp_coef_S, int exp_L1, int exp_L2, int exp_S,
207 : double sigma_r);
208 :
209 : /*!
210 : * \brief Impose damped harmonic gauge.
211 : *
212 : * \see `damped_harmonic()`
213 : */
214 1 : class DampedHarmonic final : public GaugeCondition {
215 : public:
216 : /// The width of the Gaussian for the spatial decay of the damped harmonic
217 : /// gauge.
218 1 : struct SpatialDecayWidth {
219 0 : using type = double;
220 0 : static constexpr Options::String help{
221 : "Spatial width (sigma_r) of weight function (W(x^i)) used in the "
222 : "damped harmonic gauge."};
223 : };
224 : /// The amplitudes for the L1, L2, and S terms, respectively, for the damped
225 : /// harmonic gauge.
226 1 : struct Amplitudes {
227 0 : using type = std::array<double, 3>;
228 0 : static constexpr Options::String help{
229 : "Amplitudes [A_{L1}, A_{L2}, A_{S}] for the damped harmonic gauge."};
230 : };
231 : /// The exponents for the L1, L2, and S terms, respectively, for the damped
232 : /// harmonic gauge.
233 1 : struct Exponents {
234 0 : using type = std::array<int, 3>;
235 0 : static constexpr Options::String help{
236 : "Exponents [e_{L1}, e_{L2}, e_{S}] for the damped harmonic gauge."};
237 : };
238 :
239 0 : static constexpr Options::String help{
240 : "Apply damped harmonic/damped wave gauge."};
241 :
242 0 : using options = tmpl::list<SpatialDecayWidth, Amplitudes, Exponents>;
243 :
244 0 : DampedHarmonic(double width, const std::array<double, 3>& amps,
245 : const std::array<int, 3>& exps);
246 :
247 0 : DampedHarmonic() = default;
248 0 : DampedHarmonic(const DampedHarmonic&) = default;
249 0 : DampedHarmonic& operator=(const DampedHarmonic&) = default;
250 0 : DampedHarmonic(DampedHarmonic&&) = default;
251 0 : DampedHarmonic& operator=(DampedHarmonic&&) = default;
252 0 : ~DampedHarmonic() override = default;
253 :
254 : /// \cond
255 : explicit DampedHarmonic(CkMigrateMessage* msg);
256 : using PUP::able::register_constructor;
257 : WRAPPED_PUPable_decl_template(DampedHarmonic); // NOLINT
258 : /// \endcond
259 :
260 0 : std::unique_ptr<GaugeCondition> get_clone() const override;
261 :
262 : template <size_t SpatialDim>
263 0 : void gauge_and_spacetime_derivative(
264 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame::Inertial>*> gauge_h,
265 : gsl::not_null<tnsr::ab<DataVector, SpatialDim, Frame::Inertial>*>
266 : d4_gauge_h,
267 : const Scalar<DataVector>& lapse,
268 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
269 : const tnsr::a<DataVector, SpatialDim, Frame::Inertial>&
270 : spacetime_unit_normal_one_form,
271 : const tnsr::A<DataVector, SpatialDim, Frame::Inertial>&
272 : spacetime_unit_normal,
273 : const Scalar<DataVector>& sqrt_det_spatial_metric,
274 : const tnsr::II<DataVector, SpatialDim, Frame::Inertial>&
275 : inverse_spatial_metric,
276 : const tnsr::abb<DataVector, SpatialDim, Frame::Inertial>&
277 : d4_spacetime_metric,
278 : const Scalar<DataVector>& half_pi_two_normals,
279 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
280 : half_phi_two_normals,
281 : const tnsr::aa<DataVector, SpatialDim, Frame::Inertial>& spacetime_metric,
282 : const tnsr::aa<DataVector, SpatialDim, Frame::Inertial>& pi,
283 : const tnsr::iaa<DataVector, SpatialDim, Frame::Inertial>& phi,
284 : double time,
285 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& inertial_coords)
286 : const;
287 :
288 : // NOLINTNEXTLINE(google-runtime-references)
289 0 : void pup(PUP::er& p) override;
290 :
291 : private:
292 0 : double spatial_decay_width_{std::numeric_limits<double>::signaling_NaN()};
293 0 : std::array<double, 3> amplitudes_{
294 : {std::numeric_limits<double>::signaling_NaN(),
295 : std::numeric_limits<double>::signaling_NaN(),
296 : std::numeric_limits<double>::signaling_NaN()}};
297 0 : std::array<int, 3> exponents_{{std::numeric_limits<int>::max(),
298 : std::numeric_limits<int>::max(),
299 : std::numeric_limits<int>::max()}};
300 : };
301 :
302 : namespace DampedHarmonicGauge_detail {
303 : // Used in the test
304 : double roll_on_function(double time, double t_start, double sigma_t);
305 :
306 : double time_deriv_of_roll_on_function(double time, double t_start,
307 : double sigma_t);
308 : } // namespace DampedHarmonicGauge_detail
309 : } // namespace gauges
310 : } // namespace gh
|
