SpECTRE Documentation Coverage Report
Current view: top level - ParallelAlgorithms/LinearSolver/Schwarz - OverlapHelpers.hpp Hit Total Coverage
Commit: edb1b5199a4a86c269aedbb831767801169f3e8a Lines: 11 28 39.3 %
Date: 2021-04-19 16:23:01
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 <boost/functional/hash.hpp>
       7             : #include <cstddef>
       8             : #include <limits>
       9             : #include <tuple>
      10             : 
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/FixedHashMap.hpp"
      13             : #include "DataStructures/Index.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "DataStructures/Variables.hpp"
      16             : #include "Domain/Structure/Direction.hpp"
      17             : #include "Domain/Structure/ElementId.hpp"
      18             : #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
      19             : #include "Utilities/ContainerHelpers.hpp"
      20             : #include "Utilities/ErrorHandling/Assert.hpp"
      21             : 
      22             : namespace LinearSolver::Schwarz {
      23             : 
      24             : /// Identifies a subdomain region that overlaps with another element
      25             : template <size_t Dim>
      26           1 : using OverlapId = std::pair<Direction<Dim>, ElementId<Dim>>;
      27             : 
      28             : /// Data structure that can store the `ValueType` on each possible overlap of an
      29             : /// element-centered subdomain with its neighbors. Overlaps are identified by
      30             : /// their `OverlapId`.
      31             : template <size_t Dim, typename ValueType>
      32           1 : using OverlapMap =
      33             :     FixedHashMap<maximum_number_of_neighbors(Dim), OverlapId<Dim>, ValueType,
      34             :                  boost::hash<OverlapId<Dim>>>;
      35             : 
      36             : /*!
      37             :  * \brief The number of points that an overlap extends into the `volume_extent`
      38             :  *
      39             :  * In a dimension where an element has `volume_extent` points, the overlap
      40             :  * extent is the largest number under these constraints:
      41             :  *
      42             :  * - It is at most `max_overlap`.
      43             :  * - It is smaller than the `volume_extent`.
      44             :  *
      45             :  * This means the overlap extent is always smaller than the `volume_extent`. The
      46             :  * reason for this constraint is that we define the _width_ of the overlap as
      47             :  * the element-logical coordinate distance from the face of the element to the
      48             :  * first collocation point _outside_ the overlap extent. Therefore, even an
      49             :  * overlap region that covers the full element in width does not include the
      50             :  * collocation point on the opposite side of the element.
      51             :  *
      52             :  * Here's a few notes on the definition of the overlap extent and width:
      53             :  *
      54             :  * - A typical smooth weighting function goes to zero at the overlap width, so
      55             :  *   if the grid points located at the overlap width were included in the
      56             :  *   subdomain, their solutions would not contribute to the weighted sum of
      57             :  *   subdomain solutions.
      58             :  * - Defining the overlap width as the distance to the first point _outside_ the
      59             :  *   overlap extent makes it non-zero even for a single point of overlap into a
      60             :  *   Gauss-Lobatto grid (which has points located at the element face).
      61             :  * - Boundary contributions for many (but not all) discontinuous Galerkin
      62             :  *   schemes on Gauss-Lobatto grids are limited to the grid points on the
      63             :  *   element face, e.g. for a DG operator that is pre-multiplied by the mass
      64             :  *   matrix, or one where boundary contributions are lifted using the diagonal
      65             :  *   mass-matrix approximation. Not including the grid points facing away from
      66             :  *   the subdomain in the overlap allows to ignore that face altogether in the
      67             :  *   subdomain operator.
      68             :  */
      69           1 : size_t overlap_extent(size_t volume_extent, size_t max_overlap) noexcept;
      70             : 
      71             : /*!
      72             :  * \brief Total number of grid points in an overlap region that extends
      73             :  * `overlap_extent` points into the `volume_extents` from either side in the
      74             :  * `overlap_dimension`
      75             :  *
      76             :  * The overlap region has `overlap_extent` points in the `overlap_dimension`,
      77             :  * and `volume_extents` points in the other dimensions. The number of grid
      78             :  * points returned by this function is the product of these extents.
      79             :  */
      80             : template <size_t Dim>
      81           1 : size_t overlap_num_points(const Index<Dim>& volume_extents,
      82             :                           size_t overlap_extent,
      83             :                           size_t overlap_dimension) noexcept;
      84             : 
      85             : /*!
      86             :  * \brief Width of an overlap extending `overlap_extent` points into the
      87             :  * `collocation_points` from either side.
      88             :  *
      89             :  * The "width" of an overlap is the element-logical coordinate distance from the
      90             :  * element boundary to the first collocation point outside the overlap region in
      91             :  * the overlap dimension, i.e. the dimension perpendicular to the element face.
      92             :  * See `LinearSolver::Schwarz::overlap_extent` for details.
      93             :  *
      94             :  * This function assumes the `collocation_points` are mirrored around 0.
      95             :  */
      96           1 : double overlap_width(size_t overlap_extent,
      97             :                      const DataVector& collocation_points) noexcept;
      98             : 
      99             : /*!
     100             :  * \brief Iterate over grid points in a region that extends partially into the
     101             :  * volume
     102             :  *
     103             :  * Here's an example how to use this iterator:
     104             :  *
     105             :  * \snippet Test_OverlapHelpers.cpp overlap_iterator
     106             :  */
     107           1 : class OverlapIterator {
     108             :  public:
     109             :   template <size_t Dim>
     110           0 :   OverlapIterator(const Index<Dim>& volume_extents, size_t overlap_extent,
     111             :                   const Direction<Dim>& direction) noexcept;
     112             : 
     113           0 :   explicit operator bool() const noexcept;
     114             : 
     115           0 :   OverlapIterator& operator++();
     116             : 
     117             :   /// Offset into a DataVector that holds full volume data
     118           1 :   size_t volume_offset() const noexcept;
     119             : 
     120             :   /// Offset into a DataVector that holds data only on the overlap region
     121           1 :   size_t overlap_offset() const noexcept;
     122             : 
     123           0 :   void reset() noexcept;
     124             : 
     125             :  private:
     126           0 :   size_t size_ = std::numeric_limits<size_t>::max();
     127           0 :   size_t num_slices_ = std::numeric_limits<size_t>::max();
     128           0 :   size_t stride_ = std::numeric_limits<size_t>::max();
     129           0 :   size_t stride_count_ = std::numeric_limits<size_t>::max();
     130           0 :   size_t jump_ = std::numeric_limits<size_t>::max();
     131           0 :   size_t initial_offset_ = std::numeric_limits<size_t>::max();
     132           0 :   size_t volume_offset_ = std::numeric_limits<size_t>::max();
     133           0 :   size_t overlap_offset_ = std::numeric_limits<size_t>::max();
     134             : };
     135             : 
     136             : // @{
     137             : /// The part of the tensor data that lies within the overlap region
     138             : template <size_t Dim, typename DataType, typename... TensorStructure>
     139           1 : void data_on_overlap(const gsl::not_null<Tensor<DataType, TensorStructure...>*>
     140             :                          restricted_tensor,
     141             :                      const Tensor<DataType, TensorStructure...>& tensor,
     142             :                      const Index<Dim>& volume_extents,
     143             :                      const size_t overlap_extent,
     144             :                      const Direction<Dim>& direction) noexcept {
     145             :   for (OverlapIterator overlap_iterator{volume_extents, overlap_extent,
     146             :                                         direction};
     147             :        overlap_iterator; ++overlap_iterator) {
     148             :     for (size_t tensor_component = 0; tensor_component < tensor.size();
     149             :          ++tensor_component) {
     150             :       (*restricted_tensor)[tensor_component][overlap_iterator
     151             :                                                  .overlap_offset()] =
     152             :           tensor[tensor_component][overlap_iterator.volume_offset()];
     153             :     }
     154             :   }
     155             : }
     156             : 
     157             : template <size_t Dim, typename DataType, typename... TensorStructure>
     158           0 : Tensor<DataType, TensorStructure...> data_on_overlap(
     159             :     const Tensor<DataType, TensorStructure...>& tensor,
     160             :     const Index<Dim>& volume_extents, const size_t overlap_extent,
     161             :     const Direction<Dim>& direction) noexcept {
     162             :   Tensor<DataType, TensorStructure...> restricted_tensor{overlap_num_points(
     163             :       volume_extents, overlap_extent, direction.dimension())};
     164             :   data_on_overlap(make_not_null(&restricted_tensor), tensor, volume_extents,
     165             :                   overlap_extent, direction);
     166             :   return restricted_tensor;
     167             : }
     168             : 
     169             : namespace detail {
     170             : template <size_t Dim>
     171             : void data_on_overlap_impl(double* overlap_data, const double* volume_data,
     172             :                           size_t num_components,
     173             :                           const Index<Dim>& volume_extents,
     174             :                           size_t overlap_extent,
     175             :                           const Direction<Dim>& direction) noexcept;
     176             : }  // namespace detail
     177             : 
     178             : template <size_t Dim, typename OverlapTagsList, typename VolumeTagsList>
     179           0 : void data_on_overlap(
     180             :     const gsl::not_null<Variables<OverlapTagsList>*> overlap_data,
     181             :     const Variables<VolumeTagsList>& volume_data,
     182             :     const Index<Dim>& volume_extents, const size_t overlap_extent,
     183             :     const Direction<Dim>& direction) noexcept {
     184             :   constexpr size_t num_components =
     185             :       Variables<VolumeTagsList>::number_of_independent_components;
     186             :   ASSERT(volume_data.number_of_grid_points() == volume_extents.product(),
     187             :          "volume_data has wrong number of grid points.  Expected "
     188             :              << volume_extents.product() << ", got "
     189             :              << volume_data.number_of_grid_points());
     190             :   ASSERT(overlap_data->number_of_grid_points() ==
     191             :              overlap_num_points(volume_extents, overlap_extent,
     192             :                                 direction.dimension()),
     193             :          "overlap_data has wrong number of grid points.  Expected "
     194             :              << overlap_num_points(volume_extents, overlap_extent,
     195             :                                    direction.dimension())
     196             :              << ", got " << overlap_data->number_of_grid_points());
     197             :   detail::data_on_overlap_impl(overlap_data->data(), volume_data.data(),
     198             :                                num_components, volume_extents, overlap_extent,
     199             :                                direction);
     200             : }
     201             : 
     202             : template <size_t Dim, typename TagsList>
     203           0 : Variables<TagsList> data_on_overlap(const Variables<TagsList>& volume_data,
     204             :                                     const Index<Dim>& volume_extents,
     205             :                                     const size_t overlap_extent,
     206             :                                     const Direction<Dim>& direction) noexcept {
     207             :   Variables<TagsList> overlap_data{overlap_num_points(
     208             :       volume_extents, overlap_extent, direction.dimension())};
     209             :   data_on_overlap(make_not_null(&overlap_data), volume_data, volume_extents,
     210             :                   overlap_extent, direction);
     211             :   return overlap_data;
     212             : }
     213             : // @}
     214             : 
     215             : namespace detail {
     216             : template <size_t Dim>
     217             : void add_overlap_data_impl(double* volume_data, const double* overlap_data,
     218             :                            size_t num_components,
     219             :                            const Index<Dim>& volume_extents,
     220             :                            size_t overlap_extent,
     221             :                            const Direction<Dim>& direction) noexcept;
     222             : }  // namespace detail
     223             : 
     224             : /// Add the `overlap_data` to the `volume_data`
     225             : template <size_t Dim, typename VolumeTagsList, typename OverlapTagsList>
     226           1 : void add_overlap_data(
     227             :     const gsl::not_null<Variables<VolumeTagsList>*> volume_data,
     228             :     const Variables<OverlapTagsList>& overlap_data,
     229             :     const Index<Dim>& volume_extents, const size_t overlap_extent,
     230             :     const Direction<Dim>& direction) noexcept {
     231             :   constexpr size_t num_components =
     232             :       Variables<VolumeTagsList>::number_of_independent_components;
     233             :   ASSERT(volume_data->number_of_grid_points() == volume_extents.product(),
     234             :          "volume_data has wrong number of grid points.  Expected "
     235             :              << volume_extents.product() << ", got "
     236             :              << volume_data->number_of_grid_points());
     237             :   ASSERT(overlap_data.number_of_grid_points() ==
     238             :              overlap_num_points(volume_extents, overlap_extent,
     239             :                                 direction.dimension()),
     240             :          "overlap_data has wrong number of grid points.  Expected "
     241             :              << overlap_num_points(volume_extents, overlap_extent,
     242             :                                    direction.dimension())
     243             :              << ", got " << overlap_data.number_of_grid_points());
     244             :   detail::add_overlap_data_impl(volume_data->data(), overlap_data.data(),
     245             :                                 num_components, volume_extents, overlap_extent,
     246             :                                 direction);
     247             : }
     248             : 
     249             : // @{
     250             : /// Extend the overlap data to the full mesh by filling it with zeros outside
     251             : /// the overlap region
     252             : template <size_t Dim, typename ExtendedTagsList, typename OverlapTagsList>
     253           1 : void extended_overlap_data(
     254             :     const gsl::not_null<Variables<ExtendedTagsList>*> extended_data,
     255             :     const Variables<OverlapTagsList>& overlap_data,
     256             :     const Index<Dim>& volume_extents, const size_t overlap_extent,
     257             :     const Direction<Dim>& direction) noexcept {
     258             :   *extended_data = Variables<ExtendedTagsList>{volume_extents.product(), 0.};
     259             :   add_overlap_data(extended_data, overlap_data, volume_extents, overlap_extent,
     260             :                    direction);
     261             : }
     262             : 
     263             : template <size_t Dim, typename TagsList>
     264           0 : Variables<TagsList> extended_overlap_data(
     265             :     const Variables<TagsList>& overlap_data, const Index<Dim>& volume_extents,
     266             :     const size_t overlap_extent, const Direction<Dim>& direction) noexcept {
     267             :   Variables<TagsList> extended_data{volume_extents.product()};
     268             :   extended_overlap_data(make_not_null(&extended_data), overlap_data,
     269             :                         volume_extents, overlap_extent, direction);
     270             :   return extended_data;
     271             : }
     272             : // @}
     273             : 
     274             : }  // namespace LinearSolver::Schwarz

Generated by: LCOV version 1.14