line_search_strategy.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ \.
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
 
 
#if !defined(KRATOS_LINE_SEARCH_STRATEGY )
#define  KRATOS_LINE_SEARCH_STRATEGY
 
// System includes
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "includes/model_part.h"
#include "solving_strategies/strategies/residualbased_newton_raphson_strategy.h"
 
 
namespace Kratos
{
 
 
 
 
 
 
 
 
 
 
 
 
 
 
//URL[Example of use html]{ extended_documentation/no_ex_of_use.html}
 
//URL[Example of use pdf]{ extended_documentation/no_ex_of_use.pdf}
 
//URL[Example of use doc]{ extended_documentation/no_ex_of_use.doc}
 
//URL[Example of use ps]{ extended_documentation/no_ex_of_use.ps}
 
 
//URL[Extended documentation html]{ extended_documentation/no_ext_doc.html}
 
//URL[Extended documentation pdf]{ extended_documentation/no_ext_doc.pdf}
 
//URL[Extended documentation doc]{ extended_documentation/no_ext_doc.doc}
 
//URL[Extended documentation ps]{ extended_documentation/no_ext_doc.ps}
 
 
 
template<class TSparseSpace,
         class TDenseSpace, // = DenseSpace<double>,
         class TLinearSolver //= LinearSolver<TSparseSpace,TDenseSpace>
         >
class LineSearchStrategy
    : public ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
{
public:
    typedef ConvergenceCriteria<TSparseSpace, TDenseSpace> TConvergenceCriteriaType;
 
    // Counted pointer of ClassName
    KRATOS_CLASS_POINTER_DEFINITION(LineSearchStrategy);
 
    typedef SolvingStrategy<TSparseSpace, TDenseSpace> SolvingStrategyType;
    typedef ResidualBasedNewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver> BaseType;
    typedef LineSearchStrategy<TSparseSpace,TDenseSpace,TLinearSolver> ClassType;
 
    typedef typename BaseType::TBuilderAndSolverType TBuilderAndSolverType;
 
    typedef typename BaseType::TDataType TDataType;
 
    typedef TSparseSpace SparseSpaceType;
 
    typedef typename BaseType::TSchemeType TSchemeType;
 
    //typedef typename BaseType::DofSetType DofSetType;
 
    typedef typename BaseType::DofsArrayType DofsArrayType;
 
    typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
 
    typedef typename BaseType::TSystemVectorType TSystemVectorType;
 
    typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
 
    typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
 
    typedef typename BaseType::TSystemMatrixPointerType TSystemMatrixPointerType;
    typedef typename BaseType::TSystemVectorPointerType TSystemVectorPointerType;
 
 
 
    explicit LineSearchStrategy() : BaseType()
    {
    }
 
    explicit LineSearchStrategy(ModelPart& rModelPart, Parameters ThisParameters)
        : BaseType(rModelPart)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    explicit LineSearchStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        int MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false
        ) : BaseType(rModelPart, pScheme, pNewLinearSolver,pNewConvergenceCriteria,MaxIterations,CalculateReactions,ReformDofSetAtEachStep, MoveMeshFlag)
    {
        Parameters default_settings = this->GetDefaultParameters();
        OverrideDefaultSettingsWithParameters(default_settings, MaxIterations, ReformDofSetAtEachStep, CalculateReactions);
        this->AssignSettings(default_settings);
    }
 
