solver_fire.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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/lac/diagonal_matrix.h>
#include <deal.II/lac/solver.h>

#include <functional>


DEAL_II_NAMESPACE_OPEN



template <typename VectorType = Vector<double>>
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());

  virtual ~SolverFIRE();

  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>
  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>
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>
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>
SolverFIRE<VectorType>::SolverFIRE(SolverControl &       solver_control,
                                   const AdditionalData &data)
  : SolverBase<VectorType>(solver_control)
  , additional_data(data)
{}



template <typename VectorType>
SolverFIRE<VectorType>::~SolverFIRE()
{}



template <typename 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>
template <typename MatrixType, typename PreconditionerType>
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>
void
SolverFIRE<VectorType>::print_vectors(const unsigned int,
                                      const VectorType &,
                                      const VectorType &,
                                      const VectorType &) const
{}



#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif