newton.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef PACXX_PROJECTSEMINAR_2019_SRC_NEWTON_HH
#define PACXX_PROJECTSEMINAR_2019_SRC_NEWTON_HH

// Shameless modified copy of PDELab::Newton

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <memory>

#include <dune/common/exceptions.hh>
#include <dune/common/ios_state.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/timer.hh>

#include <dune/istl/operators.hh>
#include <dune/istl/solver.hh>

#include "linearsolver.hh"

namespace PPS
{

  // Exception classes used in NewtonSolver
  class NewtonError : public Dune::Exception {};
  class NewtonDefectError : public NewtonError {};
  class NewtonLinearSolverError : public NewtonError {};
  class NewtonLineSearchError : public NewtonError {};
  class NewtonNotConverged : public NewtonError {};

  // Status information of Newton's method
  struct NewtonResult : Dune::InverseOperatorResult
  {
    double first_defect;       // the first defect
    double 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::field_type RFType;
    using Operator = Dune::LinearOperator<TrlV, TstV>;

  public:
    // export result type
    typedef NewtonResult Result;

    void setVerbosityLevel(unsigned int verbosity_level)
    {
      verbosity_level_ = verbosity_level;
    }

    //! Set whether the jacobian matrix should be kept across calls to apply().
    void setKeepMatrix(bool b)
    {
      keep_matrix_ = b;
    }

    //! Return whether the jacobian matrix is kept across calls to apply().
    bool keepMatrix() const
    {
      return keep_matrix_;
    }

    //! Discard the stored Jacobian matrix.
    void discardMatrix()
    {
      if(op_)
        op_.reset();
    }

  protected:
    const GridOperator& gridoperator_;
    TrialVector *u_;
    std::shared_ptr<TrialVector> z_;
    std::shared_ptr<TestVector> r_;
    std::shared_ptr<Operator> op_;
    Result res_;
    unsigned int verbosity_level_;
    RFType prev_defect_;
    RFType linear_reduction_;
    bool reassembled_;
    RFType reduction_;
    RFType abs_limit_;
    bool keep_matrix_;
    Dune::ParameterTree operator_params_;

    NewtonBase(const GridOperator& go, TrialVector& u)
      : gridoperator_(go)
      , u_(&u)
      , verbosity_level_(1)
      , keep_matrix_(true)
    { }

    NewtonBase(const GridOperator& go)
      : gridoperator_(go)
      , u_(0)
      , verbosity_level_(1)
      , keep_matrix_(true)
    { }

    virtual ~NewtonBase() { }

