Reference documentation for deal.II version Git d745029578 2021-12-01 09:26:47 +0100
\(\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 - 2020 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 
29 
32 
59 template <class VectorType = Vector<double>>
60 class SolverRichardson : public SolverBase<VectorType>
61 {
62 public:
67  {
71  explicit AdditionalData(const double omega = 1,
72  const bool use_preconditioned_residual = false);
73 
77  double omega;
78 
83  };
84 
90  const AdditionalData & data = AdditionalData());
91 
97  const AdditionalData &data = AdditionalData());
98 
102  virtual ~SolverRichardson() override = default;
103 
107  template <typename MatrixType, typename PreconditionerType>
108  void
109  solve(const MatrixType & A,
110  VectorType & x,
111  const VectorType & b,
112  const PreconditionerType &preconditioner);
113 
117  template <typename MatrixType, typename PreconditionerType>
118  void
119  Tsolve(const MatrixType & A,
120  VectorType & x,
121  const VectorType & b,
122  const PreconditionerType &preconditioner);
123 
127  void
128  set_omega(const double om = 1.);
129 
135  virtual void
136  print_vectors(const unsigned int step,
137  const VectorType & x,
138  const VectorType & r,
139  const VectorType & d) const;
140 
141 protected:
148  virtual typename VectorType::value_type
149  criterion(const VectorType &r, const VectorType &d) const;
150 
155 };
156 
158 /*----------------- Implementation of the Richardson Method ------------------*/
159 
160 #ifndef DOXYGEN
161 
162 template <class VectorType>
164  const double omega,
165  const bool use_preconditioned_residual)
166  : omega(omega)
167  , use_preconditioned_residual(use_preconditioned_residual)
168 {}
169 
170 
171 template <class VectorType>
174  const AdditionalData & data)
175  : SolverBase<VectorType>(cn, mem)
176  , additional_data(data)
177 {}
178 
179 
180 
181 template <class VectorType>
183  const AdditionalData &data)
185  , additional_data(data)
186 {}
187 
188 
189 
190 template <class VectorType>
191 template <typename MatrixType, typename PreconditionerType>
192 void
193 SolverRichardson<VectorType>::solve(const MatrixType & A,
194  VectorType & x,
195  const VectorType & b,
196  const PreconditionerType &preconditioner)
197 {
199 
200  double last_criterion = std::numeric_limits<double>::lowest();
201 
202  unsigned int iter = 0;
203 
204  // Memory allocation.
205  // 'Vr' holds the residual, 'Vd' the preconditioned residual
206  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
207  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
208 
209  VectorType &r = *Vr;
210  r.reinit(x);
211 
212  VectorType &d = *Vd;
213  d.reinit(x);
214 
215  LogStream::Prefix prefix("Richardson");
216 
217  // Main loop
218  while (conv == SolverControl::iterate)
219  {
220  // Do not use residual,
221  // but do it in 2 steps
222  A.vmult(r, x);
223  r.sadd(-1., 1., b);
224  preconditioner.vmult(d, r);
225 
226  // get the required norm of the (possibly preconditioned)
227  // residual
228  last_criterion = criterion(r, d);
229  conv = this->iteration_status(iter, last_criterion, x);
230  if (conv != SolverControl::iterate)
231  break;
232 
233  x.add(additional_data.omega, d);
234  print_vectors(iter, x, r, d);
235 
236  ++iter;
237  }
238 
239  // in case of failure: throw exception
240  if (conv != SolverControl::success)
241  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
242  // otherwise exit as normal
243 }
244 
245 
246 
247 template <class VectorType>
248 template <typename MatrixType, typename PreconditionerType>
249 void
250 SolverRichardson<VectorType>::Tsolve(const MatrixType & A,
251  VectorType & x,
252  const VectorType & b,
253  const PreconditionerType &preconditioner)
254 {
256  double last_criterion = std::numeric_limits<double>::lowest();
257 
258  unsigned int iter = 0;
259 
260  // Memory allocation.
261  // 'Vr' holds the residual, 'Vd' the preconditioned residual
262  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
263  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
264 
265  VectorType &r = *Vr;
266  r.reinit(x);
267 
268  VectorType &d = *Vd;
269  d.reinit(x);
270 
271  LogStream::Prefix prefix("RichardsonT");
272 
273  // Main loop
274  while (conv == SolverControl::iterate)
275  {
276  // Do not use Tresidual,
277  // but do it in 2 steps
278  A.Tvmult(r, x);
279  r.sadd(-1., 1., b);
280  preconditioner.Tvmult(d, r);
281 
282  last_criterion = criterion(r, d);
283  conv = this->iteration_status(iter, last_criterion, x);
284  if (conv != SolverControl::iterate)
285  break;
286 
287  x.add(additional_data.omega, d);
288  print_vectors(iter, x, r, d);
289 
290  ++iter;
291  }
292 
293  // in case of failure: throw exception
294  if (conv != SolverControl::success)
295  AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
296 
297  // otherwise exit as normal
298 }
299 
300 
301 template <class VectorType>
302 void
304  const VectorType &,
305  const VectorType &,
306  const VectorType &) const
307 {}
308 
309 
310 
311 template <class VectorType>
312 inline typename VectorType::value_type
314  const VectorType &d) const
315 {
317  return r.l2_norm();
318  else
319  return d.l2_norm();
320 }
321 
322 
323 template <class VectorType>
324 inline void
326 {
327  additional_data.omega = om;
328 }
329 
330 #endif // DOXYGEN
331 
333 
334 #endif
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
Continue iteration.
void set_omega(const double om=1.)
boost::signals2::signal< SolverControl::State(const unsigned int iteration, const double check_value, const VectorType &current_iterate), StateCombiner > iteration_status
Definition: solver.h:471
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:1571
SolverRichardson(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:404
virtual ~SolverRichardson() override=default
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)
static const char A
virtual VectorType::value_type criterion(const VectorType &r, const VectorType &d) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:403
AdditionalData(const double omega=1, const bool use_preconditioned_residual=false)
VectorMemory< VectorType > & memory
Definition: solver.h:420