    explicit LineSearchStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        int MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false
        ) : BaseType(rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag)
    {
        Parameters default_settings = this->GetDefaultParameters();
        OverrideDefaultSettingsWithParameters(default_settings, MaxIterations, ReformDofSetAtEachStep, CalculateReactions);
        this->AssignSettings(default_settings);
    }
 
    KRATOS_DEPRECATED_MESSAGE("Constructor deprecated, please use the constructor without linear solver")
    explicit LineSearchStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        int MaxIterations = 30,
        bool CalculateReactions = false,
        bool ReformDofSetAtEachStep = false,
        bool MoveMeshFlag = false
        ) : BaseType(rModelPart, pScheme, pNewLinearSolver,pNewConvergenceCriteria,pNewBuilderAndSolver,MaxIterations,CalculateReactions,ReformDofSetAtEachStep, MoveMeshFlag)
    {
        Parameters default_settings = this->GetDefaultParameters();
        OverrideDefaultSettingsWithParameters(default_settings, MaxIterations, ReformDofSetAtEachStep, CalculateReactions);
        this->AssignSettings(default_settings);
    }
 
    LineSearchStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        Parameters ThisParameters
        ): BaseType(rModelPart, pScheme, pNewLinearSolver,pNewConvergenceCriteria, BaseType::GetDefaultParameters())
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    LineSearchStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        Parameters ThisParameters
        ): BaseType(rModelPart, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, BaseType::GetDefaultParameters())
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    ~LineSearchStrategy() override
    {
    }
 
 
 
 
 
    typename SolvingStrategyType::Pointer Create(
        ModelPart& rModelPart,
        Parameters ThisParameters
        ) const override
    {
        return Kratos::make_shared<ClassType>(rModelPart, ThisParameters);
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"                       : "line_search_strategy",
            "max_line_search_iterations" : 5,
            "first_alpha_value"          : 0.5,
            "second_alpha_value"         : 1.0,
            "min_alpha"                  : 0.1,
            "max_alpha"                  : 2.0,
            "line_search_tolerance"      : 0.5
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "line_search_strategy";
    }
 
 
 
 
 
    std::string Info() const override
    {
        return "LineSearchStrategy";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
 
 
private:
 
 
 
    int mMaxLineSearchIterations;   //Maximum number of iterations line search do
    double mFirstAlphaValue;        //First alpha guess value used for the first iteration
    double mSecondAlphaValue;       //Second alpha guess value used for the first iteration
    double mMinAlpha;               //Minimum possible alpha value at the end of the algorithm
    double mMaxAlpha;               //Maximum possible alpha value at the end of the algorithm
    double mLineSearchTolerance;    //Tolerance of the line search algorithm, defined as the ratio between
                                    //maximum residual*alpha*dx and current iteration residual*alpha*dx
 
 
 
 
 
 
 
 
 
 
 
 
protected:
 
 
 
 
 
    void UpdateDatabase(
        TSystemMatrixType& A,
        TSystemVectorType& Dx,
        TSystemVectorType& b,
        const bool MoveMesh
    ) override
    {
        typename TSchemeType::Pointer pScheme = this->GetScheme();
        typename TBuilderAndSolverType::Pointer pBuilderAndSolver = this->GetBuilderAndSolver();
 
        TSystemVectorType aux(Dx);
 
        double x1 = mFirstAlphaValue;
        double x2 = mSecondAlphaValue;
 
        bool converged = false;
        int it = 0;
        double xprevious = 0.0;
 
        //Compute residual with 1 coefficient update (x1)
        //since no update was performed yet, this includes an increment wrt the previous
        //solution of x1*Dx
        TSparseSpace::Assign(aux,x1-xprevious, Dx);
        xprevious = x1;
        BaseType::UpdateDatabase(A,aux,b,MoveMesh);
        TSparseSpace::SetToZero(b);
        pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
        double r1 = TSparseSpace::Dot(aux,b);
 
        double rmax = std::abs(r1);
        while(!converged && it < mMaxLineSearchIterations) {
 
            //Compute residual with 2 coefficient update (x2)
            //since the database was initialized with x1*Dx
            //we need to apply ONLY THE INCREMENT, that is (x2-xprevious)*Dx
            TSparseSpace::Assign(aux,x2-xprevious, Dx);
            xprevious = x2;
            BaseType::UpdateDatabase(A,aux,b,MoveMesh);
            TSparseSpace::SetToZero(b);
            pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
            double r2 = TSparseSpace::Dot(aux,b);
 
            if(it == 0) {
                rmax = std::max(rmax,std::abs(r2));
            }
            double rmin = std::min(std::abs(r1),std::abs(r2));
 
            //Find optimum
            double x = 1.0;
            if(std::abs(r1 - r2) > 1e-10)
                x =  (r1*x2 - r2*x1)/(r1 - r2);
 
            if(x < mMinAlpha) {
                x = mMinAlpha;
            } else if(x > mMaxAlpha) {
                x = mMaxAlpha;
            }
 
            //Perform final update
            TSparseSpace::Assign(aux,x-xprevious, Dx);
            xprevious = x;
            BaseType::UpdateDatabase(A,aux,b,MoveMesh);
            if(rmin < mLineSearchTolerance*rmax) {
                KRATOS_INFO("LineSearchStrategy") << "LINE SEARCH it " << it << " coeff = " << x <<  " r1 = " << r1 << " r2 = " << r2 << std::endl;
                converged = true;
                TSparseSpace::Assign(aux,x, Dx);
                break;
            }
 
            //note that we compute the next residual only if it is strictly needed (we break on the line before if it is not needed)
            TSparseSpace::SetToZero(b);
            pBuilderAndSolver->BuildRHS(pScheme, BaseType::GetModelPart(), b );
            double rf = TSparseSpace::Dot(aux,b);
 
            KRATOS_INFO("LineSearchStrategy") << "LINE SEARCH it " << it << " coeff = " << x << " rf = " << rf << " r1 = " << r1 << " r2 = " << r2 << std::endl;
 
 
            if(std::abs(rf) < rmax*mLineSearchTolerance) {
                converged = true;
                TSparseSpace::Assign(aux,x, Dx);
            } else {
                if(std::abs(r1)>std::abs(r2)) {
                    r1 = rf;
                    x1 = x;
                } else {
                    r2 = r1;
                    x2 = x1;
                    r1 = rf;
                    x1 = x;
                }
                converged = false;
            }
 
 
            it++;
        }
        TSparseSpace::SetToZero(b);
    }
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
        mMaxLineSearchIterations = ThisParameters["max_line_search_iterations"].GetInt();
        mFirstAlphaValue = ThisParameters["first_alpha_value"].GetDouble();
        mSecondAlphaValue = ThisParameters["second_alpha_value"].GetDouble();
        mMinAlpha = ThisParameters["min_alpha"].GetDouble();
        mMaxAlpha = ThisParameters["max_alpha"].GetDouble();
        mLineSearchTolerance = ThisParameters["line_search_tolerance"].GetDouble();
    }
 
    void OverrideDefaultSettingsWithParameters(
        Parameters& rDefaultSettings,
        const double MaxIterations,
        const bool ReformDofSetAtEachStep,
        const bool CalculateReactions
        )
    {
        rDefaultSettings["max_iteration"].SetInt(MaxIterations);
        rDefaultSettings["reform_dofs_at_each_step"].SetBool(ReformDofSetAtEachStep);
        rDefaultSettings["compute_reactions"].SetBool(CalculateReactions);
    }
 
 
 
 
 
 
 
 
    LineSearchStrategy(const LineSearchStrategy& Other)
    {
    };
 
 
 
}; /* Class LineSearchStrategy */
 
 
 
 
 
} /* namespace Kratos. */
 
#endif /* KRATOS_LINE_SEARCH_STRATEGY  defined */