residual_based_bdf_scheme.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Vicente Mataix Ferrandiz
//
 
 
#if !defined(KRATOS_RESIDUAL_BASED_BDF_SCHEME )
#define  KRATOS_RESIDUAL_BASED_BDF_SCHEME
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "includes/checks.h"
#include "utilities/time_discretization.h"
#include "solving_strategies/schemes/residual_based_implicit_time_scheme.h"
 
namespace Kratos
{
 
template<class TSparseSpace,  class TDenseSpace>
class ResidualBasedBDFScheme
    : public ResidualBasedImplicitTimeScheme<TSparseSpace, TDenseSpace>
{
public:
    KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedBDFScheme );
 
    typedef Scheme<TSparseSpace,TDenseSpace>                                  BaseType;
 
    typedef typename BaseType::Pointer                                 BaseTypePointer;
 
    typedef ResidualBasedImplicitTimeScheme<TSparseSpace,TDenseSpace> ImplicitBaseType;
 
    typedef typename ImplicitBaseType::TDataType                             TDataType;
 
    typedef typename ImplicitBaseType::DofsArrayType                     DofsArrayType;
 
    typedef typename Element::DofsVectorType                            DofsVectorType;
 
    typedef typename ImplicitBaseType::TSystemMatrixType             TSystemMatrixType;
 
    typedef typename ImplicitBaseType::TSystemVectorType             TSystemVectorType;
 
    typedef typename ImplicitBaseType::LocalSystemVectorType     LocalSystemVectorType;
 
    typedef typename ImplicitBaseType::LocalSystemMatrixType     LocalSystemMatrixType;
 
    typedef ModelPart::NodesContainerType                               NodesArrayType;
 
    static constexpr double ZeroTolerance = std::numeric_limits<double>::epsilon();
 
 
    explicit ResidualBasedBDFScheme(const std::size_t Order = 2)
        :ImplicitBaseType(),
         mOrder(Order),
         mpBDFUtility(Kratos::make_unique<TimeDiscretization::BDF>(Order))
    {
        // Allocate auxiliary memory
        const std::size_t num_threads = ParallelUtilities::GetNumThreads();
 
        mVector.dotun0.resize(num_threads);
        mVector.dot2un0.resize(num_threads);
 
        // Doing a minimal check
        KRATOS_ERROR_IF(mOrder < 1) << "ERROR:: Not possible to compute a BDF of order less than 1" << std::endl;
 
        // We resize the BDF coefficients
        if (mBDF.size() != (mOrder + 1))
            mBDF.resize(mOrder + 1);
    }
 
    explicit ResidualBasedBDFScheme(ResidualBasedBDFScheme& rOther)
        :ImplicitBaseType(rOther)
        ,mOrder(rOther.mOrder)
        ,mBDF(rOther.mBDF)
        ,mVector(rOther.mVector)
        ,mpBDFUtility(nullptr)
    {
        Kratos::unique_ptr<TimeDiscretization::BDF> auxiliar_pointer = Kratos::make_unique<TimeDiscretization::BDF>(mOrder);
        mpBDFUtility.swap(auxiliar_pointer);
    }
 
    BaseTypePointer Clone() override
    {
        return BaseTypePointer( new ResidualBasedBDFScheme(*this) );
    }
 
    ~ResidualBasedBDFScheme
    () override {}
 
 
 
    void Update(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        ) override
    {
        KRATOS_TRY;
 
        // Update of displacement (by DOF)
        mpDofUpdater->UpdateDofs(rDofSet, rDx);
 
        UpdateDerivatives(rModelPart, rDofSet, rA, rDx, rb);
 
        KRATOS_CATCH( "" );
    }
 
    void Predict(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        ) override
    {
        KRATOS_TRY;
 
        KRATOS_ERROR << "Calling base BDF class" << std::endl;
 
        KRATOS_CATCH( "" );
    }
 
    void InitializeSolutionStep(
        ModelPart& rModelPart,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        ) override
    {
        KRATOS_TRY;
 
        ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
 
        ImplicitBaseType::InitializeSolutionStep(rModelPart, rA, rDx, rb);
 
        mpBDFUtility->ComputeAndSaveBDFCoefficients(r_current_process_info);
        mBDF = r_current_process_info[BDF_COEFFICIENTS];
 
        const double dt_0 = r_current_process_info[DELTA_TIME];
        const double dt_1 = r_current_process_info.GetPreviousTimeStepInfo(1)[DELTA_TIME];
        KRATOS_ERROR_IF(mOrder > 2 && std::abs(dt_0 - dt_1) > 1e-10*(dt_0 + dt_1))
        << "ResidualBasedBDFScheme. For higher orders than 2 the time step must be constant.\nPrevious time step : " << dt_1 << "\nCurrent time step : " << dt_0 << std::endl;
 
        KRATOS_CATCH( "" );
    }
 
