solver_control.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/logstream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/signaling_nan.h>

#include <deal.II/lac/solver_control.h>

#include <cmath>
#include <sstream>

DEAL_II_NAMESPACE_OPEN

/*----------------------- SolverControl ---------------------------------*/


SolverControl::SolverControl(const unsigned int maxiter,
                             const double       tolerance,
                             const bool         m_log_history,
                             const bool         m_log_result)
  : maxsteps(maxiter)
  , tol(tolerance)
  , lcheck(failure)
  , initial_val(numbers::signaling_nan<double>())
  , lvalue(numbers::signaling_nan<double>())
  , lstep(numbers::invalid_unsigned_int)
  , check_failure(false)
  , relative_failure_residual(0)
  , failure_residual(0)
  , m_log_history(m_log_history)
  , m_log_frequency(1)
  , m_log_result(m_log_result)
  , history_data_enabled(false)
{}



SolverControl::State
SolverControl::check(const unsigned int step, const double check_value)
{
  // if this is the first time we
  // come here, then store the
  // residual for later comparisons
  if (step == 0)
    {
      initial_val = check_value;
    }

  if (m_log_history && ((step % m_log_frequency) == 0))
    deallog << "Check " << step << "\t" << check_value << std::endl;

  lstep  = step;
  lvalue = check_value;

  if (step == 0)
    {
      if (check_failure)
        failure_residual = relative_failure_residual * check_value;

      if (m_log_result)
        deallog << "Starting value " << check_value << std::endl;
    }

  if (history_data_enabled)
    history_data.push_back(check_value);

  if (check_value <= tol)
    {
      if (m_log_result)
        deallog << "Convergence step " << step << " value " << check_value
                << std::endl;
      lcheck = success;
      return success;
    }

  if ((step >= maxsteps) || std::isnan(check_value) ||
      (check_failure && (check_value > failure_residual)))
    {
      if (m_log_result)
        deallog << "Failure step " << step << " value " << check_value
                << std::endl;
      lcheck = failure;
      return failure;
    }

  lcheck = iterate;
  return iterate;
}



SolverControl::State
SolverControl::last_check() const
{
  return lcheck;
}


double
SolverControl::initial_value() const
{
  return initial_val;
}


double
SolverControl::last_value() const
{
  return lvalue;
}


unsigned int
SolverControl::last_step() const
{
  return lstep;
}


unsigned int
SolverControl::log_frequency(unsigned int f)
{
  if (f == 0)
    f = 1;
  unsigned int old = m_log_frequency;
  m_log_frequency  = f;
  return old;
}


void
SolverControl::enable_history_data()
{
  history_data_enabled = true;
}



const std::vector<double> &
SolverControl::get_history_data() const
{
  Assert(history_data_enabled, ExcHistoryDataRequired());
  Assert(
    history_data.size() > 0,
    ExcMessage(
      "The SolverControl object was asked for the solver history "
      "data, but there is no data. Possibly you requested the data before the "
      "solver was run."));

  return history_data;
}



double
SolverControl::average_reduction() const
{
  if (lstep == 0)
    return 0.;

  Assert(history_data_enabled, ExcHistoryDataRequired());
  Assert(history_data.size() > lstep, ExcInternalError());
  Assert(history_data[0] > 0., ExcInternalError());
  Assert(history_data[lstep] > 0., ExcInternalError());

  return std::pow(history_data[lstep] / history_data[0], 1. / lstep);
}



double
SolverControl::step_reduction(unsigned int step) const
{
  Assert(history_data_enabled, ExcHistoryDataRequired());
  Assert(history_data.size() > lstep, ExcInternalError());
  Assert(step <= lstep, ExcIndexRange(step, 1, lstep + 1));
  Assert(step > 0, ExcIndexRange(step, 1, lstep + 1));

  return history_data[step] / history_data[step - 1];
}


double
SolverControl::final_reduction() const
{
  return step_reduction(lstep);
}


void
SolverControl::declare_parameters(ParameterHandler &param)
{
  param.declare_entry("Max steps", "100", Patterns::Integer());
  param.declare_entry("Tolerance", "1.e-10", Patterns::Double());
  param.declare_entry("Log history", "false", Patterns::Bool());
  param.declare_entry("Log frequency", "1", Patterns::Integer());
  param.declare_entry("Log result", "true", Patterns::Bool());
}


void
SolverControl::parse_parameters(ParameterHandler &param)
{
  set_max_steps(param.get_integer("Max steps"));
  set_tolerance(param.get_double("Tolerance"));
  log_history(param.get_bool("Log history"));
  log_result(param.get_bool("Log result"));
  log_frequency(param.get_integer("Log frequency"));
}

