Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <boost/functional/hash.hpp>
       8             : #include <cstddef>
       9             : #include <unordered_map>
      10             : #include <utility>
      11             : 
      12             : #include "DataStructures/DataVector.hpp"
      13             : #include "DataStructures/Index.hpp"
      14             : #include "DataStructures/SliceIterator.hpp"
      15             : #include "DataStructures/Tensor/Tensor.hpp"
      16             : #include "DataStructures/Variables.hpp"
      17             : #include "Domain/Structure/Direction.hpp"
      18             : #include "Domain/Structure/DirectionalId.hpp"
      19             : #include "Domain/Structure/Element.hpp"
      20             : #include "Domain/Structure/ElementId.hpp"
      21             : #include "Evolution/DgSubcell/Projection.hpp"
      22             : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
      23             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      24             : #include "Utilities/ErrorHandling/Assert.hpp"
      25             : #include "Utilities/Gsl.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : 
      28             : namespace evolution::dg::subcell {
      29             : /*!
      30             :  * \brief Project the DG package data to the subcells. Data received from a
      31             :  * neighboring element doing DG is always projected, while the data we sent to
      32             :  * our neighbors before doing a rollback from DG to subcell is only projected if
      33             :  * `OverwriteInternalMortarData` is `true`.
      34             :  *
      35             :  * In order for the hybrid DG-FD/FV scheme to be conservative between elements
      36             :  * using DG and elements using subcell, the boundary terms must be the same on
      37             :  * both elements. In practice this means the boundary corrections \f$G+D\f$ must
      38             :  * be computed on the same grid. Consider the element doing subcell which
      39             :  * receives data from an element doing DG. In this case the DG element's
      40             :  * ingredients going into \f$G+D\f$ are projected to the subcells and then
      41             :  * \f$G+D\f$ are computed on the subcells. Similarly, for strict conservation
      42             :  * the element doing DG must first project the data it sent to the neighbor to
      43             :  * the subcells, then compute \f$G+D\f$ on the subcells, and finally reconstrct
      44             :  * \f$G+D\f$ back to the DG grid before lifting \f$G+D\f$ to the volume.
      45             :  *
      46             :  * This function updates the `packaged_data` (ingredients into \f$G+D\f$)
      47             :  * received by an element doing subcell by projecting the neighbor's DG data
      48             :  * onto the subcells. Note that this is only half of what is required for strict
      49             :  * conservation, the DG element must also compute \f$G+D\f$ on the subcells.
      50             :  * Note that we currently do not perform the other half of the correction
      51             :  * needed to be strictly conservative.
      52             :  *
      53             :  * If we are retaking a time step after the DG step failed then maintaining
      54             :  * conservation requires additional care. If `OverwriteInternalMortarData` is
      55             :  * `true` then the local (the element switching from DG to subcell) ingredients
      56             :  * into \f$G+D\f$ are projected and overwrite the data computed from
      57             :  * the FD reconstruction to the interface. However, even this is insufficient to
      58             :  * guarantee conservation. To guarantee conservation (which we do not currently
      59             :  * do) the correction \f$G+D\f$ must be computed on the DG grid and then
      60             :  * projected to the subcells.
      61             :  *
      62             :  * Note that our practical experience shows that since the DG-subcell hybrid
      63             :  * scheme switches to the subcell solver _before_ the local solution contains
      64             :  * discontinuities, strict conservation is not necessary between DG and FD/FV
      65             :  * regions. This was also observed with a block-adaptive finite difference AMR
      66             :  * code \cite CHEN2016604
      67             :  *
      68             :  * The variable `variables_to_offset_in_dg_grid` is used in cases like the
      69             :  * combined generalized harmonic and GRMHD system where the DG grid uses
      70             :  * boundary corrections for both GH and GRMHD, but the subcell grid only does
      71             :  * GRMHD.
      72             :  */
      73             : template <bool OverwriteInternalMortarData, size_t Dim,
      74             :           typename DgPackageFieldTags>
      75           1 : void correct_package_data(
      76             :     const gsl::not_null<Variables<DgPackageFieldTags>*> lower_packaged_data,
      77             :     const gsl::not_null<Variables<DgPackageFieldTags>*> upper_packaged_data,
      78             :     const size_t logical_dimension_to_operate_in, const Element<Dim>& element,
      79             :     const Mesh<Dim>& subcell_volume_mesh,
      80             :     const std::unordered_map<DirectionalId<Dim>, evolution::dg::MortarData<Dim>,
      81             :                              boost::hash<DirectionalId<Dim>>>& mortar_data,
      82             :     const size_t variables_to_offset_in_dg_grid) {
      83             :   const Direction<Dim> upper_direction{logical_dimension_to_operate_in,
      84             :                                        Side::Upper};
      85             :   const Direction<Dim> lower_direction{logical_dimension_to_operate_in,
      86             :                                        Side::Lower};
      87             :   const bool has_upper_neighbor = element.neighbors().contains(upper_direction);
      88             :   const bool has_lower_neighbor = element.neighbors().contains(lower_direction);
      89             :   const DirectionalId<Dim> upper_mortar_id =
      90             :       has_upper_neighbor
      91             :           ? DirectionalId<Dim>{upper_direction,
      92             :                                *element.neighbors().at(upper_direction).begin()}
      93             :           : DirectionalId<Dim>{};
      94             :   const DirectionalId<Dim> lower_mortar_id =
      95             :       has_lower_neighbor
      96             :           ? DirectionalId<Dim>{lower_direction,
      97             :                                *element.neighbors().at(lower_direction).begin()}
      98             :           : DirectionalId<Dim>{};
      99             : 
     100             :   Index<Dim> subcell_extents_with_faces = subcell_volume_mesh.extents();
     101             :   ++subcell_extents_with_faces[logical_dimension_to_operate_in];
     102             :   const Mesh<Dim - 1>& subcell_face_mesh =
     103             :       subcell_volume_mesh.slice_away(logical_dimension_to_operate_in);
     104             : 
     105             :   const auto project_dg_data_to_subcells =
     106             :       [logical_dimension_to_operate_in, &subcell_extents_with_faces,
     107             :        &subcell_face_mesh, variables_to_offset_in_dg_grid](
     108             :           const gsl::not_null<Variables<DgPackageFieldTags>*>
     109             :               subcell_packaged_data,
     110             :           const size_t subcell_index, const Mesh<Dim - 1>& neighbor_face_mesh,
     111             :           const DataVector& neighbor_data) {
     112             :         const size_t dg_variables_offset_size =
     113             :             variables_to_offset_in_dg_grid *
     114             :             neighbor_face_mesh.number_of_grid_points();
     115             :         const double* slice_data =
     116             :             // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
     117             :             neighbor_data.data() + dg_variables_offset_size;
     118             :         // Warning: projected_data can't be inside the `if constexpr` since that
     119             :         // would lead to a dangling pointer.
     120             :         DataVector projected_data{};
     121             :         if constexpr (Dim > 1) {
     122             :           projected_data = DataVector{
     123             :               Variables<DgPackageFieldTags>::number_of_independent_components *
     124             :               subcell_face_mesh.number_of_grid_points()};
     125             :           evolution::dg::subcell::fd::detail::project_impl(
     126             :               gsl::make_span(projected_data.data(), projected_data.size()),
     127             :               gsl::make_span(
     128             :                   slice_data,
     129             :                   Variables<
     130             :                       DgPackageFieldTags>::number_of_independent_components *
     131             :                       neighbor_face_mesh.number_of_grid_points()),
     132             :               neighbor_face_mesh, subcell_face_mesh.extents());
     133             :           slice_data = projected_data.data();
     134             :         } else {
     135             :           (void)subcell_face_mesh;
     136             :           (void)projected_data;
     137             :         }
     138             :         const size_t volume_grid_points = subcell_extents_with_faces.product();
     139             :         const size_t slice_grid_points =
     140             :             subcell_extents_with_faces
     141             :                 .slice_away(logical_dimension_to_operate_in)
     142             :                 .product();
     143             :         double* const volume_data = subcell_packaged_data->data();
     144             :         for (SliceIterator si(subcell_extents_with_faces,
     145             :                               logical_dimension_to_operate_in, subcell_index);
     146             :              si; ++si) {
     147             :           for (size_t i = 0;
     148             :                i <
     149             :                Variables<DgPackageFieldTags>::number_of_independent_components;
     150             :                ++i) {
     151             :             // clang-tidy: do not use pointer arithmetic
     152             :             volume_data[si.volume_offset() +
     153             :                         i * volume_grid_points] =  // NOLINT
     154             :                 slice_data[si.slice_offset() +
     155             :                            i * slice_grid_points];  // NOLINT
     156             :           }
     157             :         }
     158             :       };
     159             : 
     160             :   // Project DG data to the subcells
     161             :   if (auto neighbor_mortar_data_it = mortar_data.find(upper_mortar_id);
     162             :       has_upper_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
     163             :     if (neighbor_mortar_data_it->second.neighbor_mortar_data().has_value()) {
     164             :       project_dg_data_to_subcells(
     165             :           upper_packaged_data,
     166             :           subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
     167             :           neighbor_mortar_data_it->second.neighbor_mortar_data()->first,
     168             :           neighbor_mortar_data_it->second.neighbor_mortar_data()->second);
     169             :     }
     170             :     if constexpr (OverwriteInternalMortarData) {
     171             :       if (neighbor_mortar_data_it->second.local_mortar_data().has_value()) {
     172             :         project_dg_data_to_subcells(
     173             :             lower_packaged_data,
     174             :             subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
     175             :             neighbor_mortar_data_it->second.local_mortar_data()->first,
     176             :             neighbor_mortar_data_it->second.local_mortar_data()->second);
     177             :       }
     178             :     }
     179             :   }
     180             :   if (auto neighbor_mortar_data_it = mortar_data.find(lower_mortar_id);
     181             :       has_lower_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
     182             :     if (neighbor_mortar_data_it->second.neighbor_mortar_data().has_value()) {
     183             :       project_dg_data_to_subcells(
     184             :           lower_packaged_data, 0,
     185             :           neighbor_mortar_data_it->second.neighbor_mortar_data()->first,
     186             :           neighbor_mortar_data_it->second.neighbor_mortar_data()->second);
     187             :     }
     188             :     if constexpr (OverwriteInternalMortarData) {
     189             :       if (neighbor_mortar_data_it->second.local_mortar_data().has_value()) {
     190             :         project_dg_data_to_subcells(
     191             :             upper_packaged_data, 0,
     192             :             neighbor_mortar_data_it->second.local_mortar_data()->first,
     193             :             neighbor_mortar_data_it->second.local_mortar_data()->second);
     194             :       }
     195             :     }
     196             :   }
     197             : }
     198             : }  // namespace evolution::dg::subcell