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 <memory>
       8             : #include <unordered_map>
       9             : #include <utility>
      10             : 
      11             : #include "DataStructures/Tags.hpp"
      12             : #include "Domain/Structure/Direction.hpp"
      13             : #include "Domain/Structure/DirectionMap.hpp"
      14             : #include "Domain/Structure/DirectionalId.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/TaggedTuple.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : template <size_t VolumeDim>
      21             : class Element;
      22             : template <size_t VolumeDim>
      23             : class ElementId;
      24             : template <size_t VolumeDim>
      25             : class Mesh;
      26             : enum class Side : uint8_t;
      27             : 
      28             : namespace boost {
      29             : template <class T>
      30             : struct hash;
      31             : }  // namespace boost
      32             : /// \endcond
      33             : 
      34             : namespace Limiters::Minmod_detail {
      35             : 
      36             : // Encodes the return status of the tvb_corrected_minmod function.
      37             : struct MinmodResult {
      38             :   const double value;
      39             :   const bool activated;
      40             : };
      41             : 
      42             : // The TVB-corrected minmod function, see e.g. Cockburn reference Eq. 2.26.
      43             : MinmodResult tvb_corrected_minmod(double a, double b, double c,
      44             :                                   double tvb_scale);
      45             : 
      46             : // Holds various optimization-related allocations for the Minmod TCI.
      47             : // There is no pup::er, because these allocations should be short-lived (i.e.,
      48             : // scoped within a single limiter invocation).
      49             : template <size_t VolumeDim>
      50             : class BufferWrapper {
      51             :  public:
      52             :   BufferWrapper() = delete;
      53             :   explicit BufferWrapper(const Mesh<VolumeDim>& mesh);
      54             : 
      55             :  private:
      56             :   std::unique_ptr<double[]> contiguous_boundary_buffer_;
      57             :   const std::pair<std::unique_ptr<std::pair<size_t, size_t>[]>,
      58             :                   std::array<std::pair<gsl::span<std::pair<size_t, size_t>>,
      59             :                                        gsl::span<std::pair<size_t, size_t>>>,
      60             :                              VolumeDim>>
      61             :       volume_and_slice_buffer_and_indices_ {};
      62             : 
      63             :  public:
      64             :   std::array<DataVector, VolumeDim> boundary_buffers{};
      65             :   const std::array<std::pair<gsl::span<std::pair<size_t, size_t>>,
      66             :                              gsl::span<std::pair<size_t, size_t>>>,
      67             :                    VolumeDim>& volume_and_slice_indices{};
      68             : };
      69             : 
      70             : // In each direction, average the size of all different neighbors in that
      71             : // direction. Note that only the component of neighor_size that is normal
      72             : // to the face is needed (and, therefore, computed).
      73             : //
      74             : // Expects type `PackagedData` to contain a variable `element_size` that is a
      75             : // `std::array<double, VolumeDim>`.
      76             : template <size_t VolumeDim, typename PackagedData>
      77             : DirectionMap<VolumeDim, double> compute_effective_neighbor_sizes(
      78             :     const Element<VolumeDim>& element,
      79             :     const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
      80             :                              boost::hash<DirectionalId<VolumeDim>>>&
      81             :         neighbor_data) {
      82             :   DirectionMap<VolumeDim, double> result;
      83             :   for (const auto& dir : Direction<VolumeDim>::all_directions()) {
      84             :     const bool neighbors_in_this_dir = element.neighbors().contains(dir);
      85             :     if (neighbors_in_this_dir) {
      86             :       const double effective_neighbor_size = [&dir, &element,
      87             :                                               &neighbor_data]() {
      88             :         const size_t dim = dir.dimension();
      89             :         const auto& neighbor_ids = element.neighbors().at(dir).ids();
      90             :         double size_accumulate = 0.;
      91             :         for (const auto& id : neighbor_ids) {
      92             :           size_accumulate += gsl::at(
      93             :               neighbor_data.at(DirectionalId<VolumeDim>{dir, id}).element_size,
      94             :               dim);
      95             :         }
      96             :         return size_accumulate / neighbor_ids.size();
      97             :       }();
      98             :       result.insert(std::make_pair(dir, effective_neighbor_size));
      99             :     }
     100             :   }
     101             :   return result;
     102             : }
     103             : 
     104             : // In each direction, average the mean of the specified tensor component over
     105             : // all different neighbors that direction. This produces one effective neighbor
     106             : // per direction.
     107             : //
     108             : // Expects type `PackagedData` to contain a variable `means` that is a
     109             : // `TaggedTuple<Tags::Mean<Tags>...>`. Tags... must contain Tag, the tag
     110             : // specifying the tensor to work with.
     111             : template <typename Tag, size_t VolumeDim, typename PackagedData>
     112             : DirectionMap<VolumeDim, double> compute_effective_neighbor_means(
     113             :     const size_t tensor_storage_index, const Element<VolumeDim>& element,
     114             :     const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
     115             :                              boost::hash<DirectionalId<VolumeDim>>>&
     116             :         neighbor_data) {
     117             :   DirectionMap<VolumeDim, double> result;
     118             :   for (const auto& dir : Direction<VolumeDim>::all_directions()) {
     119             :     const bool neighbors_in_this_dir = element.neighbors().contains(dir);
     120             :     if (neighbors_in_this_dir) {
     121             :       const double effective_neighbor_mean =
     122             :           [&dir, &element, &neighbor_data, &tensor_storage_index]() {
     123             :             const auto& neighbor_ids = element.neighbors().at(dir).ids();
     124             :             double mean_accumulate = 0.0;
     125             :             for (const auto& id : neighbor_ids) {
     126             :               mean_accumulate += tuples::get<::Tags::Mean<Tag>>(
     127             :                   neighbor_data.at(DirectionalId<VolumeDim>{dir, id})
     128             :                       .means)[tensor_storage_index];
     129             :             }
     130             :             return mean_accumulate / neighbor_ids.size();
     131             :           }();
     132             :       result.insert(std::make_pair(dir, effective_neighbor_mean));
     133             :     }
     134             :   }
     135             :   return result;
     136             : }
     137             : 
     138             : // Compute an effective element-center-to-neighbor-center distance that accounts
     139             : // for the possibility of different refinement levels or discontinuous maps
     140             : // (e.g., at Block boundaries). Treated naively, these domain features can make
     141             : // a smooth solution appear to be non-smooth in the logical coordinates, which
     142             : // could potentially lead to the limiter triggering erroneously. This effective
     143             : // distance is used to scale the difference in the means, so that a linear
     144             : // function at a refinement or Block boundary will still appear smooth to the
     145             : // limiter. The factor is normalized to be 1.0 on a uniform grid.
     146             : //
     147             : // Note that this is not "by the book" Minmod, but an attempt to
     148             : // generalize Minmod to work on non-uniform grids.
     149             : template <size_t VolumeDim>
     150             : double effective_difference_to_neighbor(
     151             :     double u_mean, const Element<VolumeDim>& element,
     152             :     const std::array<double, VolumeDim>& element_size, size_t dim,
     153             :     const Side& side,
     154             :     const DirectionMap<VolumeDim, double>& effective_neighbor_means,
     155             :     const DirectionMap<VolumeDim, double>& effective_neighbor_sizes);
     156             : 
     157             : }  // namespace Limiters::Minmod_detail