Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
solver_richardson.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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>
23 #include <deal.II/base/signaling_nan.h>
24 
25 #include <deal.II/lac/solver.h>
26 #include <deal.II/lac/solver_control.h>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
32 
62 template <class VectorType = Vector<double>>
63 class SolverRichardson : public SolverBase<VectorType>
64 {
65 public:
70  {
74  explicit AdditionalData(const double omega = 1,
75  const bool use_preconditioned_residual = false);
76 
80  double omega;
81 
86  };
87 
93  const AdditionalData & data = AdditionalData());
94 
100  const AdditionalData &data = AdditionalData());
101 
105  virtual ~SolverRichardson() override = default;
106 
110  template <typename MatrixType, typename PreconditionerType>
111  void
112  solve(const MatrixType & A,
113  VectorType & x,
114  const VectorType & b,
115  const PreconditionerType &preconditioner);
116 
120  template <typename MatrixType, typename PreconditionerType>
121  void
122  Tsolve(const MatrixType & A,
123  VectorType & x,
124  const VectorType & b,
125  const PreconditionerType &preconditioner);
126 
130  void
131  set_omega(const double om = 1.);
132 
138  virtual void
139  print_vectors(const unsigned int step,
140  const VectorType & x,
141  const VectorType & r,
142  const VectorType & d) const;
143 
144 protected:
151  virtual typename VectorType::value_type
152  criterion(const VectorType &r, const VectorType &d) const;
153 
158 };
159 
161 /*----------------- Implementation of the Richardson Method ------------------*/
162 
163 #ifndef DOXYGEN
164 
165 template <class VectorType>
167  const double omega,
168  const bool use_preconditioned_residual)
169  : omega(omega)
170  , use_preconditioned_residual(use_preconditioned_residual)
171 {}
172 
173 
174 template <class VectorType>
177  const AdditionalData & data)
178  : SolverBase<VectorType>(cn, mem)
179  , additional_data(data)
180 {}
181 
182 
183 
184 template <class VectorType>
186  const AdditionalData &data)
187  : SolverBase<VectorType>(cn)
188  , additional_data(data)
189 {}
190 
191 
192 
193 template <class VectorType>
194 template <typename MatrixType, typename PreconditionerType>
195 void
196 SolverRichardson<VectorType>::solve(const MatrixType & A,
197  VectorType & x,
198  const VectorType & b,
199  const PreconditionerType &preconditioner)
200 {
202 
203  double last_criterion = -std::numeric_limits<double>::max();
204 
205  unsigned int iter = 0;
206 
207  // Memory allocation.
208  // 'Vr' holds the residual, 'Vd' the preconditioned residual
209  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
210  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
211 
212  VectorType &r = *Vr;
213  r.reinit(x);
214 
215  VectorType &d = *Vd;
216  d.reinit(x);
217 
218  LogStream::Prefix prefix("Richardson");
219 
220  // Main loop
221  while (conv == SolverControl::iterate)
222  {
223  // Do not use residual,
224  // but do it in 2 steps
225  A.vmult(r, x);
226  r.sadd(-1., 1., b);
227  preconditioner.vmult(d, r);
228 
229  // get the required norm of the (possibly preconditioned)
230  // residual
231  last_criterion = criterion(r, d);
232  conv = this->iteration_status(iter, last_criterion, x);
233  if (conv != SolverControl::iterate)
234  break;
235 
236  x.add(additional_data.omega, d);
237  print_vectors(iter, x, r, d);
238 
239  ++iter;
240  }
241 
242  // in case of failure: throw exception
243  if (conv != SolverControl::success)
244  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
245  // otherwise exit as normal
246 }
247 
248 
249 
250 template <class VectorType>
251 template <typename MatrixType, typename PreconditionerType>
252 void
253 SolverRichardson<VectorType>::Tsolve(const MatrixType & A,
254  VectorType & x,
255  const VectorType & b,
256  const PreconditionerType &preconditioner)
257 {
259  double last_criterion = -std::numeric_limits<double>::max();
260 
261  unsigned int iter = 0;
262 
263  // Memory allocation.
264  // 'Vr' holds the residual, 'Vd' the preconditioned residual
265  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
266  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
267 
268  VectorType &r = *Vr;
269  r.reinit(x);
270 
271  VectorType &d = *Vd;
272  d.reinit(x);
273 
274  LogStream::Prefix prefix("RichardsonT");
275 
276  // Main loop
277  while (conv == SolverControl::iterate)
278  {
279  // Do not use Tresidual,
280  // but do it in 2 steps
281  A.Tvmult(r, x);
282  r.sadd(-1., 1., b);
283  preconditioner.Tvmult(d, r);
284 
285  last_criterion = criterion(r, d);
286  conv = this->iteration_status(iter, last_criterion, x);
287  if (conv != SolverControl::iterate)
288  break;
289 
290  x.add(additional_data.omega, d);
291  print_vectors(iter, x, r, d);
292 
293  ++iter;
294  }
295 
296  // in case of failure: throw exception
297  if (conv != SolverControl::success)
298  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
299 
300  // otherwise exit as normal
301 }
302 
303 
304 template <class VectorType>
305 void
307  const VectorType &,
308  const VectorType &,
309  const VectorType &) const
310 {}
311 
312 
313 
314 template <class VectorType>
315 inline typename VectorType::value_type
316 SolverRichardson<VectorType>::criterion(const VectorType &r,
317  const VectorType &d) const
318 {
319  if (!additional_data.use_preconditioned_residual)
320  return r.l2_norm();
321  else
322  return d.l2_norm();
323 }
324 
325 
326 template <class VectorType>
327 inline void
329 {
330  additional_data.omega = om;
331 }
332 
333 #endif // DOXYGEN
334 
335 DEAL_II_NAMESPACE_CLOSE
336 
337 #endif
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Continue iteration.
void set_omega(const double om=1.)
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
AdditionalData additional_data
void Tsolve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
SolverRichardson(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.
virtual ~SolverRichardson() override=default
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
virtual VectorType::value_type criterion(const VectorType &r, const VectorType &d) const
AdditionalData(const double omega=1, const bool use_preconditioned_residual=false)