Reference documentation for deal.II version GIT c545eda070 2023-01-27 00:25:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_richardson.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_solver_richardson_h
17 #define dealii_solver_richardson_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/logstream.h>
24 
25 #include <deal.II/lac/solver.h>
27 
28 #include <limits>
29 
31 
63 template <class VectorType = Vector<double>>
64 class SolverRichardson : public SolverBase<VectorType>
65 {
66 public:
71  {
75  explicit AdditionalData(const double omega = 1,
76  const bool use_preconditioned_residual = false);
77 
81  double omega;
82 
87  };
88 
94  const AdditionalData & data = AdditionalData());
95 
101  const AdditionalData &data = AdditionalData());
102 
106  virtual ~SolverRichardson() override = default;
107 
111  template <typename MatrixType, typename PreconditionerType>
112  void
113  solve(const MatrixType & A,
114  VectorType & x,
115  const VectorType & b,
116  const PreconditionerType &preconditioner);
117 
121  template <typename MatrixType, typename PreconditionerType>
122  void
123  Tsolve(const MatrixType & A,
124  VectorType & x,
125  const VectorType & b,
126  const PreconditionerType &preconditioner);
127 
131  void
132  set_omega(const double om = 1.);
133 
139  virtual void
140  print_vectors(const unsigned int step,
141  const VectorType & x,
142  const VectorType & r,
143  const VectorType & d) const;
144 
145 protected:
152  virtual typename VectorType::value_type
153  criterion(const VectorType &r, const VectorType &d) const;
154 
159 };
160 
162 /*----------------- Implementation of the Richardson Method ------------------*/
163 
164 #ifndef DOXYGEN
165 
166 template <class VectorType>
168  const double omega,
169  const bool use_preconditioned_residual)
170  : omega(omega)
171  , use_preconditioned_residual(use_preconditioned_residual)
172 {}
173 
174 
175 template <class VectorType>
178  const AdditionalData & data)
179  : SolverBase<VectorType>(cn, mem)
180  , additional_data(data)
181 {}
182 
183 
184 
185 template <class VectorType>
187  const AdditionalData &data)
188  : SolverBase<VectorType>(cn)
189  , additional_data(data)
190 {}
191 
192 
193 
194 template <class VectorType>
195 template <typename MatrixType, typename PreconditionerType>
196 void
197 SolverRichardson<VectorType>::solve(const MatrixType & A,
198  VectorType & x,
199  const VectorType & b,
200  const PreconditionerType &preconditioner)
201 {
203 
204  double last_criterion = std::numeric_limits<double>::lowest();
205 
206  unsigned int iter = 0;
207 
208  // Memory allocation.
209  // 'Vr' holds the residual, 'Vd' the preconditioned residual
210  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
211  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
212 
213  VectorType &r = *Vr;
214  r.reinit(x);
215 
216  VectorType &d = *Vd;
217  d.reinit(x);
218 
219  LogStream::Prefix prefix("Richardson");
220 
221  // Main loop
222  while (conv == SolverControl::iterate)
223  {
224  // Do not use residual,
225  // but do it in 2 steps
226  A.vmult(r, x);
227  r.sadd(-1., 1., b);
228  preconditioner.vmult(d, r);
229 
230  // get the required norm of the (possibly preconditioned)
231  // residual
232  last_criterion = criterion(r, d);
233  conv = this->iteration_status(iter, last_criterion, x);
234  if (conv != SolverControl::iterate)
235  break;
236 
237  x.add(additional_data.omega, d);
238  print_vectors(iter, x, r, d);
239 
240  ++iter;
241  }
242 
243  // in case of failure: throw exception
244  if (conv != SolverControl::success)
245  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
246  // otherwise exit as normal
247 }
248 
249 
250 
251 template <class VectorType>
252 template <typename MatrixType, typename PreconditionerType>
253 void
254 SolverRichardson<VectorType>::Tsolve(const MatrixType & A,
255  VectorType & x,
256  const VectorType & b,
257  const PreconditionerType &preconditioner)
258 {
260  double last_criterion = std::numeric_limits<double>::lowest();
261 
262  unsigned int iter = 0;
263 
264  // Memory allocation.
265  // 'Vr' holds the residual, 'Vd' the preconditioned residual
266  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
267  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
268 
269  VectorType &r = *Vr;
270  r.reinit(x);
271 
272  VectorType &d = *Vd;
273  d.reinit(x);
274 
275  LogStream::Prefix prefix("RichardsonT");
276 
277  // Main loop
278  while (conv == SolverControl::iterate)
279  {
280  // Do not use Tresidual,
281  // but do it in 2 steps
282  A.Tvmult(r, x);
283  r.sadd(-1., 1., b);
284  preconditioner.Tvmult(d, r);
285 
286  last_criterion = criterion(r, d);
287  conv = this->iteration_status(iter, last_criterion, x);
288  if (conv != SolverControl::iterate)
289  break;
290 
291  x.add(additional_data.omega, d);
292  print_vectors(iter, x, r, d);
293 
294  ++iter;
295  }
296 
297  // in case of failure: throw exception
298  if (conv != SolverControl::success)
299  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
300 
301  // otherwise exit as normal
302 }
303 
304 
305 template <class VectorType>
306 void
308  const VectorType &,
309  const VectorType &,
310  const VectorType &) const
311 {}
312 
313 
314 
315 template <class VectorType>
316 inline typename VectorType::value_type
317 SolverRichardson<VectorType>::criterion(const VectorType &r,
318  const VectorType &d) const
319 {
320  if (!additional_data.use_preconditioned_residual)
321  return r.l2_norm();
322  else
323  return d.l2_norm();
324 }
325 
326 
327 template <class VectorType>
328 inline void
330 {
331  additional_data.omega = om;
332 }
333 
334 #endif // DOXYGEN
335 
337 
338 #endif
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
void set_omega(const double om=1.)
SolverRichardson(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void Tsolve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
virtual ~SolverRichardson() override=default
virtual VectorType::value_type criterion(const VectorType &r, const VectorType &d) const
AdditionalData additional_data
SolverRichardson(SolverControl &cn, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
#define AssertThrow(cond, exc)
Definition: exceptions.h:1619
static const char A
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
AdditionalData(const double omega=1, const bool use_preconditioned_residual=false)