DistributedLinearSolverAlgorithmTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <algorithm>
9 #include <array>
10 #include <cstddef>
11 #include <string>
12 #include <tuple>
13 #include <vector>
14 
15 #include "AlgorithmArray.hpp"
17 #include "DataStructures/DataBox/PrefixHelpers.hpp"
19 #include "DataStructures/DataBox/Tag.hpp"
20 #include "DataStructures/DataVector.hpp"
26 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.hpp"
27 #include "ErrorHandling/Error.hpp"
29 #include "Helpers/ParallelAlgorithms/LinearSolver/LinearSolverAlgorithmTestHelpers.hpp"
30 #include "IO/Observer/Actions/RegisterWithObservers.hpp"
31 #include "IO/Observer/ObserverComponent.hpp"
32 #include "IO/Observer/Tags.hpp"
33 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
34 #include "Parallel/Actions/SetupDataBox.hpp"
35 #include "Parallel/Actions/TerminatePhase.hpp"
36 #include "Parallel/GlobalCache.hpp"
37 #include "Parallel/Info.hpp"
38 #include "Parallel/InitializationFunctions.hpp"
39 #include "Parallel/Invoke.hpp"
40 #include "Parallel/Main.hpp"
41 #include "Parallel/ParallelComponentHelpers.hpp"
42 #include "Parallel/PhaseDependentActionList.hpp"
43 #include "Parallel/Reduction.hpp"
44 #include "ParallelAlgorithms/DiscontinuousGalerkin/InitializeDomain.hpp"
45 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
46 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
48 #include "Utilities/Blas.hpp"
49 #include "Utilities/Gsl.hpp"
50 #include "Utilities/Requires.hpp"
51 #include "Utilities/TMPL.hpp"
52 
53 namespace helpers = LinearSolverAlgorithmTestHelpers;
54 
55 /// Functionality to test parallel linear solvers on multiple elements
57 
58 namespace OptionTags {
59 // This option expects a list of N matrices that each have N*M rows and M
60 // columns, where N is the number of elements and M is a nonzero integer.
61 // Therefore, this option specifies a (N*M,N*M) matrix that has its columns
62 // split over all elements. In a context where the linear operator represents a
63 // DG discretization, M is the number of collocation points per element.
65  static constexpr Options::String help = "The linear operator A to invert.";
67 };
68 // Both of the following options expect a list of N vectors that have a size of
69 // M each, so that they constitute a vector of total size N*M (see above).
70 struct Source {
71  static constexpr Options::String help = "The source b in the equation Ax=b.";
73 };
75  static constexpr Options::String help = "The solution x in the equation Ax=b";
77 };
78 } // namespace OptionTags
79 
80 // [array_allocation_tag]
83  using option_tags = tmpl::list<OptionTags::LinearOperator>;
84 
85  static constexpr bool pass_metavariables = false;
87  create_from_options(
89  linear_operator) noexcept {
90  return linear_operator;
91  }
92 };
93 // [array_allocation_tag]
94 
97  using option_tags = tmpl::list<OptionTags::Source>;
98 
99  static constexpr bool pass_metavariables = false;
100  static std::vector<DenseVector<double>> create_from_options(
101  const std::vector<DenseVector<double>>& source) noexcept {
102  return source;
103  }
104 };
105 
108  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
109 
110  static constexpr bool pass_metavariables = false;
111  static std::vector<DenseVector<double>> create_from_options(
112  const std::vector<DenseVector<double>>& expected_result) noexcept {
113  return expected_result;
114  }
115 };
116 
117 // The vector `x` we want to solve for
119  using type = Scalar<DataVector>;
120  static std::string name() noexcept { return "ScalarField"; }
121 };
122 
125 using operator_applied_to_fields_tag =
127 
128 // We are working on a one-dimension domain where each element corresponds to an
129 // entry of the vectors provided in the input file. This function translates the
130 // element to an index into these vectors. We assume the domain is composed of a
131 // single block.
132 size_t get_index(const ElementId<1>& element_id) noexcept {
133  return element_id.segment_ids()[0].index();
134 }
135 
136 // In the following `ComputeOperatorAction` and `CollectOperatorAction` actions
137 // we compute A(p)=sum_elements(A_element(p_element)) in a global reduction and
138 // then broadcast the global A(p) back to the elements so that they can extract
139 // their A_element(p). This is horribly inefficient parallelism but allows us to
140 // just provide a global matrix A (represented by the `LinearOperator` tag) in
141 // an input file.
142 
143 // Forward declare to keep these actions in the order they are used
144 template <typename OperandTag>
146 
147 template <typename OperandTag>
149  using const_global_cache_tags = tmpl::list<LinearOperator>;
150 
151  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
152  typename ActionList, typename ParallelComponent>
153  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
154  db::DataBox<DbTagsList>& box,
155  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
157  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
158  const ElementId<1>& element_id,
159  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
160  const ActionList /*meta*/,
161  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
162  const ParallelComponent* const /*meta*/) noexcept {
163  const size_t element_index = get_index(element_id);
164  const auto& operator_matrices = get<LinearOperator>(box);
165  const auto number_of_elements = operator_matrices.size();
166  const auto& linear_operator = gsl::at(operator_matrices, element_index);
167  const auto number_of_grid_points = linear_operator.columns();
168  const auto& operand = get<OperandTag>(box);
169 
170  typename OperandTag::type operator_applied_to_operand{
171  number_of_grid_points * number_of_elements};
172  dgemv_('N', linear_operator.rows(), linear_operator.columns(), 1,
173  linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
174  0, operator_applied_to_operand.data(), 1);
175 
176  Parallel::contribute_to_reduction<CollectOperatorAction<OperandTag>>(
177  Parallel::ReductionData<
179  operator_applied_to_operand},
180  Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
181  Parallel::get_parallel_component<ParallelComponent>(cache));
182 
183  // Terminate algorithm for now. The reduction will be broadcast to the
184  // next action which is responsible for restarting the algorithm.
185  return {std::move(box), true};
186  }
187 };
188 
189 template <typename OperandTag>
190 struct CollectOperatorAction {
191  using local_operator_applied_to_operand_tag =
193 
194  template <typename ParallelComponent, typename DbTagsList,
195  typename Metavariables, typename ScalarFieldOperandTag,
196  Requires<tmpl::list_contains_v<
197  DbTagsList, local_operator_applied_to_operand_tag>> = nullptr>
198  static void apply(db::DataBox<DbTagsList>& box,
200  const ElementId<1>& element_id,
201  const Variables<tmpl::list<ScalarFieldOperandTag>>&
202  Ap_global_data) noexcept {
203  const size_t element_index = get_index(element_id);
204  // This could be generalized to work on the Variables instead of the
205  // Scalar, but it's only for the purpose of this test.
206  const auto number_of_grid_points = get<LinearOperator>(box)[0].columns();
207  const auto& Ap_global = get<ScalarFieldOperandTag>(Ap_global_data).get();
208  DataVector Ap_local{number_of_grid_points};
209  std::copy(Ap_global.begin() +
210  static_cast<int>(element_index * number_of_grid_points),
211  Ap_global.begin() +
212  static_cast<int>((element_index + 1) * number_of_grid_points),
213  Ap_local.begin());
214  db::mutate<local_operator_applied_to_operand_tag>(
215  make_not_null(&box),
216  [&Ap_local, &number_of_grid_points](auto Ap) noexcept {
217  *Ap = typename local_operator_applied_to_operand_tag::type{
218  number_of_grid_points};
220  *Ap)) = Ap_local;
221  });
222  // Proceed with algorithm
223  Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
224  .perform_algorithm(true);
225  }
226 };
227 
228 // Checks for the correct solution after the algorithm has terminated.
229 template <typename OptionsGroup>
230 struct TestResult {
231  using const_global_cache_tags =
232  tmpl::list<ExpectedResult, helpers::ExpectedConvergenceReason>;
233 
234  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
235  typename ActionList, typename ParallelComponent>
236  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
237  db::DataBox<DbTagsList>& box,
238  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
239  const Parallel::GlobalCache<Metavariables>& /*cache*/,
240  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
241  const ElementId<1>& element_id,
242  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
243  const ActionList /*meta*/,
244  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
245  const ParallelComponent* const /*meta*/) noexcept {
246  const size_t element_index = get_index(element_id);
247  const auto& has_converged =
248  get<Convergence::Tags::HasConverged<OptionsGroup>>(box);
249  SPECTRE_PARALLEL_REQUIRE(has_converged);
250  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
251  get<helpers::ExpectedConvergenceReason>(box));
252  const auto& expected_result =
253  gsl::at(get<ExpectedResult>(box), element_index);
254  const auto& result = get<ScalarFieldTag>(box).get();
255  for (size_t i = 0; i < expected_result.size(); i++) {
256  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
257  }
258  return {std::move(box), true};
259  }
260 };
261 
263  using const_global_cache_tags = tmpl::list<Source>;
264 
265  using simple_tags = tmpl::list<fields_tag, sources_tag>;
266  using compute_tags = tmpl::list<>;
267  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
268  typename ActionList, typename ParallelComponent>
269  static auto apply(db::DataBox<DbTagsList>& box,
270  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
271  const Parallel::GlobalCache<Metavariables>& /*cache*/,
272  const ElementId<1>& element_id, const ActionList /*meta*/,
273  const ParallelComponent* const /*meta*/) noexcept {
274  const size_t element_index = get_index(element_id);
275  const auto& source = gsl::at(get<Source>(box), element_index);
276  const size_t num_points = source.size();
277  ::Initialization::mutate_assign<simple_tags>(
278  make_not_null(&box), typename fields_tag::type{num_points, 0.},
279  typename sources_tag::type{source});
280 
281  return std::make_tuple(std::move(box));
282  }
283 };
284 
285 namespace detail {
286 
287 template <typename Preconditioner>
288 struct run_preconditioner_impl {
289  using type =
290  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
291  typename Preconditioner::template solve<ComputeOperatorAction<
292  typename Preconditioner::operand_tag>>>;
293 };
294 
295 template <>
296 struct run_preconditioner_impl<void> {
297  using type = tmpl::list<>;
298 };
299 
300 template <typename Preconditioner>
301 using run_preconditioner =
302  typename run_preconditioner_impl<Preconditioner>::type;
303 
304 } // namespace detail
305 
306 template <typename Metavariables,
307  typename LinearSolverType = typename Metavariables::linear_solver,
308  typename PreconditionerType = typename Metavariables::preconditioner>
309 using initialization_actions =
310  tmpl::list<Actions::SetupDataBox, dg::Actions::InitializeDomain<1>,
311  InitializeElement, typename LinearSolverType::initialize_element,
312  ComputeOperatorAction<fields_tag>,
313  helpers::detail::init_preconditioner<PreconditionerType>,
315 
316 template <typename Metavariables,
317  typename LinearSolverType = typename Metavariables::linear_solver,
318  typename PreconditionerType = typename Metavariables::preconditioner>
319 using register_actions =
320  tmpl::list<typename LinearSolverType::register_element,
321  helpers::detail::register_preconditioner<PreconditionerType>,
323 
324 template <typename Metavariables,
325  typename LinearSolverType = typename Metavariables::linear_solver,
326  typename PreconditionerType = typename Metavariables::preconditioner>
327 using solve_actions = tmpl::list<
328  typename LinearSolverType::template solve<tmpl::list<
329  detail::run_preconditioner<PreconditionerType>,
330  ComputeOperatorAction<typename LinearSolverType::operand_tag>>>,
332 
333 template <typename Metavariables,
334  typename LinearSolverType = typename Metavariables::linear_solver,
335  typename PreconditionerType = typename Metavariables::preconditioner>
336 using test_actions =
337  tmpl::list<TestResult<typename LinearSolverType::options_group>>;
338 
339 template <typename Metavariables>
340 using ElementArray = elliptic::DgElementArray<
341  Metavariables,
342  tmpl::list<
343  Parallel::PhaseActions<typename Metavariables::Phase,
344  Metavariables::Phase::Initialization,
345  initialization_actions<Metavariables>>,
346  Parallel::PhaseActions<typename Metavariables::Phase,
347  Metavariables::Phase::RegisterWithObserver,
348  register_actions<Metavariables>>,
349  Parallel::PhaseActions<typename Metavariables::Phase,
350  Metavariables::Phase::PerformLinearSolve,
351  solve_actions<Metavariables>>,
352  Parallel::PhaseActions<typename Metavariables::Phase,
353  Metavariables::Phase::TestResult,
354  test_actions<Metavariables>>>>;
355 
356 template <typename Metavariables>
357 using component_list =
358  tmpl::push_back<tmpl::append<helpers::detail::get_component_list<
359  typename Metavariables::linear_solver>,
360  helpers::detail::get_component_list<
361  typename Metavariables::preconditioner>>,
362  ElementArray<Metavariables>,
366 
367 } // namespace DistributedLinearSolverAlgorithmTestHelpers
FloatingPointExceptions.hpp
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
DistributedLinearSolverAlgorithmTestHelpers::TestResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:230
std::string
Main.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:63
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:54
GlobalCache.hpp
vector
Error.hpp
TestingFramework.hpp
DenseVector.hpp
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
DenseMatrix.hpp
Tags::Variables
Definition: VariablesTag.hpp:21
Initialization::Actions::RemoveOptionsAndTerminatePhase
Definition: RemoveOptionsAndTerminatePhase.hpp:27
tuple
DistributedLinearSolverAlgorithmTestHelpers::ExpectedResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:106
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:64
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 >
Info.hpp
algorithm
Tags.hpp
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:70
DenseMatrix< double, blaze::columnMajor >
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.
DistributedLinearSolverAlgorithmTestHelpers
Functionality to test parallel linear solvers on multiple elements.
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:56
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:49
ElementId.hpp
DistributedLinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:148
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
elliptic::DgElementArray
The parallel component responsible for managing the DG elements that compose the computational domain...
Definition: DgElementArray.hpp:37
DataBox.hpp
cstddef
DistributedLinearSolverAlgorithmTestHelpers::ScalarFieldTag
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:118
array
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
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:346
DistributedLinearSolverAlgorithmTestHelpers::InitializeElement
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:262
Variables.hpp
DistributedLinearSolverAlgorithmTestHelpers::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:81
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
DistributedLinearSolverAlgorithmTestHelpers::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:95
DistributedLinearSolverAlgorithmTestHelpers::CollectOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:145
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Tensor.hpp
Requires.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
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
Blas.hpp
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::ExpectedResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:74
Prefixes.hpp
TMPL.hpp
dgemv_
void dgemv_(const char &TRANS, const size_t &M, const size_t &N, const double &ALPHA, const double *A, const size_t &LDA, const double *X, const size_t &INCX, const double &BETA, double *Y, const size_t &INCY)
Perform a matrix-vector multiplication.
Definition: Blas.hpp:172
string