    virtual bool terminate() = 0;
    virtual void prepare_step(Operator& op, 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::field_type RFType;
    using Operator = Dune::LinearOperator<TrlV, TstV>;

  public:
    typedef NewtonResult Result;

    NewtonSolver(const GridOperator& go, TrialVector& u_, Solver& solver)
      : NewtonBase<GOS,TrlV,TstV>(go,u_)
      , solver_(solver)
      , result_valid_(false)
    {}

    NewtonSolver(const GridOperator& go, Solver& solver)
      : NewtonBase<GOS,TrlV,TstV>(go)
      , solver_(solver)
      , result_valid_(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_;
    }

  protected:
    virtual void defect(TestVector& r) override
    {
      r = 0.0;
      this->gridoperator_.residual(*this->u_, r);
      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(Operator& op, 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.
      assert(this->reassembled_);
      this->solver_.apply(op, z, r, this->linear_reduction_);

      Dune::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_;
  }; // 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;
    Dune::Timer timer;

    try
      {
        if(!this->r_) {
          // std::cout << "=== Setting up residual vector ..." << std::endl;
          this->r_ = std::make_shared<TstV>(this->gridoperator_.rangeSize());
        }
        // 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
            Dune::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->op_) {
          // std::cout << "==== Setting up jacobian matrix ... " << std::endl;
          this->op_ = makeOperator<TrlV, TstV>(&this->gridoperator_,
                                               this->operator_params_);
        }
        if(!this->z_) {
          // std::cout << "==== Setting up correction vector ... " << std::endl;
          this->z_ = std::make_shared<TrlV>(this->gridoperator_.domainSize());
        }

        while (!this->terminate())
          {
            if (this->verbosity_level_ >= 3)
              std::cout << "  Newton iteration " << this->res_.iterations
                        << " --------------------------------" << std::endl;

            Dune::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->op_,*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;

            Dune::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->op_, *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
            Dune::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();

    Dune::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->op_.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::field_type 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_;
    }

    void setOperatorParams(Dune::ParameterTree operator_params)
    {
      this->operator_params_ = std::move(operator_params);
    }

    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 >= int(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::field_type RFType;
    using Operator = Dune::LinearOperator<TrlV, TstV>;

  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)
    {}

    /**\brief set the minimal reduction in the linear solver

       \note with min_linear_reduction > 0, the linear reduction will be
       determined as mininum of the min_linear_reduction and the
       linear_reduction needed to achieve second order
       Newton convergence. */
    void setMinLinearReduction(RFType min_linear_reduction)
    {
      min_linear_reduction_ = min_linear_reduction;
    }

    /**\brief set a fixed reduction in the linear solver (overwrites setMinLinearReduction)

       \note with fixed_linear_reduction > 0, the linear reduction
       rate will always be fixed to min_linear_reduction. */
    void setFixedLinearReduction(bool fixed_linear_reduction)
    {
      fixed_linear_reduction_ = fixed_linear_reduction;
    }

    /**\brief set a threshold, when the linear operator is reassembled

       We allow to keep the linear operator over several newton
       iterations. If the reduction in the newton drops below a
       given threshold the linear operator is reassembled to ensure
       convergence.
     */
    void setReassembleThreshold(RFType reassemble_threshold)
    {
      reassemble_threshold_ = reassemble_threshold;
    }

    virtual void prepare_step(Operator& op, TstV& ) override
    {
      this->reassembled_ = false;
      if (this->res_.defect/this->prev_defect_ > reassemble_threshold_)
        {
          if (this->verbosity_level_ >= 3)
            std::cout << "      Reassembling matrix..." << std::endl;
          dynamic_cast<LinearizedOperatorMixin<TrlV>&>(*this->op_).
            linearizeAt(*this->u_);
          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;

      Dune::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_;
  }; // 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::field_type RFType;

  public:
    enum Strategy {
      /** \brief don't do any line search or damping */
      noLineSearch,
      /** \brief perform a linear search for the optimal damping parameter with multiples of damping

       the strategy was described in <a href="http://dx.doi.org/10.1007/BF01406516">[Hackbusch and Reusken, 1989]</a> */
      hackbuschReusken,
      /** \brief same as hackbuschReusken, but doesn't fail if the best update is still not good enough */
      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;
      Dune::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:
    /** helper function to get the different strategies from their name */
    Strategy strategyFromName(const std::string & s) {
      if (s == "noLineSearch")
        return noLineSearch;
      if (s == "hackbuschReusken")
        return hackbuschReusken;
      if (s == "hackbuschReuskenAcceptBest")
        return hackbuschReuskenAcceptBest;
      DUNE_THROW(Dune::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_)
      : 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_)
      : 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)
    {}
    //! interpret a parameter tree as a set of options for the newton solver
    /**

       example configuration:

       \code
       [NewtonParameters]

       ReassembleThreshold = 0.1
       LineSearchMaxIterations = 10
       MaxIterations = 7
       AbsoluteLimit = 1e-6
       Reduction = 1e-4
       MinLinearReduction = 1e-3
       LineSearchDamping  = 0.9
       \endcode

       and invocation in the code:
       \code
       newton.setParameters(param.sub("NewtonParameters"));
       \endcode
    */
    void setParameters(const Dune::ParameterTree & param)
    {
      typedef typename TstV::field_type 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.hasSub("Operator"))
        this->setOperatorParams(param.sub("Operator"));
    }
  }; // end class Newton
} // end namespace PPS

#endif // PACXX_PROJECTSEMINAR_2019_SRC_NEWTON_HH