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.hpp"
24 #include "IO/Observer/Helpers.hpp"
25 #include "IO/Observer/ObserverComponent.hpp"
26 #include "IO/Observer/RegisterObservers.hpp"
27 #include "IO/Observer/Tags.hpp"
28 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
29 #include "NumericalAlgorithms/Convergence/Reason.hpp"
30 #include "Options/Options.hpp"
31 #include "Parallel/Actions/Goto.hpp"
32 #include "Parallel/Actions/TerminatePhase.hpp"
34 #include "Parallel/InitializationFunctions.hpp"
35 #include "Parallel/Invoke.hpp"
36 #include "Parallel/Main.hpp"
37 #include "Parallel/ParallelComponentHelpers.hpp"
38 #include "Parallel/PhaseDependentActionList.hpp"
39 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
40 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.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 OptionString help = "The linear operator A to invert.";
53  using type = DenseMatrix<double>;
54 };
55 struct Source {
56  static constexpr OptionString help = "The source b in the equation Ax=b.";
57  using type = DenseVector<double>;
58 };
59 struct InitialGuess {
60  static constexpr OptionString help = "The initial guess for the vector x.";
61  using type = DenseVector<double>;
62 };
64  static constexpr OptionString 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 OptionString 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(
141  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
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(
172  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
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<LinearSolver::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  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
196  typename ActionList, typename ParallelComponent>
197  static auto apply(db::DataBox<DbTagsList>& box,
198  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
200  const int /*array_index*/, const ActionList /*meta*/,
201  const ParallelComponent* const /*meta*/) noexcept {
202  return std::make_tuple(
206  std::move(box), get<InitialGuess>(box), get<Source>(box)));
207  }
208 };
209 
210 namespace detail {
211 
212 template <typename Preconditioner>
213 struct init_preconditioner_impl {
214  using type = typename Preconditioner::initialize_element;
215 };
216 
217 template <>
218 struct init_preconditioner_impl<void> {
219  using type = tmpl::list<>;
220 };
221 
222 template <typename Preconditioner>
223 using init_preconditioner =
224  typename init_preconditioner_impl<Preconditioner>::type;
225 
226 template <typename Preconditioner>
227 struct run_preconditioner_impl {
228  using type =
229  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
230  typename Preconditioner::prepare_solve,
231  ::Actions::RepeatUntil<
233  typename Preconditioner::options_group>,
234  tmpl::list<typename Preconditioner::prepare_step,
235  ComputeOperatorAction<
236  typename Preconditioner::operand_tag>,
237  typename Preconditioner::perform_step>>>;
238 };
239 
240 template <>
241 struct run_preconditioner_impl<void> {
242  using type = tmpl::list<>;
243 };
244 
245 template <typename Preconditioner>
246 using run_preconditioner =
247  typename run_preconditioner_impl<Preconditioner>::type;
248 
249 } // namespace detail
250 
251 template <typename Metavariables>
252 struct ElementArray {
253  using chare_type = Parallel::Algorithms::Array;
254  using array_index = int;
255  using metavariables = Metavariables;
256  using linear_solver = typename Metavariables::linear_solver;
257  using preconditioner = typename Metavariables::preconditioner;
258 
259  // In each step of the algorithm we must provide A(p). The linear solver then
260  // takes care of updating x and p, as well as the internal variables r, its
261  // magnitude and the iteration step number.
262  /// [action_list]
263  using phase_dependent_action_list = tmpl::list<
265  typename Metavariables::Phase, Metavariables::Phase::Initialization,
266  tmpl::list<InitializeElement,
267  typename linear_solver::initialize_element,
269  detail::init_preconditioner<preconditioner>,
272  typename Metavariables::Phase,
273  Metavariables::Phase::RegisterWithObserver,
274  tmpl::list<typename linear_solver::register_element,
275  typename linear_solver::prepare_solve,
278  typename Metavariables::Phase,
279  Metavariables::Phase::PerformLinearSolve,
281  typename linear_solver::options_group>,
282  typename linear_solver::prepare_step,
283  detail::run_preconditioner<preconditioner>,
285  typename linear_solver::perform_step>>,
286 
288  typename Metavariables::Phase, Metavariables::Phase::TestResult,
289  tmpl::list<TestResult<typename linear_solver::options_group>>>>;
290  /// [action_list]
293  using const_global_cache_tags =
294  tmpl::list<LinearOperator, Source, InitialGuess, ExpectedResult>;
295 
296  static void allocate_array(
297  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache,
298  const tuples::tagged_tuple_from_typelist<initialization_tags>&
299  initialization_items) noexcept {
300  auto& local_component = Parallel::get_parallel_component<ElementArray>(
301  *(global_cache.ckLocalBranch()));
302  local_component[0].insert(global_cache, initialization_items, 0);
303  local_component.doneInserting();
304  }
305 
306  static void execute_next_phase(
307  const typename Metavariables::Phase next_phase,
308  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
309  auto& local_component = Parallel::get_parallel_component<ElementArray>(
310  *(global_cache.ckLocalBranch()));
311  local_component.start_phase(next_phase);
312  }
313 };
314 
315 // After the algorithm completes we perform a cleanup phase that checks the
316 // expected output file was written and deletes it.
317 template <bool CheckExpectedOutput>
318 struct CleanOutput {
319  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
320  typename ArrayIndex, typename ActionList,
321  typename ParallelComponent>
322  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
324  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
326  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
327  const ParallelComponent* const /*meta*/) noexcept {
328  const auto& reductions_file_name =
329  get<observers::Tags::ReductionFileName>(box) + ".h5";
330  if (file_system::check_if_file_exists(reductions_file_name)) {
331  file_system::rm(reductions_file_name, true);
332  } else if (CheckExpectedOutput) {
333  ERROR("Expected reductions file '" << reductions_file_name
334  << "' does not exist");
335  }
336  return {std::move(box), true};
337  }
338 };
339 
340 template <typename Metavariables>
342  using chare_type = Parallel::Algorithms::Singleton;
343  using metavariables = Metavariables;
344  using phase_dependent_action_list =
345  tmpl::list<Parallel::PhaseActions<typename Metavariables::Phase,
346  Metavariables::Phase::Initialization,
347  tmpl::list<CleanOutput<false>>>,
348 
349  Parallel::PhaseActions<typename Metavariables::Phase,
350  Metavariables::Phase::CleanOutput,
351  tmpl::list<CleanOutput<true>>>>;
352  using initialization_tags = Parallel::get_initialization_tags<
354 
355  static void execute_next_phase(
356  const typename Metavariables::Phase next_phase,
357  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
358  auto& local_component = Parallel::get_parallel_component<OutputCleaner>(
359  *(global_cache.ckLocalBranch()));
360  local_component.start_phase(next_phase);
361  }
362 };
363 
364 enum class Phase {
365  Initialization,
366  RegisterWithObserver,
367  PerformLinearSolve,
368  TestResult,
369  CleanOutput,
370  Exit
371 };
372 
373 template <typename Metavariables>
374 Phase determine_next_phase(const Phase& current_phase,
375  const Parallel::CProxy_ConstGlobalCache<
376  Metavariables>& /*cache_proxy*/) noexcept {
377  switch (current_phase) {
378  case Phase::Initialization:
379  return Phase::RegisterWithObserver;
380  case Phase::RegisterWithObserver:
381  return Phase::PerformLinearSolve;
382  case Phase::PerformLinearSolve:
383  return Phase::TestResult;
384  case Phase::TestResult:
385  return Phase::CleanOutput;
386  default:
387  return Phase::Exit;
388  }
389 }
390 
391 namespace detail {
392 
393 template <typename LinearSolver>
394 struct get_component_list_impl {
395  using type = typename LinearSolver::component_list;
396 };
397 
398 template <>
399 struct get_component_list_impl<void> {
400  using type = tmpl::list<>;
401 };
402 
403 template <typename LinearSolver>
404 using get_component_list = typename get_component_list_impl<LinearSolver>::type;
405 
406 } // namespace detail
407 
408 template <typename Metavariables>
409 using component_list = tmpl::push_back<
410  tmpl::append<
411  detail::get_component_list<typename Metavariables::linear_solver>,
412  detail::get_component_list<typename Metavariables::preconditioner>>,
413  ElementArray<Metavariables>, observers::Observer<Metavariables>,
414  observers::ObserverWriter<Metavariables>, OutputCleaner<Metavariables>>;
415 
416 template <typename Metavariables>
417 using observed_reduction_data_tags =
418  observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
419  typename Metavariables::linear_solver,
420  tmpl::conditional_t<
421  std::is_same_v<typename Metavariables::preconditioner, void>,
422  tmpl::list<>, typename Metavariables::preconditioner>>>>;
423 
424 } // namespace LinearSolverAlgorithmTestHelpers
Parallel::ConstGlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
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:107
LinearSolverAlgorithmTestHelpers::CleanOutput
Definition: LinearSolverAlgorithmTestHelpers.hpp:318
LinearSolverAlgorithmTestHelpers::VectorTag
Definition: LinearSolverAlgorithmTestHelpers.hpp:129
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:292
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: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: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: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:48
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:21
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
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:341
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
LinearSolver::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the linear solver has converged,...
Definition: Tags.hpp:242
db::AddSimpleTags
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1150
Gsl.hpp
LinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: LinearSolverAlgorithmTestHelpers.hpp:51
Requires.hpp
LinearSolverAlgorithmTestHelpers::ElementArray
Definition: LinearSolverAlgorithmTestHelpers.hpp:252
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:164
Prefixes.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, 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:289
LinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: LinearSolverAlgorithmTestHelpers.hpp:136
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
LinearSolverAlgorithmTestHelpers::Source
Definition: LinearSolverAlgorithmTestHelpers.hpp:85
OptionString
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
TMPL.hpp
ConstGlobalCache.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string