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 
16 #include "DataStructures/DataBox/PrefixHelpers.hpp"
18 #include "DataStructures/DataBox/Tag.hpp"
19 #include "DataStructures/DataVector.hpp"
24 #include "Domain/CreateInitialElement.hpp"
26 #include "Domain/Structure/CreateInitialMesh.hpp"
28 #include "Elliptic/DiscontinuousGalerkin/DgElementArray.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"
35 #include "Parallel/Actions/SetupDataBox.hpp"
36 #include "Parallel/Actions/TerminatePhase.hpp"
38 #include "Parallel/GlobalCache.hpp"
39 #include "Parallel/InitializationFunctions.hpp"
40 #include "Parallel/Invoke.hpp"
41 #include "Parallel/Main.hpp"
42 #include "Parallel/ParallelComponentHelpers.hpp"
43 #include "Parallel/PhaseDependentActionList.hpp"
44 #include "Parallel/Reduction.hpp"
45 #include "ParallelAlgorithms/Initialization/Actions/RemoveOptionsAndTerminatePhase.hpp"
46 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
47 #include "ParallelAlgorithms/LinearSolver/Actions/MakeIdentityIfSkipped.hpp"
49 #include "Utilities/Blas.hpp"
52 #include "Utilities/Gsl.hpp"
53 #include "Utilities/Requires.hpp"
54 #include "Utilities/TMPL.hpp"
55 
56 namespace helpers = LinearSolverAlgorithmTestHelpers;
57 
58 /// Functionality to test parallel linear solvers on multiple elements
60 
61 namespace OptionTags {
62 // This option expects a list of N matrices that each have N*M rows and M
63 // columns, where N is the number of elements and M is a nonzero integer.
64 // Therefore, this option specifies a (N*M,N*M) matrix that has its columns
65 // split over all elements. In a context where the linear operator represents a
66 // DG discretization, M is the number of collocation points per element.
68  static constexpr Options::String help = "The linear operator A to invert.";
70 };
71 // Both of the following options expect a list of N vectors that have a size of
72 // M each, so that they constitute a vector of total size N*M (see above).
73 struct Source {
74  static constexpr Options::String help = "The source b in the equation Ax=b.";
76 };
78  static constexpr Options::String help = "The solution x in the equation Ax=b";
80 };
81 } // namespace OptionTags
82 
83 // [array_allocation_tag]
86  using option_tags = tmpl::list<OptionTags::LinearOperator>;
87 
88  static constexpr bool pass_metavariables = false;
90  create_from_options(
92  linear_operator) noexcept {
93  return linear_operator;
94  }
95 };
96 // [array_allocation_tag]
97 
100  using option_tags = tmpl::list<OptionTags::Source>;
101 
102  static constexpr bool pass_metavariables = false;
103  static std::vector<DenseVector<double>> create_from_options(
104  const std::vector<DenseVector<double>>& source) noexcept {
105  return source;
106  }
107 };
108 
111  using option_tags = tmpl::list<OptionTags::ExpectedResult>;
112 
113  static constexpr bool pass_metavariables = false;
114  static std::vector<DenseVector<double>> create_from_options(
115  const std::vector<DenseVector<double>>& expected_result) noexcept {
116  return expected_result;
117  }
118 };
119 
120 // The vector `x` we want to solve for
122  using type = Scalar<DataVector>;
123  static std::string name() noexcept { return "ScalarField"; }
124 };
125 
128 using operator_applied_to_fields_tag =
130 
131 // We are working on a one-dimension domain where each element corresponds to an
132 // entry of the vectors provided in the input file. This function translates the
133 // element to an index into these vectors. We assume the domain is composed of a
134 // single block.
135 size_t get_index(const ElementId<1>& element_id) noexcept {
136  return element_id.segment_id(0).index();
137 }
138 
139 // In the following `ComputeOperatorAction` and `CollectOperatorAction` actions
140 // we compute A(p)=sum_elements(A_element(p_element)) in a global reduction and
141 // then broadcast the global A(p) back to the elements so that they can extract
142 // their A_element(p). This is horribly inefficient parallelism but allows us to
143 // just provide a global matrix A (represented by the `LinearOperator` tag) in
144 // an input file.
145 
146 // Forward declare to keep these actions in the order they are used
147 template <typename OperandTag>
149 
150 template <typename OperandTag>
152  using const_global_cache_tags = tmpl::list<LinearOperator>;
153 
154  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
155  typename ActionList, typename ParallelComponent>
156  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
157  db::DataBox<DbTagsList>& box,
158  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
160  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
161  const ElementId<1>& element_id,
162  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
163  const ActionList /*meta*/,
164  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
165  const ParallelComponent* const /*meta*/) noexcept {
166  const size_t element_index = get_index(element_id);
167  const auto& operator_matrices = get<LinearOperator>(box);
168  const auto number_of_elements = operator_matrices.size();
169  const auto& linear_operator = gsl::at(operator_matrices, element_index);
170  const auto number_of_grid_points = linear_operator.columns();
171  const auto& operand = get<OperandTag>(box);
172 
173  typename OperandTag::type operator_applied_to_operand{
174  number_of_grid_points * number_of_elements};
175  dgemv_('N', linear_operator.rows(), linear_operator.columns(), 1,
176  linear_operator.data(), linear_operator.spacing(), operand.data(), 1,
177  0, operator_applied_to_operand.data(), 1);
178 
179  Parallel::contribute_to_reduction<CollectOperatorAction<OperandTag>>(
180  Parallel::ReductionData<
182  operator_applied_to_operand},
183  Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
184  Parallel::get_parallel_component<ParallelComponent>(cache));
185 
186  // Terminate algorithm for now. The reduction will be broadcast to the
187  // next action which is responsible for restarting the algorithm.
188  return {std::move(box), true};
189  }
190 };
191 
192 template <typename OperandTag>
193 struct CollectOperatorAction {
194  using local_operator_applied_to_operand_tag =
196 
197  template <typename ParallelComponent, typename DbTagsList,
198  typename Metavariables, typename ScalarFieldOperandTag,
199  Requires<tmpl::list_contains_v<
200  DbTagsList, local_operator_applied_to_operand_tag>> = nullptr>
201  static void apply(db::DataBox<DbTagsList>& box,
203  const ElementId<1>& element_id,
204  const Variables<tmpl::list<ScalarFieldOperandTag>>&
205  Ap_global_data) noexcept {
206  const size_t element_index = get_index(element_id);
207  // This could be generalized to work on the Variables instead of the
208  // Scalar, but it's only for the purpose of this test.
209  const auto number_of_grid_points = get<LinearOperator>(box)[0].columns();
210  const auto& Ap_global = get<ScalarFieldOperandTag>(Ap_global_data).get();
211  DataVector Ap_local{number_of_grid_points};
212  std::copy(Ap_global.begin() +
213  static_cast<int>(element_index * number_of_grid_points),
214  Ap_global.begin() +
215  static_cast<int>((element_index + 1) * number_of_grid_points),
216  Ap_local.begin());
217  db::mutate<local_operator_applied_to_operand_tag>(
218  make_not_null(&box),
219  [&Ap_local, &number_of_grid_points](auto Ap) noexcept {
220  *Ap = typename local_operator_applied_to_operand_tag::type{
221  number_of_grid_points};
223  *Ap)) = Ap_local;
224  });
225  // Proceed with algorithm
226  Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
227  .perform_algorithm(true);
228  }
229 };
230 
231 // Checks for the correct solution after the algorithm has terminated.
232 template <typename OptionsGroup>
233 struct TestResult {
234  using const_global_cache_tags =
235  tmpl::list<ExpectedResult, helpers::ExpectedConvergenceReason>;
236 
237  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
238  typename ActionList, typename ParallelComponent>
239  static std::tuple<db::DataBox<DbTagsList>&&, bool> apply(
240  db::DataBox<DbTagsList>& box,
241  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
242  const Parallel::GlobalCache<Metavariables>& /*cache*/,
243  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
244  const ElementId<1>& element_id,
245  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
246  const ActionList /*meta*/,
247  // NOLINTNEXTLINE(readability-avoid-const-params-in-decls)
248  const ParallelComponent* const /*meta*/) noexcept {
249  const size_t element_index = get_index(element_id);
250  const auto& has_converged =
251  get<Convergence::Tags::HasConverged<OptionsGroup>>(box);
252  SPECTRE_PARALLEL_REQUIRE(has_converged);
253  SPECTRE_PARALLEL_REQUIRE(has_converged.reason() ==
254  get<helpers::ExpectedConvergenceReason>(box));
255  const auto& expected_result =
256  gsl::at(get<ExpectedResult>(box), element_index);
257  const auto& result = get<ScalarFieldTag>(box).get();
258  for (size_t i = 0; i < expected_result.size(); i++) {
259  SPECTRE_PARALLEL_REQUIRE(result[i] == approx(expected_result[i]));
260  }
261  return {std::move(box), true};
262  }
263 };
264 
266  using initialization_tags =
267  tmpl::list<domain::Tags::InitialExtents<1>,
269  using const_global_cache_tags = tmpl::list<Source>;
270  using simple_tags =
271  tmpl::list<domain::Tags::Mesh<1>, domain::Tags::Element<1>,
273  sources_tag>;
274  using compute_tags = tmpl::list<>;
275  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
276  typename ActionList, typename ParallelComponent>
277  static auto apply(db::DataBox<DbTagsList>& box,
278  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
279  const Parallel::GlobalCache<Metavariables>& /*cache*/,
280  const ElementId<1>& element_id, const ActionList /*meta*/,
281  const ParallelComponent* const /*meta*/) noexcept {
282  // Domain geometry
283  const auto& domain = db::get<domain::Tags::Domain<1>>(box);
284  const auto& initial_extents = db::get<domain::Tags::InitialExtents<1>>(box);
285  const auto& initial_refinement =
286  db::get<domain::Tags::InitialRefinementLevels<1>>(box);
288  initial_extents, element_id, Spectral::Quadrature::GaussLobatto);
289  const auto& block = domain.blocks()[element_id.block_id()];
291  element_id, block, initial_refinement);
292  auto logical_coords = logical_coordinates(mesh);
293  // Element data
294  const size_t element_index = get_index(element_id);
295  const auto& source = gsl::at(get<Source>(box), element_index);
296  const size_t num_points = source.size();
297  ::Initialization::mutate_assign<simple_tags>(
298  make_not_null(&box), std::move(mesh), std::move(element),
299  std::move(logical_coords), typename fields_tag::type{num_points, 0.},
300  typename sources_tag::type{source});
301  return std::make_tuple(std::move(box));
302  }
303 };
304 
305 namespace detail {
306 
307 template <typename Preconditioner>
308 struct run_preconditioner_impl {
309  using type =
310  tmpl::list<ComputeOperatorAction<typename Preconditioner::fields_tag>,
311  typename Preconditioner::template solve<ComputeOperatorAction<
312  typename Preconditioner::operand_tag>>,
314 };
315 
316 template <>
317 struct run_preconditioner_impl<void> {
318  using type = tmpl::list<>;
319 };
320 
321 template <typename Preconditioner>
322 using run_preconditioner =
323  typename run_preconditioner_impl<Preconditioner>::type;
324 
325 } // namespace detail
326 
327 template <typename Metavariables,
328  typename LinearSolverType = typename Metavariables::linear_solver,
329  typename PreconditionerType = typename Metavariables::preconditioner>
330 using initialization_actions =
331  tmpl::list<Actions::SetupDataBox, InitializeElement,
332  typename LinearSolverType::initialize_element,
333  ComputeOperatorAction<fields_tag>,
334  helpers::detail::init_preconditioner<PreconditionerType>,
336 
337 template <typename Metavariables,
338  typename LinearSolverType = typename Metavariables::linear_solver,
339  typename PreconditionerType = typename Metavariables::preconditioner>
340 using register_actions =
341  tmpl::list<typename LinearSolverType::register_element,
342  helpers::detail::register_preconditioner<PreconditionerType>,
344 
345 template <typename Metavariables,
346  typename LinearSolverType = typename Metavariables::linear_solver,
347  typename PreconditionerType = typename Metavariables::preconditioner>
348 using solve_actions = tmpl::list<
349  typename LinearSolverType::template solve<tmpl::list<
350  detail::run_preconditioner<PreconditionerType>,
351  ComputeOperatorAction<typename LinearSolverType::operand_tag>>>,
353 
354 template <typename Metavariables,
355  typename LinearSolverType = typename Metavariables::linear_solver,
356  typename PreconditionerType = typename Metavariables::preconditioner>
357 using test_actions =
358  tmpl::list<TestResult<typename LinearSolverType::options_group>>;
359 
360 template <typename Metavariables>
361 using ElementArray = elliptic::DgElementArray<
362  Metavariables,
363  tmpl::list<
364  Parallel::PhaseActions<typename Metavariables::Phase,
365  Metavariables::Phase::Initialization,
366  initialization_actions<Metavariables>>,
367  Parallel::PhaseActions<typename Metavariables::Phase,
368  Metavariables::Phase::RegisterWithObserver,
369  register_actions<Metavariables>>,
370  Parallel::PhaseActions<typename Metavariables::Phase,
371  Metavariables::Phase::PerformLinearSolve,
372  solve_actions<Metavariables>>,
373  Parallel::PhaseActions<typename Metavariables::Phase,
374  Metavariables::Phase::TestResult,
375  test_actions<Metavariables>>>>;
376 
377 template <typename Metavariables>
378 using component_list =
379  tmpl::push_back<tmpl::append<helpers::detail::get_component_list<
380  typename Metavariables::linear_solver>,
381  helpers::detail::get_component_list<
382  typename Metavariables::preconditioner>>,
383  ElementArray<Metavariables>,
387 
388 } // namespace DistributedLinearSolverAlgorithmTestHelpers
FloatingPointExceptions.hpp
domain::Initialization::create_initial_element
Element< VolumeDim > create_initial_element(const ElementId< VolumeDim > &element_id, const Block< VolumeDim > &block, const std::vector< std::array< size_t, VolumeDim >> &initial_refinement_levels) noexcept
Creates an initial element of a Block.
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
Actions::SetupDataBox
Add into the DataBox default constructed items for the collection of tags requested by any of the act...
Definition: SetupDataBox.hpp:102
DistributedLinearSolverAlgorithmTestHelpers::TestResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:233
std::string
Main.hpp
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
domain::Tags::Coordinates
Definition: Tags.hpp:130
Parallel::ReductionDatum
The data to be reduced, and invokables to be called whenever two reduction messages are combined and ...
Definition: Reduction.hpp:65
LinearSolver::Tags::OperatorAppliedTo
The linear operator applied to the data in Tag
Definition: Tags.hpp:53
GlobalCache.hpp
vector
Error.hpp
domain::Tags::Element
Definition: Tags.hpp:97
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
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
DistributedLinearSolverAlgorithmTestHelpers::ExpectedResult
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:109
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:67
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
DistributedLinearSolverAlgorithmTestHelpers::OptionTags::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:73
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:59
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
ElementId.hpp
DistributedLinearSolverAlgorithmTestHelpers::ComputeOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:151
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
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:121
array
LogicalCoordinates.hpp
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:46
LinearSolverAlgorithmTestHelpers::OutputCleaner
Definition: LinearSolverAlgorithmTestHelpers.hpp:351
DistributedLinearSolverAlgorithmTestHelpers::InitializeElement
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:265
Variables.hpp
AlgorithmArray.hpp
DistributedLinearSolverAlgorithmTestHelpers::LinearOperator
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:84
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
logical_coordinates
void logical_coordinates(gsl::not_null< tnsr::I< DataVector, VolumeDim, Frame::Logical > * > logical_coords, const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
domain::Initialization::create_initial_mesh
Mesh< Dim > create_initial_mesh(const std::vector< std::array< size_t, Dim >> &initial_extents, const ElementId< Dim > &element_id, Spectral::Quadrature quadrature, const OrientationMap< Dim > &orientation={}) noexcept
Construct the initial Mesh of an Element.
DistributedLinearSolverAlgorithmTestHelpers::Source
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:98
DistributedLinearSolverAlgorithmTestHelpers::CollectOperatorAction
Definition: DistributedLinearSolverAlgorithmTestHelpers.hpp:148
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
domain::Tags::InitialRefinementLevels
Definition: Tags.hpp:84
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:77
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
Mesh.hpp
string