Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
solver_richardson.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2022 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
24
25#include <deal.II/lac/solver.h>
27
28#include <limits>
29
31
63template <class VectorType = Vector<double>>
64class SolverRichardson : public SolverBase<VectorType>
65{
66public:
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
145protected:
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
166template <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
175template <class VectorType>
178 const AdditionalData & data)
179 : SolverBase<VectorType>(cn, mem)
180 , additional_data(data)
181{}
182
183
184
185template <class VectorType>
187 const AdditionalData &data)
188 : SolverBase<VectorType>(cn)
189 , additional_data(data)
190{}
191
192
193
194template <class VectorType>
195template <typename MatrixType, typename PreconditionerType>
196void
197SolverRichardson<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
251template <class VectorType>
252template <typename MatrixType, typename PreconditionerType>
253void
254SolverRichardson<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
305template <class VectorType>
306void
308 const VectorType &,
309 const VectorType &,
310 const VectorType &) const
311{}
312
313
314
315template <class VectorType>
316inline typename VectorType::value_type
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
327template <class VectorType>
328inline 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define AssertThrow(cond, exc)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
AdditionalData(const double omega=1, const bool use_preconditioned_residual=false)