LinearSolverAlgorithmTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <algorithm>
9 #include <cstddef>
10 #include <string>
11 #include <tuple>
12 #include <vector>
13 
14 #include "AlgorithmArray.hpp"
15 #include "AlgorithmSingleton.hpp"
18 #include "DataStructures/DataBox/Tag.hpp"
21 #include "ErrorHandling/Error.hpp"
23 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
24 #include "IO/Observer/Helpers.hpp"
25 #include "IO/Observer/ObserverComponent.hpp"
26 #include "IO/Observer/Tags.hpp"
27 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
28 #include "NumericalAlgorithms/Convergence/Reason.hpp"
29 #include "NumericalAlgorithms/Convergence/Tags.hpp"
30 #include "Options/Options.hpp"
31 #include "Parallel/Actions/Goto.hpp"
32 #include "Parallel/Actions/SetupDataBox.hpp"
33 #include "Parallel/Actions/TerminatePhase.hpp"
34 #include "Parallel/GlobalCache.hpp"
35 #include "Parallel/InitializationFunctions.hpp"
36 #include "Parallel/Invoke.hpp"
37 #include "Parallel/Main.hpp"
38 #include "Parallel/ParallelComponentHelpers.hpp"
39 #include "Parallel/PhaseDependentActionList.hpp"
40 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
42 #include "Utilities/FileSystem.hpp"
43 #include "Utilities/Gsl.hpp"
44 #include "Utilities/Requires.hpp"
45 #include "Utilities/TMPL.hpp"
46 #include "Utilities/TaggedTuple.hpp"
47 
48 namespace LinearSolverAlgorithmTestHelpers {
49 
50 namespace OptionTags {
52  static constexpr Options::String help = "The linear operator A to invert.";
53  using type = DenseMatrix<double>;
54 };
55 struct Source {
56  static constexpr Options::String help = "The source b in the equation Ax=b.";
57  using type = DenseVector<double>;
58 };
59 struct InitialGuess {
60  static constexpr Options::String help = "The initial guess for the vector x.";
61  using type = DenseVector<double>;
62 };
64  static constexpr Options::String help = "The solution x in the equation Ax=b";
65  using type = DenseVector<double>;
66 };
68  static std::string name() noexcept { return "ConvergenceReason"; }
69  static constexpr Options::String help = "The expected convergence reason";
70  using type = Convergence::Reason;
71 };
72 } // namespace OptionTags
73 
75  using type = DenseMatrix<double>;
76  using option_tags = tmpl::list<OptionTags::LinearOperator>;
77 
78  static constexpr bool pass_metavariables = false;
79  static DenseMatrix<double> create_from_options(
80  const DenseMatrix<double>& linear_operator) noexcept {
81  return linear_operator;
82  }
83 };
84 
86  using type = DenseVector<double>;
87  using option_tags = tmpl::list<OptionTags::Source>;
88 
89  static constexpr bool pass_metavariables = false;
90  static DenseVector<double> create_from_options(
91  const DenseVector<double>& source) noexcept {
92  return source;
93  }
94 };
95 
97  using type = DenseVector<double>;
98  using option_tags = tmpl::list<OptionTags::InitialGuess>;
99 
100  static constexpr bool pass_metavariables = false;
101  static DenseVector<double> create_from_options(
102  const DenseVector<double>& initial_guess) noexcept {
103  return initial_guess;
104  }
105 };
106 
108  using type = DenseVector<double>;
109  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
110 
111  static constexpr bool pass_metavariables = false;
112  static DenseVector<double> create_from_options(
113  const DenseVector<double>& expected_result) noexcept {
114  return expected_result;
115  }
116 };
117 
119  using type = Convergence::Reason;
120  using option_tags = tmpl::list<OptionTags::ExpectedConvergenceReason>;
121 
122  static constexpr bool pass_metavariables = false;
123  static type create_from_options(const type& option_value) noexcept {
124  return option_value;
125  }
126 };
127 
128 // The vector `x` we want to solve for
130  using type = DenseVector<double>;
131 };
132 
133 using fields_tag = VectorTag;
134 
135 template <typename OperandTag>
137  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
138  typename ActionList, typename ParallelComponent>
139  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
140  db::DataBox<DbTagsList>& box,
141  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
142  const Parallel::GlobalCache<Metavariables>& /*cache*/,
143  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
144  const int /*array_index*/,
145  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
146  const ActionList /*meta*/,
147  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
148  const ParallelComponent* const /*meta*/) noexcept {
149  db::mutate<LinearSolver::Tags::OperatorAppliedTo<OperandTag>>(
150  make_not_null(&box),
152  operator_applied_to_operand,
153  const DenseMatrix<double>& linear_operator,
154  const DenseVector<double>& operand) noexcept {
155  *operator_applied_to_operand = linear_operator * operand;
156  },
157  get<LinearOperator>(box), get<OperandTag>(box));
158  return {std::move(box), false};
159  }
160 };
161 
162 // Checks for the correct solution after the algorithm has terminated.
163 template <typename OptionsGroup>
164 struct TestResult {
165  using const_global_cache_tags =
166  tmpl::list<ExpectedResult, ExpectedConvergenceReason>;
167 
168  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
169  typename ActionList, typename ParallelComponent>
170  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
171  db::DataBox<DbTagsList>& box,
172  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
173  const Parallel::GlobalCache<Metavariables>& /*cache*/,
174  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
175  const int /*array_index*/,
176  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
177  const ActionList /*meta*/,
178  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
179  const ParallelComponent* const /*meta*/) noexcept {
180  const auto& has_converged =
181  get<Convergence::Tags::HasConverged<OptionsGroup>>(box);
182  SPECTRE_PARALLEL_REQUIRE(has_converged);
183  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
184  get<ExpectedConvergenceReason>(box));
185  const auto& result = get<VectorTag>(box);
186  const auto& expected_result = get<ExpectedResult>(box);
187  for (size_t i = 0; i < expected_result.size(); i++) {
188  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
189  }
190  return {std::move(box), true};
191  }
192 };
193 
195  using simple_tags = tmpl::list<VectorTag, ::Tags::FixedSource<VectorTag>>;
196  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
197  typename ActionList, typename ParallelComponent>
198  static auto apply(db::DataBox<DbTagsList>& box,
199  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
200  const Parallel::GlobalCache<Metavariables>& /*cache*/,
201  const int /*array_index*/, const ActionList /*meta*/,
202  const ParallelComponent* const /*meta*/) noexcept {
203  Initialization::mutate_assign<simple_tags>(
204  make_not_null(&box), get<InitialGuess>(box), get<Source>(box));
205  return std::make_tuple(std::move(box));
206  }
207 };
208 
209 namespace detail {
210 
211 template <typename Preconditioner>
212 struct init_preconditioner_impl {
213  using type = typename Preconditioner::initialize_element;
214 };
215 
216 template <>
217 struct init_preconditioner_impl<void> {
218  using type = tmpl::list<>;
219 };
220 
221 template <typename Preconditioner>
222 using init_preconditioner =
223  typename init_preconditioner_impl<Preconditioner>::type;
224 
225 template <typename Preconditioner>
226 struct register_preconditioner_impl {
227  using type = typename Preconditioner::register_element;
228 };
229 
230 template <>
231 struct register_preconditioner_impl<void> {
232  using type = tmpl::list<>;
233 };
234 
235 template <typename Preconditioner>
236 using register_preconditioner =
237  typename register_preconditioner_impl<Preconditioner>::type;
238 
239 template <typename Preconditioner>
240 struct run_preconditioner_impl {
241  using type =
242  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
243  typename Preconditioner::template solve<ComputeOperatorAction<
244  typename Preconditioner::operand_tag>>>;
245 };
246 
247 template <>
248 struct run_preconditioner_impl<void> {
249  using type = tmpl::list<>;
250 };
251 
252 template <typename Preconditioner>
253 using run_preconditioner =
254  typename run_preconditioner_impl<Preconditioner>::type;
255 
256 } // namespace detail
257 
258 template <typename Metavariables>
259 struct ElementArray {
260  using chare_type = Parallel::Algorithms::Array;
261  using array_index = int;
262  using metavariables = Metavariables;
263  using linear_solver = typename Metavariables::linear_solver;
264  using preconditioner = typename Metavariables::preconditioner;
265 
266  // In each step of the algorithm we must provide A(p). The linear solver then
267  // takes care of updating x and p, as well as the internal variables r, its
268  // magnitude and the iteration step number.
269  /// [action_list]
270  using phase_dependent_action_list = tmpl::list<
272  typename Metavariables::Phase, Metavariables::Phase::Initialization,
274  typename linear_solver::initialize_element,
276  detail::init_preconditioner<preconditioner>,
279  typename Metavariables::Phase,
280  Metavariables::Phase::RegisterWithObserver,
281  tmpl::list<typename linear_solver::register_element,
282  detail::register_preconditioner<preconditioner>,
285  typename Metavariables::Phase,
286  Metavariables::Phase::PerformLinearSolve,
287  tmpl::list<
288  typename linear_solver::template solve<tmpl::list<
289  detail::run_preconditioner<preconditioner>,
293  typename Metavariables::Phase, Metavariables::Phase::TestResult,
294  tmpl::list<TestResult<typename linear_solver::options_group>>>>;
295  /// [action_list]
298  using const_global_cache_tags =
299  tmpl::list<LinearOperator, Source, InitialGuess, ExpectedResult>;
300 
301  static void allocate_array(
302  Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
303  const tuples::tagged_tuple_from_typelist<initialization_tags>&
304  initialization_items) noexcept {
305  auto& local_component = Parallel::get_parallel_component<ElementArray>(
306  *(global_cache.ckLocalBranch()));
307  local_component[0].insert(global_cache, initialization_items, 0);
308  local_component.doneInserting();
309  }
310 
311  static void execute_next_phase(
312  const typename Metavariables::Phase next_phase,
313  Parallel::CProxy_GlobalCache<Metavariables>& global_cache) noexcept {
314  auto& local_component = Parallel::get_parallel_component<ElementArray>(
315  *(global_cache.ckLocalBranch()));
316  local_component.start_phase(next_phase);
317  }
318 };
319 
320 // After the algorithm completes we perform a cleanup phase that checks the
321 // expected output file was written and deletes it.
322 template <bool CheckExpectedOutput>
323 struct CleanOutput {
324  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
325  typename ArrayIndex, typename ActionList,
326  typename ParallelComponent>
327  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
328  db::DataBox<DbTagsList>& box,
329  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
330  const Parallel::GlobalCache<Metavariables>& /*cache*/,
331  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
332  const ParallelComponent* const /*meta*/) noexcept {
333  const auto& reductions_file_name =
334  get<observers::Tags::ReductionFileName>(box) + ".h5";
335  if (file_system::check_if_file_exists(reductions_file_name)) {
336  file_system::rm(reductions_file_name, true);
337  } else if (CheckExpectedOutput) {
338  ERROR("Expected reductions file '" << reductions_file_name
339  << "' does not exist");
340  }
341  return {std::move(box), true};
342  }
343 };
344 
345 template <typename Metavariables>
347  using chare_type = Parallel::Algorithms::Singleton;
348  using metavariables = Metavariables;
349  using phase_dependent_action_list =
350  tmpl::list<Parallel::PhaseActions<typename Metavariables::Phase,
351  Metavariables::Phase::Initialization,
352  tmpl::list<CleanOutput<false>>>,
353 
354  Parallel::PhaseActions<typename Metavariables::Phase,
355  Metavariables::Phase::CleanOutput,
356  tmpl::list<CleanOutput<true>>>>;
357  using initialization_tags = Parallel::get_initialization_tags<
359 
360  static void execute_next_phase(
361  const typename Metavariables::Phase next_phase,
362  Parallel::CProxy_GlobalCache<Metavariables>& global_cache) noexcept {
363  auto& local_component = Parallel::get_parallel_component<OutputCleaner>(
364  *(global_cache.ckLocalBranch()));
365  local_component.start_phase(next_phase);
366  }
367 };
368 
369 enum class Phase {
370  Initialization,
371  RegisterWithObserver,
372  PerformLinearSolve,
373  TestResult,
374  CleanOutput,
375  Exit
376 };
377 
378 template <typename Metavariables>
379 Phase determine_next_phase(const Phase& current_phase,
380  const Parallel::CProxy_GlobalCache<
381  Metavariables>& /*cache_proxy*/) noexcept {
382  switch (current_phase) {
383  case Phase::Initialization:
384  return Phase::RegisterWithObserver;
385  case Phase::RegisterWithObserver:
386  return Phase::PerformLinearSolve;
387  case Phase::PerformLinearSolve:
388  return Phase::TestResult;
389  case Phase::TestResult:
390  return Phase::CleanOutput;
391  default:
392  return Phase::Exit;
393  }
394 }
395 
396 namespace detail {
397 
398 template <typename LinearSolver>
399 struct get_component_list_impl {
400  using type = typename LinearSolver::component_list;
401 };
402 
403 template <>
404 struct get_component_list_impl<void> {
405  using type = tmpl::list<>;
406 };
407 
408 template <typename LinearSolver>
409 using get_component_list = typename get_component_list_impl<LinearSolver>::type;
410 
411 } // namespace detail
412 
413 template <typename Metavariables>
414 using component_list = tmpl::push_back<
415  tmpl::append<
416  detail::get_component_list<typename Metavariables::linear_solver>,
417  detail::get_component_list<typename Metavariables::preconditioner>>,
418  ElementArray<Metavariables>, observers::Observer<Metavariables>,
419  observers::ObserverWriter<Metavariables>, OutputCleaner<Metavariables>>;
420 
421 template <typename Metavariables>
422 using observed_reduction_data_tags =
423  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
424  typename Metavariables::linear_solver,
425  tmpl::conditional_t<
426  std::is_same_v<typename Metavariables::preconditioner, void>,
427  tmpl::list<>, typename Metavariables::preconditioner>>>>;
428 
429 } // namespace LinearSolverAlgorithmTestHelpers
FloatingPointExceptions.hpp
Actions::SetupDataBox
Add into the DataBox default constructed items for the collection of tags requested by any of the act...
Definition: SetupDataBox.hpp:81
std::string
Main.hpp
LinearSolverAlgorithmTestHelpers::ExpectedResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:107
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
LinearSolverAlgorithmTestHelpers::CleanOutput
Definition: LinearSolverAlgorithmTestHelpers.hpp:323
LinearSolverAlgorithmTestHelpers::VectorTag
Definition: LinearSolverAlgorithmTestHelpers.hpp:129
GlobalCache.hpp
Options.hpp
vector
Error.hpp
TestingFramework.hpp
DenseVector.hpp
DenseMatrix.hpp
tuple
LinearSolverAlgorithmTestHelpers::ElementArray::initialization_tags
Parallel::get_initialization_tags< Parallel::get_initialization_actions_list< phase_dependent_action_list > > initialization_tags
[action_list]
Definition: LinearSolverAlgorithmTestHelpers.hpp:297
Parallel::Actions::TerminatePhase
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
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 >
algorithm
Tags.hpp
FileSystem.hpp
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:252
LinearSolverAlgorithmTestHelpers::OptionTags::ExpectedResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:63
file_system::check_if_file_exists
bool check_if_file_exists(const std::string &file)
Returns true if the regular file or link to the regular file exists.
Definition: FileSystem.cpp:159
LinearSolverAlgorithmTestHelpers::InitialGuess
Definition: LinearSolverAlgorithmTestHelpers.hpp:96
DenseMatrix< double >
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:216
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.
LinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: LinearSolverAlgorithmTestHelpers.hpp:55
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
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:51
LinearSolverAlgorithmTestHelpers::OptionTags::InitialGuess
Definition: LinearSolverAlgorithmTestHelpers.hpp:59
LinearSolverAlgorithmTestHelpers::OptionTags::ExpectedConvergenceReason
Definition: LinearSolverAlgorithmTestHelpers.hpp:67
DataBox.hpp
LinearSolverAlgorithmTestHelpers::ExpectedConvergenceReason
Definition: LinearSolverAlgorithmTestHelpers.hpp:118
LinearSolverAlgorithmTestHelpers::LinearOperator
Definition: LinearSolverAlgorithmTestHelpers.hpp:74
cstddef
Convergence::Reason
Reason
The reason the algorithm has converged.
Definition: Reason.hpp:23
LinearSolverAlgorithmTestHelpers::ElementArray::phase_dependent_action_list
tmpl::list< Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::Initialization, tmpl::list< Actions::SetupDataBox, InitializeElement, typename linear_solver::initialize_element, ComputeOperatorAction< fields_tag >, detail::init_preconditioner< preconditioner >, Parallel::Actions::TerminatePhase > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::RegisterWithObserver, tmpl::list< typename linear_solver::register_element, detail::register_preconditioner< preconditioner >, Parallel::Actions::TerminatePhase > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::PerformLinearSolve, tmpl::list< typename linear_solver::template solve< tmpl::list< detail::run_preconditioner< preconditioner >, ComputeOperatorAction< typename linear_solver::operand_tag > >>, Parallel::Actions::TerminatePhase > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::TestResult, tmpl::list< TestResult< typename linear_solver::options_group > >> > phase_dependent_action_list
[action_list]
Definition: LinearSolverAlgorithmTestHelpers.hpp:294
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:29
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:346
LinearSolverAlgorithmTestHelpers::InitializeElement
Definition: LinearSolverAlgorithmTestHelpers.hpp:194
file_system::rm
void rm(const std::string &path, bool recursive)
Deletes a file or directory.
Definition: FileSystem.cpp:234
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
LinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: LinearSolverAlgorithmTestHelpers.hpp:51
Requires.hpp
LinearSolverAlgorithmTestHelpers::ElementArray
Definition: LinearSolverAlgorithmTestHelpers.hpp:259
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
LinearSolverAlgorithmTestHelpers::TestResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:164
Prefixes.hpp
LinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: LinearSolverAlgorithmTestHelpers.hpp:136
LinearSolverAlgorithmTestHelpers::Source
Definition: LinearSolverAlgorithmTestHelpers.hpp:85
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string