Reference documentation for deal.II version 9.6.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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_richardson_h
16#define dealii_solver_richardson_h
17
18
19#include <deal.II/base/config.h>
20
24
25#include <deal.II/lac/solver.h>
27
28#include <limits>
29
31
63template <typename VectorType = Vector<double>>
65class SolverRichardson : public SolverBase<VectorType>
66{
67public:
72 {
76 explicit AdditionalData(const double omega = 1,
77 const bool use_preconditioned_residual = false);
78
82 double omega;
83
88 };
89
95 const AdditionalData &data = AdditionalData());
96
102 const AdditionalData &data = AdditionalData());
103
107 virtual ~SolverRichardson() override = default;
108
112 template <typename MatrixType, typename PreconditionerType>
116 void solve(const MatrixType &A,
117 VectorType &x,
118 const VectorType &b,
119 const PreconditionerType &preconditioner);
120
124 template <typename MatrixType, typename PreconditionerType>
126 (concepts::is_transpose_linear_operator_on<MatrixType, VectorType> &&
127 concepts::is_transpose_linear_operator_on<PreconditionerType, VectorType>))
128 void Tsolve(const MatrixType &A,
129 VectorType &x,
130 const VectorType &b,
131 const PreconditionerType &preconditioner);
132
136 void
137 set_omega(const double om = 1.);
138
144 virtual void
145 print_vectors(const unsigned int step,
146 const VectorType &x,
147 const VectorType &r,
148 const VectorType &d) const;
149
150protected:
157 virtual typename VectorType::value_type
158 criterion(const VectorType &r, const VectorType &d) const;
159
163 AdditionalData additional_data;
164};
165
167/*----------------- Implementation of the Richardson Method ------------------*/
168
169#ifndef DOXYGEN
170
171template <typename VectorType>
174 const double omega,
175 const bool use_preconditioned_residual)
176 : omega(omega)
177 , use_preconditioned_residual(use_preconditioned_residual)
178{}
179
180
181template <typename VectorType>
185 const AdditionalData &data)
186 : SolverBase<VectorType>(cn, mem)
187 , additional_data(data)
188{}
189
190
191
192template <typename VectorType>
195 const AdditionalData &data)
196 : SolverBase<VectorType>(cn)
197 , additional_data(data)
198{}
199
200
201
202template <typename VectorType>
204template <typename MatrixType, typename PreconditionerType>
209 const MatrixType &A,
210 VectorType &x,
211 const VectorType &b,
212 const PreconditionerType &preconditioner)
213{
215
216 double last_criterion = std::numeric_limits<double>::lowest();
217
218 unsigned int iter = 0;
219
220 // Memory allocation.
221 // 'Vr' holds the residual, 'Vd' the preconditioned residual
222 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
223 typename VectorMemory<VectorType>::Pointer Vd(this->memory);
224
225 VectorType &r = *Vr;
226 r.reinit(x);
227
228 VectorType &d = *Vd;
229 d.reinit(x);
230
231 LogStream::Prefix prefix("Richardson");
232
233 // Main loop
234 while (conv == SolverControl::iterate)
235 {
236 // Do not use residual,
237 // but do it in 2 steps
238 A.vmult(r, x);
239 r.sadd(-1., 1., b);
240 preconditioner.vmult(d, r);
241
242 // get the required norm of the (possibly preconditioned)
243 // residual
244 last_criterion = criterion(r, d);
245 conv = this->iteration_status(iter, last_criterion, x);
246 if (conv != SolverControl::iterate)
247 break;
248
249 x.add(additional_data.omega, d);
250 print_vectors(iter, x, r, d);
251
252 ++iter;
253 }
254
255 // in case of failure: throw exception
256 if (conv != SolverControl::success)
257 AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
258 // otherwise exit as normal
259}
260
261
262
263template <typename VectorType>
265template <typename MatrixType, typename PreconditionerType>
270 const MatrixType &A,
271 VectorType &x,
272 const VectorType &b,
273 const PreconditionerType &preconditioner)
274{
276 double last_criterion = std::numeric_limits<double>::lowest();
277
278 unsigned int iter = 0;
279
280 // Memory allocation.
281 // 'Vr' holds the residual, 'Vd' the preconditioned residual
282 typename VectorMemory<VectorType>::Pointer Vr(this->memory);
283 typename VectorMemory<VectorType>::Pointer Vd(this->memory);
284
285 VectorType &r = *Vr;
286 r.reinit(x);
287
288 VectorType &d = *Vd;
289 d.reinit(x);
290
291 LogStream::Prefix prefix("RichardsonT");
292
293 // Main loop
294 while (conv == SolverControl::iterate)
295 {
296 // Do not use Tresidual,
297 // but do it in 2 steps
298 A.Tvmult(r, x);
299 r.sadd(-1., 1., b);
300 preconditioner.Tvmult(d, r);
301
302 last_criterion = criterion(r, d);
303 conv = this->iteration_status(iter, last_criterion, x);
304 if (conv != SolverControl::iterate)
305 break;
306
307 x.add(additional_data.omega, d);
308 print_vectors(iter, x, r, d);
309
310 ++iter;
311 }
312
313 // in case of failure: throw exception
314 if (conv != SolverControl::success)
315 AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
316
317 // otherwise exit as normal
318}
319
320
321
322template <typename VectorType>
324void SolverRichardson<VectorType>::print_vectors(const unsigned int,
325 const VectorType &,
326 const VectorType &,
327 const VectorType &) const
328{}
329
330
331
332template <typename VectorType>
334inline typename VectorType::value_type
335 SolverRichardson<VectorType>::criterion(const VectorType &r,
336 const VectorType &d) const
337{
338 if (!additional_data.use_preconditioned_residual)
339 return r.l2_norm();
340 else
341 return d.l2_norm();
342}
343
344
345template <typename VectorType>
347inline void SolverRichardson<VectorType>::set_omega(const double om)
348{
349 additional_data.omega = om;
350}
351
352#endif // DOXYGEN
353
355
356#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
SolverRichardson(SolverControl &cn, const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:177
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#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)