SpECTRE Documentation Coverage Report
Current view: top level - ParallelAlgorithms/LinearSolver/Schwarz - ElementCenteredSubdomainData.hpp Hit Total Coverage
Commit: edb1b5199a4a86c269aedbb831767801169f3e8a Lines: 4 62 6.5 %
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 <cstddef>
       7             : #include <ostream>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <tuple>
      11             : 
      12             : #include "DataStructures/Variables.hpp"
      13             : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
      14             : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
      15             : #include "Utilities/MakeWithValue.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : namespace LinearSolver::Schwarz {
      19             : 
      20             : /// \cond
      21             : template <bool Const, size_t Dim, typename TagsList>
      22             : struct ElementCenteredSubdomainDataIterator;
      23             : /// \endcond
      24             : 
      25             : /*!
      26             :  * \brief Data on an element-centered subdomain
      27             :  *
      28             :  * An element-centered subdomain consists of a central element and overlap
      29             :  * regions with all neighboring elements. This class holds data on such a
      30             :  * subdomain. It supports vector space operations (addition and scalar
      31             :  * multiplication) and an inner product, which allows the use of this data type
      32             :  * with linear solvers (see e.g. `LinearSolver::Serial::Gmres`).
      33             :  */
      34             : template <size_t Dim, typename TagsList>
      35           1 : struct ElementCenteredSubdomainData {
      36           0 :   static constexpr size_t volume_dim = Dim;
      37           0 :   using ElementData = Variables<TagsList>;
      38           0 :   using OverlapData = ElementData;
      39           0 :   using iterator = ElementCenteredSubdomainDataIterator<false, Dim, TagsList>;
      40           0 :   using const_iterator =
      41             :       ElementCenteredSubdomainDataIterator<true, Dim, TagsList>;
      42             : 
      43           0 :   ElementCenteredSubdomainData() = default;
      44           0 :   ElementCenteredSubdomainData(const ElementCenteredSubdomainData&) = default;
      45           0 :   ElementCenteredSubdomainData& operator=(
      46             :       const ElementCenteredSubdomainData&) noexcept = default;
      47           0 :   ElementCenteredSubdomainData(ElementCenteredSubdomainData&&) noexcept =
      48             :       default;
      49           0 :   ElementCenteredSubdomainData& operator=(
      50             :       ElementCenteredSubdomainData&&) noexcept = default;
      51           0 :   ~ElementCenteredSubdomainData() noexcept = default;
      52             : 
      53           0 :   explicit ElementCenteredSubdomainData(
      54             :       const size_t element_num_points) noexcept
      55             :       : element_data{element_num_points} {}
      56             : 
      57             :   template <typename UsedForSizeTagsList>
      58           0 :   void destructive_resize(
      59             :       const ElementCenteredSubdomainData<Dim, UsedForSizeTagsList>&
      60             :           used_for_size) noexcept {
      61             :     if (UNLIKELY(element_data.number_of_grid_points() !=
      62             :                  used_for_size.element_data.number_of_grid_points())) {
      63             :       element_data.initialize(
      64             :           used_for_size.element_data.number_of_grid_points());
      65             :     }
      66             :     for (const auto& [overlap_id, used_for_overlap_size] :
      67             :          used_for_size.overlap_data) {
      68             :       if (UNLIKELY(overlap_data[overlap_id].number_of_grid_points() !=
      69             :                    used_for_overlap_size.number_of_grid_points())) {
      70             :         overlap_data[overlap_id].initialize(
      71             :             used_for_overlap_size.number_of_grid_points());
      72             :       }
      73             :     }
      74             :   }
      75             : 
      76           0 :   ElementCenteredSubdomainData(
      77             :       Variables<TagsList> local_element_data,
      78             :       OverlapMap<Dim, Variables<TagsList>> local_overlap_data) noexcept
      79             :       : element_data{std::move(local_element_data)},
      80             :         overlap_data{std::move(local_overlap_data)} {}
      81             : 
      82           0 :   size_t size() const noexcept {
      83             :     return std::accumulate(
      84             :         overlap_data.begin(), overlap_data.end(), element_data.size(),
      85             :         [](const size_t size, const auto& overlap_id_and_data) noexcept {
      86             :           return size + overlap_id_and_data.second.size();
      87             :         });
      88             :   }
      89           0 :   iterator begin() noexcept { return {this}; }
      90           0 :   iterator end() noexcept { return {}; }
      91           0 :   const_iterator begin() const noexcept { return {this}; }
      92           0 :   const_iterator end() const noexcept { return {}; }
      93           0 :   const_iterator cbegin() const noexcept { return begin(); }
      94           0 :   const_iterator cend() const noexcept { return end(); }
      95             : 
      96           0 :   void pup(PUP::er& p) noexcept {  // NOLINT
      97             :     p | element_data;
      98             :     p | overlap_data;
      99             :   }
     100             : 
     101             :   template <typename RhsTagsList>
     102           0 :   ElementCenteredSubdomainData& operator+=(
     103             :       const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
     104             :     element_data += rhs.element_data;
     105             :     for (auto& [overlap_id, data] : overlap_data) {
     106             :       data += rhs.overlap_data.at(overlap_id);
     107             :     }
     108             :     return *this;
     109             :   }
     110             : 
     111             :   template <typename RhsTagsList>
     112           0 :   ElementCenteredSubdomainData& operator-=(
     113             :       const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
     114             :     element_data -= rhs.element_data;
     115             :     for (auto& [overlap_id, data] : overlap_data) {
     116             :       data -= rhs.overlap_data.at(overlap_id);
     117             :     }
     118             :     return *this;
     119             :   }
     120             : 
     121           0 :   ElementCenteredSubdomainData& operator*=(const double scalar) noexcept {
     122             :     element_data *= scalar;
     123             :     for (auto& [overlap_id, data] : overlap_data) {
     124             :       data *= scalar;
     125             :       // Silence unused-variable warning on GCC 7
     126             :       (void)overlap_id;
     127             :     }
     128             :     return *this;
     129             :   }
     130             : 
     131           0 :   ElementCenteredSubdomainData& operator/=(const double scalar) noexcept {
     132             :     element_data /= scalar;
     133             :     for (auto& [overlap_id, data] : overlap_data) {
     134             :       data /= scalar;
     135             :       // Silence unused-variable warning on GCC 7
     136             :       (void)overlap_id;
     137             :     }
     138             :     return *this;
     139             :   }
     140             : 
     141           0 :   ElementData element_data{};
     142           0 :   OverlapMap<Dim, OverlapData> overlap_data{};
     143             : };
     144             : 
     145             : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
     146           0 : decltype(auto) operator+(
     147             :     ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
     148             :     const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
     149             :   lhs += rhs;
     150             :   return lhs;
     151             : }
     152             : 
     153             : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
     154           0 : decltype(auto) operator-(
     155             :     ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
     156             :     const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
     157             :   lhs -= rhs;
     158             :   return lhs;
     159             : }
     160             : 
     161             : template <size_t Dim, typename TagsList>
     162           0 : decltype(auto) operator*(
     163             :     const double scalar,
     164             :     ElementCenteredSubdomainData<Dim, TagsList> data) noexcept {
     165             :   data *= scalar;
     166             :   return data;
     167             : }
     168             : 
     169             : template <size_t Dim, typename TagsList>
     170           0 : decltype(auto) operator*(ElementCenteredSubdomainData<Dim, TagsList> data,
     171             :                          const double scalar) noexcept {
     172             :   data *= scalar;
     173             :   return data;
     174             : }
     175             : 
     176             : template <size_t Dim, typename TagsList>
     177           0 : decltype(auto) operator/(ElementCenteredSubdomainData<Dim, TagsList> data,
     178             :                          const double scalar) noexcept {
     179             :   data /= scalar;
     180             :   return data;
     181             : }
     182             : 
     183             : template <size_t Dim, typename TagsList>
     184           0 : std::ostream& operator<<(std::ostream& os,
     185             :                          const ElementCenteredSubdomainData<Dim, TagsList>&
     186             :                              subdomain_data) noexcept {
     187             :   os << "Element data:\n"
     188             :      << subdomain_data.element_data << "\nOverlap data:\n"
     189             :      << subdomain_data.overlap_data;
     190             :   return os;
     191             : }
     192             : 
     193             : template <size_t Dim, typename TagsList>
     194           0 : bool operator==(
     195             :     const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
     196             :     const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
     197             :   return lhs.element_data == rhs.element_data and
     198             :          lhs.overlap_data == rhs.overlap_data;
     199             : }
     200             : 
     201             : template <size_t Dim, typename TagsList>
     202           0 : bool operator!=(
     203             :     const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
     204             :     const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
     205             :   return not(lhs == rhs);
     206             : }
     207             : 
     208             : namespace detail {
     209             : // Defines a consistent ordering of overlap IDs
     210             : template <size_t Dim, typename ValueType>
     211             : std::vector<OverlapId<Dim>> ordered_overlap_ids(
     212             :     const OverlapMap<Dim, ValueType>& overlap_map) noexcept {
     213             :   std::vector<OverlapId<Dim>> overlap_ids{};
     214             :   overlap_ids.reserve(overlap_map.size());
     215             :   std::transform(overlap_map.begin(), overlap_map.end(),
     216             :                  std::back_inserter(overlap_ids),
     217             :                  [](const auto& overlap_id_and_value) noexcept {
     218             :                    return overlap_id_and_value.first;
     219             :                  });
     220             :   std::sort(overlap_ids.begin(), overlap_ids.end(),
     221             :             [](const OverlapId<Dim>& lhs, const OverlapId<Dim>& rhs) noexcept {
     222             :               if (lhs.first.axis() != rhs.first.axis()) {
     223             :                 return lhs.first.axis() < rhs.first.axis();
     224             :               }
     225             :               if (lhs.first.side() != rhs.first.side()) {
     226             :                 return lhs.first.side() < rhs.first.side();
     227             :               }
     228             :               if (lhs.second.block_id() != rhs.second.block_id()) {
     229             :                 return lhs.second.block_id() < rhs.second.block_id();
     230             :               }
     231             :               for (size_t d = 0; d < Dim; ++d) {
     232             :                 const auto& lhs_segment_id = lhs.second.segment_ids().at(d);
     233             :                 const auto& rhs_segment_id = rhs.second.segment_ids().at(d);
     234             :                 if (lhs_segment_id.refinement_level() !=
     235             :                     rhs_segment_id.refinement_level()) {
     236             :                   return lhs_segment_id.refinement_level() <
     237             :                          rhs_segment_id.refinement_level();
     238             :                 }
     239             :                 if (lhs_segment_id.index() != rhs_segment_id.index()) {
     240             :                   return lhs_segment_id.index() < rhs_segment_id.index();
     241             :                 }
     242             :               }
     243             :               return false;
     244             :             });
     245             :   return overlap_ids;
     246             : }
     247             : }  // namespace detail
     248             : 
     249             : /*!
     250             :  * \brief Iterate over `LinearSolver::Schwarz::ElementCenteredSubdomainData`
     251             :  *
     252             :  * This iterator guarantees that it steps through the data in the same order as
     253             :  * long as these conditions are satisfied:
     254             :  *
     255             :  * - The set of overlap IDs in the `overlap_data` doesn't change
     256             :  * - The extents of the `element_data` and the `overlap_data doesn't change
     257             :  *
     258             :  * Iterating requires sorting the overlap IDs. If you find this impacts
     259             :  * performance, be advised to implement the internal data storage in
     260             :  * `ElementCenteredSubdomainData` so it stores its data contiguously, e.g. by
     261             :  * implementing non-owning variables.
     262             :  */
     263             : template <bool Const, size_t Dim, typename TagsList>
     264           1 : struct ElementCenteredSubdomainDataIterator {
     265             :  private:
     266           0 :   using PtrType =
     267             :       tmpl::conditional_t<Const,
     268             :                           const ElementCenteredSubdomainData<Dim, TagsList>*,
     269             :                           ElementCenteredSubdomainData<Dim, TagsList>*>;
     270             : 
     271             :  public:
     272           0 :   using difference_type = ptrdiff_t;
     273           0 :   using value_type = double;
     274           0 :   using pointer = value_type*;
     275           0 :   using reference = value_type&;
     276           0 :   using iterator_category = std::forward_iterator_tag;
     277             : 
     278             :   /// Construct begin state
     279           1 :   ElementCenteredSubdomainDataIterator(PtrType data) noexcept : data_(data) {
     280             :     overlap_ids_ = detail::ordered_overlap_ids(data_->overlap_data);
     281             :     reset();
     282             :   }
     283             : 
     284           0 :   void reset() noexcept {
     285             :     overlap_index_ = (data_->element_data.size() == 0 and overlap_ids_.empty())
     286             :                          ? std::numeric_limits<size_t>::max()
     287             :                          : 0;
     288             :     data_index_ = 0;
     289             :   }
     290             : 
     291             :   /// Construct end state
     292           1 :   ElementCenteredSubdomainDataIterator() noexcept {
     293             :     overlap_index_ = std::numeric_limits<size_t>::max();
     294             :     data_index_ = 0;
     295             :   }
     296             : 
     297           0 :   ElementCenteredSubdomainDataIterator& operator++() noexcept {
     298             :     ++data_index_;
     299             :     if (data_index_ ==
     300             :         (overlap_index_ == 0
     301             :              ? data_->element_data
     302             :              : data_->overlap_data.at(overlap_ids_[overlap_index_ - 1]))
     303             :             .size()) {
     304             :       ++overlap_index_;
     305             :       data_index_ = 0;
     306             :     }
     307             :     if (overlap_index_ == overlap_ids_.size() + 1) {
     308             :       overlap_index_ = std::numeric_limits<size_t>::max();
     309             :     }
     310             :     return *this;
     311             :   }
     312             : 
     313           0 :   tmpl::conditional_t<Const, double, double&> operator*() const noexcept {
     314             :     if (overlap_index_ == 0) {
     315             :       return data_->element_data.data()[data_index_];
     316             :     } else {
     317             :       return data_->overlap_data.at(overlap_ids_[overlap_index_ - 1])
     318             :           .data()[data_index_];
     319             :     }
     320             :   }
     321             : 
     322             :  private:
     323           0 :   friend bool operator==(
     324             :       const ElementCenteredSubdomainDataIterator& lhs,
     325             :       const ElementCenteredSubdomainDataIterator& rhs) noexcept {
     326             :     return lhs.overlap_index_ == rhs.overlap_index_ and
     327             :            lhs.data_index_ == rhs.data_index_;
     328             :   }
     329             : 
     330           0 :   friend bool operator!=(
     331             :       const ElementCenteredSubdomainDataIterator& lhs,
     332             :       const ElementCenteredSubdomainDataIterator& rhs) noexcept {
     333             :     return not(lhs == rhs);
     334             :   }
     335             : 
     336           0 :   PtrType data_;
     337           0 :   std::vector<OverlapId<Dim>> overlap_ids_;
     338           0 :   size_t overlap_index_;
     339           0 :   size_t data_index_;
     340             : };
     341             : 
     342             : }  // namespace LinearSolver::Schwarz
     343             : 
     344             : namespace LinearSolver::InnerProductImpls {
     345             : 
     346             : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
     347             : struct InnerProductImpl<
     348             :     Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>,
     349           0 :     Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>> {
     350           0 :   static double apply(
     351             :       const Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>& lhs,
     352             :       const Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>&
     353             :           rhs) noexcept {
     354             :     double result = inner_product(lhs.element_data, rhs.element_data);
     355             :     for (const auto& [overlap_id, lhs_data] : lhs.overlap_data) {
     356             :       result += inner_product(lhs_data, rhs.overlap_data.at(overlap_id));
     357             :     }
     358             :     return result;
     359             :   }
     360             : };
     361             : 
     362             : }  // namespace LinearSolver::InnerProductImpls
     363             : 
     364             : namespace MakeWithValueImpls {
     365             : 
     366             : template <size_t Dim, typename TagsListOut, typename TagsListIn>
     367             : struct MakeWithValueImpl<
     368             :     LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>,
     369           0 :     LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>> {
     370           0 :   using SubdomainDataIn =
     371             :       LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>;
     372           0 :   using SubdomainDataOut =
     373             :       LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>;
     374             :   static SPECTRE_ALWAYS_INLINE SubdomainDataOut
     375           0 :   apply(const SubdomainDataIn& input, const double value) noexcept {
     376             :     SubdomainDataOut output{};
     377             :     output.element_data =
     378             :         make_with_value<typename SubdomainDataOut::ElementData>(
     379             :             input.element_data, value);
     380             :     for (const auto& [overlap_id, input_data] : input.overlap_data) {
     381             :       output.overlap_data.emplace(
     382             :           overlap_id, make_with_value<typename SubdomainDataOut::OverlapData>(
     383             :                           input_data, value));
     384             :     }
     385             :     return output;
     386             :   }
     387             : };
     388             : 
     389             : }  // namespace MakeWithValueImpls

Generated by: LCOV version 1.14