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 <boost/functional/hash.hpp>  // IWYU pragma: keep
       8             : #include <cstddef>
       9             : #include <limits>
      10             : #include <string>
      11             : #include <type_traits>
      12             : #include <unordered_map>
      13             : #include <utility>
      14             : 
      15             : #include "DataStructures/Tensor/Tensor.hpp"
      16             : #include "DataStructures/Variables.hpp"
      17             : #include "Domain/Structure/Direction.hpp"
      18             : #include "Domain/Structure/Element.hpp"
      19             : #include "Domain/Structure/ElementId.hpp"
      20             : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
      21             : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
      22             : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
      23             : #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
      24             : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
      25             : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      26             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      27             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      28             : #include "Utilities/Gsl.hpp"
      29             : 
      30             : namespace Limiters::Weno_detail {
      31             : 
      32             : // Compute the Simple WENO solution for one tensor component
      33             : //
      34             : // This interface is intended for use in limiters that check the troubled-cell
      35             : // indicator independently for each tensor component. These limiters generally
      36             : // need to limit only a subset of the tensor components.
      37             : //
      38             : // When calling `simple_weno_impl`,
      39             : // - `interpolator_buffer` may be empty
      40             : // - `modified_neighbor_solution_buffer` should contain one DataVector for each
      41             : //   neighboring element (i.e. for each entry in `neighbor_data`)
      42             : template <typename Tag, size_t VolumeDim, typename PackagedData>
      43             : void simple_weno_impl(
      44             :     const gsl::not_null<std::unordered_map<
      45             :         DirectionalId<VolumeDim>, intrp::RegularGrid<VolumeDim>,
      46             :         boost::hash<DirectionalId<VolumeDim>>>*>
      47             :         interpolator_buffer,
      48             :     const gsl::not_null<
      49             :         std::unordered_map<DirectionalId<VolumeDim>, DataVector,
      50             :                            boost::hash<DirectionalId<VolumeDim>>>*>
      51             :         modified_neighbor_solution_buffer,
      52             :     const gsl::not_null<typename Tag::type*> tensor,
      53             :     const double neighbor_linear_weight, const size_t tensor_storage_index,
      54             :     const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
      55             :     const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
      56             :                              boost::hash<DirectionalId<VolumeDim>>>&
      57             :         neighbor_data) {
      58             :   // Check that basis is LGL or LG
      59             :   // Note that the SimpleWeno implementation should generalize well to other
      60             :   // bases beyond LGL or LG, once the oscillation_indicator function has been
      61             :   // generalized.
      62             :   ASSERT(mesh.basis() == make_array<VolumeDim>(Spectral::Basis::Legendre),
      63             :          "Unsupported basis: " << mesh);
      64             :   ASSERT(mesh.quadrature() ==
      65             :                  make_array<VolumeDim>(Spectral::Quadrature::GaussLobatto) or
      66             :              mesh.quadrature() ==
      67             :                  make_array<VolumeDim>(Spectral::Quadrature::Gauss),
      68             :          "Unsupported quadrature: " << mesh);
      69             : 
      70             :   ASSERT(
      71             :       modified_neighbor_solution_buffer->size() == neighbor_data.size(),
      72             :       "modified_neighbor_solution_buffer->size() = "
      73             :           << modified_neighbor_solution_buffer->size()
      74             :           << "\nneighbor_data.size() = " << neighbor_data.size()
      75             :           << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
      76             :              "before calling simple_weno_impl.");
      77             : 
      78             :   // Compute the modified neighbor solutions.
      79             :   // First extrapolate neighbor data onto local grid points, then shift the
      80             :   // extrapolated data so its mean matches the local mean.
      81             :   DataVector& component_to_limit = (*tensor)[tensor_storage_index];
      82             :   const double local_mean = mean_value(component_to_limit, mesh);
      83             :   for (const auto& neighbor_and_data : neighbor_data) {
      84             :     const auto& neighbor = neighbor_and_data.first;
      85             :     const auto& data = neighbor_and_data.second;
      86             : 
      87             :     if (interpolator_buffer->find(neighbor) == interpolator_buffer->end()) {
      88             :       // No interpolator found => create one
      89             :       const auto& direction = neighbor.direction;
      90             :       const auto& source_mesh = data.mesh;
      91             :       const auto target_1d_logical_coords =
      92             :           Weno_detail::local_grid_points_in_neighbor_logical_coords(
      93             :               mesh, source_mesh, element, direction);
      94             :       interpolator_buffer->insert(std::make_pair(
      95             :           neighbor, intrp::RegularGrid<VolumeDim>(source_mesh, mesh,
      96             :                                                   target_1d_logical_coords)));
      97             :     }
      98             : 
      99             :     // Avoid allocations by working directly in the preallocated buffer
     100             :     DataVector& buffer = modified_neighbor_solution_buffer->at(neighbor);
     101             : 
     102             :     interpolator_buffer->at(neighbor).interpolate(
     103             :         make_not_null(&buffer),
     104             :         get<Tag>(data.volume_data)[tensor_storage_index]);
     105             :     const double neighbor_mean = mean_value(buffer, mesh);
     106             :     buffer += (local_mean - neighbor_mean);
     107             :   }
     108             : 
     109             :   // Sum local and modified neighbor polynomials for the WENO reconstruction
     110             :   Weno_detail::reconstruct_from_weighted_sum(
     111             :       make_not_null(&component_to_limit), neighbor_linear_weight,
     112             :       Weno_detail::DerivativeWeight::PowTwoEll, mesh,
     113             :       *modified_neighbor_solution_buffer);
     114             : }
     115             : 
     116             : // Compute the Simple WENO solution for one tensor
     117             : //
     118             : // This interface is intended for use in limiters that check the troubled-cell
     119             : // indicator for the whole cell, and apply the limiter to all fields.
     120             : //
     121             : // When calling `simple_weno_impl`,
     122             : // - `interpolator_buffer` may be empty
     123             : // - `modified_neighbor_solution_buffer` should contain one DataVector for each
     124             : //   neighboring element (i.e. for each entry in `neighbor_data`)
     125             : template <typename Tag, size_t VolumeDim, typename PackagedData>
     126             : void simple_weno_impl(
     127             :     const gsl::not_null<std::unordered_map<
     128             :         DirectionalId<VolumeDim>, intrp::RegularGrid<VolumeDim>,
     129             :         boost::hash<DirectionalId<VolumeDim>>>*>
     130             :         interpolator_buffer,
     131             :     const gsl::not_null<
     132             :         std::unordered_map<DirectionalId<VolumeDim>, DataVector,
     133             :                            boost::hash<DirectionalId<VolumeDim>>>*>
     134             :         modified_neighbor_solution_buffer,
     135             :     const gsl::not_null<typename Tag::type*> tensor,
     136             :     const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
     137             :     const Element<VolumeDim>& element,
     138             :     const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
     139             :                              boost::hash<DirectionalId<VolumeDim>>>&
     140             :         neighbor_data) {
     141             :   for (size_t tensor_storage_index = 0; tensor_storage_index < tensor->size();
     142             :        ++tensor_storage_index) {
     143             :     simple_weno_impl<Tag>(interpolator_buffer,
     144             :                           modified_neighbor_solution_buffer, tensor,
     145             :                           neighbor_linear_weight, tensor_storage_index, mesh,
     146             :                           element, neighbor_data);
     147             :   }
     148             : }
     149             : 
     150             : }  // namespace Limiters::Weno_detail