newton.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PDELAB_SOLVER_NEWTON_HH
#define DUNE_PDELAB_SOLVER_NEWTON_HH
 
#include <dune/common/exceptions.hh>
#include <dune/common/ios_state.hh>
 
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/solver/newtonerrors.hh>
#include <dune/pdelab/solver/linesearch.hh>
#include <dune/pdelab/solver/terminate.hh>
#include <dune/pdelab/solver/utility.hh>
 
namespace Dune::PDELab
{
  namespace Impl
  {
    // Some SFinae magic to execute setReuse(bool) on a backend
    template<typename T1, typename = void>
    struct HasSetReuse
      : std::false_type
    {};
 
    template<typename T>
    struct HasSetReuse<T, decltype(std::declval<T>().setReuse(true), void())>
      : std::true_type
    {};
 
    template<typename T>
    inline void setLinearSystemReuse(T& solver_backend, bool reuse, std::true_type)
    {
      if (!solver_backend.getReuse() && reuse)
        dwarn << "WARNING: Newton needed to override your choice to reuse the linear system in order to work!" << std::endl;
    solver_backend.setReuse(reuse);
    }
 
    template<typename T>
    inline void setLinearSystemReuse(T&, bool, std::false_type)
    {}
 
    template<typename T>
    inline void setLinearSystemReuse(T& solver_backend, bool reuse)
    {
      setLinearSystemReuse(solver_backend, reuse, HasSetReuse<T>());
    }
  }
 
 
  template <typename GridOperator_, typename LinearSolver_>
  class NewtonMethod
  {
  public:
    using GridOperator = GridOperator_;
 
    using LinearSolver = LinearSolver_;
 
    using Domain = typename GridOperator::Traits::Domain;
 
    using Range = typename GridOperator::Traits::Range;
 
    using Jacobian = typename GridOperator::Traits::Jacobian;
 
    using Real = typename Dune::FieldTraits<typename Domain::ElementType>::real_type;
 
    using Result = PDESolverResult<Real>;
 
    using LineSearch = LineSearchInterface<Domain>;
 
    Result& result()
    {
      if (not _resultValid)
        DUNE_THROW(NewtonError, "NewtonMethod::result() called before NewtonMethod::solve()");
      return _result;
    }
 
    virtual void prepareStep(Domain& solution)
    {
      _reassembled = false;
      if (_result.defect/_previousDefect > _reassembleThreshold){
        if (_hangingNodeModifications){
          auto dirichletValues = solution;
          // Set all non dirichlet values to zero
          Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
          // Set all constrained DOFs to zero in solution
          Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
          // Copy correct Dirichlet values back into solution vector
          Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
          // Interpolate periodic constraints / hanging nodes
          _gridOperator.localAssembler().backtransform(solution);
        }
        if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
          if (_verbosity>=3)
            std::cout << "      Reassembling matrix..." << std::endl;
          *_jacobian = Real(0.0);
          _gridOperator.jacobian(solution, *_jacobian);
          _reassembled = true;
        }
      }
 
      _linearReduction = _minLinearReduction;
      // corner case: ForceIteration==true. Use _minLinearReduction when initial defect is less than AbsoluteLimit.
      if (not _fixedLinearReduction && !(_result.iterations==0 && _result.first_defect<_absoluteLimit)){
        // Determine maximum defect, where Newton is converged.
        using std::min;
        using std::max;
        auto stop_defect = max(_result.first_defect*_reduction, _absoluteLimit);
 
        // 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.
        _linearReduction =
          max( stop_defect/(10*_result.defect),
            min(_minLinearReduction, _result.defect*_result.defect/(_previousDefect*_previousDefect)) );
        }
 
        if (_verbosity >= 3)
          std::cout << "      requested linear reduction:       "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _linearReduction << std::endl;
    }
 
    virtual void linearSolve()
    {
      if (_verbosity >= 4)
        std::cout << "      Solving linear system..." << std::endl;
 
      if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
        // If the jacobian was not reassembled we might save some work in the solver backend
        Impl::setLinearSystemReuse(_linearSolver, not _reassembled);
      }
 
      // Solve the linear system
      _correction = 0.0;
      if constexpr (linearSolverIsMatrixFree<LinearSolver>()){
        _linearSolver.apply(_correction, _residual, _linearReduction);
      }
      else{
        _linearSolver.apply(*_jacobian, _correction, _residual, _linearReduction);
      }
 
