Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <type_traits>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/TaggedContainers.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "DataStructures/VectorImpl.hpp"
19 : #include "Domain/Structure/Element.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "Domain/TagsTimeDependent.hpp"
22 : #include "Evolution/BoundaryCorrectionTags.hpp"
23 : #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
24 : #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
25 : #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
26 : #include "Evolution/DgSubcell/Projection.hpp"
27 : #include "Evolution/DgSubcell/ReconstructionOrder.hpp"
28 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
29 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
30 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
31 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
32 : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
34 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
36 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp"
37 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/BoundaryConditionGhostData.hpp"
38 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Derivatives.hpp"
39 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/FilterOptions.hpp"
40 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Filters.hpp"
41 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Reconstructor.hpp"
42 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/FiniteDifference/Tag.hpp"
43 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/StressEnergy.hpp"
44 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/System.hpp"
45 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/Tags.hpp"
46 : #include "Evolution/Systems/GrMhd/GhValenciaDivClean/TimeDerivativeTerms.hpp"
47 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp"
48 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp"
49 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp"
50 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp"
51 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TimeDerivativeTerms.hpp"
52 : #include "NumericalAlgorithms/FiniteDifference/PartialDerivatives.hpp"
53 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
54 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/DerivSpatialMetric.hpp"
55 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/ExtrinsicCurvature.hpp"
56 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/SpatialDerivOfLapse.hpp"
57 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/SpatialDerivOfShift.hpp"
58 : #include "PointwiseFunctions/Hydro/Tags.hpp"
59 : #include "Utilities/CallWithDynamicType.hpp"
60 : #include "Utilities/ErrorHandling/Assert.hpp"
61 : #include "Utilities/Gsl.hpp"
62 : #include "Utilities/TMPL.hpp"
63 :
64 : /// \cond
65 : namespace Tags {
66 : struct Time;
67 : } // namespace Tags
68 : /// \endcond
69 :
70 : namespace grmhd::GhValenciaDivClean::subcell {
71 : namespace detail {
72 : template <class GhDtTagsList, class GhTemporariesList, class GhGradientTagsList,
73 : class GhExtraTagsList, class GrmhdDtTagsList,
74 : class GrmhdSourceTagsList, class GrmhdArgumentSourceTagsList>
75 : struct ComputeTimeDerivImpl;
76 :
77 : template <class... GhDtTags, class... GhTemporaries, class... GhGradientTags,
78 : class... GhExtraTags, class... GrmhdDtTags, class... GrmhdSourceTags,
79 : class... GrmhdArgumentSourceTags>
80 : struct ComputeTimeDerivImpl<
81 : tmpl::list<GhDtTags...>, tmpl::list<GhTemporaries...>,
82 : tmpl::list<GhGradientTags...>, tmpl::list<GhExtraTags...>,
83 : tmpl::list<GrmhdDtTags...>, tmpl::list<GrmhdSourceTags...>,
84 : tmpl::list<GrmhdArgumentSourceTags...>> {
85 : template <class DbTagsList>
86 : static void apply(
87 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
88 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
89 : const Scalar<DataVector>& cell_centered_det_inv_jacobian,
90 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
91 : Frame::Inertial>&
92 : cell_centered_logical_to_inertial_inv_jacobian,
93 : const std::array<double, 3>& one_over_delta_xi,
94 : const std::array<Variables<tmpl::list<GrmhdDtTags...>>, 3>&
95 : boundary_corrections,
96 : const Variables<
97 : db::wrap_tags_in<::Tags::deriv, typename System::gradients_tags,
98 : tmpl::size_t<3>, Frame::Inertial>>& gh_derivs) {
99 : const Mesh<3>& subcell_mesh =
100 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
101 : const size_t number_of_points = subcell_mesh.number_of_grid_points();
102 : // Note: GH+GRMHD tags are always GH,GRMHD
103 : using deriv_lapse = ::Tags::deriv<gr::Tags::Lapse<DataVector>,
104 : tmpl::size_t<3>, Frame::Inertial>;
105 : using deriv_shift = ::Tags::deriv<gr::Tags::Shift<DataVector, 3>,
106 : tmpl::size_t<3>, Frame::Inertial>;
107 : using deriv_spatial_metric =
108 : ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
109 : Frame::Inertial>;
110 : using extra_tags_for_grmhd =
111 : tmpl::list<deriv_lapse, deriv_shift, deriv_spatial_metric,
112 : gr::Tags::ExtrinsicCurvature<DataVector, 3>>;
113 : using temporary_tags = tmpl::remove_duplicates<tmpl::append<
114 : typename gh::TimeDerivative<3_st>::temporary_tags,
115 : tmpl::push_front<typename grmhd::ValenciaDivClean::TimeDerivativeTerms::
116 : temporary_tags,
117 : ::gh::ConstraintDamping::Tags::ConstraintGamma0>,
118 : extra_tags_for_grmhd,
119 : tmpl::list<Tags::TraceReversedStressEnergy, Tags::FourVelocityOneForm,
120 : Tags::ComovingMagneticFieldOneForm>>>;
121 : Variables<temporary_tags> temp_tags{subcell_mesh.number_of_grid_points()};
122 : const auto temp_tags_ptr = make_not_null(&temp_tags);
123 :
124 : // Compute constraint damping terms.
125 : const double time = db::get<::Tags::Time>(*box);
126 : const auto& functions_of_time =
127 : db::get<::domain::Tags::FunctionsOfTime>(*box);
128 : const auto& grid_coords =
129 : db::get<evolution::dg::subcell::Tags::Coordinates<3, Frame::Grid>>(
130 : *box);
131 : db::get<
132 : gh::ConstraintDamping::Tags::DampingFunctionGamma0<3, Frame::Grid>> (
133 : *box)(get<gh::ConstraintDamping::Tags::ConstraintGamma0>(temp_tags_ptr),
134 : grid_coords, time, functions_of_time);
135 : db::get<
136 : gh::ConstraintDamping::Tags::DampingFunctionGamma1<3, Frame::Grid>> (
137 : *box)(get<gh::ConstraintDamping::Tags::ConstraintGamma1>(temp_tags_ptr),
138 : grid_coords, time, functions_of_time);
139 : db::get<
140 : gh::ConstraintDamping::Tags::DampingFunctionGamma2<3, Frame::Grid>> (
141 : *box)(get<gh::ConstraintDamping::Tags::ConstraintGamma2>(temp_tags_ptr),
142 : grid_coords, time, functions_of_time);
143 :
144 : using variables_tag = typename System::variables_tag;
145 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
146 : const gsl::not_null<typename dt_variables_tag::type*> dt_vars_ptr =
147 : db::mutate<dt_variables_tag>(
148 : [](const auto local_dt_vars_ptr) { return local_dt_vars_ptr; },
149 : box);
150 : dt_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
151 :
152 : using primitives_tag = typename System::primitive_variables_tag;
153 : using evolved_vars_tag =
154 : typename grmhd::GhValenciaDivClean::System::variables_tag;
155 :
156 : const auto& primitive_vars = db::get<primitives_tag>(*box);
157 : const auto& evolved_vars = db::get<evolved_vars_tag>(*box);
158 :
159 : // Velocity of the moving mesh, if applicable. We project the value
160 : // stored on the DG grid onto the subcell grid.
161 : const Mesh<3>& dg_mesh = db::get<domain::Tags::Mesh<3>>(*box);
162 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
163 : mesh_velocity_dg = db::get<domain::Tags::MeshVelocity<3>>(*box);
164 : const std::optional<Scalar<DataVector>>& div_mesh_velocity_dg =
165 : db::get<domain::Tags::DivMeshVelocity>(*box);
166 : std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>
167 : mesh_velocity_subcell = {};
168 : if (mesh_velocity_dg.has_value()) {
169 : mesh_velocity_subcell = tnsr::I<DataVector, 3, Frame::Inertial>{
170 : subcell_mesh.number_of_grid_points()};
171 : for (size_t i = 0; i < 3; i++) {
172 : mesh_velocity_subcell.value().get(i) =
173 : evolution::dg::subcell::fd::project(mesh_velocity_dg.value().get(i),
174 : dg_mesh,
175 : subcell_mesh.extents());
176 : }
177 : }
178 :
179 : gh::TimeDerivative<3_st>::apply(
180 : get<::Tags::dt<GhDtTags>>(dt_vars_ptr)...,
181 : get<GhTemporaries>(temp_tags_ptr)...,
182 : get<::Tags::deriv<GhGradientTags, tmpl::size_t<3>, Frame::Inertial>>(
183 : gh_derivs)...,
184 : get<GhExtraTags>(evolved_vars, temp_tags)...,
185 :
186 : db::get<::gh::gauges::Tags::GaugeCondition>(*box),
187 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box), time,
188 : inertial_coords, cell_centered_logical_to_inertial_inv_jacobian,
189 : mesh_velocity_subcell);
190 : if (get<gh::gauges::Tags::GaugeCondition>(*box).is_harmonic()) {
191 : get(get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(*temp_tags_ptr)) =
192 : sqrt(
193 : get(get<gr::Tags::DetSpatialMetric<DataVector>>(*temp_tags_ptr)));
194 : }
195 :
196 : // Add source terms from moving mesh
197 : if (mesh_velocity_dg.has_value()) {
198 : tmpl::for_each<tmpl::list<GhDtTags...>>([&dt_vars_ptr,
199 : &mesh_velocity_subcell,
200 : &gh_derivs](
201 : auto evolved_var_tag_v) {
202 : using evolved_var_tag = tmpl::type_from<decltype(evolved_var_tag_v)>;
203 : using dt_tag = ::Tags::dt<evolved_var_tag>;
204 : using grad_tag =
205 : ::Tags::deriv<evolved_var_tag, tmpl::size_t<3>, Frame::Inertial>;
206 : // Flux and gradients use the same indexing conventions,
207 : // replacing the direction of the face with the direction
208 : // of the derivative.
209 : using FluxTensor = typename grad_tag::type;
210 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
211 : const auto& grad_var = get<grad_tag>(gh_derivs);
212 : for (size_t i = 0; i < dt_var.size(); ++i) {
213 : const auto tensor_index = dt_var.get_tensor_index(i);
214 : for (size_t j = 0; j < 3; j++) {
215 : const auto grad_index =
216 : FluxTensor::get_storage_index(prepend(tensor_index, j));
217 : // Add (mesh_velocity)^j grad_j (var[i])
218 : dt_var[i] +=
219 : mesh_velocity_subcell.value().get(j) * grad_var[grad_index];
220 : }
221 : }
222 : });
223 : }
224 :
225 : {
226 : // Set extra tags needed for GRMHD source terms. We compute these from
227 : // quantities already computed inside the GH RHS computation to minimize
228 : // FLOPs.
229 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(temp_tags);
230 : const auto& half_phi_two_normals =
231 : get<gh::Tags::HalfPhiTwoNormals<3>>(temp_tags);
232 : const auto& phi = get<gh::Tags::Phi<DataVector, 3>>(evolved_vars);
233 : const auto& phi_one_normal = get<gh::Tags::PhiOneNormal<3>>(temp_tags);
234 : const auto& spacetime_normal_vector =
235 : get<gr::Tags::SpacetimeNormalVector<DataVector, 3>>(temp_tags);
236 : const auto& inverse_spacetime_metric =
237 : get<gr::Tags::InverseSpacetimeMetric<DataVector, 3>>(temp_tags);
238 :
239 : auto& spatial_deriv_lapse = get<deriv_lapse>(temp_tags);
240 : auto& spatial_deriv_shift = get<deriv_shift>(temp_tags);
241 : // Compute d_i beta^i
242 : for (size_t i = 0; i < 3; ++i) {
243 : // Use spatial_deriv_lapse as temp buffer to reduce number of 2*
244 : // operations.
245 : const auto& phi_two_normals_i = spatial_deriv_lapse.get(i) =
246 : 2.0 * half_phi_two_normals.get(i);
247 : for (size_t j = 0; j < 3; ++j) {
248 : spatial_deriv_shift.get(i, j) =
249 : spacetime_normal_vector.get(j + 1) * phi_two_normals_i;
250 : for (size_t a = 0; a < 4; ++a) {
251 : spatial_deriv_shift.get(i, j) +=
252 : inverse_spacetime_metric.get(j + 1, a) *
253 : phi_one_normal.get(i, a);
254 : }
255 : spatial_deriv_shift.get(i, j) *= get(lapse);
256 : }
257 : }
258 :
259 : // Compute d_i lapse
260 : for (size_t i = 0; i < 3; ++i) {
261 : spatial_deriv_lapse.get(i) = -get(lapse) * half_phi_two_normals.get(i);
262 : }
263 : // Extract d_i \gamma_{ij}
264 : for (size_t k = 0; k < 3; ++k) {
265 : for (size_t i = 0; i < 3; ++i) {
266 : for (size_t j = i; j < 3; ++j) {
267 : get<deriv_spatial_metric>(temp_tags).get(k, i, j) =
268 : phi.get(k, i + 1, j + 1);
269 : }
270 : }
271 : }
272 :
273 : // Compute extrinsic curvature
274 : const auto& pi = get<gh::Tags::Pi<DataVector, 3>>(evolved_vars);
275 : for (size_t i = 0; i < 3; ++i) {
276 : for (size_t j = i; j < 3; ++j) {
277 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(temp_tags).get(i,
278 : j) =
279 : 0.5 * (pi.get(i + 1, j + 1) + phi_one_normal.get(i, j + 1) +
280 : phi_one_normal.get(j, i + 1));
281 : }
282 : }
283 : } // End scope for computing metric terms in GRMHD source terms.
284 :
285 : grmhd::ValenciaDivClean::ComputeSources::apply(
286 : get<::Tags::dt<GrmhdSourceTags>>(dt_vars_ptr)...,
287 : get<GrmhdArgumentSourceTags>(temp_tags, primitive_vars, evolved_vars,
288 : *box)...);
289 :
290 : // Zero GRMHD tags that don't have sources.
291 : tmpl::for_each<tmpl::list<GrmhdDtTags...>>([&dt_vars_ptr](
292 : auto evolved_var_tag_v) {
293 : using evolved_var_tag = tmpl::type_from<decltype(evolved_var_tag_v)>;
294 : using dt_tag = ::Tags::dt<evolved_var_tag>;
295 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
296 : for (size_t i = 0; i < dt_var.size(); ++i) {
297 : if constexpr (not tmpl::list_contains_v<tmpl::list<GrmhdSourceTags...>,
298 : evolved_var_tag>) {
299 : // Zero the GRMHD dt(u) for variables that do not have a source term .
300 : // This is necessary to avoid `+=` to a `NaN` (debug mode) or random
301 : // garbage (release mode) when adding to dt_var below.
302 : dt_var[i] = 0.0;
303 : }
304 : }
305 : });
306 : // Correction to source terms due to moving mesh
307 : if (div_mesh_velocity_dg.has_value()) {
308 : const DataVector div_mesh_velocity_subcell =
309 : evolution::dg::subcell::fd::project(
310 : div_mesh_velocity_dg.value().get(), dg_mesh,
311 : subcell_mesh.extents());
312 : tmpl::for_each<tmpl::list<GrmhdDtTags...>>(
313 : [&dt_vars_ptr, &div_mesh_velocity_subcell,
314 : &evolved_vars](auto evolved_var_tag_v) {
315 : using evolved_var_tag =
316 : tmpl::type_from<decltype(evolved_var_tag_v)>;
317 : using dt_tag = ::Tags::dt<evolved_var_tag>;
318 : auto& dt_var = get<dt_tag>(*dt_vars_ptr);
319 : const auto& evolved_var = get<evolved_var_tag>(evolved_vars);
320 : for (size_t i = 0; i < dt_var.size(); ++i) {
321 : dt_var[i] -= div_mesh_velocity_subcell * evolved_var[i];
322 : }
323 : });
324 : }
325 :
326 : const tnsr::ii<DataVector, 3> spatial_metric{};
327 : for (size_t i = 0; i < 3; ++i) {
328 : for (size_t j = i; j < 3; ++j) {
329 : make_const_view(
330 : make_not_null(&spatial_metric.get(i, j)),
331 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(evolved_vars)
332 : .get(i + 1, j + 1),
333 : 0, number_of_points);
334 : }
335 : }
336 :
337 : tenex::evaluate<ti::i>(get<hydro::Tags::SpatialVelocityOneForm<
338 : DataVector, 3, Frame::Inertial>>(temp_tags_ptr),
339 : get<hydro::Tags::SpatialVelocity<DataVector, 3>>(
340 : primitive_vars)(ti::J) *
341 : spatial_metric(ti::i, ti::j));
342 :
343 : tenex::evaluate<ti::i>(
344 : get<hydro::Tags::MagneticFieldOneForm<DataVector, 3, Frame::Inertial>>(
345 : temp_tags_ptr),
346 : get<hydro::Tags::MagneticField<DataVector, 3>>(primitive_vars)(ti::J) *
347 : spatial_metric(ti::i, ti::j));
348 :
349 : tenex::evaluate(
350 : get<hydro::Tags::MagneticFieldSquared<DataVector>>(temp_tags_ptr),
351 : get<hydro::Tags::MagneticField<DataVector, 3>>(primitive_vars)(ti::J) *
352 : get<hydro::Tags::MagneticFieldOneForm<DataVector, 3>>(temp_tags)(
353 : ti::j));
354 :
355 : tenex::evaluate(
356 : get<hydro::Tags::MagneticFieldDotSpatialVelocity<DataVector>>(
357 : temp_tags_ptr),
358 : get<hydro::Tags::SpatialVelocity<DataVector, 3>>(primitive_vars)(
359 : ti::J) *
360 : get<hydro::Tags::MagneticFieldOneForm<DataVector, 3>>(temp_tags)(
361 : ti::j));
362 :
363 : tenex::evaluate(get<typename ValenciaDivClean::TimeDerivativeTerms::
364 : OneOverLorentzFactorSquared>(temp_tags_ptr),
365 : 1.0 / (square(get<hydro::Tags::LorentzFactor<DataVector>>(
366 : primitive_vars)())));
367 :
368 : trace_reversed_stress_energy(
369 : get<Tags::TraceReversedStressEnergy>(temp_tags_ptr),
370 : get<Tags::FourVelocityOneForm>(temp_tags_ptr),
371 : get<Tags::ComovingMagneticFieldOneForm>(temp_tags_ptr),
372 :
373 : get<hydro::Tags::RestMassDensity<DataVector>>(evolved_vars, temp_tags,
374 : primitive_vars),
375 : get<hydro::Tags::SpatialVelocityOneForm<DataVector, 3,
376 : Frame::Inertial>>(
377 : evolved_vars, temp_tags, primitive_vars),
378 :
379 : get<hydro::Tags::MagneticFieldOneForm<DataVector, 3, Frame::Inertial>>(
380 : evolved_vars, temp_tags, primitive_vars),
381 :
382 : get<hydro::Tags::MagneticFieldSquared<DataVector>>(
383 : evolved_vars, temp_tags, primitive_vars),
384 :
385 : get<hydro::Tags::MagneticFieldDotSpatialVelocity<DataVector>>(
386 : evolved_vars, temp_tags, primitive_vars),
387 : get<hydro::Tags::LorentzFactor<DataVector>>(evolved_vars, temp_tags,
388 : primitive_vars),
389 : get<typename ValenciaDivClean::TimeDerivativeTerms::
390 : OneOverLorentzFactorSquared>(evolved_vars, temp_tags,
391 : primitive_vars),
392 : get<hydro::Tags::Pressure<DataVector>>(evolved_vars, temp_tags,
393 : primitive_vars),
394 : get<hydro::Tags::SpecificInternalEnergy<DataVector>>(
395 : evolved_vars, temp_tags, primitive_vars),
396 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(evolved_vars, temp_tags,
397 : primitive_vars),
398 : get<gr::Tags::Shift<DataVector, 3>>(evolved_vars, temp_tags,
399 : primitive_vars),
400 : get<gr::Tags::Lapse<DataVector>>(evolved_vars, temp_tags,
401 : primitive_vars));
402 :
403 : add_stress_energy_term_to_dt_pi(
404 : get<::Tags::dt<gh::Tags::Pi<DataVector, 3>>>(dt_vars_ptr),
405 : get<Tags::TraceReversedStressEnergy>(temp_tags),
406 : get<gr::Tags::Lapse<DataVector>>(temp_tags));
407 :
408 : for (size_t dim = 0; dim < 3; ++dim) {
409 : const auto& boundary_correction_in_axis =
410 : gsl::at(boundary_corrections, dim);
411 : const double inverse_delta = gsl::at(one_over_delta_xi, dim);
412 : EXPAND_PACK_LEFT_TO_RIGHT([&dt_vars_ptr, &boundary_correction_in_axis,
413 : &cell_centered_det_inv_jacobian, dim,
414 : inverse_delta, &subcell_mesh]() {
415 : auto& dt_var = *get<::Tags::dt<GrmhdDtTags>>(dt_vars_ptr);
416 : const auto& var_correction =
417 : get<GrmhdDtTags>(boundary_correction_in_axis);
418 : for (size_t i = 0; i < dt_var.size(); ++i) {
419 : evolution::dg::subcell::add_cartesian_flux_divergence(
420 : make_not_null(&dt_var[i]), inverse_delta,
421 : get(cell_centered_det_inv_jacobian), var_correction[i],
422 : subcell_mesh.extents(), dim);
423 : }
424 : }());
425 : }
426 : }
427 : };
428 : } // namespace detail
429 :
430 : /*!
431 : * \brief Compute the time derivative on the subcell grid using FD
432 : * reconstruction.
433 : *
434 : * The code makes the following unchecked assumptions:
435 : * - Assumes Cartesian coordinates with a diagonal Jacobian matrix
436 : * from the logical to the inertial frame
437 : */
438 1 : struct TimeDerivative {
439 : template <typename DbTagsList>
440 0 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
441 : using evolved_vars_tag =
442 : typename grmhd::GhValenciaDivClean::System::variables_tag;
443 : using evolved_vars_tags = typename evolved_vars_tag::tags_list;
444 : using grmhd_evolved_vars_tag =
445 : typename grmhd::ValenciaDivClean::System::variables_tag;
446 : using grmhd_evolved_vars_tags = typename grmhd_evolved_vars_tag::tags_list;
447 : using fluxes_tags =
448 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
449 : tmpl::size_t<3>, Frame::Inertial>;
450 : using prim_tag = typename System::primitive_variables_tag;
451 : using prim_tags = typename prim_tag::tags_list;
452 : using recons_prim_tags = tmpl::push_front<tmpl::push_back<
453 : prim_tags,
454 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>>>;
455 : using gradients_tags = typename System::gradients_tags;
456 :
457 : const Mesh<3>& dg_mesh = db::get<domain::Tags::Mesh<3>>(*box);
458 : const Mesh<3>& subcell_mesh =
459 : db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
460 : ASSERT(
461 : subcell_mesh == Mesh<3>(subcell_mesh.extents(0), subcell_mesh.basis(0),
462 : subcell_mesh.quadrature(0)),
463 : "The subcell/FD mesh must be isotropic for the FD time derivative but "
464 : "got "
465 : << subcell_mesh);
466 : const size_t num_pts = subcell_mesh.number_of_grid_points();
467 : const size_t reconstructed_num_pts =
468 : (subcell_mesh.extents(0) + 1) *
469 : subcell_mesh.extents().slice_away(0).product();
470 :
471 : const tnsr::I<DataVector, 3, Frame::ElementLogical>&
472 : cell_centered_logical_coords =
473 : db::get<evolution::dg::subcell::Tags::Coordinates<
474 : 3, Frame::ElementLogical>>(*box);
475 : std::array<double, 3> one_over_delta_xi{};
476 : for (size_t i = 0; i < 3; ++i) {
477 : // Note: assumes isotropic extents
478 : gsl::at(one_over_delta_xi, i) =
479 : 1.0 / (get<0>(cell_centered_logical_coords)[1] -
480 : get<0>(cell_centered_logical_coords)[0]);
481 : }
482 : const auto& cell_centered_logical_to_inertial_inv_jacobian = db::get<
483 : evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToInertial<3>>(
484 : *box);
485 : const auto& inertial_coords =
486 : db::get<evolution::dg::subcell::Tags::Coordinates<3, Frame::Inertial>>(
487 : *box);
488 :
489 : const Element<3>& element = db::get<domain::Tags::Element<3>>(*box);
490 : const bool element_is_interior = element.external_boundaries().empty();
491 : constexpr bool subcell_enabled_at_external_boundary =
492 : std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
493 : *box))>::SubcellOptions::subcell_enabled_at_external_boundary;
494 :
495 : ASSERT(element_is_interior or subcell_enabled_at_external_boundary,
496 : "Subcell time derivative is called at a boundary element while "
497 : "using subcell is disabled at external boundaries."
498 : "ElementID "
499 : << element.id());
500 :
501 : const fd::Reconstructor& recons = db::get<fd::Tags::Reconstructor>(*box);
502 : // If the element has external boundaries and subcell is enabled for
503 : // boundary elements, compute FD ghost data with a given boundary condition.
504 : if constexpr (subcell_enabled_at_external_boundary) {
505 : if (not element_is_interior) {
506 : fd::BoundaryConditionGhostData::apply(box, element, recons);
507 : }
508 : }
509 : std::optional<std::array<gsl::span<std::uint8_t>, 3>>
510 : reconstruction_order{};
511 :
512 : if (const auto& filter_options =
513 : db::get<grmhd::GhValenciaDivClean::fd::Tags::FilterOptions>(*box);
514 : filter_options.spacetime_dissipation.has_value()) {
515 : db::mutate<evolved_vars_tag>(
516 : [&filter_options, &recons, &subcell_mesh](const auto evolved_vars_ptr,
517 : const auto& ghost_data) {
518 : typename evolved_vars_tag::type filtered_vars = *evolved_vars_ptr;
519 : // $(recons.ghost_zone_size() - 1) * 2 + 1$ => always use highest
520 : // order dissipation filter possible.
521 : grmhd::GhValenciaDivClean::fd::spacetime_kreiss_oliger_filter(
522 : make_not_null(&filtered_vars), *evolved_vars_ptr, ghost_data,
523 : subcell_mesh, 2 * recons.ghost_zone_size(),
524 : filter_options.spacetime_dissipation.value());
525 : *evolved_vars_ptr = filtered_vars;
526 : },
527 : box,
528 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(
529 : *box));
530 : }
531 :
532 : // Velocity of the moving mesh on the dg grid, if applicable.
533 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
534 : mesh_velocity_dg = db::get<domain::Tags::MeshVelocity<3>>(*box);
535 : // Inverse jacobian, to be projected on faces
536 : const auto& inv_jacobian_dg =
537 : db::get<domain::Tags::InverseJacobian<3, Frame::ElementLogical,
538 : Frame::Inertial>>(*box);
539 : const auto& det_inv_jacobian_dg = db::get<
540 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
541 : *box);
542 :
543 : // GH+GRMHD is a bit different.
544 : // 1. Compute GH time derivative, since this will also give us lapse, shift,
545 : // etc. that we need to reconstruct.
546 : // 2. Compute d_t Pi_{ab} source terms from MHD (or do we wait until post
547 : // MHD source terms?)
548 : // 3. Reconstruct MHD+spacetime vars to interfaces
549 : // 4. Compute MHD time derivatives.
550 : //
551 : // Compute FD GH derivatives with neighbor data
552 : // Use highest possible FD order for number of GZ, 2 * (ghost_zone_size)
553 : const auto& evolved_vars = db::get<evolved_vars_tag>(*box);
554 : Variables<db::wrap_tags_in<::Tags::deriv, gradients_tags, tmpl::size_t<3>,
555 : Frame::Inertial>>
556 : cell_centered_gh_derivs{num_pts};
557 : grmhd::GhValenciaDivClean::fd::spacetime_derivatives(
558 : make_not_null(&cell_centered_gh_derivs), evolved_vars,
559 : db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<3>>(
560 : *box),
561 : recons.ghost_zone_size()*2, subcell_mesh,
562 : cell_centered_logical_to_inertial_inv_jacobian);
563 :
564 : // Now package the data and compute the correction
565 : //
566 : // Note: Assumes a the GH and GRMHD corrections can be invoked separately.
567 : // This is reasonable since the systems are a tensor product system.
568 : const auto& base_boundary_correction =
569 : db::get<evolution::Tags::BoundaryCorrection<System>>(*box);
570 : using derived_boundary_corrections = typename std::decay_t<
571 : decltype(base_boundary_correction)>::creatable_classes;
572 : std::array<Variables<grmhd_evolved_vars_tags>, 3> boundary_corrections{};
573 : call_with_dynamic_type<void, derived_boundary_corrections>(
574 : &base_boundary_correction, [&](const auto* gh_grmhd_correction) {
575 : // Need the GH packaged tags to avoid projecting them.
576 : using gh_dg_package_field_tags = typename std::decay_t<
577 : decltype(gh_grmhd_correction
578 : ->gh_correction())>::dg_package_field_tags;
579 : // Only apply correction to GRMHD variables.
580 : const auto& boundary_correction =
581 : gh_grmhd_correction->valencia_correction();
582 : using DerivedCorrection = std::decay_t<decltype(boundary_correction)>;
583 : using dg_package_data_temporary_tags =
584 : typename DerivedCorrection::dg_package_data_temporary_tags;
585 :
586 : using dg_package_data_argument_tags = tmpl::append<
587 : evolved_vars_tags, recons_prim_tags, fluxes_tags,
588 : tmpl::remove_duplicates<tmpl::push_back<
589 : dg_package_data_temporary_tags,
590 : gr::Tags::SpatialMetric<DataVector, 3>,
591 : gr::Tags::SqrtDetSpatialMetric<DataVector>,
592 : gr::Tags::InverseSpatialMetric<DataVector, 3>,
593 : evolution::dg::Actions::detail::NormalVector<3>>>>;
594 :
595 : // Computed prims and cons on face via reconstruction
596 : auto package_data_argvars_lower_face = make_array<3>(
597 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
598 : auto package_data_argvars_upper_face = make_array<3>(
599 : Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
600 :
601 : // Reconstruct data to the face
602 : call_with_dynamic_type<void, typename grmhd::GhValenciaDivClean::fd::
603 : Reconstructor::creatable_classes>(
604 : &recons, [&box, &package_data_argvars_lower_face,
605 : &package_data_argvars_upper_face,
606 : &reconstruction_order](const auto& reconstructor) {
607 : using ReconstructorType =
608 : std::decay_t<decltype(*reconstructor)>;
609 : db::apply<
610 : typename ReconstructorType::reconstruction_argument_tags>(
611 : [&package_data_argvars_lower_face,
612 : &package_data_argvars_upper_face, &reconstructor,
613 : &reconstruction_order](const auto&... args) {
614 : if constexpr (ReconstructorType::use_adaptive_order) {
615 : reconstructor->reconstruct(
616 : make_not_null(&package_data_argvars_lower_face),
617 : make_not_null(&package_data_argvars_upper_face),
618 : make_not_null(&reconstruction_order), args...);
619 : } else {
620 : (void)reconstruction_order;
621 : reconstructor->reconstruct(
622 : make_not_null(&package_data_argvars_lower_face),
623 : make_not_null(&package_data_argvars_upper_face),
624 : args...);
625 : }
626 : },
627 : *box);
628 : });
629 :
630 : using dg_package_field_tags =
631 : typename DerivedCorrection::dg_package_field_tags;
632 : // Allocated outside for loop to reduce allocations
633 : Variables<dg_package_field_tags> upper_packaged_data{
634 : reconstructed_num_pts};
635 : Variables<dg_package_field_tags> lower_packaged_data{
636 : reconstructed_num_pts};
637 :
638 : // Compute fluxes on faces
639 : for (size_t i = 0; i < 3; ++i) {
640 : auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, i);
641 : auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, i);
642 : grmhd::ValenciaDivClean::subcell::compute_fluxes(
643 : make_not_null(&vars_upper_face));
644 : grmhd::ValenciaDivClean::subcell::compute_fluxes(
645 : make_not_null(&vars_lower_face));
646 :
647 : // Build extents of mesh shifted by half a grid cell in direction i
648 : const unsigned long& num_subcells_1d = subcell_mesh.extents(0);
649 : Index<3> face_mesh_extents(std::array<size_t, 3>{
650 : num_subcells_1d, num_subcells_1d, num_subcells_1d});
651 : face_mesh_extents[i] = num_subcells_1d + 1;
652 : // Add moving mesh corrections to the fluxes, if needed
653 : std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>
654 : mesh_velocity_on_face = {};
655 : if (mesh_velocity_dg.has_value()) {
656 : // Project mesh velocity on face mesh.
657 : // Can we get away with only doing the normal component? It
658 : // is also used in the packaged data...
659 : mesh_velocity_on_face = tnsr::I<DataVector, 3, Frame::Inertial>{
660 : reconstructed_num_pts};
661 : for (size_t j = 0; j < 3; j++) {
662 : // j^th component of the velocity on the i^th directed face
663 : mesh_velocity_on_face.value().get(j) =
664 : evolution::dg::subcell::fd::project_to_faces(
665 : mesh_velocity_dg.value().get(j), dg_mesh,
666 : face_mesh_extents, i);
667 : }
668 :
669 : tmpl::for_each<grmhd_evolved_vars_tags>(
670 : [&vars_upper_face, &vars_lower_face,
671 : &mesh_velocity_on_face](auto tag_v) {
672 : using tag = tmpl::type_from<decltype(tag_v)>;
673 : using flux_tag =
674 : ::Tags::Flux<tag, tmpl::size_t<3>, Frame::Inertial>;
675 : using FluxTensor = typename flux_tag::type;
676 : const auto& var_upper = get<tag>(vars_upper_face);
677 : const auto& var_lower = get<tag>(vars_lower_face);
678 : auto& flux_upper = get<flux_tag>(vars_upper_face);
679 : auto& flux_lower = get<flux_tag>(vars_lower_face);
680 : for (size_t storage_index = 0;
681 : storage_index < var_upper.size(); ++storage_index) {
682 : const auto tensor_index =
683 : var_upper.get_tensor_index(storage_index);
684 : for (size_t j = 0; j < 3; j++) {
685 : const auto flux_storage_index =
686 : FluxTensor::get_storage_index(
687 : prepend(tensor_index, j));
688 : flux_upper[flux_storage_index] -=
689 : mesh_velocity_on_face.value().get(j) *
690 : var_upper[storage_index];
691 : flux_lower[flux_storage_index] -=
692 : mesh_velocity_on_face.value().get(j) *
693 : var_lower[storage_index];
694 : }
695 : }
696 : });
697 : }
698 :
699 : // Normal vectors in curved spacetime normalized by inverse
700 : // spatial metric. Since we assume a Cartesian grid, this is
701 : // relatively easy. Note that we use the sign convention on
702 : // the normal vectors to be compatible with DG.
703 : //
704 : // Note that these normal vectors are on all faces inside the DG
705 : // element since there are a bunch of subcells. We don't use the
706 : // NormalCovectorAndMagnitude tag in the DataBox right now to avoid
707 : // conflicts with the DG solver. We can explore in the future if
708 : // it's possible to reuse that allocation.
709 : //
710 : // The unnormalized normal vector is
711 : // n_j = d \xi^{\hat i}/dx^j
712 : // with "i" the current face.
713 : tnsr::i<DataVector, 3, Frame::Inertial> lower_outward_conormal{
714 : reconstructed_num_pts, 0.0};
715 : for (size_t j = 0; j < 3; j++) {
716 : lower_outward_conormal.get(j) =
717 : evolution::dg::subcell::fd::project_to_faces(
718 : inv_jacobian_dg.get(i, j), dg_mesh, face_mesh_extents, i);
719 : }
720 : const auto det_inv_jacobian_face =
721 : evolution::dg::subcell::fd::project_to_faces(
722 : get(det_inv_jacobian_dg), dg_mesh, face_mesh_extents, i);
723 :
724 : const Scalar<DataVector> normalization{sqrt(get(
725 : dot_product(lower_outward_conormal, lower_outward_conormal,
726 : get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(
727 : vars_upper_face))))};
728 : for (size_t j = 0; j < 3; j++) {
729 : lower_outward_conormal.get(j) =
730 : lower_outward_conormal.get(j) / get(normalization);
731 : }
732 :
733 : tnsr::i<DataVector, 3, Frame::Inertial> upper_outward_conormal{
734 : reconstructed_num_pts, 0.0};
735 : for (size_t j = 0; j < 3; j++) {
736 : upper_outward_conormal.get(j) = -lower_outward_conormal.get(j);
737 : }
738 : // Note: we probably should compute the normal vector in addition to
739 : // the co-vector. Not a huge issue since we'll get an FPE right now
740 : // if it's used by a Riemann solver.
741 :
742 : // Compute the packaged data
743 : using dg_package_data_projected_tags = tmpl::append<
744 : grmhd_evolved_vars_tags, fluxes_tags,
745 : dg_package_data_temporary_tags,
746 : typename DerivedCorrection::dg_package_data_primitive_tags>;
747 : evolution::dg::Actions::detail::dg_package_data<System>(
748 : make_not_null(&upper_packaged_data),
749 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
750 : vars_upper_face, upper_outward_conormal, mesh_velocity_on_face,
751 : *box, typename DerivedCorrection::dg_package_data_volume_tags{},
752 : dg_package_data_projected_tags{});
753 :
754 : evolution::dg::Actions::detail::dg_package_data<System>(
755 : make_not_null(&lower_packaged_data),
756 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
757 : vars_lower_face, lower_outward_conormal, mesh_velocity_on_face,
758 : *box, typename DerivedCorrection::dg_package_data_volume_tags{},
759 : dg_package_data_projected_tags{});
760 :
761 : // Now need to check if any of our neighbors are doing DG,
762 : // because if so then we need to use whatever boundary data
763 : // they sent instead of what we computed locally.
764 : //
765 : // Note: We could check this beforehand to avoid the extra
766 : // work of reconstruction and flux computations at the
767 : // boundaries.
768 : evolution::dg::subcell::correct_package_data<true>(
769 : make_not_null(&lower_packaged_data),
770 : make_not_null(&upper_packaged_data), i, element, subcell_mesh,
771 : db::get<evolution::dg::Tags::MortarData<3>>(*box),
772 : Variables<gh_dg_package_field_tags>::
773 : number_of_independent_components);
774 :
775 : // Compute the corrections on the faces. We only need to
776 : // compute this once because we can just flip the normal
777 : // vectors then
778 : gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
779 : evolution::dg::subcell::compute_boundary_terms(
780 : make_not_null(&gsl::at(boundary_corrections, i)),
781 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
782 : upper_packaged_data, lower_packaged_data);
783 : // We need to multiply by the normal vector normalization
784 : gsl::at(boundary_corrections, i) *= get(normalization);
785 : // Also multiply by determinant of Jacobian, following Eq.(34)
786 : // of 2109.11645
787 : gsl::at(boundary_corrections, i) *= 1.0 / det_inv_jacobian_face;
788 : }
789 : });
790 :
791 : // Now compute the actual time derivatives.
792 : using gh_variables_tags =
793 : typename System::gh_system::variables_tag::tags_list;
794 : using gh_gradient_tags = typename TimeDerivativeTerms::gh_gradient_tags;
795 : using gh_temporary_tags = typename TimeDerivativeTerms::gh_temp_tags;
796 : using gh_extra_tags =
797 : tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
798 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>,
799 : ::gh::ConstraintDamping::Tags::ConstraintGamma0,
800 : ::gh::ConstraintDamping::Tags::ConstraintGamma1,
801 : ::gh::ConstraintDamping::Tags::ConstraintGamma2>;
802 : using grmhd_source_tags =
803 : tmpl::transform<ValenciaDivClean::ComputeSources::return_tags,
804 : tmpl::bind<db::remove_tag_prefix, tmpl::_1>>;
805 : using grmhd_source_argument_tags =
806 : ValenciaDivClean::ComputeSources::argument_tags;
807 : detail::ComputeTimeDerivImpl<gh_variables_tags, gh_temporary_tags,
808 : gh_gradient_tags, gh_extra_tags,
809 : grmhd_evolved_vars_tags, grmhd_source_tags,
810 : grmhd_source_argument_tags>::
811 : apply(box, inertial_coords,
812 : db::get<evolution::dg::subcell::fd::Tags::
813 : DetInverseJacobianLogicalToInertial>(*box),
814 : cell_centered_logical_to_inertial_inv_jacobian, one_over_delta_xi,
815 : boundary_corrections, cell_centered_gh_derivs);
816 : evolution::dg::subcell::store_reconstruction_order_in_databox(
817 : box, reconstruction_order);
818 : }
819 : };
820 : } // namespace grmhd::GhValenciaDivClean::subcell
|