aitken_convergence_accelerator.hpp

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ \.
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:          BSD License
//  Original author:  Ruben Zorrilla
//
 
#if !defined(KRATOS_AITKEN_CONVERGENCE_ACCELERATOR)
#define  KRATOS_AITKEN_CONVERGENCE_ACCELERATOR
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "solving_strategies/convergence_accelerators/convergence_accelerator.h"
#include "utilities/math_utils.h"
 
namespace Kratos
{
 
 
 
 
 
template<class TSparseSpace, class TDenseSpace>
class AitkenConvergenceAccelerator: public ConvergenceAccelerator<TSparseSpace, TDenseSpace>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( AitkenConvergenceAccelerator );
 
    typedef ConvergenceAccelerator<TSparseSpace, TDenseSpace>              BaseType;
 
    typedef typename BaseType::Pointer                              BaseTypePointer;
 
    typedef typename BaseType::VectorType                                VectorType;
    typedef typename BaseType::VectorPointerType                  VectorPointerType;
 
 
    explicit AitkenConvergenceAccelerator(Parameters rConvAcceleratorParameters)
    {
        Parameters aitken_default_parameters(R"(
        {
            "solver_type"            : "Relaxation",
            "acceleration_type"      : "Aitken",
            "w_0"                    : 0.825,
            "interface_block_newton" : false
        }
        )");
 
        rConvAcceleratorParameters.ValidateAndAssignDefaults(aitken_default_parameters);
 
        mOmega_0 = rConvAcceleratorParameters["w_0"].GetDouble();
    }
 
 
    AitkenConvergenceAccelerator(double rOmegaInitial = 0.825)
    {
        mOmega_0 = rOmegaInitial;
    }
 
    AitkenConvergenceAccelerator(const AitkenConvergenceAccelerator& rOther)
    {
        mOmega_0 = rOther.mOmega_0;
    }
 
    virtual BaseTypePointer Clone()
    {
        return BaseTypePointer( new AitkenConvergenceAccelerator(*this) );
    }
 
    virtual ~AitkenConvergenceAccelerator
    () {}
 
 
 
 
    void InitializeSolutionStep() override
    {
        KRATOS_TRY;
 
        mConvergenceAcceleratorIteration = 1;
 
        KRATOS_CATCH( "" );
    }
 
    void UpdateSolution(const VectorType& rResidualVector,
                        VectorType& rIterationGuess) override
    {
        KRATOS_TRY;
 
        VectorPointerType pAux(new VectorType(rResidualVector));
        std::swap(mpResidualVector_1, pAux);
 
        if (mConvergenceAcceleratorIteration == 1)
        {
            TSparseSpace::UnaliasedAdd(rIterationGuess, mOmega_0, *mpResidualVector_1);
        }
        else
        {
            VectorType Aux1minus0(*mpResidualVector_1);                  // Auxiliar copy of mResidualVector_1
            TSparseSpace::UnaliasedAdd(Aux1minus0, -1.0, *mpResidualVector_0); // mResidualVector_1 - mResidualVector_0
 
            double den = TSparseSpace::Dot(Aux1minus0, Aux1minus0);
            double num = TSparseSpace::Dot(*mpResidualVector_0, Aux1minus0);
 
            mOmega_1 = -mOmega_0*(num/den);
 
            TSparseSpace::UnaliasedAdd(rIterationGuess, mOmega_1, *mpResidualVector_1);
            mOmega_0 = mOmega_1;
        }
 
        KRATOS_CATCH( "" );
    }
 
    void FinalizeNonLinearIteration() override
    {
        KRATOS_TRY;
 
        // mpResidualVector_0 = mpResidualVector_1;
        std::swap(mpResidualVector_0, mpResidualVector_1);
        mConvergenceAcceleratorIteration += 1;
 
        KRATOS_CATCH( "" );
    }
 
 
 
 
 
 
protected:
 
 
 
    unsigned int mConvergenceAcceleratorIteration = 0;
 
    double mOmega_0;
    double mOmega_1 = 0.0;
 
    VectorPointerType mpResidualVector_0 = nullptr;
    VectorPointerType mpResidualVector_1 = nullptr;
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
 
}; /* Class AitkenConvergenceAccelerator */
 
 
 
 
}  /* namespace Kratos.*/
 
#endif /* KRATOS_AITKEN_CONVERGENCE_ACCELERATOR defined */