DistributedLinearSolverAlgorithmTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <algorithm>
9 #include <array>
10 #include <cstddef>
11 #include <string>
12 #include <tuple>
13 #include <vector>
14 
15 #include "AlgorithmArray.hpp"
18 #include "DataStructures/DataBox/PrefixHelpers.hpp"
20 #include "DataStructures/DataBox/Tag.hpp"
21 #include "DataStructures/DataVector.hpp"
26 #include "ErrorHandling/Error.hpp"
28 #include "Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp"
29 #include "IO/Observer/Actions.hpp"
30 #include "IO/Observer/ObserverComponent.hpp"
31 #include "IO/Observer/RegisterObservers.hpp"
32 #include "IO/Observer/Tags.hpp"
33 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
34 #include "Parallel/Actions/TerminatePhase.hpp"
36 #include "Parallel/Info.hpp"
37 #include "Parallel/InitializationFunctions.hpp"
38 #include "Parallel/Invoke.hpp"
39 #include "Parallel/Main.hpp"
40 #include "Parallel/ParallelComponentHelpers.hpp"
41 #include "Parallel/PhaseDependentActionList.hpp"
42 #include "Parallel/Reduction.hpp"
43 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
44 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
46 #include "Utilities/Blas.hpp"
47 #include "Utilities/Gsl.hpp"
48 #include "Utilities/Requires.hpp"
49 #include "Utilities/TMPL.hpp"
50 
51 namespace helpers = LinearSolverAlgorithmTestHelpers;
52 
53 namespace DistributedLinearSolverAlgorithmTestHelpers {
54 
55 namespace OptionTags {
57  static constexpr OptionString help =
58  "The number of elements to distribute work on.";
59  using type = size_t;
60 };
61 } // namespace OptionTags
62 
63 /// [array_allocation_tag]
64 namespace Initialization {
65 namespace Tags {
67  using type = int;
68  using option_tags = tmpl::list<OptionTags::NumberOfElements>;
69 
70  static constexpr bool pass_metavariables = false;
71  static int create_from_options(const size_t number_of_elements) noexcept {
72  return number_of_elements;
73  }
74 };
75 } // namespace Tags
76 } // namespace Initialization
77 /// [array_allocation_tag]
78 
79 namespace OptionTags {
80 // This option expects a list of N matrices that each have N*M rows and M
81 // columns, where N is the `NumberOfElements` and M is a nonzero integer.
82 // Therefore, this option specifies a (N*M,N*M) matrix that has its columns
83 // split over all elements. In a context where the linear operator represents a
84 // DG discretization, M is the number of collocation points per element.
86  static constexpr OptionString help = "The linear operator A to invert.";
88 };
89 // Both of the following options expect a list of N vectors that have a size of
90 // M each, so that they constitute a vector of total size N*M (see above).
91 struct Source {
92  static constexpr OptionString help = "The source b in the equation Ax=b.";
94 };
96  static constexpr OptionString help = "The solution x in the equation Ax=b";
98 };
99 } // namespace OptionTags
100 
103  using option_tags = tmpl::list<OptionTags::LinearOperator>;
104 
105  static constexpr bool pass_metavariables = false;
107  create_from_options(
109  linear_operator) noexcept {
110  return linear_operator;
111  }
112 };
113 
116  using option_tags = tmpl::list<OptionTags::Source>;
117 
118  static constexpr bool pass_metavariables = false;
119  static std::vector<DenseVector<double>> create_from_options(
120  const std::vector<DenseVector<double>>& source) noexcept {
121  return source;
122  }
123 };
124 
127  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
128 
129  static constexpr bool pass_metavariables = false;
130  static std::vector<DenseVector<double>> create_from_options(
131  const std::vector<DenseVector<double>>& expected_result) noexcept {
132  return expected_result;
133  }
134 };
135 
136 // The vector `x` we want to solve for
138  using type = Scalar<DataVector>;
139  static std::string name() noexcept { return "ScalarField"; }
140 };
141 
144 using operator_applied_to_fields_tag =
146 
147 // In the following `ComputeOperatorAction` and `CollectOperatorAction` actions
148 // we compute A(p)=sum_elements(A_element(p_element)) in a global reduction and
149 // then broadcast the global A(p) back to the elements so that they can extract
150 // their A_element(p). This is horribly inefficient parallelism but allows us to
151 // just provide a global matrix A (represented by the `LinearOperator` tag) in
152 // an input file.
153 
154 // Forward declare to keep these actions in the order they are used
155 template <typename OperandTag>
157 
158 template <typename OperandTag>
160  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
161  typename ActionList, typename ParallelComponent>
162  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
164  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
166  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
167  const int array_index,
168  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
169  const ActionList /*meta*/,
170  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
171  const ParallelComponent* const /*meta*/) noexcept {
172  const auto& operator_matrices = get<LinearOperator>(box);
173  const auto number_of_elements = operator_matrices.size();
174  const auto& linear_operator = gsl::at(operator_matrices, array_index);
175  const auto number_of_grid_points = linear_operator.columns();
176  const auto& operand = get<OperandTag>(box);
177 
178  db::item_type<OperandTag> operator_applied_to_operand{
179  number_of_grid_points * number_of_elements};
180  dgemv_('N', linear_operator.rows(), linear_operator.columns(), 1,
181  linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
182  0, operator_applied_to_operand.data(), 1);
183 
184  Parallel::contribute_to_reduction<CollectOperatorAction<OperandTag>>(
185  Parallel::ReductionData<
187  operator_applied_to_operand},
188  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
189  Parallel::get_parallel_component<ParallelComponent>(cache));
190 
191  // Terminate algorithm for now. The reduction will be broadcast to the
192  // next action which is responsible for restarting the algorithm.
193  return {std::move(box), true};
194  }
195 };
196 
197 template <typename OperandTag>
198 struct CollectOperatorAction {
199  using local_operator_applied_to_operand_tag =
201 
202  template <typename ParallelComponent, typename DbTagsList,
203  typename Metavariables, typename ScalarFieldOperandTag,
204  Requires<tmpl::list_contains_v<
205  DbTagsList, local_operator_applied_to_operand_tag>> = nullptr>
206  static void apply(db::DataBox<DbTagsList>& box,
208  const int array_index,
209  const Variables<tmpl::list<ScalarFieldOperandTag>>&
210  Ap_global_data) noexcept {
211  // This could be generalized to work on the Variables instead of the
212  // Scalar, but it's only for the purpose of this test.
213  const auto number_of_grid_points = get<LinearOperator>(box)[0].columns();
214  const auto& Ap_global = get<ScalarFieldOperandTag>(Ap_global_data).get();
215  DataVector Ap_local{number_of_grid_points};
216  std::copy(Ap_global.begin() +
217  array_index * static_cast<int>(number_of_grid_points),
218  Ap_global.begin() +
219  (array_index + 1) * static_cast<int>(number_of_grid_points),
220  Ap_local.begin());
221  db::mutate<local_operator_applied_to_operand_tag>(
222  make_not_null(&box),
223  [&Ap_local, &number_of_grid_points](auto Ap) noexcept {
225  number_of_grid_points};
227  *Ap)) = Ap_local;
228  });
229  // Proceed with algorithm
230  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
231  .perform_algorithm(true);
232  }
233 };
234 
235 // Checks for the correct solution after the algorithm has terminated.
236 template <typename OptionsGroup>
237 struct TestResult {
238  using const_global_cache_tags =
239  tmpl::list<ExpectedResult, helpers::ExpectedConvergenceReason>;
240 
241  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
242  typename ActionList, typename ParallelComponent>
243  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
245  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
247  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
248  const int array_index,
249  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
250  const ActionList /*meta*/,
251  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
252  const ParallelComponent* const /*meta*/) noexcept {
253  const auto& has_converged =
254  get<LinearSolver::Tags::HasConverged<OptionsGroup>>(box);
255  SPECTRE_PARALLEL_REQUIRE(has_converged);
257  has_converged.reason() == get<helpers::ExpectedConvergenceReason>(box));
258  const auto& expected_result =
259  gsl::at(get<ExpectedResult>(box), array_index);
260  const auto& result = get<ScalarFieldTag>(box).get();
261  for (size_t i = 0; i < expected_result.size(); i++) {
262  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
263  }
264  return {std::move(box), true};
265  }
266 };
267 
269  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
270  typename ActionList, typename ParallelComponent>
271  static auto apply(db::DataBox<DbTagsList>& box,
272  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
274  const int array_index, const ActionList /*meta*/,
275  const ParallelComponent* const /*meta*/) noexcept {
276  const auto& source = gsl::at(get<Source>(box), array_index);
277  const size_t num_points = source.size();
278 
279  return std::make_tuple(
282  std::move(box), db::item_type<fields_tag>{num_points, 0.},
283  db::item_type<sources_tag>{source}));
284  }
285 };
286 
287 namespace detail {
288 
289 template <typename Preconditioner>
290 struct run_preconditioner_impl {
291  using type =
292  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
293  typename Preconditioner::prepare_solve,
294  ::Actions::RepeatUntil<
296  typename Preconditioner::options_group>,
297  tmpl::list<typename Preconditioner::prepare_step,
299  typename Preconditioner::operand_tag>,
300  typename Preconditioner::perform_step>>>;
301 };
302 
303 template <>
304 struct run_preconditioner_impl<void> {
305  using type = tmpl::list<>;
306 };
307 
308 template <typename Preconditioner>
309 using run_preconditioner =
310  typename run_preconditioner_impl<Preconditioner>::type;
311 
312 } // namespace detail
313 
314 template <typename Metavariables>
315 struct ElementArray {
316  using chare_type = Parallel::Algorithms::Array;
317  using metavariables = Metavariables;
318  using linear_solver = typename Metavariables::linear_solver;
319  using preconditioner = typename Metavariables::preconditioner;
320  using phase_dependent_action_list = tmpl::list<
322  typename Metavariables::Phase, Metavariables::Phase::Initialization,
323  tmpl::list<InitializeElement,
324  typename linear_solver::initialize_element,
326  helpers::detail::init_preconditioner<preconditioner>,
329  typename Metavariables::Phase,
330  Metavariables::Phase::RegisterWithObserver,
331  tmpl::list<typename linear_solver::register_element,
332  typename linear_solver::prepare_solve,
335  typename Metavariables::Phase,
336  Metavariables::Phase::PerformLinearSolve,
338  typename linear_solver::options_group>,
339  typename linear_solver::prepare_step,
340  detail::run_preconditioner<preconditioner>,
342  typename linear_solver::perform_step>>,
343 
345  typename Metavariables::Phase, Metavariables::Phase::TestResult,
346  tmpl::list<TestResult<typename linear_solver::options_group>>>>;
347  using array_allocation_tags =
348  tmpl::list<Initialization::Tags::NumberOfElements>;
349  using initialization_tags = Parallel::get_initialization_tags<
351  array_allocation_tags>;
352  using const_global_cache_tags =
353  tmpl::list<LinearOperator, Source, ExpectedResult>;
354  using array_index = int;
355 
356  static void allocate_array(
357  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache,
358  const tuples::tagged_tuple_from_typelist<initialization_tags>&
359  initialization_items) noexcept {
360  auto& array_proxy = Parallel::get_parallel_component<ElementArray>(
361  *(global_cache.ckLocalBranch()));
362  for (int i = 0, which_proc = 0,
364  i < get<Initialization::Tags::NumberOfElements>(initialization_items);
365  i++) {
366  array_proxy[i].insert(global_cache, initialization_items, which_proc);
367  which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
368  }
369  array_proxy.doneInserting();
370  }
371 
372  static void execute_next_phase(
373  const typename Metavariables::Phase next_phase,
374  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
375  auto& local_component = Parallel::get_parallel_component<ElementArray>(
376  *(global_cache.ckLocalBranch()));
377  local_component.start_phase(next_phase);
378  }
379 };
380 
381 template <typename Metavariables>
382 using component_list =
383  tmpl::push_back<tmpl::append<helpers::detail::get_component_list<
384  typename Metavariables::linear_solver>,
385  helpers::detail::get_component_list<
386  typename Metavariables::preconditioner>>,
391 
392 } // namespace DistributedLinearSolverAlgorithmTestHelpers
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
FloatingPointExceptions.hpp
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
LinearSolver::Actions::TerminateIfConverged
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:33
DistributedLinearSolverAlgorithmTestHelpers::TestResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:237
std::string
DataBoxTag.hpp
Main.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:689
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:70
vector
Error.hpp
TestingFramework.hpp
DenseVector.hpp
DenseMatrix.hpp
Tags::Variables
Definition: VariablesTag.hpp:21
tuple
Parallel::number_of_procs
int number_of_procs()
Number of processing elements.
Definition: Info.hpp:16
DistributedLinearSolverAlgorithmTestHelpers::ExpectedResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:125
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:85
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::NumberOfElements
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:56
Parallel::Actions::TerminatePhase
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
db::SimpleTag
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
SPECTRE_PARALLEL_REQUIRE
#define SPECTRE_PARALLEL_REQUIRE(expr)
A similar to Catch's REQUIRE statement, but can be used in tests that spawn several chares with possi...
Definition: TestingFramework.hpp:65
DenseVector< double >
Info.hpp
db::add_tag_prefix
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Definition: PrefixHelpers.hpp:145
algorithm
Tags.hpp
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:91
Parallel::get_initialization_tags
tmpl::remove_duplicates< tmpl::flatten< tmpl::list< AllocationTagsList, tmpl::transform< InitializationActionsList, detail::get_initialization_tags_from_action< tmpl::_1 > >> >> get_initialization_tags
Given a list of initialization actions, and possibly a list of tags needed for allocation of an array...
Definition: ParallelComponentHelpers.hpp:140
DistributedLinearSolverAlgorithmTestHelpers::ElementArray
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:315
DenseMatrix< double, blaze::columnMajor >
Parallel::get_initialization_actions_list
tmpl::flatten< tmpl::transform< PhaseDepActionList, detail::get_initialization_actions_list< tmpl::_1 > >> get_initialization_actions_list
Given the phase dependent action list, return the list of actions in the Initialization phase (or an ...
Definition: ParallelComponentHelpers.hpp:104
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
DistributedLinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:159
Parallel::PhaseActions
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
observers::ObserverWriter
The nodegroup parallel component that is responsible for writing data to disk.
Definition: ObserverComponent.hpp:48
DataBox.hpp
cstddef
funcl::Plus
Functional for computing + of two objects.
Definition: Functional.hpp:238
DistributedLinearSolverAlgorithmTestHelpers::ScalarFieldTag
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:137
DistributedLinearSolverAlgorithmTestHelpers::Initialization::Tags::NumberOfElements
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:66
array
db::item_type
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:246
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
observers::Observer
The group parallel component that is responsible for reducing data to be observed.
Definition: ObserverComponent.hpp:27
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:341
DistributedLinearSolverAlgorithmTestHelpers::InitializeElement
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:268
Variables.hpp
LinearSolver::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the linear solver has converged,...
Definition: Tags.hpp:243
db::AddSimpleTags
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1150
DistributedLinearSolverAlgorithmTestHelpers::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:101
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
DistributedLinearSolverAlgorithmTestHelpers::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:114
DistributedLinearSolverAlgorithmTestHelpers::CollectOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:156
Tensor.hpp
Requires.hpp
Initialization::merge_into_databox
auto merge_into_databox(db::DataBox< DbTagsList > &&box, Args &&... args) noexcept
Add tags that are not yet in the DataBox.
Definition: MergeIntoDataBox.hpp:133
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
Blas.hpp
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::ExpectedResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:95
Prefixes.hpp
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp
dgemv_
void dgemv_(const char &TRANS, const size_t &M, const size_t &N, const double &ALPHA, const double *A, const size_t &LDA, const double *X, const size_t &INCX, const double &BETA, double *Y, const size_t &INCY)
Perform a matrix-vector multiplication.
Definition: Blas.hpp:161
ConstGlobalCache.hpp
string