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 "Options/Options.hpp"
30 #include "Parallel/Actions/Goto.hpp"
31 #include "Parallel/Actions/TerminatePhase.hpp"
32 #include "Parallel/GlobalCache.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 "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
39 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
41 #include "Utilities/FileSystem.hpp"
42 #include "Utilities/Gsl.hpp"
43 #include "Utilities/Requires.hpp"
44 #include "Utilities/TMPL.hpp"
45 #include "Utilities/TaggedTuple.hpp"
46 
47 namespace LinearSolverAlgorithmTestHelpers {
48 
49 namespace OptionTags {
51  static constexpr Options::String help = "The linear operator A to invert.";
52  using type = DenseMatrix<double>;
53 };
54 struct Source {
55  static constexpr Options::String help = "The source b in the equation Ax=b.";
56  using type = DenseVector<double>;
57 };
58 struct InitialGuess {
59  static constexpr Options::String help = "The initial guess for the vector x.";
60  using type = DenseVector<double>;
61 };
63  static constexpr Options::String help = "The solution x in the equation Ax=b";
64  using type = DenseVector<double>;
65 };
67  static std::string name() noexcept { return "ConvergenceReason"; }
68  static constexpr Options::String help = "The expected convergence reason";
69  using type = Convergence::Reason;
70 };
71 } // namespace OptionTags
72 
74  using type = DenseMatrix<double>;
75  using option_tags = tmpl::list<OptionTags::LinearOperator>;
76 
77  static constexpr bool pass_metavariables = false;
78  static DenseMatrix<double> create_from_options(
79  const DenseMatrix<double>& linear_operator) noexcept {
80  return linear_operator;
81  }
82 };
83 
85  using type = DenseVector<double>;
86  using option_tags = tmpl::list<OptionTags::Source>;
87 
88  static constexpr bool pass_metavariables = false;
89  static DenseVector<double> create_from_options(
90  const DenseVector<double>& source) noexcept {
91  return source;
92  }
93 };
94 
96  using type = DenseVector<double>;
97  using option_tags = tmpl::list<OptionTags::InitialGuess>;
98 
99  static constexpr bool pass_metavariables = false;
100  static DenseVector<double> create_from_options(
101  const DenseVector<double>& initial_guess) noexcept {
102  return initial_guess;
103  }
104 };
105 
107  using type = DenseVector<double>;
108  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
109 
110  static constexpr bool pass_metavariables = false;
111  static DenseVector<double> create_from_options(
112  const DenseVector<double>& expected_result) noexcept {
113  return expected_result;
114  }
115 };
116 
118  using type = Convergence::Reason;
119  using option_tags = tmpl::list<OptionTags::ExpectedConvergenceReason>;
120 
121  static constexpr bool pass_metavariables = false;
122  static type create_from_options(const type& option_value) noexcept {
123  return option_value;
124  }
125 };
126 
127 // The vector `x` we want to solve for
129  using type = DenseVector<double>;
130 };
131 
132 using fields_tag = VectorTag;
133 
134 template <typename OperandTag>
136  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
137  typename ActionList, typename ParallelComponent>
138  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
140  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
141  const Parallel::GlobalCache<Metavariables>& /*cache*/,
142  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
143  const int /*array_index*/,
144  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
145  const ActionList /*meta*/,
146  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
147  const ParallelComponent* const /*meta*/) noexcept {
148  db::mutate<LinearSolver::Tags::OperatorAppliedTo<OperandTag>>(
149  make_not_null(&box),
151  operator_applied_to_operand,
152  const DenseMatrix<double>& linear_operator,
153  const DenseVector<double>& operand) noexcept {
154  *operator_applied_to_operand = linear_operator * operand;
155  },
156  get<LinearOperator>(box), get<OperandTag>(box));
157  return {std::move(box), false};
158  }
159 };
160 
161 // Checks for the correct solution after the algorithm has terminated.
162 template <typename OptionsGroup>
163 struct TestResult {
164  using const_global_cache_tags =
165  tmpl::list<ExpectedResult, ExpectedConvergenceReason>;
166 
167  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
168  typename ActionList, typename ParallelComponent>
169  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
171  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
172  const Parallel::GlobalCache<Metavariables>& /*cache*/,
173  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
174  const int /*array_index*/,
175  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
176  const ActionList /*meta*/,
177  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
178  const ParallelComponent* const /*meta*/) noexcept {
179  const auto& has_converged =
180  get<LinearSolver::Tags::HasConverged<OptionsGroup>>(box);
181  SPECTRE_PARALLEL_REQUIRE(has_converged);
182  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
183  get<ExpectedConvergenceReason>(box));
184  const auto& result = get<VectorTag>(box);
185  const auto& expected_result = get<ExpectedResult>(box);
186  for (size_t i = 0; i < expected_result.size(); i++) {
187  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
188  }
189  return {std::move(box), true};
190  }
191 };
192 
194  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
195  typename ActionList, typename ParallelComponent>
196  static auto apply(db::DataBox<DbTagsList>& box,
197  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
198  const Parallel::GlobalCache<Metavariables>& /*cache*/,
199  const int /*array_index*/, const ActionList /*meta*/,
200  const ParallelComponent* const /*meta*/) noexcept {
201  return std::make_tuple(
205  std::move(box), get<InitialGuess>(box), get<Source>(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::prepare_solve,
244  ::Actions::RepeatUntil<
246  typename Preconditioner::options_group>,
247  tmpl::list<typename Preconditioner::prepare_step,
248  ComputeOperatorAction<
249  typename Preconditioner::operand_tag>,
250  typename Preconditioner::perform_step>>>;
251 };
252 
253 template <>
254 struct run_preconditioner_impl<void> {
255  using type = tmpl::list<>;
256 };
257 
258 template <typename Preconditioner>
259 using run_preconditioner =
260  typename run_preconditioner_impl<Preconditioner>::type;
261 
262 } // namespace detail
263 
264 template <typename Metavariables>
265 struct ElementArray {
266  using chare_type = Parallel::Algorithms::Array;
267  using array_index = int;
268  using metavariables = Metavariables;
269  using linear_solver = typename Metavariables::linear_solver;
270  using preconditioner = typename Metavariables::preconditioner;
271 
272  // In each step of the algorithm we must provide A(p). The linear solver then
273  // takes care of updating x and p, as well as the internal variables r, its
274  // magnitude and the iteration step number.
275  /// [action_list]
276  using phase_dependent_action_list = tmpl::list<
278  typename Metavariables::Phase, Metavariables::Phase::Initialization,
279  tmpl::list<InitializeElement,
280  typename linear_solver::initialize_element,
282  detail::init_preconditioner<preconditioner>,
285  typename Metavariables::Phase,
286  Metavariables::Phase::RegisterWithObserver,
287  tmpl::list<typename linear_solver::register_element,
288  detail::register_preconditioner<preconditioner>,
289  typename linear_solver::prepare_solve,
292  typename Metavariables::Phase,
293  Metavariables::Phase::PerformLinearSolve,
295  typename linear_solver::options_group>,
296  typename linear_solver::prepare_step,
297  detail::run_preconditioner<preconditioner>,
299  typename linear_solver::perform_step>>,
300 
302  typename Metavariables::Phase, Metavariables::Phase::TestResult,
303  tmpl::list<TestResult<typename linear_solver::options_group>>>>;
304  /// [action_list]
307  using const_global_cache_tags =
308  tmpl::list<LinearOperator, Source, InitialGuess, ExpectedResult>;
309 
310  static void allocate_array(
311  Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
312  const tuples::tagged_tuple_from_typelist<initialization_tags>&
313  initialization_items) noexcept {
314  auto& local_component = Parallel::get_parallel_component<ElementArray>(
315  *(global_cache.ckLocalBranch()));
316  local_component[0].insert(global_cache, initialization_items, 0);
317  local_component.doneInserting();
318  }
319 
320  static void execute_next_phase(
321  const typename Metavariables::Phase next_phase,
322  Parallel::CProxy_GlobalCache<Metavariables>& global_cache) noexcept {
323  auto& local_component = Parallel::get_parallel_component<ElementArray>(
324  *(global_cache.ckLocalBranch()));
325  local_component.start_phase(next_phase);
326  }
327 };
328 
329 // After the algorithm completes we perform a cleanup phase that checks the
330 // expected output file was written and deletes it.
331 template <bool CheckExpectedOutput>
332 struct CleanOutput {
333  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
334  typename ArrayIndex, typename ActionList,
335  typename ParallelComponent>
336  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
338  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
339  const Parallel::GlobalCache<Metavariables>& /*cache*/,
340  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
341  const ParallelComponent* const /*meta*/) noexcept {
342  const auto& reductions_file_name =
343  get<observers::Tags::ReductionFileName>(box) + ".h5";
344  if (file_system::check_if_file_exists(reductions_file_name)) {
345  file_system::rm(reductions_file_name, true);
346  } else if (CheckExpectedOutput) {
347  ERROR("Expected reductions file '" << reductions_file_name
348  << "' does not exist");
349  }
350  return {std::move(box), true};
351  }
352 };
353 
354 template <typename Metavariables>
356  using chare_type = Parallel::Algorithms::Singleton;
357  using metavariables = Metavariables;
358  using phase_dependent_action_list =
359  tmpl::list<Parallel::PhaseActions<typename Metavariables::Phase,
360  Metavariables::Phase::Initialization,
361  tmpl::list<CleanOutput<false>>>,
362 
363  Parallel::PhaseActions<typename Metavariables::Phase,
364  Metavariables::Phase::CleanOutput,
365  tmpl::list<CleanOutput<true>>>>;
366  using initialization_tags = Parallel::get_initialization_tags<
368 
369  static void execute_next_phase(
370  const typename Metavariables::Phase next_phase,
371  Parallel::CProxy_GlobalCache<Metavariables>& global_cache) noexcept {
372  auto& local_component = Parallel::get_parallel_component<OutputCleaner>(
373  *(global_cache.ckLocalBranch()));
374  local_component.start_phase(next_phase);
375  }
376 };
377 
378 enum class Phase {
379  Initialization,
380  RegisterWithObserver,
381  PerformLinearSolve,
382  TestResult,
383  CleanOutput,
384  Exit
385 };
386 
387 template <typename Metavariables>
388 Phase determine_next_phase(const Phase& current_phase,
389  const Parallel::CProxy_GlobalCache<
390  Metavariables>& /*cache_proxy*/) noexcept {
391  switch (current_phase) {
392  case Phase::Initialization:
393  return Phase::RegisterWithObserver;
394  case Phase::RegisterWithObserver:
395  return Phase::PerformLinearSolve;
396  case Phase::PerformLinearSolve:
397  return Phase::TestResult;
398  case Phase::TestResult:
399  return Phase::CleanOutput;
400  default:
401  return Phase::Exit;
402  }
403 }
404 
405 namespace detail {
406 
407 template <typename LinearSolver>
408 struct get_component_list_impl {
409  using type = typename LinearSolver::component_list;
410 };
411 
412 template <>
413 struct get_component_list_impl<void> {
414  using type = tmpl::list<>;
415 };
416 
417 template <typename LinearSolver>
418 using get_component_list = typename get_component_list_impl<LinearSolver>::type;
419 
420 } // namespace detail
421 
422 template <typename Metavariables>
423 using component_list = tmpl::push_back<
424  tmpl::append<
425  detail::get_component_list<typename Metavariables::linear_solver>,
426  detail::get_component_list<typename Metavariables::preconditioner>>,
427  ElementArray<Metavariables>, observers::Observer<Metavariables>,
428  observers::ObserverWriter<Metavariables>, OutputCleaner<Metavariables>>;
429 
430 template <typename Metavariables>
431 using observed_reduction_data_tags =
432  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
433  typename Metavariables::linear_solver,
434  tmpl::conditional_t<
435  std::is_same_v<typename Metavariables::preconditioner, void>,
436  tmpl::list<>, typename Metavariables::preconditioner>>>>;
437 
438 } // namespace LinearSolverAlgorithmTestHelpers
FloatingPointExceptions.hpp
LinearSolver::Actions::TerminateIfConverged
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:33
std::string
Main.hpp
LinearSolverAlgorithmTestHelpers::ExpectedResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:106
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
LinearSolverAlgorithmTestHelpers::CleanOutput
Definition: LinearSolverAlgorithmTestHelpers.hpp:332
LinearSolverAlgorithmTestHelpers::VectorTag
Definition: LinearSolverAlgorithmTestHelpers.hpp:128
GlobalCache.hpp
LinearSolverAlgorithmTestHelpers::ElementArray::phase_dependent_action_list
tmpl::list< Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::Initialization, tmpl::list< 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 >, typename linear_solver::prepare_solve, Parallel::Actions::TerminatePhase > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::PerformLinearSolve, tmpl::list< LinearSolver::Actions::TerminateIfConverged< typename linear_solver::options_group >, typename linear_solver::prepare_step, detail::run_preconditioner< preconditioner >, ComputeOperatorAction< typename linear_solver::operand_tag >, typename linear_solver::perform_step > >, 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:303
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:306
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 >
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:140
LinearSolverAlgorithmTestHelpers::OptionTags::ExpectedResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:62
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:95
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: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.
Tags::FixedSource
Prefix indicating a source term that is independent of dynamic variables.
Definition: Prefixes.hpp:75
LinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: LinearSolverAlgorithmTestHelpers.hpp:54
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:48
LinearSolverAlgorithmTestHelpers::OptionTags::InitialGuess
Definition: LinearSolverAlgorithmTestHelpers.hpp:58
LinearSolverAlgorithmTestHelpers::OptionTags::ExpectedConvergenceReason
Definition: LinearSolverAlgorithmTestHelpers.hpp:66
DataBox.hpp
LinearSolverAlgorithmTestHelpers::ExpectedConvergenceReason
Definition: LinearSolverAlgorithmTestHelpers.hpp:117
LinearSolverAlgorithmTestHelpers::LinearOperator
Definition: LinearSolverAlgorithmTestHelpers.hpp:73
cstddef
Convergence::Reason
Reason
The reason the algorithm has converged.
Definition: Reason.hpp:23
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
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:355
LinearSolverAlgorithmTestHelpers::InitializeElement
Definition: LinearSolverAlgorithmTestHelpers.hpp:193
file_system::rm
void rm(const std::string &path, bool recursive)
Deletes a file or directory.
Definition: FileSystem.cpp:234
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
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
LinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: LinearSolverAlgorithmTestHelpers.hpp:50
Requires.hpp
LinearSolverAlgorithmTestHelpers::ElementArray
Definition: LinearSolverAlgorithmTestHelpers.hpp:265
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
LinearSolverAlgorithmTestHelpers::TestResult
Definition: LinearSolverAlgorithmTestHelpers.hpp:163
Prefixes.hpp
LinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: LinearSolverAlgorithmTestHelpers.hpp:135
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
LinearSolverAlgorithmTestHelpers::Source
Definition: LinearSolverAlgorithmTestHelpers.hpp:84
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string