solver_fire.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2017 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_solver_fire_h
#define dealii_solver_fire_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/logstream.h>
#include <deal.II/base/template_constraints.h>
 
#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/solver.h>
 
#include <functional>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
template <typename VectorType = Vector<double>>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
class SolverFIRE : public SolverBase<VectorType>
{
public:
  struct AdditionalData
  {
    explicit AdditionalData(const double initial_timestep    = 0.1,
                            const double maximum_timestep    = 1,
                            const double maximum_linfty_norm = 1);
 
    const double initial_timestep;
 
    const double maximum_timestep;
 
    const double maximum_linfty_norm;
  };
 
  SolverFIRE(SolverControl            &solver_control,
             VectorMemory<VectorType> &vector_memory,
             const AdditionalData     &data = AdditionalData());
 
  SolverFIRE(SolverControl        &solver_control,
             const AdditionalData &data = AdditionalData());
 
  template <typename PreconditionerType = DiagonalMatrix<VectorType>>
  void
  solve(const std::function<double(VectorType &, const VectorType &)> &compute,
        VectorType                                                    &x,
        const PreconditionerType &inverse_mass_matrix);
 
  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);
 
protected:
  virtual void
  print_vectors(const unsigned int,
                const VectorType &x,
                const VectorType &v,
                const VectorType &g) const;
 
  const AdditionalData additional_data;
};
 
/*------------------------- Implementation ----------------------------*/
 
#ifndef DOXYGEN
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverFIRE<VectorType>::AdditionalData::AdditionalData(
  const double initial_timestep,
  const double maximum_timestep,
  const double maximum_linfty_norm)
  : initial_timestep(initial_timestep)
  , maximum_timestep(maximum_timestep)
  , maximum_linfty_norm(maximum_linfty_norm)
{
  AssertThrow(initial_timestep > 0. && maximum_timestep > 0. &&
                maximum_linfty_norm > 0.,
              ExcMessage("Expected positive values for initial_timestep, "
                         "maximum_timestep and maximum_linfty_norm but one "
                         "or more of the these values are not positive."));
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverFIRE<VectorType>::SolverFIRE(SolverControl            &solver_control,
                                   VectorMemory<VectorType> &vector_memory,
                                   const AdditionalData     &data)
  : SolverBase<VectorType>(solver_control, vector_memory)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
SolverFIRE<VectorType>::SolverFIRE(SolverControl        &solver_control,
                                   const AdditionalData &data)
  : SolverBase<VectorType>(solver_control)
  , additional_data(data)
{}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
template <typename PreconditionerType>
void SolverFIRE<VectorType>::solve(
  const std::function<double(VectorType &, const VectorType &)> &compute,
  VectorType                                                    &x,
  const PreconditionerType &inverse_mass_matrix)
{
  LogStream::Prefix prefix("FIRE");
 
  // FIRE algorithm constants
  const double DELAYSTEP       = 5;
  const double TIMESTEP_GROW   = 1.1;
  const double TIMESTEP_SHRINK = 0.5;
  const double ALPHA_0         = 0.1;
  const double ALPHA_SHRINK    = 0.99;
 
  using real_type = typename VectorType::real_type;
 
  typename VectorMemory<VectorType>::Pointer v(this->memory);
  typename VectorMemory<VectorType>::Pointer g(this->memory);
 
  // Set velocities to zero but not gradients
  // as we are going to compute them soon.
  v->reinit(x, false);
  g->reinit(x, true);
 
  // Refer to v and g with some readable names.
  VectorType &velocities = *v;
  VectorType &gradients  = *g;
 
  // Update gradients for the new x.
  compute(gradients, x);
 
  unsigned int iter = 0;
 
  SolverControl::State conv = SolverControl::iterate;
  conv = this->iteration_status(iter, gradients * gradients, x);
  if (conv != SolverControl::iterate)
    return;
 
  // Refer to additional data members with some readable names.
  const auto &maximum_timestep = additional_data.maximum_timestep;
  double      timestep         = additional_data.initial_timestep;
 
  // First scaling factor.
  double alpha = ALPHA_0;
 
  unsigned int previous_iter_with_positive_v_dot_g = 0;
 
  while (conv == SolverControl::iterate)
    {
      ++iter;
      // Euler integration step.
      x.add(timestep, velocities);                     // x += dt     * v
      inverse_mass_matrix.vmult(gradients, gradients); // g  = M^{-1} * g
      velocities.add(-timestep, gradients);            // v -= dt     * h
 
      // Compute gradients for the new x.
      compute(gradients, x);
 
      const real_type gradient_norm_squared = gradients * gradients;
      conv = this->iteration_status(iter, gradient_norm_squared, x);
      if (conv != SolverControl::iterate)
        break;
 
      // v_dot_g = V * G
      const real_type v_dot_g = velocities * gradients;
 
      if (v_dot_g < 0.)
        {
          const real_type velocities_norm_squared = velocities * velocities;
 
          // Check if we divide by zero in DEBUG mode.
          Assert(gradient_norm_squared > 0., ExcInternalError());
 
          // beta = - alpha |V|/|G|
          const real_type beta =
            -alpha * std::sqrt(velocities_norm_squared / gradient_norm_squared);
 
          // V = (1-alpha) V + beta G.
          velocities.sadd(1. - alpha, beta, gradients);
 
          if (iter - previous_iter_with_positive_v_dot_g > DELAYSTEP)
            {
              // Increase timestep and decrease alpha.
              timestep = std::min(timestep * TIMESTEP_GROW, maximum_timestep);
              alpha *= ALPHA_SHRINK;
            }
        }
      else
        {
          // Decrease timestep, reset alpha and set V = 0.
          previous_iter_with_positive_v_dot_g = iter;
          timestep *= TIMESTEP_SHRINK;
          alpha      = ALPHA_0;
          velocities = 0.;
        }
 
      real_type vmax = velocities.linfty_norm();
 
      // Change timestep if any dof would move more than maximum_linfty_norm.
      if (vmax > 0.)
        {
          const double minimal_timestep =
            additional_data.maximum_linfty_norm / vmax;
          if (minimal_timestep < timestep)
            timestep = minimal_timestep;
        }
 
      print_vectors(iter, x, velocities, gradients);
 
    } // While we need to iterate.
 
  // In the case of failure: throw exception.
  if (conv != SolverControl::success)
    AssertThrow(false,
                SolverControl::NoConvergence(iter, gradients * gradients));
}
 
 
 
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 SolverFIRE<VectorType>::solve(const MatrixType         &A,
                                   VectorType               &x,
                                   const VectorType         &b,
                                   const PreconditionerType &preconditioner)
{
  std::function<double(VectorType &, const VectorType &)> compute_func =
    [&](VectorType &g, const VectorType &x) -> double {
    // Residual of the quadratic form @f$ \frac{1}{2} xAx - xb @f$.
    // G = b - Ax
    A.residual(g, x, b);
 
    // Gradient G = Ax -b.
    g *= -1.;
 
    // The quadratic form @f$\frac{1}{2} xAx - xb @f$.
    return 0.5 * A.matrix_norm_square(x) - x * b;
  };
 
  this->solve(compute_func, x, preconditioner);
}
 
 
 
template <typename VectorType>
DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
void SolverFIRE<VectorType>::print_vectors(const unsigned int,
                                           const VectorType &,
                                           const VectorType &,
                                           const VectorType &) const
{}
 
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif