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

Generated by: LCOV version 1.14