ibqn_mvqn_convergence_accelerator.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main author:     Ruben Zorrilla
//
 
#if !defined(KRATOS_IBQN_MVQN_CONVERGENCE_ACCELERATOR)
#define  KRATOS_IBQN_MVQN_CONVERGENCE_ACCELERATOR
 
// System includes
 
// External includes
 
// Project includes
#include "includes/ublas_interface.h"
#include "utilities/math_utils.h"
#include "utilities/dense_householder_qr_decomposition.h"
#include "utilities/parallel_utilities.h"
 
// Application includes
 
namespace Kratos
{
 
 
 
 
 
template<class TSparseSpace, class TDenseSpace>
class MVQNFullJacobianConvergenceAccelerator;
 
template<class TSparseSpace, class TDenseSpace>
class IBQNMVQNConvergenceAccelerator: public ConvergenceAccelerator<TSparseSpace, TDenseSpace>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( IBQNMVQNConvergenceAccelerator );
 
    typedef ConvergenceAccelerator<TSparseSpace, TDenseSpace>              BaseType;
    typedef typename BaseType::Pointer                              BaseTypePointer;
 
    typedef typename BaseType::DenseVectorType                           VectorType;
    typedef typename BaseType::DenseVectorPointerType             VectorPointerType;
 
    typedef typename BaseType::DenseMatrixType                           MatrixType;
    typedef typename BaseType::DenseMatrixPointerType             MatrixPointerType;
 
    typedef MVQNFullJacobianConvergenceAccelerator<TSparseSpace, TDenseSpace>          MVQNType;
    typedef typename MVQNType::Pointer                                          MVQNPointerType;
 
 
    explicit IBQNMVQNConvergenceAccelerator() = default;
 
    explicit IBQNMVQNConvergenceAccelerator(Parameters rParameters)
    {
        rParameters.ValidateAndAssignDefaults(GetDefaultParameters());
        mInitialOmega = rParameters["w_0"].GetDouble();
 
        // Set the subdomains MVQN convergence accelerator pointers
        // Note that we call the simplified constructor with default zero relaxation omega and IBQN switch
        const double abs_cut_off_tol = rParameters["abs_cut_off_tol"].GetDouble();
        mpConvergenceAcceleratorLeft = Kratos::make_unique<MVQNType>(abs_cut_off_tol);
        mpConvergenceAcceleratorRight = Kratos::make_unique<MVQNType>(abs_cut_off_tol);
    }
 
    IBQNMVQNConvergenceAccelerator(const IBQNMVQNConvergenceAccelerator& rOther) = delete;
 
    virtual ~IBQNMVQNConvergenceAccelerator(){}
 
 
 
 
    void Initialize() override
    {
        KRATOS_TRY;
 
        mpConvergenceAcceleratorLeft->Initialize();
        mpConvergenceAcceleratorRight->Initialize();
 
        KRATOS_CATCH( "" );
    }
 
    void InitializeSolutionStep() override
    {
        KRATOS_TRY;
 
        mConvergenceAcceleratorIteration = 0;
 
        mpConvergenceAcceleratorLeft->InitializeSolutionStep();
        mpConvergenceAcceleratorRight->InitializeSolutionStep();
 
        KRATOS_CATCH( "" );
    }
 
    void InitializeNonLinearIteration() override
    {
        KRATOS_TRY;
 
        mpConvergenceAcceleratorLeft->InitializeNonLinearIteration();
        mpConvergenceAcceleratorRight->InitializeNonLinearIteration();
 
        KRATOS_CATCH( "" );
    }
 
    void UpdateSolution(
        const VectorType& rResidualVector,
        VectorType& rIterationGuess) override
    {
        KRATOS_TRY;
 
        KRATOS_ERROR << "The UpdateSolution() method cannot be called in the interface block Newton case. Use either \'UpdateSolutionLeft\' or \'UpdateSolutionRight\'" << std::endl;
 
        KRATOS_CATCH("");
    }
 
    //TODO: CREATE AN INTERMEDIATE IBQN BASE CLASS
    virtual void UpdateSolutionRight(
        const VectorType& rForceInputVector,
        const VectorType& rDisplacementOutputVector,
        VectorType& rIterationGuess)
    {
        KRATOS_TRY;
 
        // Save the uncorrected force vector for the displacement update
        auto p_uncor_disp_aux = Kratos::make_shared<VectorType>(rDisplacementOutputVector);
        std::swap(mpUncorrectedDisplacementVector, p_uncor_disp_aux);
 
        // Update right inverse Jacobian approximation with the current information
        mpConvergenceAcceleratorRight->UpdateInverseJacobianApproximation(rForceInputVector, rDisplacementOutputVector);
 
        VectorType right_correction(TSparseSpace::Size(rIterationGuess));
        if (!mFirstRightCorrectionPerformed) {
            // Calculate the first correction as a fixed point relaxation
            TSparseSpace::SetToZero(right_correction);
            TSparseSpace::UnaliasedAdd(right_correction, mInitialOmega, rDisplacementOutputVector - rIterationGuess);
 
            // Update the first iteration correction flag
            mFirstRightCorrectionPerformed = true;
 
        } else {
            // Calculate the right (displacement) correction
            SolveInterfaceBlockQuasiNewtonRightUpdate(rForceInputVector, rDisplacementOutputVector, rIterationGuess, right_correction);
 
        }
 
        // Update the iteration guess
        TSparseSpace::UnaliasedAdd(rIterationGuess, 1.0, right_correction);
 
        KRATOS_CATCH("");
    }
 
