ResidualMonitorActions.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
10 #include "Informer/Tags.hpp"
11 #include "Informer/Verbosity.hpp"
13 #include "NumericalAlgorithms/LinearSolver/Observe.hpp"
16 #include "Parallel/Info.hpp"
17 #include "Parallel/Invoke.hpp"
18 #include "Parallel/Printf.hpp"
19 #include "Utilities/EqualWithinRoundoff.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/Requires.hpp"
22 
23 /// \cond
24 namespace tuples {
25 template <typename...>
26 class TaggedTuple;
27 } // namespace tuples
28 namespace LinearSolver {
29 namespace gmres_detail {
30 struct NormalizeInitialOperand;
31 struct OrthogonalizeOperand;
32 struct NormalizeOperandAndUpdateField;
33 } // namespace gmres_detail
34 } // namespace LinearSolver
35 /// \endcond
36 
37 namespace LinearSolver {
38 namespace gmres_detail {
39 
40 template <typename BroadcastTarget>
41 struct InitializeResidualMagnitude {
42  template <typename... DbTags, typename... InboxTags, typename Metavariables,
43  typename ArrayIndex, typename ActionList,
44  typename ParallelComponent,
45  Requires<sizeof...(DbTags) != 0> = nullptr>
46  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
49  const ArrayIndex& /*array_index*/,
50  const ActionList /*meta*/,
51  const ParallelComponent* const /*meta*/,
52  const double residual_magnitude) noexcept {
53  using fields_tag = typename Metavariables::system::fields_tag;
54  using residual_magnitude_tag = db::add_tag_prefix<
57  using initial_residual_magnitude_tag =
59 
60  db::mutate<residual_magnitude_tag, initial_residual_magnitude_tag>(
61  make_not_null(&box), [residual_magnitude](
63  local_residual_magnitude,
65  initial_residual_magnitude) noexcept {
66  *local_residual_magnitude = *initial_residual_magnitude =
67  residual_magnitude;
68  });
69 
70  LinearSolver::observe_detail::contribute_to_reduction_observer(box, cache);
71 
72  // Determine whether the linear solver has already converged. This invokes
73  // the compute item.
74  const auto& has_converged = db::get<LinearSolver::Tags::HasConverged>(box);
75 
76  if (UNLIKELY(has_converged and
77  static_cast<int>(get<::Tags::Verbosity>(box)) >=
78  static_cast<int>(::Verbosity::Quiet))) {
79  Parallel::printf("%s", has_converged);
80  }
81 
82  Parallel::simple_action<NormalizeInitialOperand>(
83  Parallel::get_parallel_component<BroadcastTarget>(cache),
84  residual_magnitude, has_converged);
85  }
86 };
87 
88 template <typename BroadcastTarget>
89 struct StoreOrthogonalization {
90  template <typename... DbTags, typename... InboxTags, typename Metavariables,
91  typename ArrayIndex, typename ActionList,
92  typename ParallelComponent,
93  Requires<sizeof...(DbTags) != 0> = nullptr>
94  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
97  const ArrayIndex& /*array_index*/,
98  const ActionList /*meta*/,
99  const ParallelComponent* const /*meta*/,
100  const double orthogonalization) noexcept {
101  using fields_tag = typename Metavariables::system::fields_tag;
102  using orthogonalization_iteration_id_tag =
105  using orthogonalization_history_tag =
107  fields_tag>;
108 
109  db::mutate<orthogonalization_history_tag,
110  orthogonalization_iteration_id_tag>(
111  make_not_null(&box),
112  [orthogonalization](
114  orthogonalization_history,
115  const gsl::not_null<IterationId*> orthogonalization_iteration_id,
116  const IterationId& iteration_id) noexcept {
117  (*orthogonalization_history)(
118  orthogonalization_iteration_id->step_number,
119  iteration_id.step_number) = orthogonalization;
120  orthogonalization_iteration_id->step_number++;
121  },
122  get<LinearSolver::Tags::IterationId>(box));
123 
124  Parallel::simple_action<OrthogonalizeOperand>(
125  Parallel::get_parallel_component<BroadcastTarget>(cache),
126  orthogonalization);
127  }
128 };
129 
130 template <typename BroadcastTarget>
131 struct StoreFinalOrthogonalization {
132  template <typename... DbTags, typename... InboxTags, typename Metavariables,
133  typename ArrayIndex, typename ActionList,
134  typename ParallelComponent,
135  Requires<sizeof...(DbTags) != 0> = nullptr>
136  static void apply(db::DataBox<tmpl::list<DbTags...>>& box,
139  const ArrayIndex& /*array_index*/,
140  const ActionList /*meta*/,
141  const ParallelComponent* const /*meta*/,
142  const double orthogonalization) noexcept {
143  using fields_tag = typename Metavariables::system::fields_tag;
144  using residual_magnitude_tag = db::add_tag_prefix<
145  LinearSolver::Tags::Magnitude,
146  db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
147  using initial_residual_magnitude_tag =
149  using orthogonalization_iteration_id_tag =
152  using orthogonalization_history_tag =
154  fields_tag>;
155 
156  db::mutate<orthogonalization_history_tag>(
157  make_not_null(&box),
158  [orthogonalization](
160  orthogonalization_history,
161  const IterationId& iteration_id,
162  const IterationId& orthogonalization_iteration_id) noexcept {
163  (*orthogonalization_history)(
164  orthogonalization_iteration_id.step_number,
165  iteration_id.step_number) = sqrt(orthogonalization);
166  },
167  get<LinearSolver::Tags::IterationId>(box),
168  get<orthogonalization_iteration_id_tag>(box));
169 
170  // Perform a QR decomposition of the Hessenberg matrix that was built during
171  // the orthogonalization
172  const auto& orthogonalization_history =
173  get<orthogonalization_history_tag>(box);
174  const auto num_rows =
175  get<orthogonalization_iteration_id_tag>(box).step_number + 1;
176  DenseMatrix<double> qr_Q;
177  DenseMatrix<double> qr_R;
178  blaze::qr(orthogonalization_history, qr_Q, qr_R);
179  // Compute the residual vector from the QR decomposition
180  DenseVector<double> beta(num_rows, 0.);
181  beta[0] = get<initial_residual_magnitude_tag>(box);
182  const DenseVector<double> minres =
183  blaze::inv(qr_R) * blaze::trans(qr_Q) * beta;
184  const double residual_magnitude =
185  blaze::length(beta - orthogonalization_history * minres);
186 
187  // Store residual magnitude and prepare for the next iteration
188  db::mutate<residual_magnitude_tag, LinearSolver::Tags::IterationId,
189  orthogonalization_iteration_id_tag,
190  orthogonalization_history_tag>(
191  make_not_null(&box),
192  [residual_magnitude](
193  const gsl::not_null<double*> local_residual_magnitude,
194  const gsl::not_null<IterationId*> iteration_id,
195  const gsl::not_null<IterationId*> orthogonalization_iteration_id,
197  local_orthogonalization_history) noexcept {
198  *local_residual_magnitude = residual_magnitude;
199  // Prepare for the next iteration
200  iteration_id->step_number++;
201  orthogonalization_iteration_id->step_number = 0;
202  local_orthogonalization_history->resize(
203  iteration_id->step_number + 2, iteration_id->step_number + 1);
204  // Make sure the new entries are zero
205  for (size_t i = 0; i < local_orthogonalization_history->rows(); i++) {
206  (*local_orthogonalization_history)(
207  i, local_orthogonalization_history->columns() - 1) = 0.;
208  }
209  for (size_t j = 0; j < local_orthogonalization_history->columns();
210  j++) {
211  (*local_orthogonalization_history)(
212  local_orthogonalization_history->rows() - 1, j) = 0.;
213  }
214  });
215 
216  // At this point, the iteration is complete. We proceed with observing,
217  // logging and checking convergence before broadcasting back to the
218  // elements.
219 
220  LinearSolver::observe_detail::contribute_to_reduction_observer(box, cache);
221 
222  // Determine whether the linear solver has converged. This invokes the
223  // compute item.
224  const auto& has_converged = db::get<LinearSolver::Tags::HasConverged>(box);
225 
226  // Do some logging
227  if (UNLIKELY(static_cast<int>(get<::Tags::Verbosity>(box)) >=
228  static_cast<int>(::Verbosity::Verbose))) {
230  "Linear solver iteration %d done. Remaining residual: %e\n",
231  get<LinearSolver::Tags::IterationId>(box).step_number,
232  residual_magnitude);
233  }
234  if (UNLIKELY(has_converged and
235  static_cast<int>(get<::Tags::Verbosity>(box)) >=
236  static_cast<int>(::Verbosity::Quiet))) {
237  Parallel::printf("%s", has_converged);
238  }
239 
240  Parallel::simple_action<NormalizeOperandAndUpdateField>(
241  Parallel::get_parallel_component<BroadcastTarget>(cache),
242  sqrt(orthogonalization), minres, has_converged);
243  }
244 };
245 
246 } // namespace gmres_detail
247 } // namespace LinearSolver
The magnitude w.r.t. the LinearSolver::inner_product
Definition: Tags.hpp:115
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Wrap Tag in Prefix<_, Args...>, also wrapping variables tags if Tag is a Tags::Variables.
Definition: DataBoxTag.hpp:533
A Hessenberg matrix built up during an orthogonalization procedure.
Definition: Tags.hpp:156
Definition: Variables.hpp:46
Defines functions for interfacing with the parallelization framework.
Defines DataBox tags for the linear solver.
void mutate(const gsl::not_null< DataBox< TagList > *> box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1099
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
#define UNLIKELY(x)
Definition: Gsl.hpp:72
Definition: TaggedTuple.hpp:25
Define prefixes for DataBox tags.
Defines the type alias Requires.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
A dynamically sized matrix of arbitrary type.
Definition: DenseMatrix.hpp:29
Holds an IterationId that identifies a step in the linear solver algorithm.
Definition: Tags.hpp:66
T beta(T... args)
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
Defines class DenseMatrix.
Defines Parallel::printf for writing to stdout.
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
The prefix for tags related to an orthogonalization procedurce.
Definition: Tags.hpp:143
Defines class DenseVector.
Defines functions and classes from the GSL.
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:863
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
T sqrt(T... args)
Definition: SolvePoissonProblem.hpp:38
Defines class template ConstGlobalCache.
void printf(const std::string &format, Args &&... args)
Print an atomic message to stdout with C printf usage.
Definition: Printf.hpp:100
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
Defines class IterationId.