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

Generated by: LCOV version 1.14