/*----------------------- ReductionControl ---------------------------------*/


ReductionControl::ReductionControl(const unsigned int n,
                                   const double       tol,
                                   const double       red,
                                   const bool         m_log_history,
                                   const bool         m_log_result)
  : SolverControl(n, tol, m_log_history, m_log_result)
  , reduce(red)
  , reduced_tol(numbers::signaling_nan<double>())
{}


ReductionControl::ReductionControl(const SolverControl &c)
  : SolverControl(c)
  , reduce(numbers::signaling_nan<double>())
  , reduced_tol(numbers::signaling_nan<double>())
{
  set_reduction(0.);
}


ReductionControl &
ReductionControl::operator=(const SolverControl &c)
{
  SolverControl::operator=(c);
  set_reduction(0.);
  return *this;
}



SolverControl::State
ReductionControl::check(const unsigned int step, const double check_value)
{
  // if this is the first time we
  // come here, then store the
  // residual for later comparisons
  if (step == 0)
    {
      initial_val = check_value;
      reduced_tol = check_value * reduce;
    }

  // check whether desired reduction
  // has been achieved. also check
  // for equality in case initial
  // residual already was zero
  if (check_value <= reduced_tol)
    {
      if (m_log_result)
        deallog << "Convergence step " << step << " value " << check_value
                << std::endl;
      lstep  = step;
      lvalue = check_value;
      lcheck = success;

      if (history_data_enabled)
        history_data.push_back(check_value);

      return success;
    }
  else
    return SolverControl::check(step, check_value);
}



void
ReductionControl::declare_parameters(ParameterHandler &param)
{
  SolverControl::declare_parameters(param);
  param.declare_entry("Reduction", "1.e-2", Patterns::Double());
}


void
ReductionControl::parse_parameters(ParameterHandler &param)
{
  SolverControl::parse_parameters(param);
  set_reduction(param.get_double("Reduction"));
}


/*---------------------- IterationNumberControl -----------------------------*/


IterationNumberControl::IterationNumberControl(const unsigned int n,
                                               const double       tolerance,
                                               const bool         m_log_history,
                                               const bool         m_log_result)
  : SolverControl(n, tolerance, m_log_history, m_log_result)
{}



SolverControl::State
IterationNumberControl::check(const unsigned int step, const double check_value)
{
  // check whether the given number of iterations was reached, and return
  // success in that case. Otherwise, go on to the check of the base class.
  if (step >= this->maxsteps)
    {
      if (m_log_result)
        deallog << "Convergence step " << step << " value " << check_value
                << std::endl;
      lstep  = step;
      lvalue = check_value;

      lcheck = success;
      return success;
    }
  else
    return SolverControl::check(step, check_value);
}

/*------------------------ ConsecutiveControl -------------------------------*/


ConsecutiveControl::ConsecutiveControl(
  const unsigned int n,
  const double       tolerance,
  const unsigned int n_consecutive_iterations,
  const bool         m_log_history,
  const bool         m_log_result)
  : SolverControl(n, tolerance, m_log_history, m_log_result)
  , n_consecutive_iterations(n_consecutive_iterations)
  , n_converged_iterations(0)
{
  AssertThrow(n_consecutive_iterations > 0,
              ExcMessage("n_consecutive_iterations should be positive"));
}



ConsecutiveControl::ConsecutiveControl(const SolverControl &c)
  : SolverControl(c)
  , n_consecutive_iterations(1)
  , n_converged_iterations(0)
{}



ConsecutiveControl &
ConsecutiveControl::operator=(const SolverControl &c)
{
  SolverControl::operator  =(c);
  n_consecutive_iterations = 1;
  n_converged_iterations   = 0;
  return *this;
}



SolverControl::State
ConsecutiveControl::check(const unsigned int step, const double check_value)
{
  // reset the counter if ConsecutiveControl is being reused
  if (step == 0)
    n_converged_iterations = 0;
  else
    {
      // check two things:
      // (i)  steps are ascending without repetitions
      // (ii) user started from zero even when solver is being reused.
      Assert(step - 1 == lstep,
             ExcMessage("steps should be ascending integers."));
    }

  SolverControl::State state = SolverControl::check(step, check_value);
  // check if we need to override the success:
  if (state == success)
    {
      n_converged_iterations++;
      if (n_converged_iterations == n_consecutive_iterations)
        {
          return success;
        }
      else
        {
          lcheck = iterate;
          return iterate;
        }
    }
  else
    {
      n_converged_iterations = 0;
      return state;
    }
}

DEAL_II_NAMESPACE_CLOSE