    int Check(const ModelPart& rModelPart) const override
    {
        KRATOS_TRY;
 
        const int err = ImplicitBaseType::Check(rModelPart);
        if(err!=0) return err;
 
        // Check for minimum value of the buffer index
        // Verify buffer size
        KRATOS_ERROR_IF(rModelPart.GetBufferSize() < mOrder + 1) << "Insufficient buffer size. Buffer size should be greater than " << mOrder + 1 << ". Current size is " << rModelPart.GetBufferSize() << std::endl;
 
        KRATOS_CATCH( "" );
 
        return 0;
    }
 
    void Clear() override
    {
        this->mpDofUpdater->Clear();
    }
 
 
 
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"               : "base_bdf_scheme",
            "integration_order"  : 2
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = ImplicitBaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    std::string Info() const override
    {
        return "ResidualBasedBDFScheme";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
 
    struct GeneralVectors
    {
        std::vector< Vector > dotun0; 
        std::vector< Vector > dot2un0; 
    };
 
    const std::size_t mOrder;                      
    Vector mBDF;                                   
    GeneralVectors mVector;                        
    Kratos::unique_ptr<TimeDiscretization::BDF> mpBDFUtility; 
 
 
 
 
    inline void UpdateDerivatives(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        )
    {
        // Updating time derivatives (nodally for efficiency)
        const int num_nodes = static_cast<int>( rModelPart.Nodes().size() );
 
        // Getting first node iterator
        const auto it_node_begin = rModelPart.Nodes().begin();
 
        IndexPartition<std::size_t>(num_nodes).for_each([&](std::size_t Index){
            auto it_node = it_node_begin + Index;
 
            UpdateFirstDerivative(it_node);
            UpdateSecondDerivative(it_node);
        });
    }
 
    virtual inline void UpdateFirstDerivative(NodesArrayType::iterator itNode)
    {
        KRATOS_ERROR << "Calling base BDF class" << std::endl;
    }
 
    virtual inline void UpdateSecondDerivative(NodesArrayType::iterator itNode)
    {
        KRATOS_ERROR << "Calling base BDF class" << std::endl;
    }
 
    void AddDynamicsToLHS(
        LocalSystemMatrixType& rLHS_Contribution,
        LocalSystemMatrixType& rD,
        LocalSystemMatrixType& rM,
        const ProcessInfo& rCurrentProcessInfo
        ) override
    {
        // Adding mass contribution to the dynamic stiffness
        if (rM.size1() != 0) { // if M matrix declared
            noalias(rLHS_Contribution) += rM * std::pow(mBDF[0], 2);
        }
 
        // Adding  damping contribution
        if (rD.size1() != 0) { // if D matrix declared
            noalias(rLHS_Contribution) += rD * mBDF[0];
        }
    }
 
    template <class TObjectType>
    void TemplateAddDynamicsToRHS(
        TObjectType& rObject,
        LocalSystemVectorType& rRHS_Contribution,
        LocalSystemMatrixType& rD,
        LocalSystemMatrixType& rM,
        const ProcessInfo& rCurrentProcessInfo
        )
    {
        const std::size_t this_thread = OpenMPUtils::ThisThread();
        const auto& r_const_obj_ref = rObject;
 
        // Adding inertia contribution
        if (rM.size1() != 0) {
            r_const_obj_ref.GetSecondDerivativesVector(mVector.dot2un0[this_thread], 0);
            noalias(rRHS_Contribution) -= prod(rM, mVector.dot2un0[this_thread]);
        }
 
        // Adding damping contribution
        if (rD.size1() != 0) {
            r_const_obj_ref.GetFirstDerivativesVector(mVector.dotun0[this_thread], 0);
            noalias(rRHS_Contribution) -= prod(rD, mVector.dotun0[this_thread]);
        }
    }
 
    void AddDynamicsToRHS(
        Element& rElement,
        LocalSystemVectorType& rRHS_Contribution,
        LocalSystemMatrixType& rD,
        LocalSystemMatrixType& rM,
        const ProcessInfo& rCurrentProcessInfo
        ) override
    {
        TemplateAddDynamicsToRHS<Element>(rElement, rRHS_Contribution, rD, rM, rCurrentProcessInfo);
    }
 
    void AddDynamicsToRHS(
        Condition& rCondition,
        LocalSystemVectorType& rRHS_Contribution,
        LocalSystemMatrixType& rD,
        LocalSystemMatrixType& rM,
        const ProcessInfo& rCurrentProcessInfo
        ) override
    {
        TemplateAddDynamicsToRHS<Condition>(rCondition, rRHS_Contribution, rD, rM, rCurrentProcessInfo);
    }
 
 
 
 
private:
 
 
    typename TSparseSpace::DofUpdaterPointerType mpDofUpdater = TSparseSpace::CreateDofUpdater();
 
 
 
 
 
}; /* Class ResidualBasedBDFScheme */
}  /* namespace Kratos.*/
 
#endif /* KRATOS_RESIDUAL_BASED_BDF_SCHEME defined */