    //TODO: CREATE AN INTERMEDIATE IBQN BASE CLASS
    virtual void UpdateSolutionLeft(
        const VectorType& rDisplacementInputVector,
        const VectorType& rForceOutputVector,
        VectorType& rIterationGuess)
    {
        KRATOS_TRY;
 
        // Save the uncorrected force vector for the displacement update
        auto p_uncor_force_aux = Kratos::make_shared<VectorType>(rForceOutputVector);
        std::swap(mpUncorrectedForceVector, p_uncor_force_aux);
 
        // Update left inverse Jacobian approximation with the current information
        mpConvergenceAcceleratorLeft->UpdateInverseJacobianApproximation(rDisplacementInputVector, rForceOutputVector);
 
        VectorType left_correction(TSparseSpace::Size(rIterationGuess));
        if (!mFirstLeftCorrectionPerformed) {
            // Do nothing in the first traction correction
            // This means to set as correction the difference between the guess and the previous iteration
            left_correction = rForceOutputVector - rIterationGuess;
 
            // Update the first iteration correction flag
            mFirstLeftCorrectionPerformed = true;
 
        } else {
            // Calculate the left (traction) correction
            SolveInterfaceBlockQuasiNewtonLeftUpdate(rDisplacementInputVector, rForceOutputVector, rIterationGuess, left_correction);
 
        }
 
        // Update the iteration guess
        TSparseSpace::UnaliasedAdd(rIterationGuess, 1.0, left_correction);
 
        KRATOS_CATCH("");
    }
 
    void FinalizeNonLinearIteration() override
    {
        KRATOS_TRY;
 
        mConvergenceAcceleratorIteration += 1;
 
        mpConvergenceAcceleratorLeft->FinalizeNonLinearIteration();
        mpConvergenceAcceleratorRight->FinalizeNonLinearIteration();
 
        KRATOS_CATCH( "" );
    }
 
    void FinalizeSolutionStep() override
    {
        KRATOS_TRY;
 
        mpConvergenceAcceleratorLeft->FinalizeSolutionStep();
        mpConvergenceAcceleratorRight->FinalizeSolutionStep();
 
        KRATOS_CATCH( "" );
    }
 
    void Finalize() override
    {
        KRATOS_TRY;
 
        mpConvergenceAcceleratorLeft->Finalize();
        mpConvergenceAcceleratorRight->Finalize();
 
        KRATOS_CATCH( "" );
    }
 
    void Clear() override
    {
        KRATOS_TRY;
 
        mpConvergenceAcceleratorLeft->Clear();
        mpConvergenceAcceleratorRight->Clear();
 
        KRATOS_CATCH( "" );
    }
 
    bool IsBlockNewton() const override
    {
        return true;
    }
 
    virtual Parameters GetDefaultParameters() const
    {
        Parameters ibqn_mvqn_default_parameters(R"({
            "solver_type"            : "IBQN_MVQN",
            "w_0"                    : 0.825,
            "abs_cut_off_tol"        : 1e-8
        })");
 
        return ibqn_mvqn_default_parameters;
    }
 
 
 
 
 
 
 
 
 
