iterative_solver.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Pooyan Dadvand
//                   Riccardo Rossi
//
 
#if !defined(KRATOS_ITERATIVE_SOLVER_H_INCLUDED )
#define  KRATOS_ITERATIVE_SOLVER_H_INCLUDED
 
// System includes
#include <string>
#include <iostream>
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/kratos_parameters.h"
#include "linear_solvers/linear_solver.h"
#include "linear_solvers/preconditioner.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType,
         class TPreconditionerType = Preconditioner<TSparseSpaceType, TDenseSpaceType>,
         class TReordererType = Reorderer<TSparseSpaceType, TDenseSpaceType> >
class IterativeSolver : public LinearSolver<TSparseSpaceType, TDenseSpaceType, TReordererType>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(IterativeSolver);
 
    typedef LinearSolver<TSparseSpaceType, TDenseSpaceType, TReordererType> BaseType;
 
    typedef typename TSparseSpaceType::MatrixType SparseMatrixType;
 
    typedef typename TSparseSpaceType::VectorType VectorType;
 
    typedef typename TDenseSpaceType::MatrixType DenseMatrixType;
 
    typedef  TPreconditionerType PreconditionerType;
 
    typedef typename TSparseSpaceType::IndexType IndexType;
 
 
    IterativeSolver()
        : mResidualNorm(0)
        , mIterationsNumber(0)
        , mBNorm(0)
        , mpPreconditioner(new TPreconditionerType())
        , mTolerance(0)
        , mMaxIterationsNumber(0)
    {
    }
 
    IterativeSolver(double NewTolerance)
        : mResidualNorm(0)
        , mIterationsNumber(0)
        , mBNorm(0)
        , mpPreconditioner(new TPreconditionerType())
        ,   mTolerance(NewTolerance)
        , mMaxIterationsNumber(0)
    {
    }
 
    IterativeSolver(double NewTolerance, unsigned int NewMaxIterationsNumber)
        : mResidualNorm(0)
        , mIterationsNumber(0)
        , mBNorm(0)
        , mpPreconditioner(new TPreconditionerType())
        , mTolerance(NewTolerance)
        , mMaxIterationsNumber(NewMaxIterationsNumber) {}
 
    IterativeSolver(double NewTolerance, unsigned int NewMaxIterationsNumber, typename TPreconditionerType::Pointer pNewPreconditioner) :
        mResidualNorm(0), mIterationsNumber(0), mBNorm(0),
        mpPreconditioner(pNewPreconditioner),
        mTolerance(NewTolerance),
        mMaxIterationsNumber(NewMaxIterationsNumber) {}
 
    IterativeSolver(Parameters settings,
                    typename TPreconditionerType::Pointer pNewPreconditioner = Kratos::make_shared<TPreconditionerType>()
                   ):
        mResidualNorm(0),
        mIterationsNumber(0),
        mBNorm(0),
        mpPreconditioner(pNewPreconditioner)
    {
        KRATOS_TRY
 
        Parameters default_parameters( R"(
        {
        "solver_type": "IterativeSolver",
        "tolerance" : 1.0e-6,
        "max_iteration" : 200,
        "preconditioner_type": "none",
        "scaling":false
        }  )" );
 
        //now validate agains defaults -- this also ensures no type mismatch
        settings.ValidateAndAssignDefaults(default_parameters);
 
        this->SetTolerance( settings["tolerance"].GetDouble() );
        this->SetMaxIterationsNumber( settings["max_iteration"].GetInt() );
        
        
        KRATOS_CATCH("")
    }
 
    IterativeSolver(const IterativeSolver& Other) : BaseType(Other),
        mResidualNorm(Other.mResidualNorm), mIterationsNumber(Other.mIterationsNumber), mBNorm(Other.mBNorm),
        mpPreconditioner(Other.mpPreconditioner),
        mTolerance(Other.mTolerance),
        mMaxIterationsNumber(Other.mMaxIterationsNumber)
    {
 
    }
 
    ~IterativeSolver() override {}
 
 
 
    IterativeSolver& operator=(const IterativeSolver& Other)
    {
        BaseType::operator=(Other);
        mResidualNorm = Other.mResidualNorm;
        mFirstResidualNorm = Other.mFirstResidualNorm;
        mIterationsNumber = Other.mIterationsNumber;
        mBNorm = Other.mBNorm;
        return *this;
    }
 
    void InitializeSolutionStep (SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        GetPreconditioner()->InitializeSolutionStep(rA,rX,rB);
    }
 
    void FinalizeSolutionStep (SparseMatrixType& rA, VectorType& rX, VectorType& rB) override
    {
        GetPreconditioner()->FinalizeSolutionStep(rA,rX,rB);
    }
 
    void Clear() override
    {
        GetPreconditioner()->Clear();
    }
 
    bool AdditionalPhysicalDataIsNeeded() override
    {
        if (GetPreconditioner()->AdditionalPhysicalDataIsNeeded())
            return true;
        else
            return false;
    }
 
    void ProvideAdditionalData (
        SparseMatrixType& rA,
        VectorType& rX,
        VectorType& rB,
        typename ModelPart::DofsArrayType& rdof_set,
        ModelPart& r_model_part
    ) override
    {
        if (GetPreconditioner()->AdditionalPhysicalDataIsNeeded())
            GetPreconditioner()->ProvideAdditionalData(rA,rX,rB,rdof_set,r_model_part);
    }
 
 
 
 
    virtual typename TPreconditionerType::Pointer GetPreconditioner(void)
    {
        return mpPreconditioner;
    }
 
    virtual const typename TPreconditionerType::Pointer GetPreconditioner(void) const
    {
        return mpPreconditioner;
    }
 
    virtual void SetPreconditioner(typename TPreconditionerType::Pointer pNewPreconditioner)
    {
        mpPreconditioner = pNewPreconditioner;
    }
 
