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

Generated by: LCOV version 1.14