Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 : #pragma once
4 :
5 : #include <cmath>
6 : #include <cstddef>
7 :
8 : #include "DataStructures/DataBox/AsAccess.hpp"
9 : #include "DataStructures/DataBox/DataBox.hpp"
10 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/TaggedContainers.hpp"
14 : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "DataStructures/VectorImpl.hpp"
18 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
19 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
20 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
21 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
22 : #include "Evolution/Systems/Ccz4/BoundaryConditions/BoundaryCondition.hpp"
23 : #include "Evolution/Systems/Ccz4/BoundaryConditions/Sommerfeld.hpp"
24 : #include "Evolution/Systems/Ccz4/Christoffel.hpp"
25 : #include "Evolution/Systems/Ccz4/DerivChristoffel.hpp"
26 : #include "Evolution/Systems/Ccz4/DerivLapse.hpp"
27 : #include "Evolution/Systems/Ccz4/DerivZ4Constraint.hpp"
28 : #include "Evolution/Systems/Ccz4/FiniteDifference/BoundaryConditionGhostData.hpp"
29 : #include "Evolution/Systems/Ccz4/FiniteDifference/Derivatives.hpp"
30 : #include "Evolution/Systems/Ccz4/FiniteDifference/Tags.hpp"
31 : #include "Evolution/Systems/Ccz4/Ricci.hpp"
32 : #include "Evolution/Systems/Ccz4/RicciScalarPlusDivergenceZ4Constraint.hpp"
33 : #include "Evolution/Systems/Ccz4/Tags.hpp"
34 : #include "Evolution/Systems/Ccz4/TagsDeclarations.hpp"
35 : #include "Evolution/Systems/Ccz4/TempTags.hpp"
36 : #include "Evolution/Systems/Ccz4/TimeDerivative.hpp"
37 : #include "Evolution/Systems/Ccz4/Z4Constraint.hpp"
38 : #include "Utilities/CallWithDynamicType.hpp"
39 : #include "Utilities/ConstantExpressions.hpp"
40 : #include "Utilities/ContainerHelpers.hpp"
41 : #include "Utilities/ErrorHandling/Assert.hpp"
42 : #include "Utilities/ErrorHandling/Error.hpp"
43 : #include "Utilities/Gsl.hpp"
44 : #include "Utilities/TMPL.hpp"
45 :
46 : /*!
47 : * \brief The namespace for evolving the second-order Ccz4 system.
48 : * Spatial derivatives are computed using 4-th order finite
49 : * differencing. Currently this system only works in 3D.
50 : */
51 : namespace Ccz4::fd {
52 0 : const size_t Dim = 3;
53 :
54 : namespace detail {
55 : // Calculate the time derivative of the evolved variables for the second-order
56 : // Ccz4 system. There is quite some overlap between this apply() funcion
57 : // and the apply() function in the first-order Ccz4 system. However,
58 : // it is not straightforward to directly reuse the first-order
59 : // apply() function as the function signatures differ significantly.
60 : template <size_t Dim>
61 : static void apply(
62 : // LHS time derivatives of evolved variables: eq 4(a) - 4(i)
63 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> dt_conformal_spatial_metric,
64 : const gsl::not_null<Scalar<DataVector>*> dt_lapse,
65 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_shift,
66 : const gsl::not_null<Scalar<DataVector>*> dt_conformal_factor,
67 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> dt_a_tilde,
68 : const gsl::not_null<Scalar<DataVector>*> dt_trace_extrinsic_curvature,
69 : const gsl::not_null<Scalar<DataVector>*> dt_theta,
70 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_gamma_hat,
71 : const gsl::not_null<tnsr::I<DataVector, Dim>*> dt_b,
72 :
73 : // quantities we need for computing eq 4, 13-27
74 : const gsl::not_null<Scalar<DataVector>*> conformal_factor_squared,
75 : const gsl::not_null<Scalar<DataVector>*> det_conformal_spatial_metric,
76 : const gsl::not_null<tnsr::II<DataVector, Dim>*>
77 : inv_conformal_spatial_metric,
78 : const gsl::not_null<tnsr::II<DataVector, Dim>*> inv_spatial_metric,
79 : const gsl::not_null<tnsr::II<DataVector, Dim>*> inv_a_tilde,
80 : // temporary expressions
81 : const gsl::not_null<tnsr::ij<DataVector, Dim>*> a_tilde_times_field_b,
82 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
83 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
84 : const gsl::not_null<Scalar<DataVector>*> contracted_field_b,
85 : const gsl::not_null<tnsr::ijK<DataVector, Dim>*> symmetrized_d_field_b,
86 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
87 : contracted_symmetrized_d_field_b,
88 : const gsl::not_null<tnsr::i<DataVector, Dim>*> field_d_up_times_a_tilde,
89 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
90 : contracted_field_d_up, // temp for eq 18 -20
91 : const gsl::not_null<Scalar<DataVector>*>
92 : half_conformal_factor_squared, // temp for eq 25
93 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
94 : conformal_metric_times_field_b,
95 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
96 : conformal_metric_times_trace_a_tilde,
97 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
98 : inv_conformal_metric_times_d_a_tilde,
99 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
100 : gamma_hat_minus_contracted_conformal_christoffel,
101 : const gsl::not_null<tnsr::iJ<DataVector, Dim>*>
102 : d_gamma_hat_minus_contracted_conformal_christoffel,
103 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
104 : contracted_christoffel_second_kind, // temp for eq 18 -20
105 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
106 : contracted_d_conformal_christoffel_difference, // temp for eq 18 -20
107 : const gsl::not_null<Scalar<DataVector>*> k_minus_2_theta_c,
108 : const gsl::not_null<Scalar<DataVector>*> k_minus_k0_minus_2_theta_c,
109 : const gsl::not_null<tnsr::ii<DataVector, Dim>*> lapse_times_a_tilde,
110 : const gsl::not_null<tnsr::ijj<DataVector, Dim>*> lapse_times_d_a_tilde,
111 : const gsl::not_null<tnsr::i<DataVector, Dim>*> lapse_times_field_a,
112 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
113 : lapse_times_conformal_spatial_metric,
114 : const gsl::not_null<Scalar<DataVector>*> lapse_times_slicing_condition,
115 : const gsl::not_null<Scalar<DataVector>*>
116 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
117 : const gsl::not_null<tnsr::I<DataVector, Dim>*> shift_times_deriv_gamma_hat,
118 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
119 : inv_tau_times_conformal_metric,
120 : // expressions and identities needed for evolution equations: eq 13 - 27
121 : const gsl::not_null<Scalar<DataVector>*> trace_a_tilde, // eq 13
122 : const gsl::not_null<tnsr::iJJ<DataVector, Dim>*> field_d_up, // eq 14
123 : const gsl::not_null<tnsr::Ijj<DataVector, Dim>*>
124 : conformal_christoffel_second_kind, // eq 15
125 : const gsl::not_null<tnsr::iJkk<DataVector, Dim>*>
126 : d_conformal_christoffel_second_kind, // eq 16
127 : const gsl::not_null<tnsr::Ijj<DataVector, Dim>*>
128 : christoffel_second_kind, // eq 17
129 : const gsl::not_null<tnsr::ii<DataVector, Dim>*>
130 : spatial_ricci_tensor, // eq 18 - 20
131 : const gsl::not_null<tnsr::ij<DataVector, Dim>*> grad_grad_lapse, // eq 21
132 : const gsl::not_null<Scalar<DataVector>*> divergence_lapse, // eq 22
133 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
134 : contracted_conformal_christoffel_second_kind, // eq 23
135 : const gsl::not_null<tnsr::iJ<DataVector, Dim>*>
136 : d_contracted_conformal_christoffel_second_kind, // eq 24
137 : const gsl::not_null<tnsr::i<DataVector, Dim>*>
138 : spatial_z4_constraint, // eq 25
139 : const gsl::not_null<tnsr::I<DataVector, Dim>*>
140 : upper_spatial_z4_constraint, // eq 25
141 : const gsl::not_null<tnsr::ij<DataVector, Dim>*>
142 : grad_spatial_z4_constraint, // eq 26
143 : const gsl::not_null<Scalar<DataVector>*>
144 : ricci_scalar_plus_divergence_z4_constraint, // eq 27
145 :
146 : // fixed params for SO-CCZ4
147 : // c = 1.0, cleaning_speed = 1.0
148 : // one_over_relaxation_time = 0.0
149 : const double c, const double cleaning_speed, // e in the paper
150 : const double one_over_relaxation_time, // \tau^{-1}
151 :
152 : // free params for SO-CCZ4
153 : const Scalar<DataVector>& eta, const double f, const double kappa_1,
154 : const double kappa_2, const double kappa_3, const Scalar<DataVector>& k_0,
155 :
156 : // evolved variables
157 : const tnsr::ii<DataVector, Dim>& conformal_spatial_metric,
158 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
159 : const Scalar<DataVector>& conformal_factor,
160 : const tnsr::ii<DataVector, Dim>& a_tilde,
161 : const Scalar<DataVector>& trace_extrinsic_curvature,
162 : const Scalar<DataVector>& theta, const tnsr::I<DataVector, Dim>& gamma_hat,
163 : const tnsr::I<DataVector, Dim>& b,
164 :
165 : // auxilliary fields and their derivatives
166 : const tnsr::i<DataVector, Dim>& field_a,
167 : const tnsr::iJ<DataVector, Dim>& field_b,
168 : const tnsr::ijj<DataVector, Dim>& field_d,
169 : const tnsr::i<DataVector, Dim>& field_p,
170 : const tnsr::ii<DataVector, Dim>& d_field_a,
171 : const tnsr::iiJ<DataVector, Dim>& d_field_b,
172 : const tnsr::iijj<DataVector, Dim>& d_field_d,
173 : const tnsr::ii<DataVector, Dim>& d_field_p,
174 :
175 : // spatial derivatives of other evolved variables
176 : const tnsr::ijj<DataVector, Dim>& d_a_tilde,
177 : const tnsr::i<DataVector, Dim>& d_trace_extrinsic_curvature,
178 : const tnsr::i<DataVector, Dim>& d_theta,
179 : const tnsr::iJ<DataVector, Dim>& d_gamma_hat,
180 : const tnsr::iJ<DataVector, Dim>& d_b, const bool shifting_shift,
181 : const bool evolve_lapse_and_shift) {
182 : constexpr double one_third = 1.0 / 3.0;
183 : // quantities we need for computing eq 4, 13 - 27
184 :
185 : determinant_and_inverse(det_conformal_spatial_metric,
186 : inv_conformal_spatial_metric,
187 : conformal_spatial_metric);
188 :
189 : get(*conformal_factor_squared) = square(get(conformal_factor));
190 :
191 : ::tenex::evaluate<ti::I, ti::J>(
192 : inv_spatial_metric, (*conformal_factor_squared)() *
193 : (*inv_conformal_spatial_metric)(ti::I, ti::J));
194 :
195 : ::tenex::evaluate<ti::I, ti::J>(
196 : inv_a_tilde, a_tilde(ti::k, ti::l) *
197 : (*inv_conformal_spatial_metric)(ti::I, ti::K) *
198 : (*inv_conformal_spatial_metric)(ti::J, ti::L));
199 :
200 : ASSERT(min(get(lapse)) > 0.0,
201 : "The lapse must be positive when using 1+log slicing but is: "
202 : << get(lapse));
203 :
204 : // slicing_condition and d_slicing_condition is not used in SO-CCZ4
205 : // \alpha g(\alpha) == 2
206 : *lapse_times_slicing_condition =
207 : make_with_value<Scalar<DataVector>>(lapse, 2.0);
208 :
209 : // expressions and identities needed for evolution equations: eq 13 - 27
210 :
211 : // eq 13
212 : ::tenex::evaluate(
213 : trace_a_tilde,
214 : (*inv_conformal_spatial_metric)(ti::I, ti::J) * a_tilde(ti::i, ti::j));
215 :
216 : // from eq 14: field_d_up is the D_k^{ij} tensor
217 : ::tenex::evaluate<ti::k, ti::I, ti::J>(
218 : field_d_up, (*inv_conformal_spatial_metric)(ti::I, ti::N) *
219 : (*inv_conformal_spatial_metric)(ti::M, ti::J) *
220 : field_d(ti::k, ti::n, ti::m));
221 :
222 : // eq 15
223 : ::Ccz4::conformal_christoffel_second_kind(conformal_christoffel_second_kind,
224 : *inv_conformal_spatial_metric,
225 : field_d);
226 :
227 : // eq 16
228 : ::Ccz4::deriv_conformal_christoffel_second_kind(
229 : d_conformal_christoffel_second_kind, *inv_conformal_spatial_metric,
230 : field_d, d_field_d, *field_d_up);
231 :
232 : // eq 17
233 : ::Ccz4::christoffel_second_kind(christoffel_second_kind,
234 : conformal_spatial_metric,
235 : *inv_conformal_spatial_metric, field_p,
236 : *conformal_christoffel_second_kind);
237 :
238 : // temporary expressions needed for eq 18 - 20
239 : ::tenex::evaluate<ti::l>(contracted_christoffel_second_kind,
240 : (*christoffel_second_kind)(ti::M, ti::l, ti::m));
241 :
242 : // comment for the future: we should probably ensure the traces are taken
243 : // before computing the differences as the off-diagonal terms are not
244 : // needed
245 : ::tenex::evaluate<ti::i, ti::j>(
246 : contracted_d_conformal_christoffel_difference,
247 : (*d_conformal_christoffel_second_kind)(ti::m, ti::M, ti::i, ti::j) -
248 : (*d_conformal_christoffel_second_kind)(ti::j, ti::M, ti::i, ti::m));
249 :
250 : ::tenex::evaluate<ti::L>(contracted_field_d_up,
251 : (*field_d_up)(ti::m, ti::M, ti::L));
252 :
253 : // eq 18 - 20
254 : ::Ccz4::spatial_ricci_tensor(
255 : spatial_ricci_tensor, *christoffel_second_kind,
256 : *contracted_christoffel_second_kind,
257 : *contracted_d_conformal_christoffel_difference, conformal_spatial_metric,
258 : *inv_conformal_spatial_metric, field_d, *field_d_up,
259 : *contracted_field_d_up, field_p, d_field_p);
260 :
261 : // eq 21
262 : ::Ccz4::grad_grad_lapse(grad_grad_lapse, lapse, *christoffel_second_kind,
263 : field_a, d_field_a);
264 :
265 : // eq 22
266 : ::Ccz4::divergence_lapse(divergence_lapse, *conformal_factor_squared,
267 : *inv_conformal_spatial_metric, *grad_grad_lapse);
268 :
269 : // eq 23
270 : ::Ccz4::contracted_conformal_christoffel_second_kind(
271 : contracted_conformal_christoffel_second_kind,
272 : *inv_conformal_spatial_metric, *conformal_christoffel_second_kind);
273 :
274 : // eq 24
275 : ::Ccz4::deriv_contracted_conformal_christoffel_second_kind(
276 : d_contracted_conformal_christoffel_second_kind,
277 : *inv_conformal_spatial_metric, *field_d_up,
278 : *conformal_christoffel_second_kind, *d_conformal_christoffel_second_kind);
279 :
280 : // temp for eq 25
281 : ::tenex::evaluate<ti::I>(
282 : gamma_hat_minus_contracted_conformal_christoffel,
283 : gamma_hat(ti::I) -
284 : (*contracted_conformal_christoffel_second_kind)(ti::I));
285 :
286 : // eq 25
287 : ::Ccz4::spatial_z4_constraint(
288 : spatial_z4_constraint, conformal_spatial_metric,
289 : *gamma_hat_minus_contracted_conformal_christoffel);
290 :
291 : // temp for eq 25
292 : ::tenex::evaluate(half_conformal_factor_squared,
293 : 0.5 * (*conformal_factor_squared)());
294 :
295 : // eq 25
296 : ::Ccz4::upper_spatial_z4_constraint(
297 : upper_spatial_z4_constraint, *half_conformal_factor_squared,
298 : *gamma_hat_minus_contracted_conformal_christoffel);
299 :
300 : // temp for eq 26
301 : ::tenex::evaluate<ti::i, ti::L>(
302 : d_gamma_hat_minus_contracted_conformal_christoffel,
303 : d_gamma_hat(ti::i, ti::L) -
304 : (*d_contracted_conformal_christoffel_second_kind)(ti::i, ti::L));
305 :
306 : // eq 26
307 : ::Ccz4::grad_spatial_z4_constraint(
308 : grad_spatial_z4_constraint, *spatial_z4_constraint,
309 : conformal_spatial_metric, *christoffel_second_kind, field_d,
310 : *gamma_hat_minus_contracted_conformal_christoffel,
311 : *d_gamma_hat_minus_contracted_conformal_christoffel);
312 :
313 : // eq 27
314 : ::Ccz4::ricci_scalar_plus_divergence_z4_constraint(
315 : ricci_scalar_plus_divergence_z4_constraint, *conformal_factor_squared,
316 : *inv_conformal_spatial_metric, *spatial_ricci_tensor,
317 : *grad_spatial_z4_constraint);
318 :
319 : // temporary expressions not already computed above
320 :
321 : ::tenex::evaluate(contracted_field_b, field_b(ti::k, ti::K));
322 :
323 : ::tenex::evaluate<ti::k, ti::j, ti::I>(
324 : symmetrized_d_field_b,
325 : 0.5 * (d_field_b(ti::k, ti::j, ti::I) + d_field_b(ti::j, ti::k, ti::I)));
326 :
327 : ::tenex::evaluate<ti::k>(contracted_symmetrized_d_field_b,
328 : (*symmetrized_d_field_b)(ti::k, ti::i, ti::I));
329 :
330 : ::tenex::evaluate<ti::k>(
331 : field_d_up_times_a_tilde,
332 : (*field_d_up)(ti::k, ti::I, ti::J) * a_tilde(ti::i, ti::j));
333 :
334 : ::tenex::evaluate<ti::i, ti::j>(
335 : conformal_metric_times_field_b,
336 : conformal_spatial_metric(ti::k, ti::i) * field_b(ti::j, ti::K));
337 :
338 : ::tenex::evaluate<ti::i, ti::j>(
339 : conformal_metric_times_trace_a_tilde,
340 : conformal_spatial_metric(ti::i, ti::j) * (*trace_a_tilde)());
341 :
342 : ::tenex::evaluate<ti::k>(inv_conformal_metric_times_d_a_tilde,
343 : (*inv_conformal_spatial_metric)(ti::I, ti::J) *
344 : d_a_tilde(ti::k, ti::i, ti::j));
345 :
346 : ::tenex::evaluate<ti::i, ti::j>(
347 : a_tilde_times_field_b, a_tilde(ti::k, ti::i) * field_b(ti::j, ti::K));
348 :
349 : ::tenex::evaluate<ti::i, ti::j>(
350 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
351 : a_tilde(ti::i, ti::j) -
352 : one_third * (*conformal_metric_times_trace_a_tilde)(ti::i, ti::j));
353 :
354 : ::tenex::evaluate(k_minus_2_theta_c,
355 : trace_extrinsic_curvature() - 2.0 * c * theta());
356 :
357 : ::tenex::evaluate(k_minus_k0_minus_2_theta_c, (*k_minus_2_theta_c)() - k_0());
358 :
359 : ::tenex::evaluate<ti::i, ti::j>(lapse_times_a_tilde,
360 : (lapse)() * a_tilde(ti::i, ti::j));
361 :
362 : tenex::evaluate<ti::k, ti::i, ti::j>(
363 : lapse_times_d_a_tilde, (lapse)() * d_a_tilde(ti::k, ti::i, ti::j));
364 :
365 : ::tenex::evaluate<ti::k>(lapse_times_field_a, (lapse)() * field_a(ti::k));
366 :
367 : ::tenex::evaluate<ti::i, ti::j>(
368 : lapse_times_conformal_spatial_metric,
369 : (lapse)() * conformal_spatial_metric(ti::i, ti::j));
370 :
371 : ::tenex::evaluate(
372 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
373 : (lapse)() * (*ricci_scalar_plus_divergence_z4_constraint)());
374 :
375 : ::tenex::evaluate<ti::I>(shift_times_deriv_gamma_hat,
376 : shift(ti::K) * d_gamma_hat(ti::k, ti::I));
377 :
378 : ::tenex::evaluate<ti::i, ti::j>(
379 : inv_tau_times_conformal_metric,
380 : one_over_relaxation_time * conformal_spatial_metric(ti::i, ti::j));
381 :
382 : // time derivative computation: eq 4(a) - 4(i)
383 : // The way we evaluate the following may seem weird compared
384 : // to 4(a) - 4(i). This is because we try to reuse the first-order
385 : // Ccz4 codes as much as possible. Hence, the following is written
386 : // based on 12a-12i, with s=1, and red terms canceled.
387 :
388 : // eq 12a : time derivative of the conformal spatial metric
389 : // this reduces to eq 4a for our choice of \tau^{-1}=0.
390 : ::tenex::evaluate<ti::i, ti::j>(
391 : dt_conformal_spatial_metric,
392 : 2.0 * shift(ti::K) * field_d(ti::k, ti::i, ti::j) +
393 : (*conformal_metric_times_field_b)(ti::i, ti::j) +
394 : (*conformal_metric_times_field_b)(ti::j, ti::i) -
395 : 2.0 * one_third * conformal_spatial_metric(ti::i, ti::j) *
396 : (*contracted_field_b)() -
397 : 2.0 * (lapse)() *
398 : (*a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde)(
399 : ti::i, ti::j) -
400 : (*inv_tau_times_conformal_metric)(ti::i, ti::j) *
401 : ((*det_conformal_spatial_metric)() - 1.0));
402 :
403 : if (evolve_lapse_and_shift) {
404 : // time derivative of the lapse
405 : ::tenex::evaluate<>(dt_lapse, (shift(ti::K) * field_a(ti::k) -
406 : (*lapse_times_slicing_condition)() *
407 : (*k_minus_k0_minus_2_theta_c)()) *
408 : lapse());
409 : // time derivative of the shift
410 : ::tenex::evaluate<ti::I>(dt_shift, f * b(ti::I));
411 : if (shifting_shift) {
412 : ::tenex::update<ti::I>(
413 : dt_shift, (*dt_shift)(ti::I) + shift(ti::K) * field_b(ti::k, ti::I));
414 : }
415 : } else {
416 : *dt_lapse = make_with_value<Scalar<DataVector>>(lapse, 0.0);
417 : *dt_shift = make_with_value<tnsr::I<DataVector, Dim>>(shift, 0.0);
418 : }
419 :
420 : // time derivative of the the conformal factor
421 : ::tenex::evaluate(dt_conformal_factor,
422 : (shift(ti::K) * field_p(ti::k) +
423 : one_third * ((lapse)() * trace_extrinsic_curvature() -
424 : (*contracted_field_b)())) *
425 : conformal_factor());
426 :
427 : // time derivative of the trace-free part of the extrinsic curvature
428 : ::tenex::evaluate<ti::i, ti::j>(
429 : dt_a_tilde,
430 : shift(ti::K) * d_a_tilde(ti::k, ti::i, ti::j) +
431 : (*conformal_factor_squared)() *
432 : ((lapse)() * ((*spatial_ricci_tensor)(ti::i, ti::j) +
433 : (*grad_spatial_z4_constraint)(ti::i, ti::j) +
434 : (*grad_spatial_z4_constraint)(ti::j, ti::i)) -
435 : (*grad_grad_lapse)(ti::i, ti::j)) -
436 : one_third * conformal_spatial_metric(ti::i, ti::j) *
437 : ((*lapse_times_ricci_scalar_plus_divergence_z4_constraint)() -
438 : (*divergence_lapse)()) +
439 : (*a_tilde_times_field_b)(ti::i, ti::j) +
440 : (*a_tilde_times_field_b)(ti::j, ti::i) -
441 : 2.0 * one_third * a_tilde(ti::i, ti::j) * (*contracted_field_b)() +
442 : (*lapse_times_a_tilde)(ti::i, ti::j) * (*k_minus_2_theta_c)() -
443 : 2.0 * (*lapse_times_a_tilde)(ti::i, ti::l) *
444 : (*inv_conformal_spatial_metric)(ti::L, ti::M) *
445 : a_tilde(ti::m, ti::j) -
446 : (*inv_tau_times_conformal_metric)(ti::i, ti::j) * (*trace_a_tilde)());
447 :
448 : // time derivative of the trace of the extrinsic curvature
449 : ::tenex::evaluate(
450 : dt_trace_extrinsic_curvature,
451 : shift(ti::K) * d_trace_extrinsic_curvature(ti::k) -
452 : (*divergence_lapse)() +
453 : (*lapse_times_ricci_scalar_plus_divergence_z4_constraint)() +
454 : (lapse)() * (trace_extrinsic_curvature() * (*k_minus_2_theta_c)() -
455 : 3.0 * kappa_1 * (1.0 + kappa_2) * theta()));
456 :
457 : // time derivative of the projection of the Z4 four-vector along
458 : // the normal direction
459 : ::tenex::evaluate(
460 : dt_theta,
461 : shift(ti::K) * d_theta(ti::k) +
462 : (lapse)() *
463 : (0.5 * square(cleaning_speed) *
464 : ((*ricci_scalar_plus_divergence_z4_constraint)() +
465 : 2.0 * one_third * square(trace_extrinsic_curvature()) -
466 : a_tilde(ti::i, ti::j) * (*inv_a_tilde)(ti::I, ti::J)) -
467 : c * theta() * trace_extrinsic_curvature() -
468 : (*upper_spatial_z4_constraint)(ti::I)*field_a(ti::i) -
469 : kappa_1 * (2.0 + kappa_2) * theta()));
470 :
471 : // time derivative \hat{\Gamma}^i
472 : // first, compute terms without s
473 : ::tenex::evaluate<ti::I>(
474 : dt_gamma_hat,
475 : // terms without lapse nor s
476 : (*shift_times_deriv_gamma_hat)(ti::I) +
477 : 2.0 * one_third *
478 : (*contracted_conformal_christoffel_second_kind)(ti::I) *
479 : (*contracted_field_b)() -
480 : (*contracted_conformal_christoffel_second_kind)(ti::K)*field_b(
481 : ti::k, ti::I) +
482 : 2.0 * kappa_3 * (*spatial_z4_constraint)(ti::j) *
483 : (2.0 * one_third * (*inv_conformal_spatial_metric)(ti::I, ti::J) *
484 : (*contracted_field_b)() -
485 : (*inv_conformal_spatial_metric)(ti::J, ti::K) *
486 : field_b(ti::k, ti::I)) +
487 : // terms with lapse but not s
488 : 2.0 * (lapse)() *
489 : (-2.0 * one_third *
490 : (*inv_conformal_spatial_metric)(ti::I, ti::J) *
491 : d_trace_extrinsic_curvature(ti::j) +
492 : (*inv_conformal_spatial_metric)(ti::K, ti::I) * d_theta(ti::k) +
493 : (*conformal_christoffel_second_kind)(ti::I, ti::j, ti::k) *
494 : (*inv_a_tilde)(ti::J, ti::K) -
495 : 3.0 * (*inv_a_tilde)(ti::I, ti::J) * field_p(ti::j) -
496 : (*inv_conformal_spatial_metric)(ti::K, ti::I) *
497 : (theta() * field_a(ti::k) +
498 : 2.0 * one_third * trace_extrinsic_curvature() *
499 : (*spatial_z4_constraint)(ti::k)) -
500 : (*inv_a_tilde)(ti::I, ti::J) * field_a(ti::j) -
501 : kappa_1 * (*inv_conformal_spatial_metric)(ti::I, ti::J) *
502 : (*spatial_z4_constraint)(ti::j)));
503 : // We add the following since s=1 (Gamma-driver gauge) is assumed in SoCcz4
504 : ::tenex::update<ti::I>(dt_gamma_hat,
505 : (*dt_gamma_hat)(ti::I) +
506 : // red terms should cancel
507 : // terms with s but not lapse
508 : (*inv_conformal_spatial_metric)(ti::K, ti::L) *
509 : (*symmetrized_d_field_b)(ti::k, ti::l, ti::I) +
510 : one_third *
511 : (*inv_conformal_spatial_metric)(ti::I, ti::K) *
512 : (*contracted_symmetrized_d_field_b)(ti::k));
513 :
514 : // time derivative b^i
515 : if (evolve_lapse_and_shift) {
516 : ::tenex::evaluate<ti::I>(dt_b, (*dt_gamma_hat)(ti::I)-eta() * b(ti::I));
517 : if (shifting_shift) {
518 : ::tenex::update<ti::I>(
519 : dt_b, (*dt_b)(ti::I) + shift(ti::K) * (d_b(ti::k, ti::I) -
520 : d_gamma_hat(ti::k, ti::I)));
521 : }
522 : } else {
523 : (*dt_b).get(0) = 0.0;
524 : (*dt_b).get(1) = 0.0;
525 : (*dt_b).get(2) = 0.0;
526 : }
527 :
528 : // Note that, we do not need to evolve the auxiliary variables in SO-CCZ4.
529 : }
530 : } // namespace detail
531 :
532 : /*!
533 : * \brief Compute the RHS of the second-order CCZ4 formulation of Einstein's
534 : * equations \cite Dumbser2017okk with finite differencing. Also
535 : * update the time derivative at outermost interior points in a sphere domain
536 : * if Sommerfeld boundary conditions are applied.
537 : *
538 : * \details The evolution equations are equations 4(a) - 4(i) of
539 : * \cite Dumbser2017okk Equations 13 - 27 define identities
540 : * used in the evolution equations.
541 : *
542 : * \note Different from the first-order system, the lapse
543 : * and the conformal factor are evolved instead of their natural logs.
544 : * The four auxiliary varialbels \f$A_i\f$, \f$B_k{}^{i}\f$,
545 : * \f$D_{kij}\f$, and \f$P_i\f$ are NOT evolved.
546 : */
547 1 : struct SoTimeDerivative {
548 : template <typename DbTagsList>
549 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
550 : // compute the first spatial derivatives of the evolved variables
551 : using evolved_vars_tag = typename System::variables_tag;
552 : using gradients_tags = typename System::gradients_tags;
553 :
554 : // only 4-th order accurate second derivatives have been implemented
555 : // To keep the same order of accuracy we use the same order for
556 : // both first and second derivatives
557 : constexpr size_t fd_order = 4;
558 : const auto& evolved_vars = db::get<evolved_vars_tag>(*box);
559 : const Mesh<Dim>& subcell_mesh =
560 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
561 : const size_t num_pts = subcell_mesh.number_of_grid_points();
562 : using deriv_var_tag = db::wrap_tags_in<::Tags::deriv, gradients_tags,
563 : tmpl::size_t<Dim>, Frame::Inertial>;
564 : Variables<deriv_var_tag> cell_centered_Ccz4_derivs{num_pts};
565 : const auto& cell_centered_logical_to_inertial_inv_jacobian =
566 : db::get<evolution::dg::subcell::fd::Tags::
567 : InverseJacobianLogicalToInertial<Dim>>(*box);
568 : const auto& cell_centered_logical_to_inertial_inv_hessian = db::get<
569 : evolution::dg::subcell::fd::Tags::InverseHessianLogicalToInertial<Dim>>(
570 : *box);
571 :
572 : constexpr bool subcell_enabled_at_external_boundary =
573 : std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
574 : *box))>::SubcellOptions::subcell_enabled_at_external_boundary;
575 :
576 : const Element<3>& element = db::get<domain::Tags::Element<3>>(*box);
577 :
578 : const Ccz4::fd::Reconstructor& recons =
579 : db::get<Ccz4::fd::Tags::Reconstructor>(*box);
580 :
581 : // If the element has external boundaries and subcell is enabled for
582 : // boundary elements, compute FD ghost data with a given boundary condition.
583 : if constexpr (subcell_enabled_at_external_boundary) {
584 : if (not element.external_boundaries().empty()) {
585 : fd::BoundaryConditionGhostData::apply(box, element, recons);
586 : }
587 : }
588 :
589 : const bool evolve_lapse_and_shift =
590 : get<::Ccz4::fd::Tags::EvolveLapseAndShift>(*box);
591 :
592 : ::Ccz4::fd::spacetime_derivatives(
593 : make_not_null(&cell_centered_Ccz4_derivs), evolved_vars,
594 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>>(
595 : *box),
596 : fd_order, subcell_mesh, cell_centered_logical_to_inertial_inv_jacobian);
597 :
598 : // calculate the four auxiliary fields in eq. 6
599 : // auxiliary variables NOT evolved in SO-CCZ4
600 : const auto& d_lapse =
601 : get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<Dim>,
602 : Frame::Inertial>>(cell_centered_Ccz4_derivs);
603 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(evolved_vars);
604 : const auto field_a = ::tenex::evaluate<ti::i>(d_lapse(ti::i) / lapse());
605 :
606 : const auto& field_b =
607 : get<::Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
608 : Frame::Inertial>>(cell_centered_Ccz4_derivs);
609 :
610 : const auto& d_spatial_conformal_metric =
611 : get<::Tags::deriv<::Ccz4::Tags::ConformalMetric<DataVector, Dim>,
612 : tmpl::size_t<Dim>, Frame::Inertial>>(
613 : cell_centered_Ccz4_derivs);
614 : tnsr::ijj<DataVector, Dim> field_d;
615 : ::tenex::evaluate<ti::i, ti::j, ti::k>(
616 : make_not_null(&field_d),
617 : 0.5 * d_spatial_conformal_metric(ti::i, ti::j, ti::k));
618 :
619 : const auto& conformal_factor =
620 : get<::Ccz4::Tags::ConformalFactor<DataVector>>(evolved_vars);
621 : const auto& d_conformal_factor =
622 : get<::Tags::deriv<::Ccz4::Tags::ConformalFactor<DataVector>,
623 : tmpl::size_t<Dim>, Frame::Inertial>>(
624 : cell_centered_Ccz4_derivs);
625 : const auto field_p = ::tenex::evaluate<ti::i>(d_conformal_factor(ti::i) /
626 : conformal_factor());
627 :
628 : // compute second derivatives of the evolved variables
629 : Variables<db::wrap_tags_in<::Tags::second_deriv, gradients_tags,
630 : tmpl::size_t<Dim>, Frame::Inertial>>
631 : cell_centered_Ccz4_second_derivs{num_pts};
632 :
633 : Ccz4::fd::second_spacetime_derivatives(
634 : make_not_null(&cell_centered_Ccz4_second_derivs), evolved_vars,
635 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>>(
636 : *box),
637 : fd_order, subcell_mesh, cell_centered_logical_to_inertial_inv_jacobian,
638 : cell_centered_logical_to_inertial_inv_hessian);
639 :
640 : // compute spatial derivative of the four auxiliary fields
641 : const auto& d_d_lapse =
642 : get<::Tags::second_deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<Dim>,
643 : Frame::Inertial>>(
644 : cell_centered_Ccz4_second_derivs);
645 : tnsr::ii<DataVector, Dim> d_field_a;
646 : ::tenex::evaluate<ti::i, ti::j>(
647 : make_not_null(&d_field_a),
648 : (d_d_lapse(ti::i, ti::j) - d_lapse(ti::i) * d_lapse(ti::j) / lapse()) /
649 : lapse());
650 :
651 : const auto& d_field_b =
652 : get<::Tags::second_deriv<gr::Tags::Shift<DataVector, Dim>,
653 : tmpl::size_t<Dim>, Frame::Inertial>>(
654 : cell_centered_Ccz4_second_derivs);
655 :
656 : const auto& d_d_conformal_metric =
657 : get<::Tags::second_deriv<::Ccz4::Tags::ConformalMetric<DataVector, Dim>,
658 : tmpl::size_t<Dim>, Frame::Inertial>>(
659 : cell_centered_Ccz4_second_derivs);
660 :
661 : tnsr::iijj<DataVector, Dim> d_field_d;
662 : ::tenex::evaluate<ti::i, ti::j, ti::k, ti::l>(
663 : make_not_null(&d_field_d),
664 : 0.5 * d_d_conformal_metric(ti::i, ti::j, ti::k, ti::l));
665 :
666 : const auto& d_d_conformal_factor =
667 : get<::Tags::second_deriv<::Ccz4::Tags::ConformalFactor<DataVector>,
668 : tmpl::size_t<Dim>, Frame::Inertial>>(
669 : cell_centered_Ccz4_second_derivs);
670 :
671 : tnsr::ii<DataVector, Dim> d_field_p;
672 : ::tenex::evaluate<ti::i, ti::j>(
673 : make_not_null(&d_field_p),
674 : (d_d_conformal_factor(ti::i, ti::j) - d_conformal_factor(ti::i) *
675 : d_conformal_factor(ti::j) /
676 : conformal_factor()) /
677 : conformal_factor());
678 :
679 : // intialize containers to be supplied in the SO-CCZ4 TimeDerivative.cpp
680 : // apply() function quantities we need for computing eq 4, 13 - 27
681 : using TempVars = Variables<tmpl::list<
682 : ::Ccz4::Tags::ConformalFactorSquared<DataVector>,
683 : ::Ccz4::Tags::DetConformalSpatialMetric<DataVector>,
684 : ::Ccz4::Tags::InverseConformalMetric<DataVector, Dim>,
685 : gr::Tags::InverseSpatialMetric<DataVector, Dim>,
686 : ::Ccz4::Tags::InvATilde<DataVector, Dim>,
687 : ::Ccz4::Tags::ATildeTimesFieldB<DataVector, Dim>,
688 : ::Ccz4::Tags::ATildeMinusOneThirdConformalMetricTimesTraceATilde<
689 : DataVector, Dim>,
690 : ::Ccz4::Tags::ContractedFieldB<DataVector>,
691 : ::Ccz4::Tags::SymmetrizedDerivFieldB<DataVector, Dim>,
692 : ::Ccz4::Tags::ContractedSymmetrizedDerivFieldB<DataVector, Dim>,
693 : ::Ccz4::Tags::FieldDUpTimesATilde<DataVector, Dim>,
694 : ::Ccz4::Tags::ContractedFieldDUp<DataVector, Dim>,
695 : ::Ccz4::Tags::HalfConformalFactorSquared<DataVector>,
696 : ::Ccz4::Tags::ConformalMetricTimesFieldB<DataVector, Dim>,
697 : ::Ccz4::Tags::ConformalMetricTimesTraceATilde<DataVector, Dim>,
698 : ::Ccz4::Tags::InverseConformalMetricTimesDerivATilde<DataVector, Dim>,
699 : ::Ccz4::Tags::GammaHatMinusContractedConformalChristoffel<DataVector,
700 : Dim>,
701 : ::Ccz4::Tags::DerivGammaHatMinusContractedConformalChristoffel<
702 : DataVector, Dim>,
703 : ::Ccz4::Tags::ContractedChristoffelSecondKind<DataVector, Dim>,
704 : ::Ccz4::Tags::ContractedDerivConformalChristoffelDifference<DataVector,
705 : Dim>,
706 : ::Ccz4::Tags::KMinus2ThetaC<DataVector>,
707 : ::Ccz4::Tags::KMinusK0Minus2ThetaC<DataVector>,
708 : ::Ccz4::Tags::LapseTimesATilde<DataVector, Dim>,
709 : ::Ccz4::Tags::LapseTimesDerivATilde<DataVector, Dim>,
710 : ::Ccz4::Tags::LapseTimesFieldA<DataVector, Dim>,
711 : ::Ccz4::Tags::LapseTimesConformalMetric<DataVector, Dim>,
712 : ::Ccz4::Tags::LapseTimesSlicingCondition<DataVector>,
713 : ::Ccz4::Tags::LapseTimesRicciScalarPlus2DivergenceZ4Constraint<
714 : DataVector>,
715 : ::Ccz4::Tags::ShiftTimesDerivGammaHat<DataVector, Dim>,
716 : ::Ccz4::Tags::InverseTauTimesConformalMetric<DataVector, Dim>,
717 : ::Ccz4::Tags::TraceATilde<DataVector>,
718 : ::Ccz4::Tags::FieldDUp<DataVector, Dim>,
719 : ::Ccz4::Tags::ConformalChristoffelSecondKind<DataVector, Dim>,
720 : ::Ccz4::Tags::DerivConformalChristoffelSecondKind<DataVector, Dim>,
721 : ::Ccz4::Tags::ChristoffelSecondKind<DataVector, Dim>,
722 : ::Ccz4::Tags::SpatialRicciTensor<DataVector, Dim>,
723 : ::Ccz4::Tags::GradGradLapse<DataVector, Dim>,
724 : ::Ccz4::Tags::DivergenceLapse<DataVector>,
725 : ::Ccz4::Tags::ContractedConformalChristoffelSecondKind<DataVector, Dim>,
726 : ::Ccz4::Tags::DerivContractedConformalChristoffelSecondKind<DataVector,
727 : Dim>,
728 : ::Ccz4::Tags::SpatialZ4Constraint<DataVector, Dim>,
729 : ::Ccz4::Tags::GradSpatialZ4Constraint<DataVector, Dim>,
730 : ::Ccz4::Tags::RicciScalarPlusDivergenceZ4Constraint<DataVector>>>;
731 :
732 : TempVars temp_vars(num_pts);
733 : tnsr::I<DataVector, Dim> upper_spatial_z4_constraint(num_pts);
734 :
735 : // free params
736 : const double c = 1.0; // c = 1.0 in SO-CCZ4
737 : const double cleaning_speed = 1.0; // e in the paper; e = 1.0 for SO-CCZ4
738 : const Scalar<DataVector>& eta = get<::Ccz4::Tags::Eta<DataVector>>(*box);
739 : const double f = Ccz4::fd::System::f;
740 : const Scalar<DataVector>& k_0 = get<::Ccz4::Tags::K0<DataVector>>(*box);
741 : const double kappa_1 = get<::Ccz4::Tags::Kappa1>(*box);
742 : const double kappa_2 = get<::Ccz4::Tags::Kappa2>(*box);
743 : const double kappa_3 = get<::Ccz4::Tags::Kappa3>(*box);
744 : const double one_over_relaxation_time = 0.0; // \tau^{-1} = 0 in SO-CCZ4
745 : const bool shifting_shift = Ccz4::fd::System::shifting_shift;
746 :
747 : // we assume the databox already has tags corresponding to dt of the evolved
748 : // variables
749 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, evolved_vars_tag>;
750 :
751 : db::mutate<dt_variables_tag,
752 : ::Ccz4::Tags::SpatialZ4ConstraintUp<DataVector, Dim>>(
753 : [&](const auto dt_vars_ptr,
754 : const auto upper_spatial_z4_constraint_ptr) {
755 : // resize here
756 : dt_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
757 : auto& [conformal_factor_squared, det_conformal_spatial_metric,
758 : inv_conformal_spatial_metric, inv_spatial_metric, inv_a_tilde,
759 : a_tilde_times_field_b,
760 : a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde,
761 : contracted_field_b, symmetrized_d_field_b,
762 : contracted_symmetrized_d_field_b, field_d_up_times_a_tilde,
763 : contracted_field_d_up, half_conformal_factor_squared,
764 : conformal_metric_times_field_b,
765 : conformal_metric_times_trace_a_tilde,
766 : inv_conformal_metric_times_d_a_tilde,
767 : gamma_hat_minus_contracted_conformal_christoffel,
768 : d_gamma_hat_minus_contracted_conformal_christoffel,
769 : contracted_christoffel_second_kind,
770 : contracted_d_conformal_christoffel_difference,
771 : k_minus_2_theta_c, k_minus_k0_minus_2_theta_c,
772 : lapse_times_a_tilde, lapse_times_d_a_tilde,
773 : lapse_times_field_a, lapse_times_conformal_spatial_metric,
774 : lapse_times_slicing_condition,
775 : lapse_times_ricci_scalar_plus_divergence_z4_constraint,
776 : shift_times_deriv_gamma_hat, inv_tau_times_conformal_metric,
777 : trace_a_tilde, field_d_up, conformal_christoffel_second_kind,
778 : d_conformal_christoffel_second_kind, christoffel_second_kind,
779 : spatial_ricci_tensor, grad_grad_lapse, divergence_lapse,
780 : contracted_conformal_christoffel_second_kind,
781 : d_contracted_conformal_christoffel_second_kind,
782 : spatial_z4_constraint, grad_spatial_z4_constraint,
783 : ricci_scalar_plus_divergence_z4_constraint] = temp_vars;
784 : detail::apply(
785 : // LHS time derivatives of evolved variables: eq 4a - 4i
786 : make_not_null(
787 : &get<::Tags::dt<
788 : ::Ccz4::Tags::ConformalMetric<DataVector, Dim>>>(
789 : *dt_vars_ptr)), // eq 4a
790 : make_not_null(&get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(
791 : *dt_vars_ptr)), // eq 4g
792 : make_not_null(&get<::Tags::dt<gr::Tags::Shift<DataVector, Dim>>>(
793 : *dt_vars_ptr)), // eq 4h
794 : make_not_null(
795 : &get<::Tags::dt<::Ccz4::Tags::ConformalFactor<DataVector>>>(
796 : *dt_vars_ptr)), // eq 4c
797 : make_not_null(
798 : &get<::Tags::dt<::Ccz4::Tags::ATilde<DataVector, Dim>>>(
799 : *dt_vars_ptr)), // eq 4b
800 : make_not_null(&get<::Tags::dt<
801 : gr::Tags::TraceExtrinsicCurvature<DataVector>>>(
802 : *dt_vars_ptr)), // eq 4d
803 : make_not_null(&get<::Tags::dt<::Ccz4::Tags::Theta<DataVector>>>(
804 : *dt_vars_ptr)), // eq 4e
805 : make_not_null(
806 : &get<::Tags::dt<::Ccz4::Tags::GammaHat<DataVector, Dim>>>(
807 : *dt_vars_ptr)), // eq 4f
808 : make_not_null(
809 : &get<::Tags::dt<
810 : ::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>>>(
811 : *dt_vars_ptr)), // eq 4i
812 :
813 : // quantities we need for computing eq 4, 13 - 27
814 : make_not_null(&conformal_factor_squared),
815 : make_not_null(&det_conformal_spatial_metric),
816 : make_not_null(&inv_conformal_spatial_metric),
817 : make_not_null(&inv_spatial_metric), make_not_null(&inv_a_tilde),
818 : // temporary expressions
819 : make_not_null(&a_tilde_times_field_b),
820 : make_not_null(
821 : &a_tilde_minus_one_third_conformal_metric_times_trace_a_tilde),
822 : make_not_null(&contracted_field_b),
823 : make_not_null(&symmetrized_d_field_b),
824 : make_not_null(&contracted_symmetrized_d_field_b),
825 : make_not_null(&field_d_up_times_a_tilde),
826 : make_not_null(&contracted_field_d_up), // temp for eq 18 -20
827 : make_not_null(&half_conformal_factor_squared), // temp for eq 25
828 : make_not_null(&conformal_metric_times_field_b),
829 : make_not_null(&conformal_metric_times_trace_a_tilde),
830 : make_not_null(&inv_conformal_metric_times_d_a_tilde),
831 : make_not_null(&gamma_hat_minus_contracted_conformal_christoffel),
832 : make_not_null(
833 : &d_gamma_hat_minus_contracted_conformal_christoffel),
834 : make_not_null(
835 : &contracted_christoffel_second_kind), // temp for eq 18 -20
836 : make_not_null(
837 : &contracted_d_conformal_christoffel_difference), // temp for
838 : // eq 18 -20
839 : make_not_null(&k_minus_2_theta_c),
840 : make_not_null(&k_minus_k0_minus_2_theta_c),
841 : make_not_null(&lapse_times_a_tilde),
842 : make_not_null(&lapse_times_d_a_tilde),
843 : make_not_null(&lapse_times_field_a),
844 : make_not_null(&lapse_times_conformal_spatial_metric),
845 : make_not_null(&lapse_times_slicing_condition),
846 : make_not_null(
847 : &lapse_times_ricci_scalar_plus_divergence_z4_constraint),
848 : make_not_null(&shift_times_deriv_gamma_hat),
849 : make_not_null(&inv_tau_times_conformal_metric),
850 : // expressions and identities needed for evolution equations: eq
851 : // 13
852 : // - 27
853 : make_not_null(&trace_a_tilde), // eq 13
854 : make_not_null(&field_d_up), // eq 14
855 : make_not_null(&conformal_christoffel_second_kind), // eq 15
856 : make_not_null(&d_conformal_christoffel_second_kind), // eq 16
857 : make_not_null(&christoffel_second_kind), // eq 17
858 : make_not_null(&spatial_ricci_tensor), // eq 18 - 20
859 : make_not_null(&grad_grad_lapse), // eq 21
860 : make_not_null(&divergence_lapse), // eq 22
861 : make_not_null(
862 : &contracted_conformal_christoffel_second_kind), // eq 23
863 : make_not_null(
864 : &d_contracted_conformal_christoffel_second_kind), // eq 24
865 : make_not_null(&spatial_z4_constraint), // eq 25
866 : make_not_null(&upper_spatial_z4_constraint), // eq 25
867 : make_not_null(&grad_spatial_z4_constraint), // eq 26
868 : make_not_null(
869 : &ricci_scalar_plus_divergence_z4_constraint), // eq 27
870 : // fixed params
871 : c, cleaning_speed, // e in the paper
872 : one_over_relaxation_time, // \tau^{-1}
873 : // free params
874 : eta, f, kappa_1, kappa_2, kappa_3, k_0,
875 : // evolved variables
876 : get<::Ccz4::Tags::ConformalMetric<DataVector, Dim>>(evolved_vars),
877 : get<gr::Tags::Lapse<DataVector>>(evolved_vars),
878 : get<gr::Tags::Shift<DataVector, Dim>>(evolved_vars),
879 : get<::Ccz4::Tags::ConformalFactor<DataVector>>(evolved_vars),
880 : get<::Ccz4::Tags::ATilde<DataVector, Dim>>(evolved_vars),
881 : get<gr::Tags::TraceExtrinsicCurvature<DataVector>>(evolved_vars),
882 : get<::Ccz4::Tags::Theta<DataVector>>(evolved_vars),
883 : get<::Ccz4::Tags::GammaHat<DataVector, Dim>>(evolved_vars),
884 : get<::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>>(evolved_vars),
885 :
886 : field_a, // auxiliary variables NOT evolved in SO-CCZ4
887 : field_b, field_d, field_p,
888 : d_field_a, // spatial derivative of auxiliary variables
889 : d_field_b, d_field_d, d_field_p,
890 :
891 : // spatial derivatives of other evolved variables
892 : get<::Tags::deriv<::Ccz4::Tags::ATilde<DataVector, Dim>,
893 : tmpl::size_t<Dim>, Frame::Inertial>>(
894 : cell_centered_Ccz4_derivs),
895 : get<::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataVector>,
896 : tmpl::size_t<Dim>, Frame::Inertial>>(
897 : cell_centered_Ccz4_derivs),
898 : get<::Tags::deriv<::Ccz4::Tags::Theta<DataVector>,
899 : tmpl::size_t<Dim>, Frame::Inertial>>(
900 : cell_centered_Ccz4_derivs),
901 : get<::Tags::deriv<::Ccz4::Tags::GammaHat<DataVector, Dim>,
902 : tmpl::size_t<Dim>, Frame::Inertial>>(
903 : cell_centered_Ccz4_derivs),
904 : get<::Tags::deriv<::Ccz4::Tags::AuxiliaryShiftB<DataVector, Dim>,
905 : tmpl::size_t<Dim>, Frame::Inertial>>(
906 : cell_centered_Ccz4_derivs),
907 : shifting_shift, evolve_lapse_and_shift);
908 :
909 : *upper_spatial_z4_constraint_ptr =
910 : std::move(upper_spatial_z4_constraint);
911 : },
912 : box);
913 :
914 : // handling Sommerfeld BCs
915 : // WARNING: We assume a complete sphere domain (all wedges)
916 : if constexpr (not subcell_enabled_at_external_boundary) {
917 : return;
918 : }
919 : if (element.external_boundaries().empty()) {
920 : return;
921 : }
922 :
923 : const auto& external_bcs =
924 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(*box).at(
925 : element.id().block_id());
926 : for (const auto& direction : element.external_boundaries()) {
927 : ASSERT(direction == Direction<Dim>::upper_zeta(),
928 : "external bc direction is not upper_zeta but " << direction);
929 : // this is potentially expensive; we should add a boolean flag for
930 : // time-dependent bc. This should be done in a future PR.
931 : if (dynamic_cast<const Ccz4::BoundaryConditions::Sommerfeld*>(
932 : external_bcs.at(direction).get()) == nullptr) {
933 : return; // no need to check more, as we do not allow mixed BCs
934 : }
935 : }
936 :
937 : db::mutate<dt_variables_tag>(
938 : [&](const auto dt_vars_ptr, const auto& inertial_coords,
939 : const auto& evolved_var) {
940 : const size_t num_face_pts =
941 : subcell_mesh.extents(0) * subcell_mesh.extents(1);
942 :
943 : ASSERT(inertial_coords.get(0).size() ==
944 : subcell_mesh.number_of_grid_points() and
945 : inertial_coords.get(1).size() ==
946 : subcell_mesh.number_of_grid_points() and
947 : inertial_coords.get(2).size() ==
948 : subcell_mesh.number_of_grid_points(),
949 : "Inertial coords size ["
950 : << inertial_coords.get(0).size()
951 : << ", " << inertial_coords.get(1).size()
952 : << ", " << inertial_coords.get(2).size()
953 : << "] does not match number of subcell points "
954 : << subcell_mesh.number_of_grid_points());
955 :
956 : std::array<DataVector, Dim> outermost_inertial_coords{};
957 : for (size_t i = 0; i < Dim; ++i) {
958 : make_const_view<DataVector>(
959 : make_not_null(&outermost_inertial_coords.at(i)),
960 : inertial_coords.get(i),
961 : (subcell_mesh.extents(2) - 1) * num_face_pts, num_face_pts);
962 : }
963 : const DataVector outermost_radial_coords =
964 : sqrt(square(outermost_inertial_coords[0]) +
965 : square(outermost_inertial_coords[1]) +
966 : square(outermost_inertial_coords[2]));
967 :
968 : // modify the time derivatives at the outermost interior points
969 : // per Sommerfeld BCs
970 : tmpl::for_each<Ccz4::fd::System::variables_tag_list>(
971 : [&]<typename Tag>(tmpl::type_<Tag> /*meta*/) {
972 : const auto& var = get<Tag>(evolved_var);
973 : const auto& d_var =
974 : get<::Tags::deriv<Tag, tmpl::size_t<Dim>, Frame::Inertial>>(
975 : cell_centered_Ccz4_derivs);
976 : auto& dt_var = get<::Tags::dt<Tag>>(*dt_vars_ptr);
977 :
978 : for (size_t tensor_index = 0; tensor_index < Tag::type::size();
979 : ++tensor_index) {
980 : DataVector outermost_var{};
981 : make_const_view<DataVector>(
982 : make_not_null(&outermost_var), var[tensor_index],
983 : (subcell_mesh.extents(2) - 1) * num_face_pts,
984 : num_face_pts);
985 :
986 : std::array<DataVector, Dim> outermost_d_var{};
987 : for (size_t i = 0; i < Dim; ++i) {
988 : make_const_view<DataVector>(
989 : make_not_null(&outermost_d_var.at(i)),
990 : d_var[Dim * tensor_index + i],
991 : (subcell_mesh.extents(2) - 1) * num_face_pts,
992 : num_face_pts);
993 : }
994 :
995 : const DataVector outermost_dr_var =
996 : (outermost_inertial_coords[0] * outermost_d_var[0] +
997 : outermost_inertial_coords[1] * outermost_d_var[1] +
998 : outermost_inertial_coords[2] * outermost_d_var[2]) /
999 : outermost_radial_coords;
1000 : DataVector outermost_dt_var{};
1001 : outermost_dt_var.set_data_ref(
1002 : dt_var[tensor_index].data() +
1003 : (subcell_mesh.extents(2) - 1) * num_face_pts,
1004 : num_face_pts);
1005 : outermost_dt_var =
1006 : outermost_var * (-1.0 / outermost_radial_coords) -
1007 : outermost_dr_var;
1008 : }
1009 : });
1010 : },
1011 : box,
1012 : db::get<evolution::dg::subcell::Tags::Coordinates<
1013 : Dim, Frame::Inertial>>(*box),
1014 : db::get<evolved_vars_tag>(*box));
1015 : }
1016 : };
1017 : } // namespace Ccz4::fd
|
