SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell - TimeDerivative.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 1 3 33.3 %
Date: 2024-03-29 00:33:31
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14