//       virtual typename TStopCriteriaType::Pointer GetStopCriteria(void)
//  {
//    return mpStopCriteria;
//  }
 
//       virtual void SetStopCriteria(typename TStopCriteriaType::Pointer pNewStopCriteria)
//  {
//    mpStopCriteria = pNewStopCriteria;
//  }
 
    virtual void SetMaxIterationsNumber(unsigned int NewMaxIterationsNumber)
    {
        mMaxIterationsNumber = NewMaxIterationsNumber;
    }
 
    virtual IndexType GetMaxIterationsNumber()
    {
        return mMaxIterationsNumber;
    }
 
    virtual void SetIterationsNumber(unsigned int NewIterationNumber)
    {
        mIterationsNumber = NewIterationNumber;
    }
 
    IndexType GetIterationsNumber() override
    {
        return mIterationsNumber;
    }
 
    void SetTolerance(double NewTolerance) override
    {
        mTolerance = NewTolerance;
    }
 
    double GetTolerance() override
    {
        return mTolerance;
    }
 
    virtual void SetResidualNorm(double NewResidualNorm)
    {
        if (mIterationsNumber == 1)
            mFirstResidualNorm = NewResidualNorm;
        mResidualNorm = NewResidualNorm;
    }
 
    virtual double GetResidualNorm()
    {
        return mResidualNorm;
    }
 
 
    virtual bool IterationNeeded()
    {
        return (mIterationsNumber < mMaxIterationsNumber) && (mResidualNorm > mTolerance * mBNorm);
    }
 
    virtual bool IsConverged()
    {
        return (mResidualNorm <= mTolerance * mBNorm);
    }
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "Iterative solver with " << GetPreconditioner()->Info();
        return  buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        if (mBNorm == 0.00)
            if (mResidualNorm != 0.00)
                rOStream << "    Residual ratio : infinite" << std::endl;
            else
                rOStream << "    Residual ratio : 0" << std::endl;
        else
        {
            rOStream << "    Initial Residual ratio : " << mBNorm << std::endl;
            rOStream << "    Final Residual ratio : " << mResidualNorm << std::endl;
            rOStream << "    Residual ratio : " << mResidualNorm / mBNorm << std::endl;
            rOStream << "    Slope : " << (mResidualNorm - mBNorm) / mIterationsNumber << std::endl;
        }
 
        rOStream << "    Tolerance : " << mTolerance << std::endl;
        rOStream << "    Number of iterations : " << mIterationsNumber << std::endl;
        rOStream << "    Maximum number of iterations : " << mMaxIterationsNumber;
        if (mMaxIterationsNumber == mIterationsNumber)
            rOStream << std::endl << "!!!!!!!!!!!! ITERATIVE SOLVER NON CONVERGED !!!!!!!!!!!!" << mMaxIterationsNumber;
    }
 
 
 
 
 
protected:
 
 
 
    double mResidualNorm;
 
    double mFirstResidualNorm;
 
    IndexType mIterationsNumber;
 
    double mBNorm;
 
 
 
 
    void PreconditionedMult(SparseMatrixType& rA, VectorType& rX, VectorType& rY)
    {
        GetPreconditioner()->Mult(rA, rX, rY);
    }
 
    void PreconditionedTransposeMult(SparseMatrixType& rA, VectorType& rX, VectorType& rY)
    {
        GetPreconditioner()->TransposeMult(rA, rX, rY);
    }
 
 
 
 
 
 
 
 
private:
 
 
 
    typename TPreconditionerType::Pointer mpPreconditioner;
 
    //      typename TStopCriteriaType::Pointer mpStopCriteria;
 
    double mTolerance;
 
    IndexType mMaxIterationsNumber;
 
 
 
 
 
 
 
 
 
 
 
 
}; // Class IterativeSolver
 
 
 
 
 
 
template<class TSparseSpaceType, class TDenseSpaceType, class TPreconditionerType,
         class TReordererType>
inline std::istream& operator >> (std::istream& IStream,
                                  IterativeSolver<TSparseSpaceType, TDenseSpaceType,
                                  TPreconditionerType, TReordererType>& rThis)
{
    return IStream;
}
 
template<class TSparseSpaceType, class TDenseSpaceType, class TPreconditionerType,
         class TReordererType>
inline std::ostream& operator << (std::ostream& OStream,
                                  const IterativeSolver<TSparseSpaceType, TDenseSpaceType,
                                  TPreconditionerType, TReordererType>& rThis)
{
    rThis.PrintInfo(OStream);
    OStream << std::endl;
    rThis.PrintData(OStream);
 
    return OStream;
}
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_ITERATIVE_SOLVER_H_INCLUDED  defined