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/RegisterWithObservers.hpp"
30 #include "IO/Observer/ObserverComponent.hpp"
31 #include "IO/Observer/Tags.hpp"
32 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
33 #include "Parallel/Actions/TerminatePhase.hpp"
34 #include "Parallel/GlobalCache.hpp"
35 #include "Parallel/Info.hpp"
36 #include "Parallel/InitializationFunctions.hpp"
37 #include "Parallel/Invoke.hpp"
38 #include "Parallel/Main.hpp"
39 #include "Parallel/ParallelComponentHelpers.hpp"
40 #include "Parallel/PhaseDependentActionList.hpp"
41 #include "Parallel/Reduction.hpp"
42 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
43 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
45 #include "Utilities/Blas.hpp"
46 #include "Utilities/Gsl.hpp"
47 #include "Utilities/Requires.hpp"
48 #include "Utilities/TMPL.hpp"
49 
50 namespace helpers = LinearSolverAlgorithmTestHelpers;
51 
52 namespace DistributedLinearSolverAlgorithmTestHelpers {
53 
54 namespace OptionTags {
56  static constexpr Options::String help =
57  "The number of elements to distribute work on.";
58  using type = size_t;
59 };
60 } // namespace OptionTags
61 
62 /// [array_allocation_tag]
63 namespace Initialization {
64 namespace Tags {
66  using type = int;
67  using option_tags = tmpl::list<OptionTags::NumberOfElements>;
68 
69  static constexpr bool pass_metavariables = false;
70  static int create_from_options(const size_t number_of_elements) noexcept {
71  return number_of_elements;
72  }
73 };
74 } // namespace Tags
75 } // namespace Initialization
76 /// [array_allocation_tag]
77 
78 namespace OptionTags {
79 // This option expects a list of N matrices that each have N*M rows and M
80 // columns, where N is the `NumberOfElements` and M is a nonzero integer.
81 // Therefore, this option specifies a (N*M,N*M) matrix that has its columns
82 // split over all elements. In a context where the linear operator represents a
83 // DG discretization, M is the number of collocation points per element.
85  static constexpr Options::String help = "The linear operator A to invert.";
87 };
88 // Both of the following options expect a list of N vectors that have a size of
89 // M each, so that they constitute a vector of total size N*M (see above).
90 struct Source {
91  static constexpr Options::String help = "The source b in the equation Ax=b.";
93 };
95  static constexpr Options::String help = "The solution x in the equation Ax=b";
97 };
98 } // namespace OptionTags
99 
102  using option_tags = tmpl::list<OptionTags::LinearOperator>;
103 
104  static constexpr bool pass_metavariables = false;
106  create_from_options(
108  linear_operator) noexcept {
109  return linear_operator;
110  }
111 };
112 
115  using option_tags = tmpl::list<OptionTags::Source>;
116 
117  static constexpr bool pass_metavariables = false;
118  static std::vector<DenseVector<double>> create_from_options(
119  const std::vector<DenseVector<double>>& source) noexcept {
120  return source;
121  }
122 };
123 
126  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
127 
128  static constexpr bool pass_metavariables = false;
129  static std::vector<DenseVector<double>> create_from_options(
130  const std::vector<DenseVector<double>>& expected_result) noexcept {
131  return expected_result;
132  }
133 };
134 
135 // The vector `x` we want to solve for
137  using type = Scalar<DataVector>;
138  static std::string name() noexcept { return "ScalarField"; }
139 };
140 
143 using operator_applied_to_fields_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 template <typename OperandTag>
156 
157 template <typename OperandTag>
159  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
160  typename ActionList, typename ParallelComponent>
161  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
163  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
165  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
166  const int array_index,
167  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
168  const ActionList /*meta*/,
169  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
170  const ParallelComponent* const /*meta*/) noexcept {
171  const auto& operator_matrices = get<LinearOperator>(box);
172  const auto number_of_elements = operator_matrices.size();
173  const auto& linear_operator = gsl::at(operator_matrices, array_index);
174  const auto number_of_grid_points = linear_operator.columns();
175  const auto& operand = get<OperandTag>(box);
176 
177  typename OperandTag::type operator_applied_to_operand{
178  number_of_grid_points * number_of_elements};
179  dgemv_('N', linear_operator.rows(), linear_operator.columns(), 1,
180  linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
181  0, operator_applied_to_operand.data(), 1);
182 
183  Parallel::contribute_to_reduction<CollectOperatorAction<OperandTag>>(
184  Parallel::ReductionData<
186  operator_applied_to_operand},
187  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
188  Parallel::get_parallel_component<ParallelComponent>(cache));
189 
190  // Terminate algorithm for now. The reduction will be broadcast to the
191  // next action which is responsible for restarting the algorithm.
192  return {std::move(box), true};
193  }
194 };
195 
196 template <typename OperandTag>
197 struct CollectOperatorAction {
198  using local_operator_applied_to_operand_tag =
200 
201  template <typename ParallelComponent, typename DbTagsList,
202  typename Metavariables, typename ScalarFieldOperandTag,
203  Requires<tmpl::list_contains_v<
204  DbTagsList, local_operator_applied_to_operand_tag>> = nullptr>
205  static void apply(db::DataBox<DbTagsList>& box,
207  const int array_index,
208  const Variables<tmpl::list<ScalarFieldOperandTag>>&
209  Ap_global_data) noexcept {
210  // This could be generalized to work on the Variables instead of the
211  // Scalar, but it's only for the purpose of this test.
212  const auto number_of_grid_points = get<LinearOperator>(box)[0].columns();
213  const auto& Ap_global = get<ScalarFieldOperandTag>(Ap_global_data).get();
214  DataVector Ap_local{number_of_grid_points};
215  std::copy(Ap_global.begin() +
216  array_index * static_cast<int>(number_of_grid_points),
217  Ap_global.begin() +
218  (array_index + 1) * static_cast<int>(number_of_grid_points),
219  Ap_local.begin());
220  db::mutate<local_operator_applied_to_operand_tag>(
221  make_not_null(&box),
222  [&Ap_local, &number_of_grid_points ](auto Ap) noexcept {
223  *Ap = typename local_operator_applied_to_operand_tag::type{
224  number_of_grid_points};
226  *Ap)) = Ap_local;
227  });
228  // Proceed with algorithm
229  Parallel::get_parallel_component<ParallelComponent>(cache)[array_index]
230  .perform_algorithm(true);
231  }
232 };
233 
234 // Checks for the correct solution after the algorithm has terminated.
235 template <typename OptionsGroup>
236 struct TestResult {
237  using const_global_cache_tags =
238  tmpl::list<ExpectedResult, helpers::ExpectedConvergenceReason>;
239 
240  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
241  typename ActionList, typename ParallelComponent>
242  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
244  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
245  const Parallel::GlobalCache<Metavariables>& /*cache*/,
246  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
247  const int array_index,
248  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
249  const ActionList /*meta*/,
250  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
251  const ParallelComponent* const /*meta*/) noexcept {
252  const auto& has_converged =
253  get<LinearSolver::Tags::HasConverged<OptionsGroup>>(box);
254  SPECTRE_PARALLEL_REQUIRE(has_converged);
256  has_converged.reason() == get<helpers::ExpectedConvergenceReason>(box));
257  const auto& expected_result =
258  gsl::at(get<ExpectedResult>(box), array_index);
259  const auto& result = get<ScalarFieldTag>(box).get();
260  for (size_t i = 0; i < expected_result.size(); i++) {
261  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
262  }
263  return {std::move(box), true};
264  }
265 };
266 
268  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
269  typename ActionList, typename ParallelComponent>
270  static auto apply(db::DataBox<DbTagsList>& box,
271  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
272  const Parallel::GlobalCache<Metavariables>& /*cache*/,
273  const int array_index, const ActionList /*meta*/,
274  const ParallelComponent* const /*meta*/) noexcept {
275  const auto& source = gsl::at(get<Source>(box), array_index);
276  const size_t num_points = source.size();
277 
278  return std::make_tuple(
281  std::move(box), typename fields_tag::type{num_points, 0.},
282  typename sources_tag::type{source}));
283  }
284 };
285 
286 namespace detail {
287 
288 template <typename Preconditioner>
289 struct run_preconditioner_impl {
290  using type =
291  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
292  typename Preconditioner::prepare_solve,
293  ::Actions::RepeatUntil<
295  typename Preconditioner::options_group>,
296  tmpl::list<typename Preconditioner::prepare_step,
298  typename Preconditioner::operand_tag>,
299  typename Preconditioner::perform_step>>>;
300 };
301 
302 template <>
303 struct run_preconditioner_impl<void> {
304  using type = tmpl::list<>;
305 };
306 
307 template <typename Preconditioner>
308 using run_preconditioner =
309  typename run_preconditioner_impl<Preconditioner>::type;
310 
311 } // namespace detail
312 
313 template <typename Metavariables>
314 struct ElementArray {
315  using chare_type = Parallel::Algorithms::Array;
316  using metavariables = Metavariables;
317  using linear_solver = typename Metavariables::linear_solver;
318  using preconditioner = typename Metavariables::preconditioner;
319  using phase_dependent_action_list = tmpl::list<
321  typename Metavariables::Phase, Metavariables::Phase::Initialization,
322  tmpl::list<InitializeElement,
323  typename linear_solver::initialize_element,
325  helpers::detail::init_preconditioner<preconditioner>,
328  typename Metavariables::Phase,
329  Metavariables::Phase::RegisterWithObserver,
330  tmpl::list<typename linear_solver::register_element,
331  helpers::detail::register_preconditioner<preconditioner>,
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_GlobalCache<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_GlobalCache<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
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:236
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:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:63
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:70
GlobalCache.hpp
vector
Error.hpp
TestingFramework.hpp
DenseVector.hpp
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:52
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:124
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:84
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::NumberOfElements
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:55
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
algorithm
Tags.hpp
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:90
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:314
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:158
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
DistributedLinearSolverAlgorithmTestHelpers::ScalarFieldTag
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:136
DistributedLinearSolverAlgorithmTestHelpers::Initialization::Tags::NumberOfElements
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:65
array
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
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:355
DistributedLinearSolverAlgorithmTestHelpers::InitializeElement
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:267
Variables.hpp
LinearSolver::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the linear solver has converged,...
Definition: Tags.hpp:234
db::AddSimpleTags
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1132
DistributedLinearSolverAlgorithmTestHelpers::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:100
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
DistributedLinearSolverAlgorithmTestHelpers::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:113
DistributedLinearSolverAlgorithmTestHelpers::CollectOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:155
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
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:94
Prefixes.hpp
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
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:172
string