nonlinear_matrix_solver.h

Go to the documentation of this file.

// This file is part of HermesCommon
//
// Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
// Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
//
// Hermes2D is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation; either version 2 of the License,
// or (at your option) any later version.
//
// Hermes2D is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#ifndef __HERMES_COMMON_NONLINEAR_MATRIX_SOLVER_H_
#define __HERMES_COMMON_NONLINEAR_MATRIX_SOLVER_H_

#include "precond.h"
#include "exceptions.h"
#include "algebra/cs_matrix.h"
#include "algebra/vector.h"
#include "mixins.h"
#include "algebra/algebra_mixins.h"
#include "solvers/nonlinear_convergence_measurement.h"
#include "solvers/matrix_solver.h"

namespace Hermes
{
  namespace Solvers
  {
    const int NonlinearConvergenceMeasurementTypeCount = 6;

    enum NonlinearConvergenceMeasurementType
    {
      ResidualNormRelativeToPrevious = 0x0001,
      ResidualNormRatioToInitial = 0x0002,
      ResidualNormRatioToPrevious = 0x0004,
      ResidualNormAbsolute = 0x0008,
      SolutionChangeAbsolute = 0x0010,
      SolutionChangeRelative = 0x0020
    };

    enum NonlinearConvergenceState
    {
      Converged,
      NotConverged,
      BelowMinDampingCoeff,
      AboveMaxAllowedResidualNorm,
      AboveMaxIterations,
      Error
    };

    template<typename Scalar> class HERMES_API NonlinearConvergenceMeasurement;

    template <typename Scalar>
    class NonlinearMatrixSolver :
      public virtual Hermes::Solvers::MatrixSolver<Scalar>,
      public virtual Hermes::Mixins::OutputAttachable,
      public virtual Hermes::Mixins::StateQueryable,
      public virtual Hermes::Mixins::TimeMeasurable
    {
    public:
      NonlinearMatrixSolver();
      virtual ~NonlinearMatrixSolver();

      virtual void solve(Scalar* coeff_vec);

      void set_max_allowed_iterations(int max_allowed_iterations);

      virtual void clear_tolerances();

      void set_max_allowed_residual_norm(double max_allowed_residual_norm_to_set);

      void set_tolerance(double newton_tol, NonlinearConvergenceMeasurementType toleranceType, bool handleMultipleTolerancesAnd = false);

      int get_num_iters() const;

#pragma region damping-public
      void set_min_allowed_damping_coeff(double min_allowed_damping_coeff_to_set);

      void set_manual_damping_coeff(bool onOff, double coeff);

      void set_initial_auto_damping_coeff(double coeff);

      void set_auto_damping_ratio(double ratio);

      void set_sufficient_improvement_factor(double ratio);

      void set_necessary_successful_steps_to_increase(unsigned int steps);
#pragma endregion

#pragma region jacobian_recalculation-public
      void set_sufficient_improvement_factor_jacobian(double ratio);

      void set_max_steps_with_reused_jacobian(unsigned int steps);
#pragma endregion

      virtual void free();

    protected:

#pragma region damping-private
      bool manual_damping;

      double manual_damping_factor;

      double auto_damping_ratio;
      double initial_auto_damping_factor;
      double sufficient_improvement_factor;
      unsigned int necessary_successful_steps_to_increase;
      double min_allowed_damping_coeff;
#pragma endregion

#pragma region jacobian_recalculation-private
      bool force_reuse_jacobian_values(unsigned int& successful_steps_with_reused_jacobian);
      bool jacobian_reused_okay(unsigned int& successful_steps_with_reused_jacobian);

      double sufficient_improvement_factor_jacobian;
      unsigned int max_steps_with_reused_jacobian;

      Vector<Scalar>* residual_back;
#pragma endregion

      virtual void assemble_residual(bool store_previous_residual) = 0;
      virtual bool assemble_jacobian(bool store_previous_jacobian) = 0;
      virtual bool assemble(bool store_previous_jacobian, bool store_previous_residual) = 0;

      virtual void on_damping_factor_updated();
      virtual void on_reused_jacobian_step_begin();
      virtual void on_reused_jacobian_step_end();

      virtual bool handle_convergence_state_return_finished(NonlinearConvergenceState state);

      virtual void init_solving(Scalar* coeff_vec);

      virtual NonlinearConvergenceState get_convergence_state();

      bool do_initial_step_return_finished();

      virtual void solve_linear_system();

      virtual double update_solution_return_change_norm(Scalar* linear_system_solution) = 0;

      void finalize_solving();

      virtual void deinit_solving();

      bool calculate_damping_factor(unsigned int& successful_steps);

      virtual bool damping_factor_condition();

      int get_current_iteration_number();

      void step_info();

      virtual double calculate_residual_norm();

      virtual bool isOkay() const;
      inline std::string getClassName() const { return "NonlinearMatrixSolver"; }

      void init_nonlinear();

      double max_allowed_residual_norm;

      int max_allowed_iterations;

      double tolerance[NonlinearConvergenceMeasurementTypeCount];

      bool tolerance_set[NonlinearConvergenceMeasurementTypeCount];

      bool handleMultipleTolerancesAnd;

      int num_iters;

      SparseMatrix<Scalar>* previous_jacobian;
      Vector<Scalar>* previous_residual;

#pragma region OutputAttachable
      // For derived classes - read-only access.
      const OutputParameterUnsignedInt& iteration() const { return this->p_iteration; };
      const OutputParameterDoubleVector& residual_norms() const { return this->p_residual_norms; };
      const OutputParameterDoubleVector& solution_norms() const { return this->p_solution_norms; };
      const OutputParameterDoubleVector& solution_change_norms() const { return this->p_solution_change_norms; };
      const OutputParameterUnsignedInt& successful_steps_damping() const { return this->p_successful_steps_damping; };
      const OutputParameterUnsignedInt& successful_steps_jacobian() const { return this->p_successful_steps_jacobian; };
      const OutputParameterDoubleVector& damping_factors() const { return this->p_damping_factors; };
      const OutputParameterBool& residual_norm_drop() const { return this->p_residual_norm_drop; };
      const OutputParameterBoolVector& iterations_with_recalculated_jacobian() const { return this->p_iterations_with_recalculated_jacobian; };

      OutputParameterDoubleVector p_residual_norms;
      OutputParameterDoubleVector p_solution_norms;
      OutputParameterDoubleVector p_solution_change_norms;
      OutputParameterBoolVector p_iterations_with_recalculated_jacobian;
      OutputParameterUnsignedInt p_successful_steps_damping;
      OutputParameterUnsignedInt p_successful_steps_jacobian;
      OutputParameterDoubleVector p_damping_factors;
      OutputParameterBool p_residual_norm_drop;
      OutputParameterUnsignedInt p_iteration;
#pragma endregion

      Scalar* previous_sln_vector;
      bool use_initial_guess_for_iterative_solvers;
      friend class NonlinearConvergenceMeasurement < Scalar > ;
    };
  }

  namespace Exceptions
  {
    class HERMES_API NonlinearException : public Hermes::Exceptions::Exception
    {
    public:
      NonlinearException(Hermes::Solvers::NonlinearConvergenceState convergenceState);

      Hermes::Solvers::NonlinearConvergenceState get_exception_state();

    protected:
      Hermes::Solvers::NonlinearConvergenceState convergenceState;
    };
  }
}
#endif