newton.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
//
// Note: This Newton is deprecated and will be removed for PDELab 2.8. Please
// use the Newton from <dune/pdelab/solver/newton.hh>
//
#ifndef DUNE_PDELAB_NEWTON_NEWTON_HH
#define DUNE_PDELAB_NEWTON_NEWTON_HH
 
#include <iostream>
#include <iomanip>
#include <cmath>
#include <memory>
 
#include <type_traits>
 
#include <math.h>
 
#include <dune/common/stdstreams.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/ios_state.hh>
#include <dune/common/timer.hh>
#include <dune/common/parametertree.hh>
 
#include <dune/pdelab/backend/solver.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/solver/newton.hh>
 
 
// Note: Before introducing the new Newton there was no doxygen documentation
// for Newton at all. This means that it doesn't make any sense to have doxygen
// documentation for the old Newton now.
#ifndef DOXYGEN
 
namespace Dune
{
  namespace PDELab
  {
 
    class NewtonLineSearchError : public NewtonError {};
    class NewtonNotConverged : public NewtonError {};
 
    // Status information of Newton's method
    template<class RFType>
    struct NewtonResult : LinearSolverResult<RFType>
    {
      RFType first_defect;       // the first defect
      RFType defect;             // the final defect
      double assembler_time;     // Cumulative time for matrix assembly
      double linear_solver_time; // Cumulative time for linear solver
      int linear_solver_iterations; // Total number of linear iterations
 
      NewtonResult() :
        first_defect(0.0), defect(0.0), assembler_time(0.0), linear_solver_time(0.0),
        linear_solver_iterations(0) {}
    };
 
    template<class GOS, class TrlV, class TstV>
    class NewtonBase
    {
      typedef GOS GridOperator;
      typedef TrlV TrialVector;
      typedef TstV TestVector;
 
      typedef typename TestVector::ElementType RFType;
      typedef typename GOS::Traits::Jacobian Matrix;
 
 
    public:
      // export result type
      typedef NewtonResult<RFType> Result;
 
      void setVerbosityLevel(unsigned int verbosity_level)
      {
        if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
          verbosity_level_ = 0;
        else
          verbosity_level_ = verbosity_level;
      }
 
      void setKeepMatrix(bool b)
      {
        keep_matrix_ = b;
      }
 
      bool keepMatrix() const
      {
        return keep_matrix_;
      }
 
      void discardMatrix()
      {
        if(A_)
          A_.reset();
      }
 
    protected:
      const GridOperator& gridoperator_;
      TrialVector *u_;
      std::shared_ptr<TrialVector> z_;
      std::shared_ptr<TestVector> r_;
      std::shared_ptr<Matrix> A_;
      Result res_;
      unsigned int verbosity_level_;
      RFType prev_defect_;
      RFType linear_reduction_;
      bool reassembled_;
      RFType reduction_;
      RFType abs_limit_;
      bool keep_matrix_;
 
      NewtonBase(const GridOperator& go, TrialVector& u)
        : gridoperator_(go)
        , u_(&u)
        , verbosity_level_(1)
        , keep_matrix_(true)
      {
        if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
          verbosity_level_ = 0;
      }
 
      NewtonBase(const GridOperator& go)
        : gridoperator_(go)
        , u_(0)
        , verbosity_level_(1)
        , keep_matrix_(true)
      {
        if (gridoperator_.trialGridFunctionSpace().gridView().comm().rank()>0)
          verbosity_level_ = 0;
      }
 
      virtual ~NewtonBase() { }
 
      virtual bool terminate() = 0;
      virtual void prepare_step(Matrix& A, TestVector& r) = 0;
      virtual void line_search(TrialVector& z, TestVector& r) = 0;
      virtual void defect(TestVector& r) = 0;
    }; // end class NewtonBase
 
    template<class GOS, class S, class TrlV, class TstV>
    class NewtonSolver : public virtual NewtonBase<GOS,TrlV,TstV>
    {
      typedef S Solver;
      typedef GOS GridOperator;
      typedef TrlV TrialVector;
      typedef TstV TestVector;
 
