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