SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/NewtonianEuler/Subcell - NeighborPackagedData.hpp Hit Total Coverage
Commit: 3528f39684ab2ee5d689cee48331779e729b0a07 Lines: 1 3 33.3 %
Date: 2024-02-27 07:22:14
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 <cstddef>
       7             : #include <optional>
       8             : #include <type_traits>
       9             : #include <utility>
      10             : #include <vector>
      11             : 
      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/Index.hpp"
      17             : #include "DataStructures/Tensor/Tensor.hpp"
      18             : #include "DataStructures/Variables.hpp"
      19             : #include "DataStructures/VariablesTag.hpp"
      20             : #include "Domain/Structure/Direction.hpp"
      21             : #include "Domain/Structure/DirectionalIdMap.hpp"
      22             : #include "Domain/Structure/Element.hpp"
      23             : #include "Domain/Structure/ElementId.hpp"
      24             : #include "Domain/Tags.hpp"
      25             : #include "Domain/TagsTimeDependent.hpp"
      26             : #include "Evolution/BoundaryCorrectionTags.hpp"
      27             : #include "Evolution/DgSubcell/Projection.hpp"
      28             : #include "Evolution/DgSubcell/Reconstruction.hpp"
      29             : #include "Evolution/DgSubcell/ReconstructionMethod.hpp"
      30             : #include "Evolution/DgSubcell/SubcellOptions.hpp"
      31             : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
      32             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      33             : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
      34             : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
      35             : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
      36             : #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/BoundaryCorrection.hpp"
      37             : #include "Evolution/Systems/NewtonianEuler/BoundaryCorrections/Factory.hpp"
      38             : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
      39             : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Tag.hpp"
      40             : #include "Evolution/Systems/NewtonianEuler/Subcell/ComputeFluxes.hpp"
      41             : #include "Evolution/Systems/NewtonianEuler/System.hpp"
      42             : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
      43             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      44             : #include "Parallel/Tags/Metavariables.hpp"
      45             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      46             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      47             : #include "Utilities/CallWithDynamicType.hpp"
      48             : #include "Utilities/Gsl.hpp"
      49             : #include "Utilities/TMPL.hpp"
      50             : 
      51             : namespace NewtonianEuler::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 <size_t Dim, typename DbTagsList>
      68           0 :   static DirectionalIdMap<Dim, DataVector> apply(
      69             :       const db::DataBox<DbTagsList>& box,
      70             :       const std::vector<DirectionalId<Dim>>& mortars_to_reconstruct_to) {
      71             :     using system = typename std::decay_t<decltype(
      72             :         db::get<Parallel::Tags::Metavariables>(box))>::system;
      73             :     using evolved_vars_tag = typename system::variables_tag;
      74             :     using evolved_vars_tags = typename evolved_vars_tag::tags_list;
      75             :     using prim_tags = typename system::primitive_variables_tag::tags_list;
      76             :     using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
      77             :                                          tmpl::size_t<Dim>, Frame::Inertial>;
      78             : 
      79             :     ASSERT(not db::get<domain::Tags::MeshVelocity<Dim>>(box).has_value(),
      80             :            "Haven't yet added support for moving mesh to DG-subcell. This "
      81             :            "should be easy to generalize, but we will want to consider "
      82             :            "storing the mesh velocity on the faces instead of "
      83             :            "re-slicing/projecting.");
      84             : 
      85             :     DirectionalIdMap<Dim, DataVector> neighbor_package_data{};
      86             :     if (mortars_to_reconstruct_to.empty()) {
      87             :       return neighbor_package_data;
      88             :     }
      89             : 
      90             :     const auto& ghost_subcell_data =
      91             :         db::get<evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>>(
      92             :             box);
      93             :     const Mesh<Dim>& subcell_mesh =
      94             :         db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box);
      95             :     const Mesh<Dim>& dg_mesh = db::get<domain::Tags::Mesh<Dim>>(box);
      96             :     const auto& subcell_options =
      97             :         db::get<evolution::dg::subcell::Tags::SubcellOptions<Dim>>(box);
      98             : 
      99             :     // Note: we need to compare if projecting the entire mesh or only ghost
     100             :     // zones needed is faster. This probably depends on the number of neighbors
     101             :     // we have doing FD.
     102             :     const auto volume_prims = evolution::dg::subcell::fd::project(
     103             :         db::get<typename system::primitive_variables_tag>(box), dg_mesh,
     104             :         subcell_mesh.extents());
     105             : 
     106             :     const auto& recons =
     107             :         db::get<NewtonianEuler::fd::Tags::Reconstructor<Dim>>(box);
     108             :     const auto& boundary_correction =
     109             :         db::get<evolution::Tags::BoundaryCorrection<system>>(box);
     110             :     using derived_boundary_corrections =
     111             :         typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
     112             :     tmpl::for_each<
     113             :         derived_boundary_corrections>([&box, &boundary_correction, &dg_mesh,
     114             :                                        &mortars_to_reconstruct_to,
     115             :                                        &neighbor_package_data,
     116             :                                        &ghost_subcell_data, &recons,
     117             :                                        &subcell_mesh, &subcell_options,
     118             :                                        &volume_prims](
     119             :                                           auto derived_correction_v) {
     120             :       using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
     121             :       if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
     122             :         using dg_package_data_temporary_tags =
     123             :             typename DerivedCorrection::dg_package_data_temporary_tags;
     124             :         using dg_package_data_argument_tags =
     125             :             tmpl::append<evolved_vars_tags, prim_tags, fluxes_tags,
     126             :                          dg_package_data_temporary_tags>;
     127             : 
     128             :         const auto& element = db::get<domain::Tags::Element<Dim>>(box);
     129             :         const auto& eos = get<hydro::Tags::EquationOfStateBase>(box);
     130             : 
     131             :         using dg_package_field_tags =
     132             :             typename DerivedCorrection::dg_package_field_tags;
     133             :         Variables<dg_package_data_argument_tags> vars_on_face;
     134             :         Variables<dg_package_field_tags> packaged_data;
     135             :         for (const auto& mortar_id : mortars_to_reconstruct_to) {
     136             :           const Direction<Dim>& direction = mortar_id.direction;
     137             : 
     138             :           Index<Dim> extents = subcell_mesh.extents();
     139             :           // Switch to face-centered instead of cell-centered points on the FD.
     140             :           // There are num_cell_centered+1 face-centered points.
     141             :           ++extents[direction.dimension()];
     142             : 
     143             :           // Computed prims and cons on face via reconstruction
     144             :           const size_t num_face_pts = subcell_mesh.extents()
     145             :                                           .slice_away(direction.dimension())
     146             :                                           .product();
     147             :           vars_on_face.initialize(num_face_pts);
     148             : 
     149             :           call_with_dynamic_type<void,
     150             :                                  typename NewtonianEuler::fd::Reconstructor<
     151             :                                      Dim>::creatable_classes>(
     152             :               &recons,
     153             :               [&element, &eos, &mortar_id, &ghost_subcell_data, &subcell_mesh,
     154             :                &vars_on_face, &volume_prims](const auto& reconstructor) {
     155             :                 reconstructor->reconstruct_fd_neighbor(
     156             :                     make_not_null(&vars_on_face), volume_prims, eos, element,
     157             :                     ghost_subcell_data, subcell_mesh, mortar_id.direction);
     158             :               });
     159             : 
     160             :           NewtonianEuler::subcell::compute_fluxes<Dim>(
     161             :               make_not_null(&vars_on_face));
     162             : 
     163             :           tnsr::i<DataVector, Dim, Frame::Inertial> normal_covector = get<
     164             :               evolution::dg::Tags::NormalCovector<Dim>>(
     165             :               *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(
     166             :                    box)
     167             :                    .at(mortar_id.direction));
     168             :           for (auto& t : normal_covector) {
     169             :             t *= -1.0;
     170             :           }
     171             :           if constexpr (Dim > 1) {
     172             :             const auto dg_normal_covector = normal_covector;
     173             :             for (size_t i = 0; i < Dim; ++i) {
     174             :               normal_covector.get(i) = evolution::dg::subcell::fd::project(
     175             :                   dg_normal_covector.get(i),
     176             :                   dg_mesh.slice_away(mortar_id.direction.dimension()),
     177             :                   subcell_mesh.extents().slice_away(
     178             :                       mortar_id.direction.dimension()));
     179             :             }
     180             :           }
     181             : 
     182             :           // Compute the packaged data
     183             :           packaged_data.initialize(num_face_pts);
     184             :           using dg_package_data_projected_tags = tmpl::append<
     185             :               evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
     186             :               typename DerivedCorrection::dg_package_data_primitive_tags>;
     187             :           evolution::dg::Actions::detail::dg_package_data<system>(
     188             :               make_not_null(&packaged_data),
     189             :               dynamic_cast<const DerivedCorrection&>(boundary_correction),
     190             :               vars_on_face, normal_covector, {std::nullopt}, box,
     191             :               typename DerivedCorrection::dg_package_data_volume_tags{},
     192             :               dg_package_data_projected_tags{});
     193             : 
     194             :           if constexpr (Dim == 1) {
     195             :             (void)dg_mesh;
     196             :             (void)subcell_options;
     197             :             // Make a view so we can use iterators with std::copy
     198             :             DataVector packaged_data_view{packaged_data.data(),
     199             :                                           packaged_data.size()};
     200             :             neighbor_package_data[mortar_id] = DataVector{packaged_data.size()};
     201             :             std::copy(packaged_data_view.begin(), packaged_data_view.end(),
     202             :                       neighbor_package_data[mortar_id].begin());
     203             :           } else {
     204             :             // Reconstruct the DG solution.
     205             :             // Really we should be solving the boundary correction and
     206             :             // then reconstructing, but away from a shock this doesn't
     207             :             // matter.
     208             :             auto dg_packaged_data = evolution::dg::subcell::fd::reconstruct(
     209             :                 packaged_data,
     210             :                 dg_mesh.slice_away(mortar_id.direction.dimension()),
     211             :                 subcell_mesh.extents().slice_away(
     212             :                     mortar_id.direction.dimension()),
     213             :                 subcell_options.reconstruction_method());
     214             :             // Make a view so we can use iterators with std::copy
     215             :             DataVector dg_packaged_data_view{dg_packaged_data.data(),
     216             :                                              dg_packaged_data.size()};
     217             :             neighbor_package_data[mortar_id] =
     218             :                 DataVector{dg_packaged_data.size()};
     219             :             std::copy(dg_packaged_data_view.begin(),
     220             :                       dg_packaged_data_view.end(),
     221             :                       neighbor_package_data[mortar_id].begin());
     222             :           }
     223             :         }
     224             :       }
     225             :     });
     226             : 
     227             :     return neighbor_package_data;
     228             :   }
     229             : };
     230             : }  // namespace NewtonianEuler::subcell

Generated by: LCOV version 1.14