SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/NewtonianEuler/Subcell - TimeDerivative.hpp Hit Total Coverage
Commit: 27a3bc9d135e816b53755f48a5e499ff1f2c3d56 Lines: 1 4 25.0 %
Date: 2024-11-02 01:14:36
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/AsAccess.hpp"
      12             : #include "DataStructures/DataBox/DataBox.hpp"
      13             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      14             : #include "DataStructures/DataBox/Prefixes.hpp"
      15             : #include "DataStructures/DataVector.hpp"
      16             : #include "DataStructures/Tensor/Tensor.hpp"
      17             : #include "DataStructures/Variables.hpp"
      18             : #include "Domain/CoordinateMaps/Tags.hpp"
      19             : #include "Domain/Structure/Element.hpp"
      20             : #include "Domain/Tags.hpp"
      21             : #include "Evolution/BoundaryCorrectionTags.hpp"
      22             : #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
      23             : #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
      24             : #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
      25             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      26             : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
      27             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      28             : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
      29             : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
      30             : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
      31             : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
      32             : #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp"
      33             : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
      34             : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp"
      35             : #include "Evolution/Systems/NewtonianEuler/Fluxes.hpp"
      36             : #include "Evolution/Systems/NewtonianEuler/Subcell/ComputeFluxes.hpp"
      37             : #include "Evolution/Systems/NewtonianEuler/System.hpp"
      38             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      39             : #include "Parallel/Tags/Metavariables.hpp"
      40             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      41             : #include "Utilities/CallWithDynamicType.hpp"
      42             : #include "Utilities/ErrorHandling/Assert.hpp"
      43             : #include "Utilities/Gsl.hpp"
      44             : #include "Utilities/TMPL.hpp"
      45             : 
      46             : namespace NewtonianEuler::subcell {
      47             : /*!
      48             :  * \brief Compute the time derivative on the subcell grid using FD
      49             :  * reconstruction.
      50             :  *
      51             :  * The code makes the following unchecked assumptions:
      52             :  * - Assumes Cartesian coordinates with a diagonal Jacobian matrix
      53             :  * from the logical to the inertial frame
      54             :  * - Assumes the mesh is not moving (grid and inertial frame are the same)
      55             :  */
      56             : template <size_t Dim>
      57           1 : struct TimeDerivative {
      58             :   template <typename DbTagsList>
      59           0 :   static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box) {
      60             :     using metavariables = typename std::decay_t<decltype(
      61             :         db::get<Parallel::Tags::Metavariables>(*box))>;
      62             :     using system = typename metavariables::system;
      63             :     using evolved_vars_tag = typename system::variables_tag;
      64             :     using evolved_vars_tags = typename evolved_vars_tag::tags_list;
      65             :     using prim_tags = typename system::primitive_variables_tag::tags_list;
      66             :     using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
      67             :                                          tmpl::size_t<Dim>, Frame::Inertial>;
      68             : 
      69             :     ASSERT((db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
      70             :                 Dim, Frame::Grid, Frame::Inertial>>(*box))
      71             :                .is_identity(),
      72             :            "Do not yet support moving mesh with DG-subcell.");
      73             : 
      74             :     // The copy of Mesh is intentional to avoid a GCC-7 internal compiler error.
      75             :     const Mesh<Dim> subcell_mesh =
      76             :         db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
      77             :     ASSERT(
      78             :         subcell_mesh == Mesh<Dim>(subcell_mesh.extents(0),
      79             :                                   subcell_mesh.basis(0),
      80             :                                   subcell_mesh.quadrature(0)),
      81             :         "The subcell/FD mesh must be isotropic for the FD time derivative but "
      82             :         "got "
      83             :             << subcell_mesh);
      84             :     const size_t num_pts = subcell_mesh.number_of_grid_points();
      85             :     const size_t reconstructed_num_pts =
      86             :         (subcell_mesh.extents(0) + 1) *
      87             :         subcell_mesh.extents().slice_away(0).product();
      88             : 
      89             :     const tnsr::I<DataVector, Dim, Frame::ElementLogical>&
      90             :         cell_centered_logical_coords =
      91             :             db::get<evolution::dg::subcell::Tags::Coordinates<
      92             :                 Dim, Frame::ElementLogical>>(*box);
      93             :     std::array<double, Dim> one_over_delta_xi{};
      94             :     for (size_t i = 0; i < Dim; ++i) {
      95             :       // Note: assumes isotropic extents
      96             :       gsl::at(one_over_delta_xi, i) =
      97             :           1.0 / (get<0>(cell_centered_logical_coords)[1] -
      98             :                  get<0>(cell_centered_logical_coords)[0]);
      99             :     }
     100             : 
     101             :     const NewtonianEuler::fd::Reconstructor<Dim>& recons =
     102             :         db::get<NewtonianEuler::fd::Tags::Reconstructor<Dim>>(*box);
     103             : 
     104             :     const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
     105             :     ASSERT(element.external_boundaries().size() == 0,
     106             :            "Can't have external boundaries right now with subcell. ElementID "
     107             :                << element.id());
     108             : 
     109             :     // Now package the data and compute the correction
     110             :     const auto& boundary_correction =
     111             :         db::get<evolution::Tags::BoundaryCorrection<system>>(*box);
     112             :     using derived_boundary_corrections =
     113             :         typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
     114             :     std::array<Variables<evolved_vars_tags>, Dim> boundary_corrections{};
     115             :     tmpl::for_each<derived_boundary_corrections>([&boundary_correction,
     116             :                                                   &reconstructed_num_pts,
     117             :                                                   &recons, &box, &element,
     118             :                                                   &subcell_mesh,
     119             :                                                   &boundary_corrections](
     120             :                                                      auto
     121             :                                                          derived_correction_v) {
     122             :       using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
     123             :       if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
     124             :         using dg_package_data_temporary_tags =
     125             :             typename DerivedCorrection::dg_package_data_temporary_tags;
     126             :         using dg_package_data_argument_tags =
     127             :             tmpl::append<evolved_vars_tags, prim_tags, fluxes_tags,
     128             :                          dg_package_data_temporary_tags>;
     129             :         // Computed prims and cons on face via reconstruction
     130             :         auto package_data_argvars_lower_face = make_array<Dim>(
     131             :             Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
     132             :         auto package_data_argvars_upper_face = make_array<Dim>(
     133             :             Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
     134             : 
     135             :         // Reconstruct data to the face
     136             :         call_with_dynamic_type<void, typename NewtonianEuler::fd::Reconstructor<
     137             :                                          Dim>::creatable_classes>(
     138             :             &recons,
     139             :             [&box, &package_data_argvars_lower_face,
     140             :              &package_data_argvars_upper_face](const auto& reconstructor) {
     141             :               db::apply<typename std::decay_t<decltype(
     142             :                   *reconstructor)>::reconstruction_argument_tags>(
     143             :                   [&package_data_argvars_lower_face,
     144             :                    &package_data_argvars_upper_face,
     145             :                    &reconstructor](const auto&... args) {
     146             :                     reconstructor->reconstruct(
     147             :                         make_not_null(&package_data_argvars_lower_face),
     148             :                         make_not_null(&package_data_argvars_upper_face),
     149             :                         args...);
     150             :                   },
     151             :                   *box);
     152             :             });
     153             : 
     154             :         using dg_package_field_tags =
     155             :             typename DerivedCorrection::dg_package_field_tags;
     156             :         // Allocated outside for loop to reduce allocations
     157             :         Variables<dg_package_field_tags> upper_packaged_data{
     158             :             reconstructed_num_pts};
     159             :         Variables<dg_package_field_tags> lower_packaged_data{
     160             :             reconstructed_num_pts};
     161             : 
     162             :         // Compute fluxes on faces
     163             :         for (size_t i = 0; i < Dim; ++i) {
     164             :           auto& vars_upper_face = gsl::at(package_data_argvars_upper_face, i);
     165             :           auto& vars_lower_face = gsl::at(package_data_argvars_lower_face, i);
     166             :           NewtonianEuler::subcell::compute_fluxes<Dim>(
     167             :               make_not_null(&vars_upper_face));
     168             :           NewtonianEuler::subcell::compute_fluxes<Dim>(
     169             :               make_not_null(&vars_lower_face));
     170             : 
     171             :           tnsr::i<DataVector, Dim, Frame::Inertial> lower_outward_conormal{
     172             :               reconstructed_num_pts, 0.0};
     173             :           lower_outward_conormal.get(i) = 1.0;
     174             : 
     175             :           tnsr::i<DataVector, Dim, Frame::Inertial> upper_outward_conormal{
     176             :               reconstructed_num_pts, 0.0};
     177             :           upper_outward_conormal.get(i) = -1.0;
     178             : 
     179             :           // Compute the packaged data
     180             :           using dg_package_data_projected_tags = tmpl::append<
     181             :               evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
     182             :               typename DerivedCorrection::dg_package_data_primitive_tags>;
     183             :           evolution::dg::Actions::detail::dg_package_data<system>(
     184             :               make_not_null(&upper_packaged_data),
     185             :               dynamic_cast<const DerivedCorrection&>(boundary_correction),
     186             :               vars_upper_face, upper_outward_conormal, {std::nullopt}, *box,
     187             :               typename DerivedCorrection::dg_package_data_volume_tags{},
     188             :               dg_package_data_projected_tags{});
     189             : 
     190             :           evolution::dg::Actions::detail::dg_package_data<system>(
     191             :               make_not_null(&lower_packaged_data),
     192             :               dynamic_cast<const DerivedCorrection&>(boundary_correction),
     193             :               vars_lower_face, lower_outward_conormal, {std::nullopt}, *box,
     194             :               typename DerivedCorrection::dg_package_data_volume_tags{},
     195             :               dg_package_data_projected_tags{});
     196             : 
     197             :           // Now need to check if any of our neighbors are doing DG,
     198             :           // because if so then we need to use whatever boundary data
     199             :           // they sent instead of what we computed locally.
     200             :           //
     201             :           // Note: We could check this beforehand to avoid the extra
     202             :           // work of reconstruction and flux computations at the
     203             :           // boundaries.
     204             :           evolution::dg::subcell::correct_package_data<true>(
     205             :               make_not_null(&lower_packaged_data),
     206             :               make_not_null(&upper_packaged_data), i, element, subcell_mesh,
     207             :               db::get<evolution::dg::Tags::MortarData<Dim>>(*box), 0);
     208             : 
     209             :           // Compute the corrections on the faces. We only need to
     210             :           // compute this once because we can just flip the normal
     211             :           // vectors then
     212             :           gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
     213             :           evolution::dg::subcell::compute_boundary_terms(
     214             :               make_not_null(&gsl::at(boundary_corrections, i)),
     215             :               dynamic_cast<const DerivedCorrection&>(boundary_correction),
     216             :               upper_packaged_data, lower_packaged_data, db::as_access(*box),
     217             :               typename DerivedCorrection::dg_boundary_terms_volume_tags{});
     218             :         }
     219             :       }
     220             :     });
     221             : 
     222             :     // Now compute the actual time derivatives.
     223             :     using dt_variables_tag = db::add_tag_prefix<::Tags::dt, evolved_vars_tag>;
     224             :     using source_argument_tags = tmpl::list<
     225             :         Tags::MassDensityCons, Tags::MomentumDensity<Dim>, Tags::EnergyDensity,
     226             :         hydro::Tags::SpatialVelocity<DataVector, Dim>,
     227             :         hydro::Tags::Pressure<DataVector>,
     228             :         hydro::Tags::SpecificInternalEnergy<DataVector>,
     229             :         hydro::Tags::EquationOfState<false, 2>,
     230             :         evolution::dg::subcell::Tags::Coordinates<Dim, Frame::Inertial>,
     231             :         ::Tags::Time, NewtonianEuler::Tags::SourceTerm<Dim>>;
     232             :     db::mutate_apply<tmpl::list<dt_variables_tag>, source_argument_tags>(
     233             :         [&num_pts, &boundary_corrections, &subcell_mesh, &one_over_delta_xi,
     234             :          &cell_centered_logical_to_grid_inv_jacobian =
     235             :              db::get<evolution::dg::subcell::fd::Tags::
     236             :                          InverseJacobianLogicalToGrid<Dim>>(*box)](
     237             :             const auto dt_vars_ptr, const Scalar<DataVector>& mass_density_cons,
     238             :             const tnsr::I<DataVector, Dim>& momentum_density,
     239             :             const Scalar<DataVector>& energy_density,
     240             :             const tnsr::I<DataVector, Dim>& velocity,
     241             :             const Scalar<DataVector>& pressure,
     242             :             const Scalar<DataVector>& specific_internal_energy,
     243             :             const EquationsOfState::EquationOfState<false, 2>& eos,
     244             :             const tnsr::I<DataVector, Dim>& coords, const double time,
     245             :             const Sources::Source<Dim>& source) {
     246             :           dt_vars_ptr->initialize(num_pts, 0.0);
     247             :           using MassDensityCons = NewtonianEuler::Tags::MassDensityCons;
     248             :           using MomentumDensity = NewtonianEuler::Tags::MomentumDensity<Dim>;
     249             :           using EnergyDensity = NewtonianEuler::Tags::EnergyDensity;
     250             : 
     251             :           auto& dt_mass = get<::Tags::dt<MassDensityCons>>(*dt_vars_ptr);
     252             :           auto& dt_momentum = get<::Tags::dt<MomentumDensity>>(*dt_vars_ptr);
     253             :           auto& dt_energy = get<::Tags::dt<EnergyDensity>>(*dt_vars_ptr);
     254             : 
     255             :           const auto eos_2d = eos.promote_to_2d_eos();
     256             :           source(make_not_null(&dt_mass), make_not_null(&dt_momentum),
     257             :                  make_not_null(&dt_energy), mass_density_cons, momentum_density,
     258             :                  energy_density, velocity, pressure, specific_internal_energy,
     259             :                  *eos_2d, coords, time);
     260             : 
     261             :           for (size_t dim = 0; dim < Dim; ++dim) {
     262             :             Scalar<DataVector>& mass_density_correction =
     263             :                 get<MassDensityCons>(gsl::at(boundary_corrections, dim));
     264             :             Scalar<DataVector>& energy_density_correction =
     265             :                 get<EnergyDensity>(gsl::at(boundary_corrections, dim));
     266             :             tnsr::I<DataVector, Dim, Frame::Inertial>&
     267             :                 momentum_density_correction =
     268             :                     get<MomentumDensity>(gsl::at(boundary_corrections, dim));
     269             :             evolution::dg::subcell::add_cartesian_flux_divergence(
     270             :                 make_not_null(&get(dt_mass)), gsl::at(one_over_delta_xi, dim),
     271             :                 cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
     272             :                 get(mass_density_correction), subcell_mesh.extents(), dim);
     273             :             evolution::dg::subcell::add_cartesian_flux_divergence(
     274             :                 make_not_null(&get(dt_energy)), gsl::at(one_over_delta_xi, dim),
     275             :                 cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
     276             :                 get(energy_density_correction), subcell_mesh.extents(), dim);
     277             :             for (size_t d = 0; d < Dim; ++d) {
     278             :               evolution::dg::subcell::add_cartesian_flux_divergence(
     279             :                   make_not_null(&dt_momentum.get(d)),
     280             :                   gsl::at(one_over_delta_xi, dim),
     281             :                   cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
     282             :                   momentum_density_correction.get(d), subcell_mesh.extents(),
     283             :                   dim);
     284             :             }
     285             :           }
     286             :         },
     287             :         box);
     288             :   }
     289             : 
     290             :  private:
     291             :   template <typename SourceTerm, typename... SourceTermArgs,
     292             :             typename... SourcedVars>
     293           0 :   static void sources_impl(std::tuple<gsl::not_null<Scalar<DataVector>*>,
     294             :                                       gsl::not_null<tnsr::I<DataVector, Dim>*>,
     295             :                                       gsl::not_null<Scalar<DataVector>*>>
     296             :                                dt_vars,
     297             :                            tmpl::list<SourcedVars...> /*meta*/,
     298             :                            const SourceTerm& source,
     299             :                            const SourceTermArgs&... source_term_args) {
     300             :     using dt_vars_list =
     301             :         tmpl::list<Tags::MassDensityCons, Tags::MomentumDensity<Dim>,
     302             :                    Tags::EnergyDensity>;
     303             : 
     304             :     source.apply(
     305             :         std::get<tmpl::index_of<dt_vars_list, SourcedVars>::value>(dt_vars)...,
     306             :         source_term_args...);
     307             :   }
     308             : };
     309             : }  // namespace NewtonianEuler::subcell

Generated by: LCOV version 1.14