newtonsolver.h

Go to the documentation of this file.

#ifndef NEWTON_SOLVER_H_
#define NEWTON_SOLVER_H_

#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/numerics/vector_tools.h>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>

#include "parameterreader.h"




namespace DOpE
{
  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
    class NewtonSolver : public LINEARSOLVER
  {
  public:
    NewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader);
    ~NewtonSolver();

    static void declare_params(ParameterReader &param_reader);

    template<typename PROBLEM>
      void ReInit(PROBLEM& pde);

    template<typename PROBLEM>
      bool NonlinearSolve(PROBLEM& pde, VECTOR &solution, bool apply_boundary_values=true, 
                          bool force_matrix_build=false,
                          int priority = 5, std::string algo_level = "\t\t ");
    
  protected:
    
    inline INTEGRATOR& GetIntegrator();
   
  private:
    INTEGRATOR &integrator_;

    bool build_matrix_;

    double nonlinear_global_tol_, nonlinear_tol_, nonlinear_rho_;
    double linesearch_rho_;
    int nonlinear_maxiter_, line_maxiter_;
  };

  /**********************************Implementation*******************************************/

template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
 void NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
    ::declare_params(ParameterReader &param_reader)
  {
      param_reader.SetSubsection("newtonsolver parameters");
      param_reader.declare_entry("nonlinear_global_tol", "1.e-12",Patterns::Double(0),"global tolerance for the newton iteration");
      param_reader.declare_entry("nonlinear_tol", "1.e-10",Patterns::Double(0),"relative tolerance for the newton iteration");
      param_reader.declare_entry("nonlinear_maxiter", "10",Patterns::Integer(0),"maximal number of newton iterations");
      param_reader.declare_entry("nonlinear_rho", "0.1",Patterns::Double(0),"minimal  newton reduction, if actual reduction is less, matrix is rebuild ");
           
      param_reader.declare_entry("line_maxiter", "4",Patterns::Integer(0),"maximal number of linesearch steps");
      param_reader.declare_entry("linesearch_rho", "0.9",Patterns::Double(0),"reduction rate for the linesearch damping paramete");
 
      LINEARSOLVER::declare_params(param_reader);
  }

  /*******************************************************************************************/

  template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
    NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
    ::NewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader)
    : LINEARSOLVER(param_reader), integrator_(integrator)
    {
       param_reader.SetSubsection("newtonsolver parameters");
       nonlinear_global_tol_ = param_reader.get_double ("nonlinear_global_tol");
       nonlinear_tol_        = param_reader.get_double ("nonlinear_tol"); 
       nonlinear_maxiter_    = param_reader.get_integer ("nonlinear_maxiter"); 
       nonlinear_rho_        = param_reader.get_double ("nonlinear_rho"); 

       line_maxiter_   = param_reader.get_integer ("line_maxiter");
       linesearch_rho_ = param_reader.get_double ("linesearch_rho");
     
    }

  /*******************************************************************************************/

    template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
    NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
                   ::~NewtonSolver()
    {
    }

  /*******************************************************************************************/
 template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
    template<typename PROBLEM>
   void NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
                        ::ReInit(PROBLEM& pde)
    {
      LINEARSOLVER::ReInit(pde);
    }

  /*******************************************************************************************/
 template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
   template<typename PROBLEM>
   bool NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
   ::NonlinearSolve(PROBLEM& pde,
                    VECTOR &solution, 
                    bool apply_boundary_values, 
                    bool force_matrix_build,
                    int priority,
                    std::string algo_level)
    {
      bool build_matrix = force_matrix_build;
      VECTOR residual;
      VECTOR du;
      std::stringstream out;
      pde.GetOutputHandler()->InitNewtonOut(out);

      du.reinit(solution);
      residual.reinit(solution);
      
      if(apply_boundary_values)
      {
        GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
      }
      pde.GetDoFConstraints().distribute(solution);

      GetIntegrator().AddDomainData("last_newton_solution",&solution);
    
      GetIntegrator().ComputeNonlinearResidual(pde,residual);
      residual *= -1.;
           
      pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
      pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
  
      double res = residual.linfty_norm();
      double firstres = res;
      double lastres = res;

   
      out<< algo_level << "Newton step: " <<0<<"\t Residual (abs.): "
         <<pde.GetOutputHandler()->ZeroTolerance(res, 1.0)
         <<"\n";
        
      out<< algo_level << "Newton step: " <<0<<"\t Residual (rel.):   " << std::scientific << firstres/firstres; 
         

      pde.GetOutputHandler()->Write(out,priority);
      
      int iter=0;
      while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
      {
        iter++;
        
        if(iter > nonlinear_maxiter_)
        {
          GetIntegrator().DeleteDomainData("last_newton_solution");
          throw DOpEIterationException("Iteration count exceeded bounds!","NewtonSolver::NonlinearSolve");
        }
        
        pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");
        
        LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
        build_matrix = false;

        //Linesearch
        {
          solution += du;
          GetIntegrator().ComputeNonlinearResidual(pde,residual);
          residual *= -1.;

          pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
          pde.GetOutputHandler()->Write(du,"Update"+pde.GetType(),pde.GetDoFType());

          double newres = residual.linfty_norm();
          int lineiter=0;
          double rho = linesearch_rho_;
          double alpha=1;
          
          while(newres > res)
          {
            out<< algo_level << "Newton step: " <<iter<<"\t Residual (rel.): "
               <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0)
               << "\t LineSearch {"<<lineiter<<"} ";
            
            pde.GetOutputHandler()->Write(out,priority+1);

            lineiter++;
            if(lineiter > line_maxiter_)
            {
              GetIntegrator().DeleteDomainData("last_newton_solution");
              throw DOpEIterationException("Line-Iteration count exceeded bounds!","NewtonSolver::NonlinearSolve");
            }
            solution.add(alpha*(rho-1.),du);
            alpha*= rho;
            
            GetIntegrator().ComputeNonlinearResidual(pde,residual);
            residual *= -1.;
            pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
                
            newres = residual.linfty_norm();        

          }
          if(res/lastres > nonlinear_rho_)
          {
            build_matrix=true;
          }
          lastres=res;
          res=newres;

          out << algo_level 
              << "Newton step: " 
              <<iter
              <<"\t Residual (rel.): "
              << pde.GetOutputHandler()->ZeroTolerance(res/firstres, 1.0)
              << "\t LineSearch {"
              <<lineiter
              <<"} ";

          
          pde.GetOutputHandler()->Write(out,priority);

        }//End of Linesearch
      }
      GetIntegrator().DeleteDomainData("last_newton_solution");

      return build_matrix;
    }

  /*******************************************************************************************/
    template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
    INTEGRATOR& NewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
                               ::GetIntegrator()
    {
      return integrator_;
    }

  /*******************************************************************************************/

}
#endif