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 "IO/Observer/Tags.hpp"
29 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
30 #include "Parallel/Actions/TerminatePhase.hpp"
32 #include "Parallel/Info.hpp"
33 #include "Parallel/InitializationFunctions.hpp"
34 #include "Parallel/Invoke.hpp"
35 #include "Parallel/Main.hpp"
36 #include "Parallel/ParallelComponentHelpers.hpp"
37 #include "Parallel/PhaseDependentActionList.hpp"
38 #include "Parallel/Reduction.hpp"
39 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
40 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
42 #include "Utilities/Blas.hpp"
43 #include "Utilities/Gsl.hpp"
44 #include "Utilities/Requires.hpp"
45 #include "Utilities/TMPL.hpp"
46 
47 // IWYU pragma: no_forward_declare db::DataBox
48 
50 
51 namespace OptionTags {
53  static constexpr OptionString help =
54  "The number of elements to distribute work on.";
55  using type = size_t;
56 };
57 } // namespace OptionTags
58 
59 /// [array_allocation_tag]
60 namespace Initialization {
61 namespace Tags {
63  using type = int;
64  using option_tags = tmpl::list<OptionTags::NumberOfElements>;
65 
66  static constexpr bool pass_metavariables = false;
67  static int create_from_options(const size_t number_of_elements) noexcept {
68  return number_of_elements;
69  }
70 };
71 } // namespace Tags
72 } // namespace Initialization
73 /// [array_allocation_tag]
74 
75 namespace OptionTags {
76 // This option expects a list of N matrices that each have N*M rows and M
77 // columns, where N is the `NumberOfElements` and M is a nonzero integer.
78 // Therefore, this option specifies a (N*M,N*M) matrix that has its columns
79 // split over all elements. In a context where the linear operator represents a
80 // DG discretization, M is the number of collocation points per element.
82  static constexpr OptionString help = "The linear operator A to invert.";
84 };
85 // Both of the following options expect a list of N vectors that have a size of
86 // M each, so that they constitute a vector of total size N*M (see above).
87 struct Source {
88  static constexpr OptionString help = "The source b in the equation Ax=b.";
90 };
92  static constexpr OptionString help = "The solution x in the equation Ax=b";
94 };
95 } // namespace OptionTags
96 
99  using option_tags = tmpl::list<OptionTags::LinearOperator>;
100 
101  static constexpr bool pass_metavariables = false;
105  linear_operator) noexcept {
106  return linear_operator;
107  }
108 };
109 
112  using option_tags = tmpl::list<OptionTags::Source>;
113 
114  static constexpr bool pass_metavariables = false;
116  const std::vector<DenseVector<double>>& source) noexcept {
117  return source;
118  }
119 };
120 
123  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
124 
125  static constexpr bool pass_metavariables = false;
127  const std::vector<DenseVector<double>>& expected_result) noexcept {
128  return expected_result;
129  }
130 };
131 
132 // The vector `x` we want to solve for
134  using type = Scalar<DataVector>;
135  static std::string name() noexcept { return "ScalarField"; }
136 };
137 
140 using operator_applied_to_fields_tag =
143 using operator_applied_to_operand_tag =
145 
146 // In the following `ComputeOperatorAction` and `CollectOperatorAction` actions
147 // we compute A(p)=sum_elements(A_element(p_element)) in a global reduction and
148 // then broadcast the global A(p) back to the elements so that they can extract
149 // their A_element(p). This is horribly inefficient parallelism but allows us to
150 // just provide a global matrix A (represented by the `LinearOperator` tag) in
151 // an input file.
152 
153 // Forward declare to keep these actions in the order they are used
154 struct CollectOperatorAction;
155 
157  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
158  typename ActionList, typename ParallelComponent>
161  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
163  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
164  const int array_index,
165  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
166  const ActionList /*meta*/,
167  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
168  const ParallelComponent* const /*meta*/) noexcept {
169  const auto& operator_matrices = get<LinearOperator>(cache);
170  const auto number_of_elements = operator_matrices.size();
171  const auto& linear_operator = gsl::at(operator_matrices, array_index);
172  const auto number_of_grid_points = linear_operator.columns();
173  const auto& operand = get<operand_tag>(box);
174 
175  db::item_type<fields_tag> operator_applied_to_operand{
176  number_of_grid_points * number_of_elements};
177  dgemv_('N', linear_operator.rows(), linear_operator.columns(), 1,
178  linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
179  0, operator_applied_to_operand.data(), 1);
180 
181  Parallel::contribute_to_reduction<CollectOperatorAction>(
182  Parallel::ReductionData<
184  operator_applied_to_operand},
185  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
186  Parallel::get_parallel_component<ParallelComponent>(cache));
187 
188  // Terminate algorithm for now. The reduction will be broadcast to the
189  // next action which is responsible for restarting the algorithm.
190  return {std::move(box), true};
191  }
192 };
193 
195  template <
196  typename ParallelComponent, typename DbTagsList, typename Metavariables,
198  static void apply(db::DataBox<DbTagsList>& box,
200  const int array_index,
201  const db::item_type<fields_tag>& Ap_global_data) noexcept {
202  // This could be generalized to work on the Variables instead of the
203  // Scalar, but it's only for the purpose of this test.
204  const auto number_of_grid_points = get<LinearOperator>(cache)[0].columns();
205  const auto& Ap_global = get<ScalarFieldTag>(Ap_global_data).get();
206  DataVector Ap_local{number_of_grid_points};
207  std::copy(Ap_global.begin() +
208  array_index * static_cast<int>(number_of_grid_points),
209  Ap_global.begin() +
210  (array_index + 1) * static_cast<int>(number_of_grid_points),
211  Ap_local.begin());
214  make_not_null(&box),
215  [&Ap_local](auto Ap) noexcept { *Ap = Scalar<DataVector>(Ap_local); });
216  // Proceed with algorithm
217  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
218  .perform_algorithm(true);
219  }
220 };
221 
222 // Checks for the correct solution after the algorithm has terminated.
223 struct TestResult {
224  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
225  typename ActionList, typename ParallelComponent>
228  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
230  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
231  const int array_index,
232  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
233  const ActionList /*meta*/,
234  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
235  const ParallelComponent* const /*meta*/) noexcept {
236  const auto& has_converged = get<LinearSolver::Tags::HasConverged>(box);
237  SPECTRE_PARALLEL_REQUIRE(has_converged);
238  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
239  Convergence::Reason::AbsoluteResidual);
240  const auto& expected_result =
241  gsl::at(get<ExpectedResult>(cache), array_index);
242  const auto& result = get<ScalarFieldTag>(box).get();
243  for (size_t i = 0; i < expected_result.size(); i++) {
244  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
245  }
246  return {std::move(box), true};
247  }
248 };
249 
251  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
252  typename ActionList, typename ParallelComponent>
253  static auto apply(db::DataBox<DbTagsList>& box,
254  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
256  const int array_index, const ActionList /*meta*/,
257  const ParallelComponent* const /*meta*/) noexcept {
258  const auto& source = gsl::at(get<Source>(cache), array_index);
259  const size_t num_points = source.size();
260 
261  return std::make_tuple(
264  db::AddSimpleTags<fields_tag, sources_tag,
265  operator_applied_to_fields_tag, operand_tag,
266  operator_applied_to_operand_tag>>(
267  std::move(box), db::item_type<fields_tag>{num_points, 0.},
274  }
275 };
276 
277 template <typename Metavariables>
278 struct ElementArray {
279  using chare_type = Parallel::Algorithms::Array;
280  using metavariables = Metavariables;
281  using phase_dependent_action_list = tmpl::list<
283  typename Metavariables::Phase, Metavariables::Phase::Initialization,
284  tmpl::list<InitializeElement,
285  typename Metavariables::linear_solver::initialize_element,
287 
289  typename Metavariables::Phase,
290  Metavariables::Phase::PerformLinearSolve,
292  typename Metavariables::linear_solver::prepare_step,
294  typename Metavariables::linear_solver::perform_step>>,
295 
296  Parallel::PhaseActions<typename Metavariables::Phase,
297  Metavariables::Phase::TestResult,
298  tmpl::list<TestResult>>>;
299  using array_allocation_tags =
300  tmpl::list<Initialization::Tags::NumberOfElements>;
301  using initialization_tags = Parallel::get_initialization_tags<
303  array_allocation_tags>;
304  using const_global_cache_tags =
305  tmpl::list<LinearOperator, Source, ExpectedResult>;
306  using array_index = int;
307 
308  static void allocate_array(
309  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache,
310  const tuples::tagged_tuple_from_typelist<initialization_tags>&
311  initialization_items) noexcept {
312  auto& array_proxy = Parallel::get_parallel_component<ElementArray>(
313  *(global_cache.ckLocalBranch()));
314  for (int i = 0, which_proc = 0,
316  i < get<Initialization::Tags::NumberOfElements>(initialization_items);
317  i++) {
318  array_proxy[i].insert(global_cache, initialization_items, which_proc);
319  which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
320  }
321  array_proxy.doneInserting();
322  }
323 
324  static void execute_next_phase(
325  const typename Metavariables::Phase next_phase,
326  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
327  auto& local_component = Parallel::get_parallel_component<ElementArray>(
328  *(global_cache.ckLocalBranch()));
329  local_component.start_phase(next_phase);
330  }
331 };
332 
333 } // namespace DistributedLinearSolverAlgorithmTestHelpers
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Wrap Tag in Prefix<_, Args...>, also wrapping variables tags if Tag is a Tags::Variables.
Definition: PrefixHelpers.hpp:145
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:52
Definition: VariablesTag.hpp:21
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
Defines functions for interfacing with the parallelization framework.
Definition: ConservativeSystem.hpp:50
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:62
void mutate(const gsl::not_null< DataBox< TagList > *> box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1101
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:49
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:87
T signaling_NaN(T... args)
int number_of_procs()
Number of processing elements.
Definition: Info.hpp:16
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:91
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:121
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:156
Defines the type alias Requires.
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
auto merge_into_databox(db::DataBox< DbTagsList > &&box, Args &&... args) noexcept
Add tags that are not yet in the DataBox.
Definition: MergeIntoDataBox.hpp:133
The operand that the local linear operator is applied to.
Definition: Tags.hpp:53
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:250
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
Functional for computing + of two objects.
Definition: Functional.hpp:238
Defines class DenseMatrix.
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:32
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1231
Definition: Strahlkorper.hpp:167
Definition: InterpolationTargetWedgeSectionTorus.hpp:25
Defines class Variables.
Definition: DataBoxTag.hpp:27
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:136
Declares the interfaces for the BLAS used.
Defines classes for Tensor.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1633
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:110
Functions to enable/disable termination on floating point exceptions.
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:194
#define SPECTRE_PARALLEL_REQUIRE(expr)
A similar to Catch&#39;s REQUIRE statement, but can be used in tests that spawn several chares with possi...
Definition: TestingFramework.hpp:65
tuples::TaggedTuple< Tags... > create_from_options(const tuples::TaggedTuple< OptionTags... > &options, tmpl::list< Tags... >) noexcept
Given a list of tags and a tagged tuple containing items created from input options, return a tagged tuple of items constructed by calls to create_from_options for each tag in the list.
Definition: CreateFromOptions.hpp:37
Defines class DenseVector.
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
Stores a collection of function values.
Definition: DataVector.hpp:42
Wraps the template metaprogramming library used (brigand)
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:81
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:133
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
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
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
Code to wrap or improve the Catch testing framework used for unit tests.
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
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:278
The linear operator applied to the data in Tag
Definition: Tags.hpp:66
Definition: Test_ActionTesting.cpp:365
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Defines class template ConstGlobalCache.
Defines macro ERROR.
Defines the Charm++ mainchare.
Defines DataBox tags for the linear solver.
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:223
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
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:158
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:97
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:64
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:124