Reference documentation for deal.II version GIT relicensing-218-g1f873f3929 2024-03-28 15:00: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\}}\)
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 - 2023 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
23
24#include <deal.II/lac/solver.h>
26
27#include <limits>
28
30
62template <typename VectorType = Vector<double>>
63class SolverRichardson : public SolverBase<VectorType>
64{
65public:
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
144protected:
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
165template <typename VectorType>
167 const double omega,
168 const bool use_preconditioned_residual)
169 : omega(omega)
170 , use_preconditioned_residual(use_preconditioned_residual)
171{}
172
173
174template <typename VectorType>
177 const AdditionalData &data)
178 : SolverBase<VectorType>(cn, mem)
179 , additional_data(data)
180{}
181
182
183
184template <typename VectorType>
186 const AdditionalData &data)
187 : SolverBase<VectorType>(cn)
188 , additional_data(data)
189{}
190
191
192
193template <typename VectorType>
194template <typename MatrixType, typename PreconditionerType>
195void
196SolverRichardson<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>::lowest();
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
250template <typename VectorType>
251template <typename MatrixType, typename PreconditionerType>
252void
253SolverRichardson<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>::lowest();
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
304template <typename VectorType>
305void
307 const VectorType &,
308 const VectorType &,
309 const VectorType &) const
310{}
311
312
313
314template <typename VectorType>
315inline typename VectorType::value_type
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
326template <typename VectorType>
327inline void
329{
330 additional_data.omega = om;
331}
332
333#endif // DOXYGEN
334
336
337#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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)