instat_step_newtonsolver.h

Go to the documentation of this file.

#ifndef INSTAT_STEP_NEWTON_SOLVER_H_
#define INSTAT_STEP_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 InstatStepNewtonSolver : public LINEARSOLVER
  {
  public:
    InstatStepNewtonSolver(INTEGRATOR &integrator, ParameterReader &param_reader);
    ~InstatStepNewtonSolver();

    static void declare_params(ParameterReader &param_reader);

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

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

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

    template<typename PROBLEM>
      bool NonlinearSolve(PROBLEM& pde, const VECTOR &last_time_solution, VECTOR &solution, 
                          bool apply_boundary_values=true, 
                          bool force_matrix_build=false, int priority = 5, std::string algo_level = "\t\t ");
    /******************************************************/

    template<typename PROBLEM>
      bool NonlinearSolve_Initial(PROBLEM& pde, VECTOR &solution, bool apply_boundary_values=true, 
                        bool force_matrix_build=false, int priority = 5, std::string algo_level = "\t\t ");

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

    template<typename PROBLEM>
      void NonlinearLastTimeEvals(PROBLEM& pde, const VECTOR &last_time_solution, VECTOR &residual);
    
  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 InstatStepNewtonSolver<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>
    InstatStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
    ::InstatStepNewtonSolver(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>
    InstatStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
                   ::~InstatStepNewtonSolver()
    {
    }

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

 /*******************************************************************************************/
 template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
   template<typename PROBLEM>
   void InstatStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
   ::NonlinearLastTimeEvals(PROBLEM& pde, const VECTOR &last_time_solution, VECTOR &residual)
   { 
     VECTOR tmp_residual;
     tmp_residual.reinit(residual);
     residual =0.; 
     GetIntegrator().AddDomainData("last_newton_solution",&last_time_solution);
     GetIntegrator().AddDomainData("last_time_solution",&last_time_solution);
     pde.SetStepPart("Old");
     GetIntegrator().ComputeNonlinearLhs(pde,residual);   
     GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);
     tmp_residual *= -1;
     residual += tmp_residual;

     GetIntegrator().DeleteDomainData("last_newton_solution");
     GetIntegrator().DeleteDomainData("last_time_solution");

   }
   /*******************************************************************************************/

 template <typename INTEGRATOR, typename LINEARSOLVER, typename VECTOR>
   template<typename PROBLEM>
   bool InstatStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
   ::NonlinearSolve_Initial(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 << pde.GetOutputHandler()->ZeroTolerance(firstres/firstres,1.0);
   
   
   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!","InstatStepNewtonSolver::NonlinearSolve_Initial");
     }
     
     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!","InstatStepNewtonSolver::NonlinearSolve_Initial");
         }
         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>
   template<typename PROBLEM>
    bool InstatStepNewtonSolver<INTEGRATOR,LINEARSOLVER, VECTOR>
                        ::NonlinearSolve(PROBLEM& pde,
                                         const VECTOR &last_time_solution, 
                                         VECTOR &solution, 
                                         bool apply_boundary_values, 
                                         bool force_matrix_build,
                                         int priority, 
                                         std::string /*algo_level*/)
    {      
       
      bool build_matrix = force_matrix_build;
      VECTOR residual, time_residual, tmp_residual;
      VECTOR du;
      std::stringstream out;
      pde.GetOutputHandler()->InitNewtonOut(out);

      double res = 0.0;
      double firstres = 0.0;
      double lastres = 0.0;

      du.reinit(solution);
      residual.reinit(solution);
      time_residual.reinit(solution);
      tmp_residual.reinit(solution);
      
      //Transfer from previous timestep
      residual +=solution;
      // last_time_solution is very good starting value?
      solution = last_time_solution;
   
      if(apply_boundary_values)
      {
        GetIntegrator().ApplyInitialBoundaryValues(pde,solution);
      }
      
      // Righthandside for the current timestep f^{n+1}
      GetIntegrator().AddDomainData("last_time_solution",&last_time_solution);
      GetIntegrator().AddDomainData("last_newton_solution",&solution);
      pde.SetStepPart("New");
      GetIntegrator().ComputeNonlinearRhs(pde,tmp_residual);    
      tmp_residual *= -1;
      residual += tmp_residual;

      // Save the part of the residual which is independent of  u^{n+1}
      time_residual = residual;
      time_residual *=-1;

      // Calculate the "real" residual in the current timestep
      GetIntegrator().ComputeNonlinearLhs(pde,tmp_residual); // modi, new
      residual += tmp_residual;
               
      residual *=-1.; // due to A(U)(\psi) = - A(U)(du,\psi)
     
      pde.GetOutputHandler()->SetIterationNumber(0,"PDENewton");
      pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
  
      res = residual.linfty_norm();
      firstres = res;
      lastres = res;
      int iter=0;

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

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


          
      pde.GetOutputHandler()->Write(out,priority);
      while(res > nonlinear_global_tol_ && res > firstres * nonlinear_tol_)
      {
        iter++;
        if(iter > nonlinear_maxiter_)
        {
          throw DOpEIterationException("Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
        }
        
        pde.GetOutputHandler()->SetIterationNumber(iter,"PDENewton");   
        LINEARSOLVER::Solve(pde,GetIntegrator(),residual,du,build_matrix);
        build_matrix = false;   
        //Linesearch
        {
          solution += du;        
          GetIntegrator().ComputeNonlinearLhs(pde,residual);    
          residual -= time_residual; 
          residual *= -1.;
          pde.GetOutputHandler()->Write(residual,"Residual"+pde.GetType(),pde.GetDoFType());
          
          double newres = residual.linfty_norm();
          int lineiter=0;
          double rho = linesearch_rho_;
          double alpha=1;
          
          while(newres > res)    
          {
            out<<"\t\t\t Linesearch step: " <<lineiter<<"\t Residual (rel.): "
               <<pde.GetOutputHandler()->ZeroTolerance(newres/firstres, 1.0);
            pde.GetOutputHandler()->Write(out,priority+1);
            lineiter++;
            if(lineiter > line_maxiter_)
            {
              throw DOpEIterationException("Line-Iteration count exceeded bounds!","StatSolver::NonlinearSolve");
            }
            solution.add(alpha*(rho-1.),du);
            alpha*= rho;
            
            GetIntegrator().ComputeNonlinearLhs(pde,residual);
            residual -= time_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<<"\t\t 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_time_solution");
      GetIntegrator().DeleteDomainData("last_newton_solution");
    

      return build_matrix;        
    }

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

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

}
#endif