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 <utility>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      12             : #include "DataStructures/DataBox/Prefixes.hpp"
      13             : #include "DataStructures/DataVector.hpp"
      14             : #include "DataStructures/Index.hpp"
      15             : #include "Domain/Structure/Direction.hpp"
      16             : #include "Domain/Structure/DirectionMap.hpp"
      17             : #include "Domain/Structure/Element.hpp"
      18             : #include "Domain/Structure/ElementId.hpp"
      19             : #include "Domain/Structure/OrientationMapHelpers.hpp"
      20             : #include "Domain/Tags.hpp"
      21             : #include "Domain/Tags/NeighborMesh.hpp"
      22             : #include "Evolution/DgSubcell/Projection.hpp"
      23             : #include "Evolution/DgSubcell/RdmpTci.hpp"
      24             : #include "Evolution/DgSubcell/RdmpTciData.hpp"
      25             : #include "Evolution/DgSubcell/SliceData.hpp"
      26             : #include "Evolution/DgSubcell/SubcellOptions.hpp"
      27             : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
      28             : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
      29             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      30             : #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
      31             : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
      32             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      33             : #include "Utilities/Algorithm.hpp"
      34             : #include "Utilities/ErrorHandling/Assert.hpp"
      35             : #include "Utilities/Gsl.hpp"
      36             : #include "Utilities/TMPL.hpp"
      37             : 
      38             : namespace evolution::dg::subcell {
      39             : /*!
      40             :  * \brief Add local data for our and our neighbor's relaxed discrete maximum
      41             :  * principle troubled-cell indicator, and for reconstruction.
      42             :  *
      43             :  * The local maximum and minimum of the evolved variables is added to
      44             :  * `Tags::DataForRdmpTci` for the RDMP TCI. Then the
      45             :  * data needed by neighbor elements to do reconstruction on the FD grid is sent.
      46             :  * The data to be sent is computed in the mutator
      47             :  * `Metavariables::SubcellOptions::GhostVariables`, which returns a
      48             :  * `Variables` of the tensors to send to the neighbors. The main reason for
      49             :  * having the mutator `GhostVariables` is to allow sending primitive or
      50             :  * characteristic variables for reconstruction.
      51             :  *
      52             :  * \note If all neighbors are using DG then we send our DG volume data _without_
      53             :  * orienting it. This elides the expense of projection and slicing. If any
      54             :  * neighbors are doing FD, we project and slice to all neighbors. A future
      55             :  * optimization would be to measure the cost of slicing data, and figure out how
      56             :  * many neighbors need to be doing FD before it's worth projecting and slicing
      57             :  * vs. just projecting to the ghost cells. Another optimization is to always
      58             :  * send DG volume data to neighbors that we think are doing DG, so that we can
      59             :  * elide any slicing or projection cost in that direction.
      60             :  */
      61             : template <typename Metavariables, typename DbTagsList, size_t Dim>
      62           1 : void prepare_neighbor_data(
      63             :     const gsl::not_null<DirectionMap<Dim, DataVector>*>
      64             :         all_neighbor_data_for_reconstruction,
      65             :     const gsl::not_null<Mesh<Dim>*> ghost_data_mesh,
      66             :     const gsl::not_null<db::DataBox<DbTagsList>*> box,
      67             :     [[maybe_unused]] const Variables<db::wrap_tags_in<
      68             :         ::Tags::Flux, typename Metavariables::system::flux_variables,
      69             :         tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes) {
      70             :   const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(*box);
      71             :   const Mesh<Dim>& subcell_mesh =
      72             :       db::get<::evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
      73             :   const auto& element = db::get<::domain::Tags::Element<Dim>>(*box);
      74             : 
      75             :   const std::unordered_set<Direction<Dim>>& directions_to_slice =
      76             :       element.internal_boundaries();
      77             : 
      78             :   const auto& neighbor_meshes = db::get<domain::Tags::NeighborMesh<Dim>>(*box);
      79             :   const subcell::RdmpTciData& rdmp_data =
      80             :       db::get<subcell::Tags::DataForRdmpTci>(*box);
      81             :   const ::fd::DerivativeOrder fd_derivative_order =
      82             :       get<evolution::dg::subcell::Tags::SubcellOptions<Dim>>(*box)
      83             :           .finite_difference_derivative_order();
      84             :   const size_t rdmp_size = rdmp_data.max_variables_values.size() +
      85             :                            rdmp_data.min_variables_values.size();
      86             :   const size_t extra_size_for_ghost_data =
      87             :       (fd_derivative_order == ::fd::DerivativeOrder::Two
      88             :            ? 0
      89             :            : volume_fluxes.size());
      90             : 
      91             :   if (DataVector ghost_variables =
      92             :           [&box, &extra_size_for_ghost_data, &rdmp_size, &volume_fluxes]() {
      93             :             if (extra_size_for_ghost_data == 0) {
      94             :               return db::mutate_apply(
      95             :                   typename Metavariables::SubcellOptions::GhostVariables{}, box,
      96             :                   rdmp_size);
      97             :             } else {
      98             :               DataVector ghost_vars = db::mutate_apply(
      99             :                   typename Metavariables::SubcellOptions::GhostVariables{}, box,
     100             :                   rdmp_size + extra_size_for_ghost_data);
     101             :               std::copy(
     102             :                   volume_fluxes.data(),
     103             :                   std::next(volume_fluxes.data(),
     104             :                             static_cast<std::ptrdiff_t>(volume_fluxes.size())),
     105             :                   std::prev(ghost_vars.end(),
     106             :                             static_cast<std::ptrdiff_t>(
     107             :                                 rdmp_size + extra_size_for_ghost_data)));
     108             :               return ghost_vars;
     109             :             }
     110             :           }();
     111             :       alg::all_of(neighbor_meshes,
     112             :                   [](const auto& directional_element_id_and_mesh) {
     113             :                     ASSERT(directional_element_id_and_mesh.second.basis(0) !=
     114             :                                Spectral::Basis::Chebyshev,
     115             :                            "Don't yet support Chebyshev basis with DG-FD");
     116             :                     return directional_element_id_and_mesh.second.basis(0) ==
     117             :                            Spectral::Basis::Legendre;
     118             :                   })) {
     119             :     *ghost_data_mesh = dg_mesh;
     120             :     const size_t total_to_slice = directions_to_slice.size();
     121             :     size_t slice_count = 0;
     122             :     for (const auto& direction : directions_to_slice) {
     123             :       ++slice_count;
     124             :       // Move instead of copy on the last iteration. Elides a memory
     125             :       // allocation and copy.
     126             :       [[maybe_unused]] const auto insert_result =
     127             :           UNLIKELY(slice_count == total_to_slice)
     128             :               ? all_neighbor_data_for_reconstruction->insert_or_assign(
     129             :                     direction, std::move(ghost_variables))
     130             :               : all_neighbor_data_for_reconstruction->insert_or_assign(
     131             :                     direction, ghost_variables);
     132             :       ASSERT(all_neighbor_data_for_reconstruction->size() == total_to_slice or
     133             :                  insert_result.second,
     134             :              "Failed to insert the neighbor data in direction " << direction);
     135             :     }
     136             :   } else {
     137             :     *ghost_data_mesh = subcell_mesh;
     138             :     const size_t ghost_zone_size =
     139             :         db::get<evolution::dg::subcell::Tags::Reconstructor>(*box)
     140             :             .ghost_zone_size();
     141             : 
     142             :     const DataVector data_to_project{};
     143             :     make_const_view(make_not_null(&data_to_project), ghost_variables, 0,
     144             :                     ghost_variables.size() - rdmp_size);
     145             :     const DataVector projected_data = evolution::dg::subcell::fd::project(
     146             :         data_to_project, dg_mesh, subcell_mesh.extents());
     147             :     typename evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<
     148             :         Dim>::type directions_not_to_slice{};
     149             :     for (const auto& [directional_element_id, _] :
     150             :          db::get<evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<
     151             :              Dim>>(*box)) {
     152             :       directions_not_to_slice[directional_element_id] = std::nullopt;
     153             :     }
     154             :     *all_neighbor_data_for_reconstruction = subcell::slice_data(
     155             :         projected_data, db::get<subcell::Tags::Mesh<Dim>>(*box).extents(),
     156             :         ghost_zone_size, directions_to_slice, rdmp_size,
     157             :         directions_not_to_slice);
     158             :     for (const auto& [directional_element_id, interpolator] :
     159             :          db::get<evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<
     160             :              Dim>>(*box)) {
     161             :       // NOTE: no orienting is needed, that's done by the interpolation.
     162             :       ASSERT(interpolator.has_value(),
     163             :              "All interpolators must have values. The std::optional is used "
     164             :              "for some local optimizations.");
     165             :       DataVector& data_in_direction = all_neighbor_data_for_reconstruction->at(
     166             :           directional_element_id.direction);
     167             :       gsl::span result_in_direction{data_in_direction.data(),
     168             :                                     data_in_direction.size() - rdmp_size};
     169             :       interpolator->interpolate(
     170             :           make_not_null(&result_in_direction),
     171             :           {data_to_project.data(), data_to_project.size()});
     172             :     }
     173             :   }
     174             : 
     175             :   // Note: The RDMP TCI data must be filled by the TCIs before getting to this
     176             :   // call. That means once in the initial data and then in both the DG and FD
     177             :   // TCIs.
     178             :   for (const auto& [direction, neighbors] : element.neighbors()) {
     179             :     // Add the RDMP TCI data to what we will be sending.
     180             :     // Note that this is added _after_ the reconstruction data has been
     181             :     // re-oriented (in the case where we have neighbors doing FD).
     182             :     DataVector& neighbor_data =
     183             :         all_neighbor_data_for_reconstruction->at(direction);
     184             :     std::copy(rdmp_data.max_variables_values.begin(),
     185             :               rdmp_data.max_variables_values.end(),
     186             :               std::prev(neighbor_data.end(), static_cast<int>(rdmp_size)));
     187             :     std::copy(
     188             :         rdmp_data.min_variables_values.begin(),
     189             :         rdmp_data.min_variables_values.end(),
     190             :         std::prev(neighbor_data.end(),
     191             :                   static_cast<int>(rdmp_data.min_variables_values.size())));
     192             :   }
     193             : }
     194             : }  // namespace evolution::dg::subcell