      typedef typename TestVector::ElementType RFType;
      typedef typename GOS::Traits::Jacobian Matrix;
 
    public:
      typedef NewtonResult<RFType> Result;
 
      NewtonSolver(const GridOperator& go, TrialVector& u_, Solver& solver)
        : NewtonBase<GOS,TrlV,TstV>(go,u_)
        , solver_(solver)
        , result_valid_(false)
        , use_maxnorm_(false)
      {}
 
      NewtonSolver(const GridOperator& go, Solver& solver)
        : NewtonBase<GOS,TrlV,TstV>(go)
        , solver_(solver)
        , result_valid_(false)
        , use_maxnorm_(false)
        , hanging_node_modifications_(false)
      {}
 
      void apply();
 
      void apply(TrialVector& u_);
 
      const Result& result() const
      {
        if (!result_valid_)
          DUNE_THROW(NewtonError,
                     "NewtonSolver::result() called before NewtonSolver::solve()");
        return this->res_;
      }
 
      void setUseMaxNorm(bool b)
      {
        use_maxnorm_ = b;
      }
 
      void setHangingNodeModifications(bool b)
      {
        hanging_node_modifications_ = b;
      }
 
    protected:
      virtual void defect(TestVector& r) override
      {
        if (hanging_node_modifications_){
          auto dirichletValues = *this->u_;
          // Set all non dirichlet values to zero
          Dune::PDELab::set_shifted_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, dirichletValues);
          // Set all constrained DOFs to zero in solution
          Dune::PDELab::set_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, *this->u_);
          // Copy correct Dirichlet values back into solution vector
          Dune::PDELab::copy_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), dirichletValues, *this->u_);
          // Interpolate periodic constraints / hanging nodes
          this->gridoperator_.localAssembler().backtransform(*this->u_);
        }
 
        r = 0.0;
        this->gridoperator_.residual(*this->u_, r);
        // use maximum norm as a stopping criterion, this helps loosen the tolerance
        // when solving for stationary solutions of nonlinear time-dependent problems
        // the default is to use the Euclidean norm (in the else-block) as before
        if(use_maxnorm_)
        {
          auto rank_max = Backend::native(r).infinity_norm();
          this->res_.defect = this->gridoperator_.testGridFunctionSpace().gridView().comm().max(rank_max);
        }
        else
        {
          this->res_.defect = this->solver_.norm(r);
        }
        if (!std::isfinite(this->res_.defect))
          DUNE_THROW(NewtonDefectError,
                     "NewtonSolver::defect(): Non-linear defect is NaN or Inf");
      }
 
 
    private:
      void linearSolve(Matrix& A, TrialVector& z, TestVector& r) const
      {
        if (this->verbosity_level_ >= 4)
          std::cout << "      Solving linear system..." << std::endl;
        z = 0.0;
        // If possible tell solver backend to reuse linear system when we did not reassemble.
        Impl::setLinearSystemReuse(this->solver_, !this->reassembled_);
        this->solver_.apply(A, z, r, this->linear_reduction_);
 
        ios_base_all_saver restorer(std::cout); // store old ios flags
 
        if (!this->solver_.result().converged)
          DUNE_THROW(NewtonLinearSolverError,
                     "NewtonSolver::linearSolve(): Linear solver did not converge "
                     "in " << this->solver_.result().iterations << " iterations");
        if (this->verbosity_level_ >= 4)
          std::cout << "          linear solver iterations:     "
                    << std::setw(12) << solver_.result().iterations << std::endl
                    << "          linear defect reduction:      "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << solver_.result().reduction << std::endl;
      }
 
      Solver& solver_;
      bool result_valid_;
      bool use_maxnorm_;
      bool hanging_node_modifications_;
    }; // end class NewtonSolver
 
    template<class GOS, class S, class TrlV, class TstV>
    void NewtonSolver<GOS,S,TrlV,TstV>::apply(TrialVector& u)
    {
      this->u_ = &u;
      apply();
    }
 
    template<class GOS, class S, class TrlV, class TstV>
    void NewtonSolver<GOS,S,TrlV,TstV>::apply()
    {
      this->res_.iterations = 0;
      this->res_.converged = false;
      this->res_.reduction = 1.0;
      this->res_.conv_rate = 1.0;
      this->res_.elapsed = 0.0;
      this->res_.assembler_time = 0.0;
      this->res_.linear_solver_time = 0.0;
      this->res_.linear_solver_iterations = 0;
      result_valid_ = true;
      Timer timer;
 
      try
        {
          if(!this->r_) {
            // std::cout << "=== Setting up residual vector ..." << std::endl;
            this->r_ = std::make_shared<TestVector>(this->gridoperator_.testGridFunctionSpace());
          }
          // residual calculation in member function "defect":
          //--------------------------------------------------
          // - set residual vector to zero
          // - calculate new residual
          // - store norm of residual in "this->res_.defect"
          this->defect(*this->r_);
          this->res_.first_defect = this->res_.defect;
          this->prev_defect_ = this->res_.defect;
 
          if (this->verbosity_level_ >= 2)
            {
              // store old ios flags
              ios_base_all_saver restorer(std::cout);
              std::cout << "  Initial defect: "
                        << std::setw(12) << std::setprecision(4) << std::scientific
                        << this->res_.defect << std::endl;
            }
 
          if(!this->A_) {
            // std::cout << "==== Setting up jacobian matrix ... " << std::endl;
            this->A_ = std::make_shared<Matrix>(this->gridoperator_);
          }
          if(!this->z_) {
            // std::cout << "==== Setting up correction vector ... " << std::endl;
            this->z_ = std::make_shared<TrialVector>(this->gridoperator_.trialGridFunctionSpace());
          }
 
          while (!this->terminate())
            {
              if (this->verbosity_level_ >= 3)
                std::cout << "  Newton iteration " << this->res_.iterations
                          << " --------------------------------" << std::endl;
 
              Timer assembler_timer;
              try
                {
                  // jacobian calculation in member function "prepare_step"
                  //-------------------------------------------------------
                  // - if above reassemble threshold
                  //   - set jacobian to zero
                  //   - calculate new jacobian
                  // - set linear reduction
                  this->prepare_step(*this->A_,*this->r_);
                }
              catch (...)
                {
                  this->res_.assembler_time += assembler_timer.elapsed();
                  throw;
                }
              double assembler_time = assembler_timer.elapsed();
              this->res_.assembler_time += assembler_time;
              if (this->verbosity_level_ >= 3)
                std::cout << "      matrix assembly time:             "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << assembler_time << std::endl;
 
              Timer linear_solver_timer;
              try
                {
                  // solution of linear system in member function "linearSolve"
                  //-----------------------------------------------------------
                  // - set initial guess for correction z to zero
                  // - call linear solver
                  this->linearSolve(*this->A_, *this->z_, *this->r_);
                }
              catch (...)
                {
                  this->res_.linear_solver_time += linear_solver_timer.elapsed();
                  this->res_.linear_solver_iterations += this->solver_.result().iterations;
                  throw;
                }
              double linear_solver_time = linear_solver_timer.elapsed();
              this->res_.linear_solver_time += linear_solver_time;
              this->res_.linear_solver_iterations += this->solver_.result().iterations;
 
              try
                {
                  // line search with correction z
                  // the undamped version is also integrated in here
                  this->line_search(*this->z_, *this->r_);
                }
              catch (NewtonLineSearchError&)
                {
                  if (this->reassembled_)
                    throw;
                  if (this->verbosity_level_ >= 3)
                    std::cout << "      line search failed - trying again with reassembled matrix" << std::endl;
                  continue;
                }
 
              this->res_.reduction = this->res_.defect/this->res_.first_defect;
              this->res_.iterations++;
              this->res_.conv_rate = std::pow(this->res_.reduction, 1.0/this->res_.iterations);
 
              // store old ios flags
              ios_base_all_saver restorer(std::cout);
 
              if (this->verbosity_level_ >= 3)
                std::cout << "      linear solver time:               "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << linear_solver_time << std::endl
                          << "      defect reduction (this iteration):"
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.defect/this->prev_defect_ << std::endl
                          << "      defect reduction (total):         "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.reduction << std::endl
                          << "      new defect:                       "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.defect << std::endl;
              if (this->verbosity_level_ == 2)
                std::cout << "  Newton iteration " << std::setw(2) << this->res_.iterations
                          << ".  New defect: "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.defect
                          << ".  Reduction (this): "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.defect/this->prev_defect_
                          << ".  Reduction (total): "
                          << std::setw(12) << std::setprecision(4) << std::scientific
                          << this->res_.reduction << std::endl;
            } // end while
        } // end try
      catch(...)
        {
          this->res_.elapsed = timer.elapsed();
          throw;
        }
      this->res_.elapsed = timer.elapsed();
 
      ios_base_all_saver restorer(std::cout); // store old ios flags
 
      if (this->verbosity_level_ == 1)
        std::cout << "  Newton converged after " << std::setw(2) << this->res_.iterations
                  << " iterations.  Reduction: "
                  << std::setw(12) << std::setprecision(4) << std::scientific
                  << this->res_.reduction
                  << "   (" << std::setprecision(4) << this->res_.elapsed << "s)"
                  << std::endl;
 
      if(!this->keep_matrix_)
        this->A_.reset();
    } // end apply
 
    template<class GOS, class TrlV, class TstV>
    class NewtonTerminate : public virtual NewtonBase<GOS,TrlV,TstV>
    {
      typedef GOS GridOperator;
      typedef TrlV TrialVector;
 
      typedef typename TstV::ElementType RFType;
 
    public:
      NewtonTerminate(const GridOperator& go, TrialVector& u_)
        : NewtonBase<GOS,TrlV,TstV>(go,u_)
        , maxit_(40)
        , force_iteration_(false)
      {
        this->reduction_ = 1e-8;
        this->abs_limit_ = 1e-12;
      }
 
      NewtonTerminate(const GridOperator& go)
        : NewtonBase<GOS,TrlV,TstV>(go)
        , maxit_(40)
        , force_iteration_(false)
      {
        this->reduction_ = 1e-8;
        this->abs_limit_ = 1e-12;
      }
 
      void setReduction(RFType reduction)
      {
        this->reduction_ = reduction;
      }
 
      void setMaxIterations(unsigned int maxit)
      {
        maxit_ = maxit;
      }
 
      void setForceIteration(bool force_iteration)
      {
        force_iteration_ = force_iteration;
      }
 
      void setAbsoluteLimit(RFType abs_limit_)
      {
        this->abs_limit_ = abs_limit_;
      }
 
      virtual bool terminate() override
      {
        if (force_iteration_ && this->res_.iterations == 0)
          return false;
        this->res_.converged = this->res_.defect < this->abs_limit_
          || this->res_.defect < this->res_.first_defect * this->reduction_;
        if (this->res_.iterations >= maxit_ && !this->res_.converged)
          DUNE_THROW(NewtonNotConverged,
                     "NewtonTerminate::terminate(): Maximum iteration count reached");
        return this->res_.converged;
      }
 
    private:
      unsigned int maxit_;
      bool force_iteration_;
    }; // end class NewtonTerminate
 
    template<class GOS, class TrlV, class TstV>
    class NewtonPrepareStep : public virtual NewtonBase<GOS,TrlV,TstV>
    {
      typedef GOS GridOperator;
      typedef TrlV TrialVector;
 
      typedef typename TstV::ElementType RFType;
      typedef typename GOS::Traits::Jacobian Matrix;
 
    public:
      NewtonPrepareStep(const GridOperator& go, TrialVector& u_)
        : NewtonBase<GOS,TrlV,TstV>(go,u_)
        , min_linear_reduction_(1e-3)
        , fixed_linear_reduction_(0.0)
        , reassemble_threshold_(0.0)
      {}
 
      NewtonPrepareStep(const GridOperator& go)
        : NewtonBase<GOS,TrlV,TstV>(go)
        , min_linear_reduction_(1e-3)
        , fixed_linear_reduction_(0.0)
        , reassemble_threshold_(0.0)
        , hanging_node_modifications_(false)
      {}
 
      void setMinLinearReduction(RFType min_linear_reduction)
      {
        min_linear_reduction_ = min_linear_reduction;
      }
 
      void setFixedLinearReduction(bool fixed_linear_reduction)
      {
        fixed_linear_reduction_ = fixed_linear_reduction;
      }
 
      void setReassembleThreshold(RFType reassemble_threshold)
      {
        reassemble_threshold_ = reassemble_threshold;
      }
 
      void setHangingNodeModifications(bool b)
      {
        hanging_node_modifications_ = b;
      }
 
      virtual void prepare_step(Matrix& A, TstV& ) override
      {
        this->reassembled_ = false;
        if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
          {
            if (hanging_node_modifications_){
              auto dirichletValues = *this->u_;
              // Set all non dirichlet values to zero
              Dune::PDELab::set_shifted_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, dirichletValues);
              // Set all constrained DOFs to zero in solution
              Dune::PDELab::set_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), 0.0, *this->u_);
              // Copy correct Dirichlet values back into solution vector
              Dune::PDELab::copy_constrained_dofs(this->gridoperator_.localAssembler().trialConstraints(), dirichletValues, *this->u_);
              // Interpolate periodic constraints / hanging nodes
              this->gridoperator_.localAssembler().backtransform(*this->u_);
            }
 
            if (this->verbosity_level_ >= 3)
              std::cout << "      Reassembling matrix..." << std::endl;
            A = 0.0;
            this->gridoperator_.jacobian(*this->u_, A);
            this->reassembled_ = true;
          }
 
        if (fixed_linear_reduction_ == true)
          this->linear_reduction_ = min_linear_reduction_;
        else {
          // determine maximum defect, where Newton is converged.
          RFType stop_defect =
            std::max(this->res_.first_defect * this->reduction_,
                     this->abs_limit_);
 
          /*
            To achieve second order convergence of newton
            we need a linear reduction of at least
            current_defect^2/prev_defect^2.
            For the last newton step a linear reduction of
            1/10*end_defect/current_defect
            is sufficient for convergence.
          */
          if ( stop_defect/(10*this->res_.defect) >
               this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_) )
            this->linear_reduction_ =
              stop_defect/(10*this->res_.defect);
          else
            this->linear_reduction_ =
              std::min(min_linear_reduction_,this->res_.defect*this->res_.defect/(this->prev_defect_*this->prev_defect_));
        }
 
        this->prev_defect_ = this->res_.defect;
 
        ios_base_all_saver restorer(std::cout); // store old ios flags
 
        if (this->verbosity_level_ >= 3)
          std::cout << "      requested linear reduction:       "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << this->linear_reduction_ << std::endl;
      }
 
    private:
      RFType min_linear_reduction_;
      bool fixed_linear_reduction_;
      RFType reassemble_threshold_;
      bool hanging_node_modifications_;
    }; // end class NewtonPrepareStep
 
    template<class GOS, class TrlV, class TstV>
    class NewtonLineSearch : public virtual NewtonBase<GOS,TrlV,TstV>
    {
      typedef GOS GridOperator;
      typedef TrlV TrialVector;
      typedef TstV TestVector;
 
      typedef typename TestVector::ElementType RFType;
 
    public:
      enum Strategy {
        noLineSearch,
        hackbuschReusken,
        hackbuschReuskenAcceptBest };
 
      NewtonLineSearch(const GridOperator& go, TrialVector& u_)
        : NewtonBase<GOS,TrlV,TstV>(go,u_)
        , strategy_(hackbuschReusken)
        , maxit_(10)
        , damping_factor_(0.5)
      {}
 
      NewtonLineSearch(const GridOperator& go)
        : NewtonBase<GOS,TrlV,TstV>(go)
        , strategy_(hackbuschReusken)
        , maxit_(10)
        , damping_factor_(0.5)
      {}
 
      void setLineSearchStrategy(Strategy strategy)
      {
        strategy_ = strategy;
      }
 
      void setLineSearchStrategy(std::string strategy)
      {
        strategy_ = strategyFromName(strategy);
      }
 
      void setLineSearchMaxIterations(unsigned int maxit)
      {
        maxit_ = maxit;
      }
 
      void setLineSearchDampingFactor(RFType damping_factor)
      {
        damping_factor_ = damping_factor;
      }
 
      virtual void line_search(TrialVector& z, TestVector& r) override
      {
        if ((strategy_ == noLineSearch) || (this->res_.defect < this->abs_limit_))
          {
            this->u_->axpy(-1.0, z);
            this->defect(r);
            return;
          }
 
        if (this->verbosity_level_ >= 4)
          std::cout << "      Performing line search..." << std::endl;
        RFType lambda = 1.0;
        RFType best_lambda = 0.0;
        RFType best_defect = this->res_.defect;
        TrialVector prev_u(*this->u_);
        unsigned int i = 0;
        ios_base_all_saver restorer(std::cout); // store old ios flags
 
        while (1)
          {
            if (this->verbosity_level_ >= 4)
              std::cout << "          trying line search damping factor:   "
                        << std::setw(12) << std::setprecision(4) << std::scientific
                        << lambda
                        << std::endl;
 
            this->u_->axpy(-lambda, z);
            try {
              this->defect(r);
            }
             catch (NewtonDefectError&)
              {
                if (this->verbosity_level_ >= 4)
                  std::cout << "          NaNs detected" << std::endl;
              }       // ignore NaNs and try again with lower lambda
 
            if (this->res_.defect <= (1.0 - lambda/4) * this->prev_defect_)
              {
                if (this->verbosity_level_ >= 4)
                  std::cout << "          line search converged" << std::endl;
                break;
              }
 
            if (this->res_.defect < best_defect)
              {
                best_defect = this->res_.defect;
                best_lambda = lambda;
              }
 
            if (++i >= maxit_)
              {
                if (this->verbosity_level_ >= 4)
                  std::cout << "          max line search iterations exceeded" << std::endl;
                switch (strategy_)
                  {
                  case hackbuschReusken:
                    *this->u_ = prev_u;
                    this->defect(r);
                    DUNE_THROW(NewtonLineSearchError,
                               "NewtonLineSearch::line_search(): line search failed, "
                               "max iteration count reached, "
                               "defect did not improve enough");
                  case hackbuschReuskenAcceptBest:
                    if (best_lambda == 0.0)
                      {
                        *this->u_ = prev_u;
                        this->defect(r);
                        DUNE_THROW(NewtonLineSearchError,
                                   "NewtonLineSearch::line_search(): line search failed, "
                                   "max iteration count reached, "
                                   "defect did not improve in any of the iterations");
                      }
                    if (best_lambda != lambda)
                      {
                        *this->u_ = prev_u;
                        this->u_->axpy(-best_lambda, z);
                        this->defect(r);
                      }
                    break;
                  case noLineSearch:
                    break;
                  }
                break;
              }
 
            lambda *= damping_factor_;
            *this->u_ = prev_u;
          }
        if (this->verbosity_level_ >= 4)
          std::cout << "          line search damping factor:   "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << lambda << std::endl;
      } // end line_search
 
    protected:
      Strategy strategyFromName(const std::string & s) {
        if (s == "noLineSearch")
          return noLineSearch;
        if (s == "hackbuschReusken")
          return hackbuschReusken;
        if (s == "hackbuschReuskenAcceptBest")
          return hackbuschReuskenAcceptBest;
        DUNE_THROW(Exception, "unknown line search strategy" << s);
      }
 
    private:
      Strategy strategy_;
      unsigned int maxit_;
      RFType damping_factor_;
    }; // end class NewtonLineSearch
 
    template<class GOS, class S, class TrlV, class TstV = TrlV>
    class Newton : public NewtonSolver<GOS,S,TrlV,TstV>
                    , public NewtonTerminate<GOS,TrlV,TstV>
                    , public NewtonLineSearch<GOS,TrlV,TstV>
                    , public NewtonPrepareStep<GOS,TrlV,TstV>
    {
      typedef GOS GridOperator;
      typedef S Solver;
      typedef TrlV TrialVector;
 
    public:
      Newton(const GridOperator& go, TrialVector& u_, Solver& solver_)
        DUNE_DEPRECATED_MSG("This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")
        : NewtonBase<GOS,TrlV,TstV>(go,u_)
        , NewtonSolver<GOS,S,TrlV,TstV>(go,u_,solver_)
        , NewtonTerminate<GOS,TrlV,TstV>(go,u_)
        , NewtonLineSearch<GOS,TrlV,TstV>(go,u_)
        , NewtonPrepareStep<GOS,TrlV,TstV>(go,u_)
      {}
      Newton(const GridOperator& go, Solver& solver_)
        DUNE_DEPRECATED_MSG("This Newton is deprecated. Use NewtonMethod from dune/pdelab/solver/newton.hh instead.")
        : NewtonBase<GOS,TrlV,TstV>(go)
        , NewtonSolver<GOS,S,TrlV,TstV>(go,solver_)
        , NewtonTerminate<GOS,TrlV,TstV>(go)
        , NewtonLineSearch<GOS,TrlV,TstV>(go)
        , NewtonPrepareStep<GOS,TrlV,TstV>(go)
      {}
 
      void setParameters(const Dune::ParameterTree & param)
      {
        typedef typename TstV::ElementType RFType;
        if (param.hasKey("VerbosityLevel"))
          this->setVerbosityLevel(
            param.get<unsigned int>("VerbosityLevel"));
        if (param.hasKey("Reduction"))
          this->setReduction(
            param.get<RFType>("Reduction"));
        if (param.hasKey("MaxIterations"))
          this->setMaxIterations(
            param.get<unsigned int>("MaxIterations"));
        if (param.hasKey("ForceIteration"))
          this->setForceIteration(
            param.get<bool>("ForceIteration"));
        if (param.hasKey("AbsoluteLimit"))
          this->setAbsoluteLimit(
            param.get<RFType>("AbsoluteLimit"));
        if (param.hasKey("MinLinearReduction"))
          this->setMinLinearReduction(
            param.get<RFType>("MinLinearReduction"));
        if (param.hasKey("FixedLinearReduction"))
          this->setFixedLinearReduction(
            param.get<bool>("FixedLinearReduction"));
        if (param.hasKey("ReassembleThreshold"))
          this->setReassembleThreshold(
            param.get<RFType>("ReassembleThreshold"));
        if (param.hasKey("LineSearchStrategy"))
          this->setLineSearchStrategy(
            param.get<std::string>("LineSearchStrategy"));
        if (param.hasKey("LineSearchMaxIterations"))
          this->setLineSearchMaxIterations(
            param.get<unsigned int>("LineSearchMaxIterations"));
        if (param.hasKey("LineSearchDampingFactor"))
          this->setLineSearchDampingFactor(
            param.get<RFType>("LineSearchDampingFactor"));
        if (param.hasKey("KeepMatrix"))
          this->setKeepMatrix(
            param.get<bool>("KeepMatrix"));
        if (param.hasKey("UseMaxNorm"))
          this->setUseMaxNorm(
            param.get<bool>("UseMaxNorm"));
      }
    }; // end class Newton
  } // end namespace PDELab
} // end namespace Dune
 
#endif // DOXYGEN
 
#endif // DUNE_PDELAB_NEWTON_NEWTON_HH