SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GrMhd/ValenciaDivClean/Subcell - NeighborPackagedData.hpp Hit Total Coverage
Commit: 3f09028930c0450a2fb61ee918b22882f5d03d2b Lines: 1 3 33.3 %
Date: 2021-10-22 20:52:16
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 <boost/functional/hash.hpp>
       7             : #include <cstddef>
       8             : #include <optional>
       9             : #include <type_traits>
      10             : #include <utility>
      11             : #include <vector>
      12             : 
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      15             : #include "DataStructures/DataBox/Prefixes.hpp"
      16             : #include "DataStructures/DataVector.hpp"
      17             : #include "DataStructures/FixedHashMap.hpp"
      18             : #include "DataStructures/Index.hpp"
      19             : #include "DataStructures/Tensor/Tensor.hpp"
      20             : #include "DataStructures/Variables.hpp"
      21             : #include "DataStructures/VariablesTag.hpp"
      22             : #include "Domain/Structure/Direction.hpp"
      23             : #include "Domain/Structure/Element.hpp"
      24             : #include "Domain/Structure/ElementId.hpp"
      25             : #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
      26             : #include "Domain/Tags.hpp"
      27             : #include "Domain/TagsTimeDependent.hpp"
      28             : #include "Evolution/BoundaryCorrectionTags.hpp"
      29             : #include "Evolution/DgSubcell/NeighborData.hpp"
      30             : #include "Evolution/DgSubcell/Projection.hpp"
      31             : #include "Evolution/DgSubcell/Reconstruction.hpp"
      32             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      33             : #include "Evolution/DgSubcell/Tags/NeighborData.hpp"
      34             : #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
      35             : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
      36             : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
      37             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp"
      38             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/Factory.hpp"
      39             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
      40             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp"
      41             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp"
      42             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
      43             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      44             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      45             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      46             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      47             : #include "Utilities/FakeVirtual.hpp"
      48             : #include "Utilities/Gsl.hpp"
      49             : #include "Utilities/TMPL.hpp"
      50             : 
      51             : namespace grmhd::ValenciaDivClean::subcell {
      52             : /*!
      53             :  * \brief On elements using DG, reconstructs the interface data from a
      54             :  * neighboring element doing subcell.
      55             :  *
      56             :  * The neighbor's packaged data needed by the boundary correction is computed
      57             :  * and returned so that it can be used for solving the Riemann problem on the
      58             :  * interfaces.
      59             :  *
      60             :  * Note that for strict conservation the Riemann solve should be done on the
      61             :  * subcells, with the correction being projected back to the DG interface.
      62             :  * However, in practice such strict conservation doesn't seem to be necessary
      63             :  * and can be explained by that we only need strict conservation at shocks, and
      64             :  * if one element is doing DG, then we aren't at a shock.
      65             :  */
      66           1 : struct NeighborPackagedData {
      67             :   template <typename DbTagsList>
      68             :   static FixedHashMap<
      69             :       maximum_number_of_neighbors(3), std::pair<Direction<3>, ElementId<3>>,
      70             :       std::vector<double>, boost::hash<std::pair<Direction<3>, ElementId<3>>>>
      71           0 :   apply(const db::DataBox<DbTagsList>& box,
      72             :         const std::vector<std::pair<Direction<3>, ElementId<3>>>&
      73             :             mortars_to_reconstruct_to) {
      74             :     using evolved_vars_tag = typename System::variables_tag;
      75             :     using evolved_vars_tags = typename evolved_vars_tag::tags_list;
      76             :     using prim_tags = typename System::primitive_variables_tag::tags_list;
      77             :     using recons_prim_tags = tmpl::push_back<
      78             :         prim_tags,
      79             :         hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>>;
      80             :     using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
      81             :                                          tmpl::size_t<3>, Frame::Inertial>;
      82             : 
      83             :     ASSERT(not db::get<domain::Tags::MeshVelocity<3>>(box).has_value(),
      84             :            "Haven't yet added support for moving mesh to DG-subcell. This "
      85             :            "should be easy to generalize, but we will want to consider "
      86             :            "storing the mesh velocity on the faces instead of "
      87             :            "re-slicing/projecting.");
      88             : 
      89             :     FixedHashMap<maximum_number_of_neighbors(3),
      90             :                  std::pair<Direction<3>, ElementId<3>>, std::vector<double>,
      91             :                  boost::hash<std::pair<Direction<3>, ElementId<3>>>>
      92             :         nhbr_package_data{};
      93             : 
      94             :     const auto& nhbr_subcell_data =
      95             :         db::get<evolution::dg::subcell::Tags::
      96             :                     NeighborDataForReconstructionAndRdmpTci<3>>(box);
      97             :     const Mesh<3>& subcell_mesh =
      98             :         db::get<evolution::dg::subcell::Tags::Mesh<3>>(box);
      99             :     const Mesh<3>& dg_mesh = db::get<domain::Tags::Mesh<3>>(box);
     100             : 
     101             :     const auto volume_prims = evolution::dg::subcell::fd::project(
     102             :         db::get<typename System::primitive_variables_tag>(box), dg_mesh,
     103             :         subcell_mesh.extents());
     104             : 
     105             :     const auto& recons =
     106             :         db::get<grmhd::ValenciaDivClean::fd::Tags::Reconstructor>(box);
     107             :     const auto& boundary_correction =
     108             :         db::get<evolution::Tags::BoundaryCorrection<System>>(box);
     109             :     using derived_boundary_corrections =
     110             :         typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
     111             :     tmpl::for_each<
     112             :         derived_boundary_corrections>([&box, &boundary_correction, &dg_mesh,
     113             :                                        &mortars_to_reconstruct_to,
     114             :                                        &nhbr_package_data, &nhbr_subcell_data,
     115             :                                        &recons, &subcell_mesh, &volume_prims](
     116             :                                           auto derived_correction_v) {
     117             :       using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
     118             :       if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
     119             :         using dg_package_data_temporary_tags =
     120             :             typename DerivedCorrection::dg_package_data_temporary_tags;
     121             :         using dg_package_data_argument_tags = tmpl::append<
     122             :             evolved_vars_tags, recons_prim_tags, fluxes_tags,
     123             :             tmpl::remove_duplicates<tmpl::push_back<
     124             :                 dg_package_data_temporary_tags, gr::Tags::SpatialMetric<3>,
     125             :                 gr::Tags::SqrtDetSpatialMetric<DataVector>,
     126             :                 gr::Tags::InverseSpatialMetric<3, Frame::Inertial, DataVector>,
     127             :                 evolution::dg::Actions::detail::NormalVector<3>>>>;
     128             : 
     129             :         const auto& element = db::get<domain::Tags::Element<3>>(box);
     130             :         const auto& eos = get<hydro::Tags::EquationOfStateBase>(box);
     131             : 
     132             :         using dg_package_field_tags =
     133             :             typename DerivedCorrection::dg_package_field_tags;
     134             :         Variables<dg_package_data_argument_tags> vars_on_face{0};
     135             :         Variables<dg_package_field_tags> packaged_data{0};
     136             :         for (const auto& mortar_id : mortars_to_reconstruct_to) {
     137             :           const Direction<3>& direction = mortar_id.first;
     138             : 
     139             :           Index<3> extents = subcell_mesh.extents();
     140             :           // Switch to face-centered instead of cell-centered points on the FD.
     141             :           // There are num_cell_centered+1 face-centered points.
     142             :           ++extents[direction.dimension()];
     143             : 
     144             :           // Computed prims and cons on face via reconstruction
     145             :           const size_t num_face_pts = subcell_mesh.extents()
     146             :                                           .slice_away(direction.dimension())
     147             :                                           .product();
     148             :           vars_on_face.initialize(num_face_pts);
     149             :           // Copy spacetime vars over from volume.
     150             :           using spacetime_vars_to_copy = tmpl::list<
     151             :               gr::Tags::Lapse<DataVector>,
     152             :               gr::Tags::Shift<3, Frame::Inertial, DataVector>,
     153             :               gr::Tags::SpatialMetric<3>,
     154             :               gr::Tags::SqrtDetSpatialMetric<DataVector>,
     155             :               gr::Tags::InverseSpatialMetric<3, Frame::Inertial, DataVector>>;
     156             :           tmpl::for_each<spacetime_vars_to_copy>(
     157             :               [&direction, &extents, &vars_on_face,
     158             :                &spacetime_vars_on_faces =
     159             :                    db::get<evolution::dg::subcell::Tags::OnSubcellFaces<
     160             :                        typename System::flux_spacetime_variables_tag, 3>>(box)](
     161             :                   auto tag_v) {
     162             :                 using tag = tmpl::type_from<decltype(tag_v)>;
     163             :                 data_on_slice(make_not_null(&get<tag>(vars_on_face)),
     164             :                               get<tag>(gsl::at(spacetime_vars_on_faces,
     165             :                                                direction.dimension())),
     166             :                               extents, direction.dimension(),
     167             :                               direction.side() == Side::Lower
     168             :                                   ? 0
     169             :                                   : extents[direction.dimension()] - 1);
     170             :               });
     171             : 
     172             :           call_with_dynamic_type<void, typename grmhd::ValenciaDivClean::fd::
     173             :                                            Reconstructor::creatable_classes>(
     174             :               &recons,
     175             :               [&element, &eos, &mortar_id, &nhbr_subcell_data, &subcell_mesh,
     176             :                &vars_on_face, &volume_prims](const auto& reconstructor) {
     177             :                 reconstructor->reconstruct_fd_neighbor(
     178             :                     make_not_null(&vars_on_face), volume_prims, eos, element,
     179             :                     nhbr_subcell_data, subcell_mesh, mortar_id.first);
     180             :               });
     181             : 
     182             :           grmhd::ValenciaDivClean::subcell::compute_fluxes(
     183             :               make_not_null(&vars_on_face));
     184             : 
     185             :           // Note: since the spacetime isn't dynamical we don't need to
     186             :           // worry about different normal vectors on the different sides
     187             :           // of the element. If the spacetime is dynamical, then the normal
     188             :           // covector needs to be sent, or at least the inverse spatial metric
     189             :           // from the neighbor so that the normal covector can be computed
     190             :           // correctly.
     191             :           tnsr::i<DataVector, 3, Frame::Inertial> normal_covector =
     192             :               get<evolution::dg::Tags::NormalCovector<3>>(
     193             :                   *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<3>>(
     194             :                        box)
     195             :                        .at(mortar_id.first));
     196             :           for (auto& t : normal_covector) {
     197             :             t *= -1.0;
     198             :           }
     199             :           // Note: Only need to do the projection in 2d and 3d, but GRMHD is
     200             :           // always 3d currently.
     201             :           const auto dg_normal_covector = normal_covector;
     202             :           for (size_t i = 0; i < 3; ++i) {
     203             :             normal_covector.get(i) = evolution::dg::subcell::fd::project(
     204             :                 dg_normal_covector.get(i),
     205             :                 dg_mesh.slice_away(mortar_id.first.dimension()),
     206             :                 subcell_mesh.extents().slice_away(mortar_id.first.dimension()));
     207             :           }
     208             : 
     209             :           // Compute the packaged data
     210             :           packaged_data.initialize(num_face_pts);
     211             :           using dg_package_data_projected_tags = tmpl::append<
     212             :               evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
     213             :               typename DerivedCorrection::dg_package_data_primitive_tags>;
     214             :           evolution::dg::Actions::detail::dg_package_data<System>(
     215             :               make_not_null(&packaged_data),
     216             :               dynamic_cast<const DerivedCorrection&>(boundary_correction),
     217             :               vars_on_face, normal_covector, {std::nullopt}, box,
     218             :               typename DerivedCorrection::dg_package_data_volume_tags{},
     219             :               dg_package_data_projected_tags{});
     220             : 
     221             :           // Reconstruct the DG solution.
     222             :           // Really we should be solving the boundary correction and
     223             :           // then reconstructing, but away from a shock this doesn't
     224             :           // matter.
     225             :           auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct(
     226             :               packaged_data, dg_mesh.slice_away(mortar_id.first.dimension()),
     227             :               subcell_mesh.extents().slice_away(mortar_id.first.dimension()));
     228             :           nhbr_package_data[mortar_id] = std::vector<double>{
     229             :               dg_packaged_data.data(),
     230             :               dg_packaged_data.data() + dg_packaged_data.size()};
     231             :         }
     232             :       }
     233             :     });
     234             : 
     235             :     return nhbr_package_data;
     236             :   }
     237             : };
     238             : }  // namespace grmhd::ValenciaDivClean::subcell

Generated by: LCOV version 1.14