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

Generated by: LCOV version 1.14