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"
17 #include "DataStructures/DataBox/Tag.hpp"
20 #include "ErrorHandling/Error.hpp"
22 #include "IO/Observer/Tags.hpp"
23 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
24 #include "Options/Options.hpp"
25 #include "Parallel/Actions/TerminatePhase.hpp"
27 #include "Parallel/InitializationFunctions.hpp"
28 #include "Parallel/Invoke.hpp"
29 #include "Parallel/Main.hpp"
30 #include "Parallel/ParallelComponentHelpers.hpp"
31 #include "Parallel/PhaseDependentActionList.hpp"
32 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
33 #include "ParallelAlgorithms/LinearSolver/Actions/TerminateIfConverged.hpp"
34 #include "ParallelAlgorithms/LinearSolver/Tags.hpp" // IWYU pragma: keep
35 #include "Utilities/FileSystem.hpp"
36 #include "Utilities/Gsl.hpp"
37 #include "Utilities/Requires.hpp"
38 #include "Utilities/TMPL.hpp"
39 #include "Utilities/TaggedTuple.hpp"
40 
42 
43 namespace OptionTags {
45  static constexpr OptionString help = "The linear operator A to invert.";
46  using type = DenseMatrix<double>;
47 };
48 struct Source {
49  static constexpr OptionString help = "The source b in the equation Ax=b.";
50  using type = DenseVector<double>;
51 };
52 struct InitialGuess {
53  static constexpr OptionString help = "The initial guess for the vector x.";
54  using type = DenseVector<double>;
55 };
57  static constexpr OptionString help = "The solution x in the equation Ax=b";
58  using type = DenseVector<double>;
59 };
60 } // namespace OptionTags
61 
63  using type = DenseMatrix<double>;
64  using option_tags = tmpl::list<OptionTags::LinearOperator>;
65 
66  static constexpr bool pass_metavariables = false;
68  const DenseMatrix<double>& linear_operator) noexcept {
69  return linear_operator;
70  }
71 };
72 
74  using type = DenseVector<double>;
75  using option_tags = tmpl::list<OptionTags::Source>;
76 
77  static constexpr bool pass_metavariables = false;
78  static DenseVector<double> create_from_options(
79  const DenseVector<double>& source) noexcept {
80  return source;
81  }
82 };
83 
85  using type = DenseVector<double>;
86  using option_tags = tmpl::list<OptionTags::InitialGuess>;
87 
88  static constexpr bool pass_metavariables = false;
89  static DenseVector<double> create_from_options(
90  const DenseVector<double>& initial_guess) noexcept {
91  return initial_guess;
92  }
93 };
94 
96  using type = DenseVector<double>;
97  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
98 
99  static constexpr bool pass_metavariables = false;
100  static DenseVector<double> create_from_options(
101  const DenseVector<double>& expected_result) noexcept {
102  return expected_result;
103  }
104 };
105 
106 // The vector `x` we want to solve for
108  using type = DenseVector<double>;
109 };
110 
111 using fields_tag = VectorTag;
114 
116  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
117  typename ActionList, typename ParallelComponent>
120  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
122  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
123  const int /*array_index*/,
124  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
125  const ActionList /*meta*/,
126  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
127  const ParallelComponent* const /*meta*/) noexcept {
128  db::mutate<operator_tag>(
129  make_not_null(&box),
130  [
131  ](const gsl::not_null<DenseVector<double>*> operator_applied_to_operand,
132  const DenseMatrix<double>& linear_operator,
133  const DenseVector<double>& operand) noexcept {
134  *operator_applied_to_operand = linear_operator * operand;
135  },
136  get<LinearOperator>(cache), get<operand_tag>(box));
137  return {std::move(box), false};
138  }
139 };
140 
141 // Checks for the correct solution after the algorithm has terminated.
142 struct TestResult {
143  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
144  typename ActionList, typename ParallelComponent>
147  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
149  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
150  const int /*array_index*/,
151  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
152  const ActionList /*meta*/,
153  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
154  const ParallelComponent* const /*meta*/) noexcept {
155  const auto& has_converged = get<LinearSolver::Tags::HasConverged>(box);
156  SPECTRE_PARALLEL_REQUIRE(has_converged);
157  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
158  Convergence::Reason::AbsoluteResidual);
159  const auto& result = get<VectorTag>(box);
160  const auto& expected_result = get<ExpectedResult>(cache);
161  for (size_t i = 0; i < expected_result.size(); i++) {
162  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
163  }
164  return {std::move(box), true};
165  }
166 };
167 
169  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
170  typename ActionList, typename ParallelComponent>
171  static auto apply(db::DataBox<DbTagsList>& box,
172  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
174  const int /*array_index*/, const ActionList /*meta*/,
175  const ParallelComponent* const /*meta*/) noexcept {
176  const auto& linear_operator = get<LinearOperator>(cache);
177  const auto& source = get<Source>(cache);
178  const auto& initial_guess = get<InitialGuess>(cache);
179 
180  return std::make_tuple(
186  std::move(box), initial_guess, source,
187  DenseVector<double>{linear_operator * initial_guess},
188  make_with_value<DenseVector<double>>(
190  make_with_value<DenseVector<double>>(
191  initial_guess, std::numeric_limits<double>::signaling_NaN())));
192  }
193 }; // namespace
194 
195 template <typename Metavariables>
196 struct ElementArray {
197  using chare_type = Parallel::Algorithms::Array;
198  using array_index = int;
199  using metavariables = Metavariables;
200  // In each step of the algorithm we must provide A(p). The linear solver then
201  // takes care of updating x and p, as well as the internal variables r, its
202  // magnitude and the iteration step number.
203  /// [action_list]
204  using phase_dependent_action_list = tmpl::list<
206  typename Metavariables::Phase, Metavariables::Phase::Initialization,
207  tmpl::list<InitializeElement,
208  typename Metavariables::linear_solver::initialize_element,
210 
212  typename Metavariables::Phase,
213  Metavariables::Phase::PerformLinearSolve,
215  typename Metavariables::linear_solver::prepare_step,
217  typename Metavariables::linear_solver::perform_step>>,
218 
219  Parallel::PhaseActions<typename Metavariables::Phase,
220  Metavariables::Phase::TestResult,
221  tmpl::list<TestResult>>>;
222  /// [action_list]
225  using const_global_cache_tags =
226  tmpl::list<LinearOperator, Source, InitialGuess, ExpectedResult>;
227 
228  static void allocate_array(
229  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache,
230  const tuples::tagged_tuple_from_typelist<initialization_tags>&
231  initialization_items) noexcept {
232  auto& local_component = Parallel::get_parallel_component<ElementArray>(
233  *(global_cache.ckLocalBranch()));
234  local_component[0].insert(global_cache, initialization_items, 0);
235  local_component.doneInserting();
236  }
237 
238  static void execute_next_phase(
239  const typename Metavariables::Phase next_phase,
240  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
241  auto& local_component = Parallel::get_parallel_component<ElementArray>(
242  *(global_cache.ckLocalBranch()));
243  local_component.start_phase(next_phase);
244  }
245 };
246 
247 // After the algorithm completes we perform a cleanup phase that checks the
248 // expected output file was written and deletes it.
249 template <bool CheckExpectedOutput>
250 struct CleanOutput {
251  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
252  typename ArrayIndex, typename ActionList,
253  typename ParallelComponent>
256  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
258  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
259  const ParallelComponent* const /*meta*/) noexcept {
260  const auto& reductions_file_name =
261  get<observers::Tags::ReductionFileName>(cache) + ".h5";
262  if (file_system::check_if_file_exists(reductions_file_name)) {
263  file_system::rm(reductions_file_name, true);
264  } else if (CheckExpectedOutput) {
265  ERROR("Expected reductions file '" << reductions_file_name
266  << "' does not exist");
267  }
268  return {std::move(box), true};
269  }
270 };
271 
272 template <typename Metavariables>
274  using chare_type = Parallel::Algorithms::Singleton;
275  using metavariables = Metavariables;
276  using phase_dependent_action_list =
277  tmpl::list<Parallel::PhaseActions<typename Metavariables::Phase,
278  Metavariables::Phase::Initialization,
279  tmpl::list<CleanOutput<false>>>,
280 
281  Parallel::PhaseActions<typename Metavariables::Phase,
282  Metavariables::Phase::CleanOutput,
283  tmpl::list<CleanOutput<true>>>>;
284  using initialization_tags = Parallel::get_initialization_tags<
286 
287  static void execute_next_phase(
288  const typename Metavariables::Phase next_phase,
289  Parallel::CProxy_ConstGlobalCache<Metavariables>& global_cache) noexcept {
290  auto& local_component = Parallel::get_parallel_component<OutputCleaner>(
291  *(global_cache.ckLocalBranch()));
292  local_component.start_phase(next_phase);
293  }
294 };
295 
296 } // namespace LinearSolverAlgorithmTestHelpers
tmpl::list< Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::Initialization, tmpl::list< InitializeElement, typename Metavariables::linear_solver::initialize_element, Parallel::Actions::TerminatePhase > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::PerformLinearSolve, tmpl::list< LinearSolver::Actions::TerminateIfConverged, typename Metavariables::linear_solver::prepare_step, ComputeOperatorAction, typename Metavariables::linear_solver::perform_step > >, Parallel::PhaseActions< typename Metavariables::Phase, Metavariables::Phase::TestResult, tmpl::list< TestResult > >> phase_dependent_action_list
[action_list]
Definition: LinearSolverAlgorithmTestHelpers.hpp:221
Definition: LinearSolverAlgorithmTestHelpers.hpp:115
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
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
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
Parallel::get_initialization_tags< Parallel::get_initialization_actions_list< phase_dependent_action_list > > initialization_tags
[action_list]
Definition: LinearSolverAlgorithmTestHelpers.hpp:224
T signaling_NaN(T... args)
Definition: LinearSolverAlgorithmTestHelpers.hpp:52
Definition: LinearSolverAlgorithmTestHelpers.hpp:48
Defines classes and functions for making classes creatable from input files.
Definition: LinearSolverAlgorithmTestHelpers.hpp:41
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
Definition: LinearSolverAlgorithmTestHelpers.hpp:95
Definition: LinearSolverAlgorithmTestHelpers.hpp:84
Defines the type alias Requires.
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
Definition: LinearSolverAlgorithmTestHelpers.hpp:44
auto merge_into_databox(db::DataBox< DbTagsList > &&box, Args &&... args) noexcept
Add tags that are not yet in the DataBox.
Definition: MergeIntoDataBox.hpp:133
The operand that the local linear operator is applied to.
Definition: Tags.hpp:53
Definition: LinearSolverAlgorithmTestHelpers.hpp:142
void rm(const std::string &path, bool recursive)
Deletes a file or directory.
Definition: FileSystem.cpp:234
Definition: LinearSolverAlgorithmTestHelpers.hpp:56
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
Definition: LinearSolverAlgorithmTestHelpers.hpp:273
A dynamically sized matrix of arbitrary type.
Definition: DenseMatrix.hpp:40
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
std::remove_const_t< R > make_with_value(const T &input, const ValueType &value) noexcept
Given an object of type T, create an object of type R whose elements are initialized to value...
Definition: MakeWithValue.hpp:42
Defines class DenseMatrix.
Terminate the algorithm if the linear solver has converged.
Definition: TerminateIfConverged.hpp:32
tmpl::flatten< tmpl::list< Tags... > > AddSimpleTags
List of Tags to add to the DataBox.
Definition: DataBox.hpp:1150
Definition: LinearSolverAlgorithmTestHelpers.hpp:250
Definition: Strahlkorper.hpp:167
Definition: InterpolationTargetWedgeSectionTorus.hpp:25
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:136
Definition: LinearSolverAlgorithmTestHelpers.hpp:73
Definition: LinearSolverAlgorithmTestHelpers.hpp:62
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1444
Functions to enable/disable termination on floating point exceptions.
#define SPECTRE_PARALLEL_REQUIRE(expr)
A similar to Catch&#39;s REQUIRE statement, but can be used in tests that spawn several chares with possi...
Definition: TestingFramework.hpp:65
tuples::TaggedTuple< Tags... > create_from_options(const tuples::TaggedTuple< OptionTags... > &options, tmpl::list< Tags... >) noexcept
Given a list of tags and a tagged tuple containing items created from input options, return a tagged tuple of items constructed by calls to create_from_options for each tag in the list.
Definition: CreateFromOptions.hpp:40
Defines class DenseVector.
Terminate the algorithm to proceed to the next phase.
Definition: TerminatePhase.hpp:26
Wraps the template metaprogramming library used (brigand)
Declares functions to do file system manipulations.
Definition: LinearSolverAlgorithmTestHelpers.hpp:196
List of all the actions to be executed in the specified phase.
Definition: PhaseDependentActionList.hpp:16
Defines functions and classes from the GSL.
Definition: LinearSolverAlgorithmTestHelpers.hpp:107
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
Code to wrap or improve the Catch testing framework used for unit tests.
The linear operator applied to the data in Tag
Definition: Tags.hpp:66
Definition: Test_ActionTesting.cpp:365
Defines class template ConstGlobalCache.
Defines macro ERROR.
Defines the Charm++ mainchare.
Defines DataBox tags for the linear solver.
Prefix indicating a source term that is independent of dynamic variables.
Definition: Prefixes.hpp:100
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182
Definition: LinearSolverAlgorithmTestHelpers.hpp:168