protected:
 
 
 
    virtual void SolveInterfaceBlockQuasiNewtonLeftUpdate(
        const VectorType& rDisplacementInputVector,
        const VectorType& rForceOutputVector,
        const VectorType& rIterationGuess,
        VectorType& rLeftCorrection)
    {
        // Get the interface problem size
        std::size_t problem_size = mpConvergenceAcceleratorLeft->GetProblemSize();
 
        // Get both left and inverse Jacobian
        auto p_inv_jac_left = mpConvergenceAcceleratorLeft->pGetInverseJacobianApproximation();
        auto p_inv_jac_right = mpConvergenceAcceleratorRight->pGetInverseJacobianApproximation();
 
        // Set the residual of the update problem to be solved
        VectorType aux_RHS = rForceOutputVector - rIterationGuess;
        VectorType aux_left_onto_right(problem_size);
        VectorType displacement_iteration_update = *mpUncorrectedDisplacementVector - rDisplacementInputVector;
        TSparseSpace::Mult(*p_inv_jac_left, displacement_iteration_update, aux_left_onto_right);
        TSparseSpace::UnaliasedAdd(aux_RHS, 1.0, aux_left_onto_right);
 
        // Set the LHS of the update problem to be solved
        MatrixType aux_LHS = IdentityMatrix(problem_size, problem_size);
        IndexPartition<std::size_t>(problem_size).for_each(
            VectorType(problem_size),
            [&aux_LHS, &p_inv_jac_left, &p_inv_jac_right, &problem_size](std::size_t Col, VectorType& rAuxColumnTLS)
        {
            TDenseSpace::GetColumn(Col, *p_inv_jac_right, rAuxColumnTLS);
            for (std::size_t row = 0; row < problem_size; ++row) {
                aux_LHS(row,Col) -= TDenseSpace::RowDot(row, *p_inv_jac_left, rAuxColumnTLS);
            }
        });
 
        // Calculate the correction
        // Do the QR decomposition of (I - J_{F}J_{S}) and solve for the force update
        DenseHouseholderQRDecomposition<TDenseSpace> qr_decomposition;
        qr_decomposition.Compute(aux_LHS);
        qr_decomposition.Solve(aux_RHS, rLeftCorrection);
    }
 
    virtual void SolveInterfaceBlockQuasiNewtonRightUpdate(
        const VectorType& rForceInputVector,
        const VectorType& rDisplacementOutputVector,
        const VectorType& rIterationGuess,
        VectorType& rRightCorrection)
    {
        // Get the interface problem size
        std::size_t problem_size = mpConvergenceAcceleratorRight->GetProblemSize();
 
        // Get both left and inverse Jacobian
        auto p_inv_jac_right = mpConvergenceAcceleratorRight->pGetInverseJacobianApproximation();
        auto p_inv_jac_left = mpConvergenceAcceleratorLeft->pGetInverseJacobianApproximation();
 
        // Set the residual of the update problem to be solved
        VectorType aux_RHS = rDisplacementOutputVector - rIterationGuess;
        VectorType aux_right_onto_left(problem_size);
        VectorType force_iteration_update = *mpUncorrectedForceVector - rForceInputVector;
        TSparseSpace::Mult(*p_inv_jac_right, force_iteration_update, aux_right_onto_left);
        TSparseSpace::UnaliasedAdd(aux_RHS, 1.0, aux_right_onto_left);
 
        // Set the LHS of the update problem to be solved
        MatrixType aux_LHS = IdentityMatrix(problem_size, problem_size);
        IndexPartition<std::size_t>(problem_size).for_each(
            VectorType(problem_size),
            [&aux_LHS, &p_inv_jac_right, &p_inv_jac_left, &problem_size](std::size_t Col, VectorType& rAuxColumnTLS)
        {
            TDenseSpace::GetColumn(Col, *p_inv_jac_left, rAuxColumnTLS);
            for (std::size_t row = 0; row < problem_size; ++row) {
                aux_LHS(row,Col) -= TDenseSpace::RowDot(row, *p_inv_jac_right, rAuxColumnTLS);
            }
        });
 
        // Calculate the correction
        // Do the QR decomposition of (I - J_{S}J_{F}) and solve for the force update
        DenseHouseholderQRDecomposition<TDenseSpace> qr_decomposition;
        qr_decomposition.Compute(aux_LHS);
        qr_decomposition.Solve(aux_RHS, rRightCorrection);
    }
 
 
    void SetInitialRelaxationOmega(const double Omega)
    {
        mInitialOmega = Omega;
    }
 
    void SetLeftConvergenceAcceleratorPointer(MVQNPointerType pConvergenceAcceleratorLeft)
    {
        mpConvergenceAcceleratorLeft = pConvergenceAcceleratorLeft;
    }
 
    void SetRightConvergenceAcceleratorPointer(MVQNPointerType pConvergenceAcceleratorRight)
    {
        mpConvergenceAcceleratorRight = pConvergenceAcceleratorRight;
    }
 
    VectorPointerType pGetUncorrectedForceVector()
    {
        return mpUncorrectedForceVector;
    }
 
    VectorPointerType pGetUncorrectedDisplacementVector()
    {
        return mpUncorrectedDisplacementVector;
    }
 
    MVQNPointerType pGetConvergenceAcceleratorLeft()
    {
        return mpConvergenceAcceleratorLeft;
    }
 
    MVQNPointerType pGetConvergenceAcceleratorRight()
    {
        return mpConvergenceAcceleratorRight;
    }
 
 
 
 
 
 
 
private:
 
 
 
    double mInitialOmega;
 
    bool mFirstLeftCorrectionPerformed = false; // Indicates that the initial fixed point iteration has been already performed
    bool mFirstRightCorrectionPerformed = false; // Indicates that the initial fixed point iteration has been already performed
    unsigned int mConvergenceAcceleratorIteration = 0; // Convergence accelerator iteration counter
 
    MVQNPointerType mpConvergenceAcceleratorLeft;
    MVQNPointerType mpConvergenceAcceleratorRight;
 
    VectorPointerType mpUncorrectedForceVector;
    VectorPointerType mpUncorrectedDisplacementVector;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
}; /* Class IBQNMVQNConvergenceAccelerator */
 
 
 
 
 
}  /* namespace Kratos.*/
 
#endif /* KRATOS_IBQN_MVQN_CONVERGENCE_ACCELERATOR defined */