      if (not _linearSolver.result().converged)
        DUNE_THROW(NewtonLinearSolverError,
                   "NewtonMethod::linearSolve(): Linear solver did not converge "
                   "in " << _linearSolver.result().iterations << " iterations");
      if (_verbosity >= 4)
        std::cout << "          linear solver iterations:     "
                  << std::setw(12) << _linearSolver.result().iterations << std::endl
                  << "          linear defect reduction:      "
                  << std::setw(12) << std::setprecision(4) << std::scientific
                  << _linearSolver.result().reduction << std::endl;
    }
 
    virtual void apply(Domain& solution)
    {
      // Reset solver statistics
      _result.clear();
      _resultValid = true;
 
      // Store old ios flags (will be reset when this goes out of scope)
      ios_base_all_saver restorer(std::cout);
 
      // Prepare time measuring
      using Clock = std::chrono::steady_clock;
      using Duration = Clock::duration;
      auto assembler_time = Duration::zero();
      auto linear_solver_time = Duration::zero();
      auto line_search_time   = Duration::zero();
      auto to_seconds =
        [](Duration duration){
          return std::chrono::duration<double>(duration).count();
        };
      auto start_solve = Clock::now();
 
      //=========================
      // Calculate initial defect
      //=========================
      updateDefect(solution);
      _result.first_defect = _result.defect;
      _previousDefect = _result.defect;
 
      if (_verbosity >= 2)
        std::cout << "  Initial defect: "
                  << std::setw(12) << std::setprecision(4) << std::scientific
                  << _result.defect << std::endl;
 
      //==========================
      // Calculate Jacobian matrix
      //==========================
      if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
        if (not _jacobian)
          _jacobian = std::make_shared<Jacobian>(_gridOperator);
      }
 
      //=========================
      // Nonlinear iteration loop
      //=========================
      while (not _terminate->terminate()){
        if(_verbosity >= 3)
          std::cout << "  Newton iteration " << _result.iterations
                    << " --------------------------------" << std::endl;
 
        //=============
        // Prepare step
        //=============
        auto start = Clock::now();
        try{
          prepareStep(solution);
        }
        catch (...)
        {
          // Keep track of statistics when the method fails. We record
          // independently the time spent in non-converging attempts.
          // Check OneStepMethod to see how these data are propagated.
          auto end = Clock::now();
          assembler_time += end-start;
          _result.assembler_time = to_seconds(assembler_time);
          throw;
        }
        auto end = Clock::now();
        assembler_time += end -start;
        _result.assembler_time = to_seconds(assembler_time);
 
        // Store defect
        _previousDefect = _result.defect;
 
        //====================
        // Solve linear system
        //====================
        start = Clock::now();
 
        // Important: In the matrix-free case we need to set the current linearization point!
        if constexpr (linearSolverIsMatrixFree<LinearSolver>()){
          _linearSolver.setLinearizationPoint(solution);
        }
        try{
          linearSolve();
        }
        catch (...)
        {
          // Separately catch statistics for linear solver failures.
          end = Clock::now();
          linear_solver_time += end-start;
          _result.linear_solver_time = to_seconds(linear_solver_time);
          _result.linear_solver_iterations += _linearSolver.result().iterations;
          throw;
        }
        end = Clock::now();
        linear_solver_time += end -start;
        _result.linear_solver_time = to_seconds(linear_solver_time);
        _result.linear_solver_iterations += _linearSolver.result().iterations;
 
        //===================================
        // Do line search and update solution
        //===================================
        start = Clock::now();
        _lineSearch->lineSearch(solution, _correction);
        // lineSearch(solution);
        end = Clock::now();
        line_search_time += end -start;
 
        //========================================
        // Store statistics and create some output
        //========================================
        _result.reduction = _result.defect/_result.first_defect;
        _result.iterations++;
        _result.conv_rate = std::pow(_result.reduction, 1.0/_result.iterations);
 
        if (_verbosity >= 3)
          std::cout << "      linear solver time:               "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << to_seconds(linear_solver_time) << std::endl
                    << "      defect reduction (this iteration):"
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.defect/_previousDefect << std::endl
                    << "      defect reduction (total):         "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.reduction << std::endl
                    << "      new defect:                       "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.defect << std::endl;
        if (_verbosity == 2)
          std::cout << "  Newton iteration "
                    << std::setw(2)
                    << _result.iterations
                    << ".  New defect: "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.defect
                    << ".  Reduction (this): "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.defect/_previousDefect
                    << ".  Reduction (total): "
                    << std::setw(12) << std::setprecision(4) << std::scientific
                    << _result.reduction
                    << std::endl;
      } // end while loop of nonlinear iterations
 
      auto end_solve = Clock::now();
      auto solve_time = end_solve - start_solve;
      _result.elapsed = to_seconds(solve_time);
 
      if (_verbosity == 1)
        std::cout << "  Newton converged after "
                  << std::setw(2)
                  << _result.iterations
                  << " iterations.  Reduction: "
                  << std::setw(12) << std::setprecision(4) << std::scientific
                  << _result.reduction
                  << "   (" << std::setprecision(4) << _result.elapsed << "s)"
                  << std::endl;
 
      if constexpr (!linearSolverIsMatrixFree<LinearSolver>()){
        if (not _keepMatrix)
          _jacobian.reset();
      }
    }
 
    virtual void updateDefect(Domain& solution)
    {
      if (_hangingNodeModifications){
        auto dirichletValues = solution;
        // Set all non dirichlet values to zero
        Dune::PDELab::set_shifted_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, dirichletValues);
        // Set all constrained DOFs to zero in solution
        Dune::PDELab::set_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), 0.0, solution);
        // Copy correct Dirichlet values back into solution vector
        Dune::PDELab::copy_constrained_dofs(_gridOperator.localAssembler().trialConstraints(), dirichletValues, solution);
        // Interpolate periodic constraints / hanging nodes
        _gridOperator.localAssembler().backtransform(solution);
      }
 
      _residual = 0.0;
      _gridOperator.residual(solution, _residual);
 
      // Use the 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 (_useMaxNorm){
        auto rankMax = Backend::native(_residual).infinity_norm();
        _result.defect = _gridOperator.testGridFunctionSpace().gridView().comm().max(rankMax);
      }
      else
        _result.defect =  _linearSolver.norm(_residual);
    }
 
    void setVerbosityLevel(unsigned int verbosity)
    {
      if (_gridOperator.trialGridFunctionSpace().gridView().comm().rank()>0)
        _verbosity = 0;
      else
        _verbosity = verbosity;
    }
 
    unsigned int getVerbosityLevel() const
    {
      return _verbosity;
    }
 
    void setReduction(Real reduction)
    {
      _reduction = reduction;
    }
 
    Real getReduction() const
    {
      return _reduction;
    }
 
    void setAbsoluteLimit(Real absoluteLimit)
    {
      _absoluteLimit = absoluteLimit;
    }
 
    Real getAbsoluteLimit() const
    {
      return _absoluteLimit;
    }
 
    void setKeepMatrix(bool b)
    {
      _keepMatrix = b;
    }
 
    void setUseMaxNorm(bool b)
    {
      _useMaxNorm = b;
    }
 
    void setHangingNodeModifications(bool b)
    {
      _hangingNodeModifications = b;
    }
 
 
    bool keepMatrix() const
    {
      return _keepMatrix;
    }
 
    void discardMatrix()
    {
      if(_jacobian)
        _jacobian.reset();
    }
 
    void setMinLinearReduction(Real minLinearReduction)
    {
      _minLinearReduction = minLinearReduction;
    }
 
    void setFixedLinearReduction(bool fixedLinearReduction)
    {
      _fixedLinearReduction = fixedLinearReduction;
    }
 
    void setReassembleThreshold(Real reassembleThreshold)
    {
      _reassembleThreshold = reassembleThreshold;
    }
 
    void setParameters(const ParameterTree& parameterTree){
      _verbosity = parameterTree.get("VerbosityLevel", _verbosity);
      _reduction = parameterTree.get("Reduction", _reduction);
      _absoluteLimit = parameterTree.get("AbsoluteLimit", _absoluteLimit);
      _keepMatrix = parameterTree.get("KeepMatrix", _keepMatrix);
      _useMaxNorm = parameterTree.get("UseMaxNorm", _useMaxNorm);
      _hangingNodeModifications = parameterTree.get("HangingNodeModifications", _hangingNodeModifications);
      _minLinearReduction = parameterTree.get("MinLinearReduction", _minLinearReduction);
      _fixedLinearReduction = parameterTree.get("FixedLinearReduction", _fixedLinearReduction);
      _reassembleThreshold = parameterTree.get("ReassembleThreshold", _reassembleThreshold);
 
      // first create the linesearch, depending on the parameter
      std::string lineSearchStrategy = parameterTree.get("LineSearchStrategy","hackbuschReusken");
      auto strategy = lineSearchStrategyFromString(lineSearchStrategy);
      _lineSearch = createLineSearch(*this, strategy);
 
      // now set parameters
      if (parameterTree.hasSub("Terminate")){
        _terminate->setParameters(parameterTree.sub("Terminate"));
      }
      else{
        ParameterTree terminateTree;
        terminateTree["MaxIterations"] = std::to_string(parameterTree.get("MaxIterations", 40));
        terminateTree["ForceIteration"] = std::to_string(parameterTree.get("ForceIteration", false));
        _terminate->setParameters(terminateTree);
      }
      if (parameterTree.hasSub("LineSearch")){
        _lineSearch->setParameters(parameterTree.sub("LineSearch"));
      }
      else{
        ParameterTree lineSearchTree;
        lineSearchTree["MaxIterations"] = std::to_string(parameterTree.get("LineSearchMaxIterations", 10));
        lineSearchTree["DampingFactor"] = std::to_string(parameterTree.get("LineSearchDampingFactor", 0.5));
        lineSearchTree["AcceptBest"] = std::to_string(parameterTree.get("LineSearchAcceptBest", false));
        _lineSearch->setParameters(lineSearchTree);
      }
    }
 
    void setTerminate(std::shared_ptr<TerminateInterface> terminate)
    {
      _terminate = terminate;
    }
 
    std::shared_ptr<TerminateInterface> getTerminate() const
    {
      return _terminate;
    }
 
    void setLineSearch(std::shared_ptr<LineSearch> lineSearch)
    {
      _lineSearch = lineSearch;
    }
 
    std::shared_ptr<LineSearch> getLineSearch() const
    {
      return _lineSearch;
    }
 
    void printParameters(const std::string& _name = "NewtonMethod") const
    {
      // Change boolalpha flag to output true/false in full, not 0/1.
      // Restoring to previous setting is guaranteed.
      auto ioflags = std::cout.flags();
      std::cout.setf(std::ios_base::boolalpha);
      std::cout << _name << " parameters:\n";
      std::cout << "Verbosity............... " << _verbosity << std::endl;
      std::cout << "Reduction............... " << _reduction << std::endl;
      std::cout << "AbsoluteLimit........... " << _absoluteLimit << std::endl;
      std::cout << "KeepMatrix.............. " << _keepMatrix << std::endl;
      std::cout << "UseMaxNorm.............. " << _useMaxNorm << std::endl;
      std::cout << "MinLinearReduction...... " << _minLinearReduction << std::endl;
      std::cout << "FixedLinearReduction.... " << _fixedLinearReduction << std::endl;
      std::cout << "ReassembleThreshold..... " << _reassembleThreshold << std::endl;
      std::cout << "HangingNodeModifications " << _hangingNodeModifications << std::endl;
      try
      {
        _terminate->printParameters();
        _lineSearch->printParameters();
      }
      catch (...)
      {
        // guarantee restoring previous std::cout flags
        std::cout.flags(ioflags);
        throw;
      }
      std::cout.flags(ioflags);
    }
 
 
    NewtonMethod(
      const GridOperator& gridOperator,
      LinearSolver& linearSolver)
      : _gridOperator(gridOperator)
      , _linearSolver(linearSolver)
      , _residual(gridOperator.testGridFunctionSpace())
      , _correction(gridOperator.trialGridFunctionSpace())
    {
      _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
      _lineSearch = createLineSearch(*this, LineSearchStrategy::hackbuschReusken);
    }
 
    NewtonMethod(
      const GridOperator& gridOperator,
      LinearSolver& linearSolver,
      const ParameterTree& parameterTree)
      : _gridOperator(gridOperator)
      , _linearSolver(linearSolver)
      , _residual(gridOperator.testGridFunctionSpace())
      , _correction(gridOperator.trialGridFunctionSpace())
 
    {
      _terminate = std::make_shared<DefaultTerminate<NewtonMethod>> (*this);
      setParameters(parameterTree);
    }
 
    virtual ~NewtonMethod() {}
 
  private:
    const GridOperator& _gridOperator;
    LinearSolver& _linearSolver;
 
    // Vectors and Jacobi matrix we set up only once
    Range _residual;
    Domain _correction;
    std::shared_ptr<Jacobian> _jacobian;
 
    std::shared_ptr<TerminateInterface> _terminate;
    std::shared_ptr<LineSearch> _lineSearch;
 
    // Class for storing results
    Result _result;
    bool _resultValid = false; // result class only valid after calling apply
    Real _previousDefect = 0.0;
 
    // Remember if jacobian was reassembled in prepareStep
    bool _reassembled = true; // will be set in prepare step
    Real _linearReduction = 0.0; // will be set in prepare step
 
    // User parameters
    unsigned int _verbosity = 0;
    Real _reduction = 1e-8;
    Real _absoluteLimit = 1e-12;
    bool _keepMatrix = true;
    bool _useMaxNorm = false;
 
    // Special treatment if we have hanging nodes
    bool _hangingNodeModifications = false;
 
    // User parameters for prepareStep()
    Real _minLinearReduction = 1e-3;
    bool _fixedLinearReduction = false;
    Real _reassembleThreshold = 0.0;
  };
 
} // namespace Dune::PDELab
 
#endif