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