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