solver_richardson.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 1999 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_solver_richardson_h
#define dealii_solver_richardson_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/logstream.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/template_constraints.h>
 
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>
 
#include <limits>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename VectorType = Vector<double>>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
class SolverRichardson : public SolverBase<VectorType>
{
public:
  struct AdditionalData
  {
    explicit AdditionalData(const double omega                       = 1,
                            const bool   use_preconditioned_residual = false);
 
    double omega;
 
    bool use_preconditioned_residual;
  };
 
  SolverRichardson(SolverControl            &cn,
                   VectorMemory<VectorType> &mem,
                   const AdditionalData     &data = AdditionalData());
 
  SolverRichardson(SolverControl        &cn,
                   const AdditionalData &data = AdditionalData());
 
  virtual ~SolverRichardson() override = default;
 
  template <typename MatrixType, typename PreconditionerType>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_linear_operator_on<MatrixType, VectorType> &&
     concepts::is_linear_operator_on<PreconditionerType, VectorType>))
  void solve(const MatrixType         &A,
             VectorType               &x,
             const VectorType         &b,
             const PreconditionerType &preconditioner);
 
  template <typename MatrixType, typename PreconditionerType>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_transpose_linear_operator_on<MatrixType, VectorType> &&
     concepts::is_transpose_linear_operator_on<PreconditionerType, VectorType>))
  void Tsolve(const MatrixType         &A,
              VectorType               &x,
              const VectorType         &b,
              const PreconditionerType &preconditioner);
 
  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;
 
protected:
  virtual typename VectorType::value_type
  criterion(const VectorType &r, const VectorType &d) const;
 
  AdditionalData additional_data;
};
 
/*----------------- Implementation of the Richardson Method ------------------*/
 
#ifndef DOXYGEN
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
inline SolverRichardson<VectorType>::AdditionalData::AdditionalData(
  const double omega,
  const bool   use_preconditioned_residual)
  : omega(omega)
  , use_preconditioned_residual(use_preconditioned_residual)
{}
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverRichardson<VectorType>::SolverRichardson(SolverControl            &cn,
                                               VectorMemory<VectorType> &mem,
                                               const AdditionalData     &data)
  : SolverBase<VectorType>(cn, mem)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverRichardson<VectorType>::SolverRichardson(SolverControl        &cn,
                                               const AdditionalData &data)
  : SolverBase<VectorType>(cn)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType, typename PreconditionerType>
DEAL_II_CXX20_REQUIRES(
  (concepts::is_linear_operator_on<MatrixType, VectorType> &&
   concepts::is_linear_operator_on<PreconditionerType, VectorType>))
void SolverRichardson<VectorType>::solve(
  const MatrixType         &A,
  VectorType               &x,
  const VectorType         &b,
  const PreconditionerType &preconditioner)
{
  SolverControl::State conv = SolverControl::iterate;
 
  double last_criterion = std::numeric_limits<double>::lowest();
 
  unsigned int iter = 0;
 
  // Memory allocation.
  // 'Vr' holds the residual, 'Vd' the preconditioned residual
  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
 
  VectorType &r = *Vr;
  r.reinit(x);
 
  VectorType &d = *Vd;
  d.reinit(x);
 
  LogStream::Prefix prefix("Richardson");
 
  // Main loop
  while (conv == SolverControl::iterate)
    {
      // Compute the residual:
      A.vmult(r, x);
      r.sadd(-1., 1., b);
 
      preconditioner.vmult(d, r);
 
      // get the required norm of the (possibly preconditioned)
      // residual
      last_criterion = criterion(r, d);
      conv           = this->iteration_status(iter, last_criterion, x);
      if (conv != SolverControl::iterate)
        break;
 
      // Add the correction to the current iterate. In many cases, one runs
      // Richardson's iteration with a step length (damping factor) of one,
      // in which case we can optimize the addition.
      if (additional_data.omega != 1.0)
        x.add(additional_data.omega, d);
      else
        x += d;
 
      print_vectors(iter, x, r, d);
 
      ++iter;
    }
 
  // in case of failure: throw exception
  if (conv != SolverControl::success)
    AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
  // otherwise exit as normal
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename MatrixType, typename PreconditionerType>
DEAL_II_CXX20_REQUIRES(
  (concepts::is_transpose_linear_operator_on<MatrixType, VectorType> &&
   concepts::is_transpose_linear_operator_on<PreconditionerType, VectorType>))
void SolverRichardson<VectorType>::Tsolve(
  const MatrixType         &A,
  VectorType               &x,
  const VectorType         &b,
  const PreconditionerType &preconditioner)
{
  SolverControl::State conv           = SolverControl::iterate;
  double               last_criterion = std::numeric_limits<double>::lowest();
 
  unsigned int iter = 0;
 
  // Memory allocation.
  // 'Vr' holds the residual, 'Vd' the preconditioned residual
  typename VectorMemory<VectorType>::Pointer Vr(this->memory);
  typename VectorMemory<VectorType>::Pointer Vd(this->memory);
 
  VectorType &r = *Vr;
  r.reinit(x);
 
  VectorType &d = *Vd;
  d.reinit(x);
 
  LogStream::Prefix prefix("RichardsonT");
 
  // Main loop
  while (conv == SolverControl::iterate)
    {
      // Do not use Tresidual,
      // but do it in 2 steps
      A.Tvmult(r, x);
      r.sadd(-1., 1., b);
      preconditioner.Tvmult(d, r);
 
      last_criterion = criterion(r, d);
      conv           = this->iteration_status(iter, last_criterion, x);
      if (conv != SolverControl::iterate)
        break;
 
      x.add(additional_data.omega, d);
      print_vectors(iter, x, r, d);
 
      ++iter;
    }
 
  // in case of failure: throw exception
  if (conv != SolverControl::success)
    AssertThrow(false, SolverControl::NoConvergence(iter, last_criterion));
 
  // otherwise exit as normal
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
void SolverRichardson<VectorType>::print_vectors(const unsigned int,
                                                 const VectorType &,
                                                 const VectorType &,
                                                 const VectorType &) const
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
inline typename VectorType::value_type
  SolverRichardson<VectorType>::criterion(const VectorType &r,
                                          const VectorType &d) const
{
  if (!additional_data.use_preconditioned_residual)
    return r.l2_norm();
  else
    return d.l2_norm();
}
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
inline void SolverRichardson<VectorType>::set_omega(const double om)
{
  additional_data.omega = om;
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif