nonlinear/newtonsolver.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_NEWTON_SOLVER_HH
#define DUMUX_NEWTON_SOLVER_HH
 
#include <cmath>
#include <memory>
#include <iostream>
#include <type_traits>
 
#include <dune/common/timer.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpicollectivecommunication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/multitypeblockvector.hh>
 
#include <dumux/common/parameters.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/typetraits/vector.hh>
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/pdesolver.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
#include <dumux/assembly/partialreassembler.hh>
 
#include "newtonconvergencewriter.hh"
 
namespace Dumux {
namespace Detail {
 
struct supportsPartialReassembly
{
    template<class Assembler>
    auto operator()(Assembler&& a)
    -> decltype(a.assembleJacobianAndResidual(std::declval<const typename Assembler::ResidualType&>(),
                                              std::declval<const PartialReassembler<Assembler>*>()))
    {}
};
 
template<class Assembler>
using DetectPVSwitch = typename Assembler::GridVariables::VolumeVariables::PrimaryVariableSwitch;
 
template<class Assembler>
using GetPVSwitch = Dune::Std::detected_or<int, DetectPVSwitch, Assembler>;
 
} // end namespace Detail
 
template <class Assembler, class LinearSolver,
          class Reassembler = PartialReassembler<Assembler>,
          class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
class NewtonSolver : public PDESolver<Assembler, LinearSolver>
{
    using ParentType = PDESolver<Assembler, LinearSolver>;
    using Scalar = typename Assembler::Scalar;
    using JacobianMatrix = typename Assembler::JacobianMatrix;
    using SolutionVector = typename Assembler::ResidualType;
    using ConvergenceWriter = ConvergenceWriterInterface<SolutionVector>;
    using TimeLoop = TimeLoopBase<Scalar>;
 
    using PrimaryVariableSwitch = typename Detail::GetPVSwitch<Assembler>::type;
    using HasPriVarsSwitch = typename Detail::GetPVSwitch<Assembler>::value_t; // std::true_type or std::false_type
    static constexpr bool hasPriVarsSwitch() { return HasPriVarsSwitch::value; };
 
public:
 
    using Communication = Comm;
 
    NewtonSolver(std::shared_ptr<Assembler> assembler,
                 std::shared_ptr<LinearSolver> linearSolver,
                 const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
                 const std::string& paramGroup = "")
    : ParentType(assembler, linearSolver)
    , endIterMsgStream_(std::ostringstream::out)
    , comm_(comm)
    , paramGroup_(paramGroup)
    {
        initParams_(paramGroup);
 
        // set the linear system (matrix & residual) in the assembler
        this->assembler().setLinearSystem();
 
        // set a different default for the linear solver residual reduction
        // within the Newton the linear solver doesn't need to solve too exact
        this->linearSolver().setResidualReduction(getParamFromGroup<Scalar>(paramGroup, "LinearSolver.ResidualReduction", 1e-6));
 
        // initialize the partial reassembler
        if (enablePartialReassembly_)
            partialReassembler_ = std::make_unique<Reassembler>(this->assembler());
 
        if (hasPriVarsSwitch())
        {
            const int priVarSwitchVerbosity = getParamFromGroup<int>(paramGroup, "PrimaryVariableSwitch.Verbosity", 1);
            priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(priVarSwitchVerbosity);
        }
    }
 
    const Communication& comm() const
    { return comm_; }
 
    void setMaxRelativeShift(Scalar tolerance)
    { shiftTolerance_ = tolerance; }
 
    void setMaxAbsoluteResidual(Scalar tolerance)
    { residualTolerance_ = tolerance; }
 
    void setResidualReduction(Scalar tolerance)
    { reductionTolerance_ = tolerance; }
 
    void setTargetSteps(int targetSteps)
    { targetSteps_ = targetSteps; }
 
    void setMinSteps(int minSteps)
    { minSteps_ = minSteps; }
 
    void setMaxSteps(int maxSteps)
    { maxSteps_ = maxSteps; }
 
    [[deprecated("Use attachConvergenceWriter(convWriter) and solve(x, *timeLoop) instead")]]
    void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop,
               std::shared_ptr<ConvergenceWriter> convWriter)
    {
        if (!convergenceWriter_)
            attachConvergenceWriter(convWriter);
 
        solve(uCurrentIter, timeLoop);
    }
 
    void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop) override
    {
        if (this->assembler().isStationaryProblem())
            DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!");
 
        // try solving the non-linear system
        for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
        {
            // linearize & solve
            const bool converged = solve_(uCurrentIter);
 
            if (converged)
                return;
 
            else if (!converged && i < maxTimeStepDivisions_)
            {
                // set solution to previous solution
                uCurrentIter = this->assembler().prevSol();
 
                // reset the grid variables to the previous solution
                this->assembler().resetTimeStep(uCurrentIter);
 
                if (verbosity_ >= 1)
                    std::cout << "Newton solver did not converge with dt = "
                              << timeLoop.timeStepSize() << " seconds. Retrying with time step of "
                              << timeLoop.timeStepSize() * retryTimeStepReductionFactor_ << " seconds\n";
 
                // try again with dt = dt * retryTimeStepReductionFactor_
                timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_);
            }
 
            else
            {
                DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
                                             << maxTimeStepDivisions_ << " time-step divisions. dt="
                                             << timeLoop.timeStepSize() << '\n');
            }
        }
    }
 
    void solve(SolutionVector& uCurrentIter) override
    {
        const bool converged = solve_(uCurrentIter);
        if (!converged)
            DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
                                         << numSteps_ << " iterations.\n");
    }
 
    virtual void newtonBegin(SolutionVector& u)
    {
        numSteps_ = 0;
        initPriVarSwitch_(u, HasPriVarsSwitch{});
 
        // write the initial residual if a convergence writer was set
        if (convergenceWriter_)
        {
            this->assembler().assembleResidual(u);
            SolutionVector delta(u);
            delta = 0; // dummy vector, there is no delta before solving the linear system
            convergenceWriter_->write(u, delta, this->assembler().residual());
        }
    }
 
    virtual bool newtonProceed(const SolutionVector &uCurrentIter, bool converged)
    {
        if (numSteps_ < minSteps_)
            return true;
        else if (converged)
            return false; // we are below the desired tolerance
        else if (numSteps_ >= maxSteps_)
        {
            // We have exceeded the allowed number of steps. If the
            // maximum relative shift was reduced by a factor of at least 4,
            // we proceed even if we are above the maximum number of steps.
            if (enableShiftCriterion_)
                return shift_*4.0 < lastShift_;
            else
                return reduction_*4.0 < lastReduction_;
        }
 
        return true;
    }
 
    virtual void newtonBeginStep(const SolutionVector& u)
    {
        lastShift_ = shift_;
        if (numSteps_ == 0)
        {
            lastReduction_ = 1.0;
        }
        else
        {
            lastReduction_ = reduction_;
        }
    }
 
    virtual void assembleLinearSystem(const SolutionVector& uCurrentIter)
    {
        assembleLinearSystem_(this->assembler(), uCurrentIter);
 
        if (enablePartialReassembly_)
            partialReassembler_->report(comm_, endIterMsgStream_);
    }
 
    void solveLinearSystem(SolutionVector& deltaU)
    {
        auto& b = this->assembler().residual();
 
        try
        {
            if (numSteps_ == 0)
            {
                Scalar norm2 = b.two_norm2();
                if (comm_.size() > 1)
                    norm2 = comm_.sum(norm2);
 
                using std::sqrt;
                initialResidual_ = sqrt(norm2);
            }
 
            // solve by calling the appropriate implementation depending on whether the linear solver
            // is capable of handling MultiType matrices or not
            bool converged = solveLinearSystem_(deltaU);
 
            // make sure all processes converged
            int convergedRemote = converged;
            if (comm_.size() > 1)
                convergedRemote = comm_.min(converged);
 
            if (!converged) {
                DUNE_THROW(NumericalProblem,
                           "Linear solver did not converge");
            }
            else if (!convergedRemote) {
                DUNE_THROW(NumericalProblem,
                           "Linear solver did not converge on a remote process");
            }
        }
        catch (const Dune::Exception &e) {
            // make sure all processes converged
            int converged = 0;
            if (comm_.size() > 1)
                converged = comm_.min(converged);
 
            NumericalProblem p;
            p.message(e.what());
            throw p;
        }
    }
 
    void newtonUpdate(SolutionVector &uCurrentIter,
                      const SolutionVector &uLastIter,
                      const SolutionVector &deltaU)
    {
        if (enableShiftCriterion_ || enablePartialReassembly_)
            newtonUpdateShift_(uLastIter, deltaU);
 
        if (enablePartialReassembly_) {
            // Determine the threshold 'eps' that is used for the partial reassembly.
            // Every entity where the primary variables exhibit a relative shift
            // summed up since the last linearization above 'eps' will be colored
            // red yielding a reassembly.
            // The user can provide three parameters to influence the threshold:
            // 'minEps' by 'Newton.ReassemblyMinThreshold' (1e-1*shiftTolerance_ default)
            // 'maxEps' by 'Newton.ReassemblyMaxThreshold' (1e2*shiftTolerance_ default)
            // 'omega'  by 'Newton.ReassemblyShiftWeight'  (1e-3 default)
            // The threshold is calculated from the currently achieved maximum
            // relative shift according to the formula
            // eps = max( minEps, min(maxEps, omega*shift) ).
            // Increasing/decreasing 'minEps' leads to less/more reassembly if
            // 'omega*shift' is small, i.e., for the last Newton iterations.
            // Increasing/decreasing 'maxEps' leads to less/more reassembly if
            // 'omega*shift' is large, i.e., for the first Newton iterations.
            // Increasing/decreasing 'omega' leads to more/less first and last
            // iterations in this sense.
            using std::max;
            using std::min;
            auto reassemblyThreshold = max(reassemblyMinThreshold_,
                                           min(reassemblyMaxThreshold_,
                                               shift_*reassemblyShiftWeight_));
 
            updateDistanceFromLastLinearization_(uLastIter, deltaU);
            partialReassembler_->computeColors(this->assembler(),
                                               distanceFromLastLinearization_,
                                               reassemblyThreshold);
 
            // set the discrepancy of the red entities to zero
            for (unsigned int i = 0; i < distanceFromLastLinearization_.size(); i++)
                if (partialReassembler_->dofColor(i) == EntityColor::red)
                    distanceFromLastLinearization_[i] = 0;
        }
 
        if (useLineSearch_)
            lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
 
        else if (useChop_)
            choppedUpdate_(uCurrentIter, uLastIter, deltaU);
 
        else
        {
            uCurrentIter = uLastIter;
            uCurrentIter -= deltaU;
 
            if (enableResidualCriterion_)
                computeResidualReduction_(uCurrentIter);
 
            else
            {
                // If we get here, the convergence criterion does not require
                // additional residual evaluations. Thus, the grid variables have
                // not yet been updated to the new uCurrentIter.
                this->assembler().updateGridVariables(uCurrentIter);
            }
        }
    }
 
    virtual void newtonEndStep(SolutionVector &uCurrentIter,
                               const SolutionVector &uLastIter)
    {
        invokePriVarSwitch_(uCurrentIter, HasPriVarsSwitch{});
 
        ++numSteps_;
 
        if (verbosity_ >= 1)
        {
            if (enableDynamicOutput_)
                std::cout << '\r'; // move cursor to beginning of line
 
            auto width = std::to_string(maxSteps_).size();
            std::cout << "Newton iteration " << std::setw(width) << numSteps_ << " done";
 
            auto formatFlags = std::cout.flags();
            auto prec = std::cout.precision();
            std::cout << std::scientific << std::setprecision(3);
 
            if (enableShiftCriterion_)
                std::cout << ", maximum relative shift = " << shift_;
            if (enableResidualCriterion_ && enableAbsoluteResidualCriterion_)
                std::cout << ", residual = " << residualNorm_;
            else if (enableResidualCriterion_)
                std::cout << ", residual reduction = " << reduction_;
 
            std::cout.flags(formatFlags);
            std::cout.precision(prec);
 
            std::cout << endIterMsgStream_.str() << "\n";
        }
        endIterMsgStream_.str("");
 
        // When the Newton iterations are done: ask the model to check whether it makes sense
        // TODO: how do we realize this? -> do this here in the Newton solver
        // model_().checkPlausibility();
    }
 
    virtual void newtonEnd()  {}
 
    virtual bool newtonConverged() const
    {
        // in case the model has a priVar switch and some some primary variables
        // actually switched their state in the last iteration, enforce another iteration
        if (priVarsSwitchedInLastIteration_)
            return false;
 
        if (enableShiftCriterion_ && !enableResidualCriterion_)
        {
            return shift_ <= shiftTolerance_;
        }
        else if (!enableShiftCriterion_ && enableResidualCriterion_)
        {
            if(enableAbsoluteResidualCriterion_)
                return residualNorm_ <= residualTolerance_;
            else
                return reduction_ <= reductionTolerance_;
        }
        else if (satisfyResidualAndShiftCriterion_)
        {
            if(enableAbsoluteResidualCriterion_)
                return shift_ <= shiftTolerance_
                        && residualNorm_ <= residualTolerance_;
            else
                return shift_ <= shiftTolerance_
                        && reduction_ <= reductionTolerance_;
        }
        else if(enableShiftCriterion_ && enableResidualCriterion_)
        {
            if(enableAbsoluteResidualCriterion_)
                return shift_ <= shiftTolerance_
                        || residualNorm_ <= residualTolerance_;
            else
                return shift_ <= shiftTolerance_
                        || reduction_ <= reductionTolerance_;
        }
        else
        {
            return shift_ <= shiftTolerance_
                    || reduction_ <= reductionTolerance_
                    || residualNorm_ <= residualTolerance_;
        }
 
        return false;
    }
 
    virtual void newtonFail(SolutionVector& u) {}
 
    virtual void newtonSucceed()  {}
 
    void report(std::ostream& sout = std::cout) const
    {
        sout << '\n'
             << "Newton statistics\n"
             << "----------------------------------------------\n"
             << "-- Total Newton iterations:           " << totalWastedIter_ + totalSucceededIter_ << '\n'
             << "-- Total wasted Newton iterations:    " << totalWastedIter_ << '\n'
             << "-- Total succeeded Newton iterations: " << totalSucceededIter_ << '\n'
             << "-- Average iterations per solve:      " << std::setprecision(3) << double(totalSucceededIter_) / double(numConverged_) << '\n'
             << std::endl;
    }
 
    void resetReport()
    {
        totalWastedIter_ = 0;
        totalSucceededIter_ = 0;
        numConverged_ = 0;
    }
 
    void reportParams(std::ostream& sout = std::cout) const
    {
        sout << "\nNewton solver configured with the following options and parameters:\n";
        // options
        if (useLineSearch_) sout << " -- Newton.UseLineSearch = true\n";
        if (useChop_) sout << " -- Newton.EnableChop = true\n";
        if (enablePartialReassembly_) sout << " -- Newton.EnablePartialReassembly = true\n";
        if (enableAbsoluteResidualCriterion_) sout << " -- Newton.EnableAbsoluteResidualCriterion = true\n";
        if (enableShiftCriterion_) sout << " -- Newton.EnableShiftCriterion = true (relative shift convergence criterion)\n";
        if (enableResidualCriterion_) sout << " -- Newton.EnableResidualCriterion = true\n";
        if (satisfyResidualAndShiftCriterion_) sout << " -- Newton.SatisfyResidualAndShiftCriterion = true\n";
        // parameters
        if (enableShiftCriterion_) sout << " -- Newton.MaxRelativeShift = " << shiftTolerance_ << '\n';
        if (enableAbsoluteResidualCriterion_) sout << " -- Newton.MaxAbsoluteResidual = " << residualTolerance_ << '\n';
        if (enableResidualCriterion_) sout << " -- Newton.ResidualReduction = " << reductionTolerance_ << '\n';
        sout << " -- Newton.MinSteps = " << minSteps_ << '\n';
        sout << " -- Newton.MaxSteps = " << maxSteps_ << '\n';
        sout << " -- Newton.TargetSteps = " << targetSteps_ << '\n';
        if (enablePartialReassembly_)
        {
            sout << " -- Newton.ReassemblyMinThreshold = " << reassemblyMinThreshold_ << '\n';
            sout << " -- Newton.ReassemblyMaxThreshold = " << reassemblyMaxThreshold_ << '\n';
            sout << " -- Newton.ReassemblyShiftWeight = " << reassemblyShiftWeight_ << '\n';
        }
        sout << " -- Newton.RetryTimeStepReductionFactor = " << retryTimeStepReductionFactor_ << '\n';
        sout << " -- Newton.MaxTimeStepDivisions = " << maxTimeStepDivisions_ << '\n';
        sout << std::endl;
    }
 
    Scalar suggestTimeStepSize(Scalar oldTimeStep) const
    {
        // be aggressive reducing the time-step size but
        // conservative when increasing it. the rationale is
        // that we want to avoid failing in the next Newton
        // iteration which would require another linearization
        // of the problem.
        if (numSteps_ > targetSteps_) {
            Scalar percent = Scalar(numSteps_ - targetSteps_)/targetSteps_;
            return oldTimeStep/(1.0 + percent);
        }
 
        Scalar percent = Scalar(targetSteps_ - numSteps_)/targetSteps_;
        return oldTimeStep*(1.0 + percent/1.2);
    }
 
    [[deprecated("Has been replaced by setVerbosity(int). Will be removed after 3.1 release!")]]
    void setVerbose(bool val)
    { verbosity_ = val; }
 
    [[deprecated("Has been replaced by int verbosity(). Will be removed after 3.1 release!")]]
    bool verbose() const
    { return verbosity_ ; }
 
    void setVerbosity(int val)
    { verbosity_ = val; }
 
    int verbosity() const
    { return verbosity_ ; }
 
    const std::string& paramGroup() const
    { return paramGroup_; }
 
    void attachConvergenceWriter(std::shared_ptr<ConvergenceWriter> convWriter)
    { convergenceWriter_ = convWriter; }
 
    void detachConvergenceWriter()
    { convergenceWriter_ = nullptr; }
 
protected:
 
    void computeResidualReduction_(const SolutionVector &uCurrentIter)
    {
        residualNorm_ = this->assembler().residualNorm(uCurrentIter);
        reduction_ = residualNorm_;
        reduction_ /= initialResidual_;
    }
 
    bool enableResidualCriterion() const
    { return enableResidualCriterion_; }
 
    void initPriVarSwitch_(SolutionVector&, std::false_type) {}
 
    void initPriVarSwitch_(SolutionVector& sol, std::true_type)
    {
        priVarSwitch_->reset(sol.size());
        priVarsSwitchedInLastIteration_ = false;
 
        const auto& problem = this->assembler().problem();
        const auto& gridGeometry = this->assembler().gridGeometry();
        auto& gridVariables = this->assembler().gridVariables();
        priVarSwitch_->updateBoundary(problem, gridGeometry, gridVariables, sol);
    }
 
    void invokePriVarSwitch_(SolutionVector&, std::false_type) {}
 
    void invokePriVarSwitch_(SolutionVector& uCurrentIter, std::true_type)
    {
        // update the variable switch (returns true if the pri vars at at least one dof were switched)
        // for disabled grid variable caching
        const auto& gridGeometry = this->assembler().gridGeometry();
        const auto& problem = this->assembler().problem();
        auto& gridVariables = this->assembler().gridVariables();
 
        // invoke the primary variable switch
        priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables,
                                                                problem, gridGeometry);
 
        if (priVarsSwitchedInLastIteration_)
        {
            for (const auto& element : elements(gridGeometry.gridView()))
            {
                // if the volume variables are cached globally, we need to update those where the primary variables have been switched
                priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, uCurrentIter);
 
                // if the flux variables are cached globally, we need to update those where the primary variables have been switched
                priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, uCurrentIter);
            }
        }
    }
 
    int targetSteps_;
    int minSteps_;
    int maxSteps_;
    int numSteps_;
 
    // residual criterion variables
    Scalar reduction_;
    Scalar residualNorm_;
    Scalar lastReduction_;
    Scalar initialResidual_;
 
    // shift criterion variables
    Scalar shift_;
    Scalar lastShift_;
 
    std::ostringstream endIterMsgStream_;
 
 
private:
 
    bool solve_(SolutionVector& uCurrentIter)
    {
        // the given solution is the initial guess
        SolutionVector uLastIter(uCurrentIter);
        SolutionVector deltaU(uCurrentIter);
 
        Dune::Timer assembleTimer(false);
        Dune::Timer solveTimer(false);
        Dune::Timer updateTimer(false);
 
        if (enablePartialReassembly_)
        {
            partialReassembler_->resetColors();
            resizeDistanceFromLastLinearization_(uCurrentIter, distanceFromLastLinearization_);
        }
 
        try
        {
            newtonBegin(uCurrentIter);
 
            // execute the method as long as the solver thinks
            // that we should do another iteration
            bool converged = false;
            while (newtonProceed(uCurrentIter, converged))
            {
                // notify the solver that we're about to start
                // a new timestep
                newtonBeginStep(uCurrentIter);
 
                // make the current solution to the old one
                if (numSteps_ > 0)
                    uLastIter = uCurrentIter;
 
                if (verbosity_ >= 1 && enableDynamicOutput_)
                    std::cout << "Assemble: r(x^k) = dS/dt + div F - q;   M = grad r"
                              << std::flush;
 
                // assemble
 
                // linearize the problem at the current solution
                assembleTimer.start();
                assembleLinearSystem(uCurrentIter);
                assembleTimer.stop();
 
                // linear solve
 
                // Clear the current line using an ansi escape
                // sequence.  for an explanation see
                // http://en.wikipedia.org/wiki/ANSI_escape_code
                const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
 
                if (verbosity_ >= 1 && enableDynamicOutput_)
                    std::cout << "\rSolve: M deltax^k = r"
                              << clearRemainingLine << std::flush;
 
                // solve the resulting linear equation system
                solveTimer.start();
 
                // set the delta vector to zero before solving the linear system!
                deltaU = 0;
 
                solveLinearSystem(deltaU);
                solveTimer.stop();
 
                // update
                if (verbosity_ >= 1 && enableDynamicOutput_)
                    std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k"
                              << clearRemainingLine << std::flush;
 
                updateTimer.start();
                // update the current solution (i.e. uOld) with the delta
                // (i.e. u). The result is stored in u
                newtonUpdate(uCurrentIter, uLastIter, deltaU);
                updateTimer.stop();
 
                // tell the solver that we're done with this iteration
                newtonEndStep(uCurrentIter, uLastIter);
 
                // if a convergence writer was specified compute residual and write output
                if (convergenceWriter_)
                {
                    this->assembler().assembleResidual(uCurrentIter);
                    convergenceWriter_->write(uCurrentIter, deltaU, this->assembler().residual());
                }
 
                // detect if the method has converged
                converged = newtonConverged();
            }
 
            // tell solver we are done
            newtonEnd();
 
            // reset state if Newton failed
            if (!newtonConverged())
            {
                totalWastedIter_ += numSteps_;
                newtonFail(uCurrentIter);
                return false;
            }
 
            totalSucceededIter_ += numSteps_;
            numConverged_++;
 
            // tell solver we converged successfully
            newtonSucceed();
 
            if (verbosity_ >= 1) {
                const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
                std::cout << "Assemble/solve/update time: "
                          <<  assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
                          <<  solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
                          <<  updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
                          << "\n";
            }
            return true;
 
        }
        catch (const NumericalProblem &e)
        {
            if (verbosity_ >= 1)
                std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
 
            totalWastedIter_ += numSteps_;
            newtonFail(uCurrentIter);
            return false;
        }
    }
 
    template<class A>
    auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter)
    -> typename std::enable_if_t<decltype(isValid(Detail::supportsPartialReassembly())(assembler))::value, void>
    {
        this->assembler().assembleJacobianAndResidual(uCurrentIter, partialReassembler_.get());
    }
 
    template<class A>
    auto assembleLinearSystem_(const A& assembler, const SolutionVector& uCurrentIter)
    -> typename std::enable_if_t<!decltype(isValid(Detail::supportsPartialReassembly())(assembler))::value, void>
    {
        this->assembler().assembleJacobianAndResidual(uCurrentIter);
    }
 
    virtual void newtonUpdateShift_(const SolutionVector &uLastIter,
                                    const SolutionVector &deltaU)
    {
        shift_ = 0;
        newtonUpdateShiftImpl_(uLastIter, deltaU);
 
        if (comm_.size() > 1)
            shift_ = comm_.max(shift_);
    }
 
    template<class SolVec>
    void newtonUpdateShiftImpl_(const SolVec &uLastIter,
                                const SolVec &deltaU)
    {
        for (int i = 0; i < int(uLastIter.size()); ++i) {
            typename SolVec::block_type uNewI = uLastIter[i];
            uNewI -= deltaU[i];
 
            Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI);
            using std::max;
            shift_ = max(shift_, shiftAtDof);
        }
    }
 
    template<class ...Args>
    void newtonUpdateShiftImpl_(const Dune::MultiTypeBlockVector<Args...> &uLastIter,
                                const Dune::MultiTypeBlockVector<Args...> &deltaU)
    {
        // There seems to be a bug in g++5 which which prevents compilation when
        // passing the call to the implementation directly to Dune::Hybrid::forEach.
        // We therefore store this call in a lambda and pass it to the for loop afterwards.
        auto doUpdate = [&](const auto subVectorIdx)
        {
            this->newtonUpdateShiftImpl_(uLastIter[subVectorIdx], deltaU[subVectorIdx]);
        };
 
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(uLastIter)), doUpdate);
    }
 
    virtual void lineSearchUpdate_(SolutionVector &uCurrentIter,
                                   const SolutionVector &uLastIter,
                                   const SolutionVector &deltaU)
    {
        Scalar lambda = 1.0;
        SolutionVector tmp(uLastIter);
 
        while (true)
        {
            uCurrentIter = deltaU;
            uCurrentIter *= -lambda;
            uCurrentIter += uLastIter;
 
            computeResidualReduction_(uCurrentIter);
 
            if (reduction_ < lastReduction_ || lambda <= 0.125) {
                endIterMsgStream_ << ", residual reduction " << lastReduction_ << "->"  << reduction_ << "@lambda=" << lambda;
                return;
            }
 
            // try with a smaller update
            lambda /= 2.0;
        }
    }
 
    virtual void choppedUpdate_(SolutionVector &uCurrentIter,
                                const SolutionVector &uLastIter,
                                const SolutionVector &deltaU)
    {
        DUNE_THROW(Dune::NotImplemented,
                   "Chopped Newton update strategy not implemented.");
    }
 
    virtual bool solveLinearSystem_(SolutionVector& deltaU)
    {
        return solveLinearSystemImpl_(this->linearSolver(),
                                      this->assembler().jacobian(),
                                      deltaU,
                                      this->assembler().residual());
    }
 
    template<class V = SolutionVector>
    typename std::enable_if_t<!isMultiTypeBlockVector<V>(), bool>
    solveLinearSystemImpl_(LinearSolver& ls,
                           JacobianMatrix& A,
                           SolutionVector& x,
                           SolutionVector& b)
    {
        constexpr auto blockSize = JacobianMatrix::block_type::rows;
        using BlockType = Dune::FieldVector<Scalar, blockSize>;
        Dune::BlockVector<BlockType> xTmp; xTmp.resize(b.size());
        Dune::BlockVector<BlockType> bTmp(xTmp);
        for (unsigned int i = 0; i < b.size(); ++i)
            for (unsigned int j = 0; j < blockSize; ++j)
                bTmp[i][j] = b[i][j];
 
        const int converged = ls.solve(A, xTmp, bTmp);
 
        for (unsigned int i = 0; i < x.size(); ++i)
            for (unsigned int j = 0; j < blockSize; ++j)
                x[i][j] = xTmp[i][j];
 
        return converged;
    }
 
 
    template<class LS = LinearSolver, class V = SolutionVector>
    typename std::enable_if_t<linearSolverAcceptsMultiTypeMatrix<LS>() &&
                              isMultiTypeBlockVector<V>(), bool>
    solveLinearSystemImpl_(LinearSolver& ls,
                           JacobianMatrix& A,
                           SolutionVector& x,
                           SolutionVector& b)
    {
        // check matrix sizes
        assert(checkMatrix_(A) && "Sub blocks of MultiType matrix have wrong sizes!");
 
        // TODO: automatically derive the precondBlockLevel
        return ls.template solve</*precondBlockLevel=*/2>(A, x, b);
    }
 
    template<class LS = LinearSolver, class V = SolutionVector>
    typename std::enable_if_t<!linearSolverAcceptsMultiTypeMatrix<LS>() &&
                              isMultiTypeBlockVector<V>(), bool>
    solveLinearSystemImpl_(LinearSolver& ls,
                           JacobianMatrix& A,
                           SolutionVector& x,
                           SolutionVector& b)
    {
        // check matrix sizes
        assert(checkMatrix_(A) && "Sub blocks of MultiType matrix have wrong sizes!");
 
        // create the bcrs matrix the IterativeSolver backend can handle
        const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
 
        // get the new matrix sizes
        const std::size_t numRows = M.N();
        assert(numRows == M.M());
 
        // create the vector the IterativeSolver backend can handle
        const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
        assert(bTmp.size() == numRows);
 
        // create a blockvector to which the linear solver writes the solution
        using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
        using BlockVector = typename Dune::BlockVector<VectorBlock>;
        BlockVector y(numRows);
 
        // solve
        const bool converged = ls.solve(M, y, bTmp);
 
        // copy back the result y into x
        if(converged)
            VectorConverter<SolutionVector>::retrieveValues(x, y);
 
        return converged;
    }
 
    template<class M = JacobianMatrix>
    typename std::enable_if_t<!isBCRSMatrix<M>(), bool>
    checkMatrix_(const JacobianMatrix& A)
    {
        bool matrixHasCorrectSize = true;
        using namespace Dune::Hybrid;
        using namespace Dune::Indices;
        forEach(A, [&matrixHasCorrectSize](const auto& rowOfMultiTypeMatrix)
        {
            const auto numRowsLeftMostBlock = rowOfMultiTypeMatrix[_0].N();
 
            forEach(rowOfMultiTypeMatrix, [&matrixHasCorrectSize, &numRowsLeftMostBlock](const auto& subBlock)
            {
                if (subBlock.N() != numRowsLeftMostBlock)
                    matrixHasCorrectSize = false;
            });
        });
        return matrixHasCorrectSize;
    }
 
    void initParams_(const std::string& group = "")
    {
        useLineSearch_ = getParamFromGroup<bool>(group, "Newton.UseLineSearch");
        useChop_ = getParamFromGroup<bool>(group, "Newton.EnableChop");
        if(useLineSearch_ && useChop_)
            DUNE_THROW(Dune::InvalidStateException, "Use either linesearch OR chop!");
 
        enableAbsoluteResidualCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableAbsoluteResidualCriterion");
        enableShiftCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableShiftCriterion");
        enableResidualCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableResidualCriterion") || enableAbsoluteResidualCriterion_;
        satisfyResidualAndShiftCriterion_ = getParamFromGroup<bool>(group, "Newton.SatisfyResidualAndShiftCriterion");
        enableDynamicOutput_ = getParamFromGroup<bool>(group, "Newton.EnableDynamicOutput", true);
 
        if (!enableShiftCriterion_ && !enableResidualCriterion_)
        {
            DUNE_THROW(Dune::NotImplemented,
                       "at least one of NewtonEnableShiftCriterion or "
                       << "NewtonEnableResidualCriterion has to be set to true");
        }
 
        setMaxRelativeShift(getParamFromGroup<Scalar>(group, "Newton.MaxRelativeShift"));
        setMaxAbsoluteResidual(getParamFromGroup<Scalar>(group, "Newton.MaxAbsoluteResidual"));
        setResidualReduction(getParamFromGroup<Scalar>(group, "Newton.ResidualReduction"));
        setTargetSteps(getParamFromGroup<int>(group, "Newton.TargetSteps"));
        setMinSteps(getParamFromGroup<int>(group, "Newton.MinSteps"));
        setMaxSteps(getParamFromGroup<int>(group, "Newton.MaxSteps"));
 
        enablePartialReassembly_ = getParamFromGroup<bool>(group, "Newton.EnablePartialReassembly");
        reassemblyMinThreshold_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyMinThreshold", 1e-1*shiftTolerance_);
        reassemblyMaxThreshold_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyMaxThreshold", 1e2*shiftTolerance_);
        reassemblyShiftWeight_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyShiftWeight", 1e-3);
 
        maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, "Newton.MaxTimeStepDivisions", 10);
        retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, "Newton.RetryTimeStepReductionFactor", 0.5);
 
        verbosity_ = comm_.rank() == 0 ? getParamFromGroup<int>(group, "Newton.Verbosity", 2) : 0;
        numSteps_ = 0;
 
        // output a parameter report
        if (verbosity_ >= 2)
            reportParams();
    }
 
    template<class Sol>
    void updateDistanceFromLastLinearization_(const Sol& u, const Sol& uDelta)
    {
        for (size_t i = 0; i < u.size(); ++i) {
            const auto& currentPriVars(u[i]);
            auto nextPriVars(currentPriVars);
            nextPriVars -= uDelta[i];
 
            // add the current relative shift for this degree of freedom
            auto shift = relativeShiftAtDof_(currentPriVars, nextPriVars);
            distanceFromLastLinearization_[i] += shift;
        }
    }
 
    template<class ...Args>
    void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<Args...>& uLastIter,
                                              const Dune::MultiTypeBlockVector<Args...>& deltaU)
    {
        DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
    }
 
    template<class Sol>
    void resizeDistanceFromLastLinearization_(const Sol& u, std::vector<Scalar>& dist)
    {
        dist.assign(u.size(), 0.0);
    }
 
    template<class ...Args>
    void resizeDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector<Args...>& u,
                                              std::vector<Scalar>& dist)
    {
        DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector");
    }
 
    template<class PrimaryVariables>
    Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
                               const PrimaryVariables &priVars2) const
    {
        Scalar result = 0.0;
        using std::abs;
        using std::max;
        // iterate over all primary variables
        for (int j = 0; j < PrimaryVariables::dimension; ++j) {
            Scalar eqErr = abs(priVars1[j] - priVars2[j]);
            eqErr /= max<Scalar>(1.0,abs(priVars1[j] + priVars2[j])/2);
 
            result = max(result, eqErr);
        }
        return result;
    }
 
    Communication comm_;
 
    int verbosity_;
 
    Scalar shiftTolerance_;
    Scalar reductionTolerance_;
    Scalar residualTolerance_;
 
    // time step control
    std::size_t maxTimeStepDivisions_;
    Scalar retryTimeStepReductionFactor_;
 
    // further parameters
    bool useLineSearch_;
    bool useChop_;
    bool enableAbsoluteResidualCriterion_;
    bool enableShiftCriterion_;
    bool enableResidualCriterion_;
    bool satisfyResidualAndShiftCriterion_;
    bool enableDynamicOutput_;
 
    std::string paramGroup_;
 
    // infrastructure for partial reassembly
    bool enablePartialReassembly_;
    std::unique_ptr<Reassembler> partialReassembler_;
    std::vector<Scalar> distanceFromLastLinearization_;
    Scalar reassemblyMinThreshold_;
    Scalar reassemblyMaxThreshold_;
    Scalar reassemblyShiftWeight_;
 
    // statistics for the optional report
    std::size_t totalWastedIter_ = 0; 
    std::size_t totalSucceededIter_ = 0; 
    std::size_t numConverged_ = 0; 
 
    std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
    bool priVarsSwitchedInLastIteration_ = false;
 
    std::shared_ptr<ConvergenceWriter> convergenceWriter_ = nullptr;
};
 
} // end namespace Dumux
 
#endif