residual_based_bdf_custom_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_CUSTOM_SCHEME )
#define  KRATOS_RESIDUAL_BASED_BDF_CUSTOM_SCHEME
 
/* System includes */
 
/* External includes */
 
/* Project includes */
#include "solving_strategies/schemes/residual_based_bdf_scheme.h"
#include "includes/variables.h"
#include "includes/kratos_parameters.h"
#include "includes/checks.h"
 
namespace Kratos
{
 
template<class TSparseSpace,  class TDenseSpace>
class ResidualBasedBDFCustomScheme
    : public ResidualBasedBDFScheme<TSparseSpace, TDenseSpace>
{
public:
    KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedBDFCustomScheme );
 
    typedef Scheme<TSparseSpace,TDenseSpace>                                  BaseType;
 
    typedef ResidualBasedImplicitTimeScheme<TSparseSpace,TDenseSpace> ImplicitBaseType;
 
    typedef ResidualBasedBDFScheme<TSparseSpace,TDenseSpace>               BDFBaseType;
 
    typedef ResidualBasedBDFCustomScheme<TSparseSpace, TDenseSpace>          ClassType;
 
    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;
 
    typedef ModelPart::ElementsContainerType                         ElementsArrayType;
 
    typedef ModelPart::ConditionsContainerType                     ConditionsArrayType;
 
    typedef typename BaseType::Pointer                                 BaseTypePointer;
 
 
    explicit ResidualBasedBDFCustomScheme(Parameters ThisParameters)
        :BDFBaseType()
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
 
        // Creating variables list
        CreateVariablesList(ThisParameters);
    }
 
    explicit ResidualBasedBDFCustomScheme(
        const std::size_t Order = 2,
        Parameters ThisParameters =  Parameters(R"({})")
        )
        :BDFBaseType(Order)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
 
        // Creating variables list
        CreateVariablesList(ThisParameters);
    }
 
    explicit ResidualBasedBDFCustomScheme(ResidualBasedBDFCustomScheme& rOther)
        :BDFBaseType(rOther)
        ,mDoubleVariable(rOther.mDoubleVariable)
        ,mFirstDoubleDerivatives(rOther.mFirstDoubleDerivatives)
        ,mSecondDoubleDerivatives(rOther.mSecondDoubleDerivatives)
    {
    }
 
    BaseTypePointer Clone() override
    {
        return BaseTypePointer( new ResidualBasedBDFCustomScheme(*this) );
    }
 
    ~ResidualBasedBDFCustomScheme
    () override {}
 
 
 
    typename BaseType::Pointer Create(Parameters ThisParameters) const override
    {
        return Kratos::make_shared<ClassType>(ThisParameters);
    }
 
    void Initialize(ModelPart& rModelPart) override
    {
        KRATOS_TRY
 
        BDFBaseType::Initialize(rModelPart);
 
        // The current process info
        const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
 
        // Getting dimension
        KRATOS_WARNING_IF("ResidualBasedBDFCustomScheme", !r_current_process_info.Has(DOMAIN_SIZE)) << "DOMAIN_SIZE not defined. Please define DOMAIN_SIZE. 3D case will be assumed" << std::endl;
        const std::size_t domain_size = r_current_process_info.Has(DOMAIN_SIZE) ? r_current_process_info.GetValue(DOMAIN_SIZE) : 3;
        if (domain_size != mDomainSize) {
            const std::size_t total_number_of_variables = mDoubleVariable.size();
 
            // We remove the third component
            if (domain_size == 2) {
                const std::size_t number_variables_added = total_number_of_variables/3;
                for (std::size_t i = 0; i < number_variables_added; ++i) {
                    mDoubleVariable.erase(mDoubleVariable.begin() + (2 + 2 * i));
                    mFirstDoubleDerivatives.erase(mDoubleVariable.begin() + (2 + 2 * i));
                    mSecondDoubleDerivatives.erase(mDoubleVariable.begin() + (2 + 2 * i));
                }
            } else if (domain_size == 3) { // We need to add the third component
                const std::size_t number_variables_added = total_number_of_variables/2;
                for (std::size_t i = 0; i < number_variables_added; ++i) {
                    const std::string variable_name = ((*(mDoubleVariable.begin() + 2 * i))->GetSourceVariable()).Name();
                    const auto& r_var_z = KratosComponents<Variable<double>>::Get(variable_name + "_Z");
                    mDoubleVariable.push_back(&r_var_z);
                    mFirstDoubleDerivatives.push_back(&(r_var_z.GetTimeDerivative()));
                    mSecondDoubleDerivatives.push_back(&((r_var_z.GetTimeDerivative()).GetTimeDerivative()));
                }
            } else {
                KRATOS_ERROR << "DOMAIN_SIZE can onbly be 2 or 3. It is: " << domain_size << std::endl;
            }
 
            mDomainSize = domain_size;
        }
 
        KRATOS_CATCH("")
    }
 
    void InitializeSolutionStep(
        ModelPart& rModelPart,
        TSystemMatrixType& rA,
        TSystemVectorType& rDx,
        TSystemVectorType& rb
        ) override
    {
        KRATOS_TRY;
 
        BDFBaseType::InitializeSolutionStep(rModelPart, rA, rDx, rb);
 
        // Updating time derivatives (nodally for efficiency)
        const int num_nodes = static_cast<int>( rModelPart.Nodes().size() );
        const auto it_node_begin = rModelPart.Nodes().begin();
 
        // Auxiliar fixed value
        bool fixed = false;
 
        IndexPartition<std::size_t>(num_nodes).for_each([&](std::size_t Index){
            auto it_node = it_node_begin + Index;
 
            std::size_t counter = 0;
            for (auto p_var : mDoubleVariable) {
 
                fixed = false;
 
                // Derivatives
                const auto& dvar = *mFirstDoubleDerivatives[counter];
                const auto& d2var = *mSecondDoubleDerivatives[counter];
 
                if (it_node->HasDofFor(d2var)) {
                    if (it_node->IsFixed(d2var)) {
                        it_node->Fix(*p_var);
                        fixed = true;
                    }
                }
 
                if (it_node->HasDofFor(dvar)) {
                    if (it_node->IsFixed(dvar) && !fixed) {
                        it_node->Fix(*p_var);
                    }
                }
 
                counter++;
            }
        });
 
        KRATOS_CATCH("ResidualBasedBDFCustomScheme.InitializeSolutionStep");
    }
 
    void Predict(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        TSystemMatrixType& A,
        TSystemVectorType& Dx,
        TSystemVectorType& b
        ) override
    {
        KRATOS_TRY;
 
        // Getting process info
        const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
 
        // Getting delta time
        const double delta_time = r_current_process_info[DELTA_TIME];
 
        // 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;
 
            std::size_t counter = 0;
            for (auto p_var : mDoubleVariable) {
                // Derivatives
                const auto& dvar = *mFirstDoubleDerivatives[counter];
                const auto& d2var = *mSecondDoubleDerivatives[counter];
 
                ComputePredictComponent(it_node, *p_var, dvar, d2var, delta_time);
 
                counter++;
            }
 
            // Updating time derivatives
            UpdateFirstDerivative(it_node);
            UpdateSecondDerivative(it_node);
        });
 
        KRATOS_CATCH( "" );
    }
 
    int Check(const ModelPart& rModelPart) const override
    {
        KRATOS_TRY;
 
        const int err = BDFBaseType::Check(rModelPart);
        if(err!=0) return err;
 
        // Check that variables are correctly allocated
        for(auto& r_node : rModelPart.Nodes()) {
            for ( auto p_var : mDoubleVariable)
                KRATOS_CHECK_VARIABLE_IN_NODAL_DATA((*p_var), r_node)
            for ( auto p_var : mFirstDoubleDerivatives)
                KRATOS_CHECK_VARIABLE_IN_NODAL_DATA((*p_var), r_node)
            for ( auto p_var : mSecondDoubleDerivatives)
                KRATOS_CHECK_VARIABLE_IN_NODAL_DATA((*p_var), r_node)
 
            for ( auto p_var : mDoubleVariable)
                KRATOS_CHECK_DOF_IN_NODE((*p_var), r_node)
        }
 
        KRATOS_CATCH( "" );
 
        return 0;
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"               : "bdf_scheme",
            "domain_size"        : 3,
            "integration_order"  : 2,
            "solution_variables" : ["DISPLACEMENT"]
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BDFBaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "bdf_scheme";
    }
 
 
 
 
    std::string Info() const override
    {
        return "ResidualBasedBDFCustomScheme";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
 
    std::vector<const Variable<double>*> mDoubleVariable;          
    std::vector<const Variable<double>*> mFirstDoubleDerivatives;   
    std::vector<const Variable<double>*> mSecondDoubleDerivatives; 
 
 
 
    inline void UpdateFirstDerivative(NodesArrayType::iterator itNode) override
    {
        // DOUBLES
        std::size_t counter = 0;
        for (auto p_var : mDoubleVariable) {
            double& dotun0 = itNode->FastGetSolutionStepValue(*mFirstDoubleDerivatives[counter]);
            dotun0 = BDFBaseType::mBDF[0] * itNode->FastGetSolutionStepValue(*p_var);
            for (std::size_t i_order = 1; i_order < BDFBaseType::mOrder + 1; ++i_order)
                dotun0 += BDFBaseType::mBDF[i_order] * itNode->FastGetSolutionStepValue(*p_var, i_order);
            counter++;
        }
    }
 
    inline void UpdateSecondDerivative(NodesArrayType::iterator itNode) override
    {
        // DOUBLES
        std::size_t counter = 0;
        for (auto p_var : mFirstDoubleDerivatives) {
            double& dot2un0 = itNode->FastGetSolutionStepValue(*mSecondDoubleDerivatives[counter]);
            dot2un0 = BDFBaseType::mBDF[0] * itNode->FastGetSolutionStepValue(*p_var);
            for (std::size_t i_order = 1; i_order < BDFBaseType::mOrder + 1; ++i_order)
                dot2un0 += BDFBaseType::mBDF[i_order] * itNode->FastGetSolutionStepValue(*p_var, i_order);
            counter++;
        }
    }
 
 
 
 
private:
 
 
 
    std::size_t mDomainSize = 3; 
 
 
 
    template<class TClassVar>
    void ComputePredictComponent(
        NodesArrayType::iterator itNode,
        const TClassVar& rVariable,
        const TClassVar& rDerivedVariable,
        const TClassVar& rDerived2Variable,
        const double DeltaTime
        )
    {
        // Values
        const double dot2un1 = itNode->FastGetSolutionStepValue(rDerived2Variable, 1);
        const double dotun1 = itNode->FastGetSolutionStepValue(rDerivedVariable, 1);
        const double un1 = itNode->FastGetSolutionStepValue(rVariable, 1);
        const double dot2un0 = itNode->FastGetSolutionStepValue(rDerived2Variable);
        double& dotun0 = itNode->FastGetSolutionStepValue(rDerivedVariable);
        double& un0 = itNode->FastGetSolutionStepValue(rVariable);
 
        if (itNode->HasDofFor(rDerived2Variable) && itNode->IsFixed(rDerived2Variable)) {
            dotun0 = dot2un0;
            for (std::size_t i_order = 1; i_order < BDFBaseType::mOrder + 1; ++i_order)
                dotun0 -= BDFBaseType::mBDF[i_order] * itNode->FastGetSolutionStepValue(rDerivedVariable, i_order);
            dotun0 /= BDFBaseType::mBDF[0];
 
            un0 = dotun0;
            for (std::size_t i_order = 1; i_order < BDFBaseType::mOrder + 1; ++i_order)
                un0 -= BDFBaseType::mBDF[i_order] * itNode->FastGetSolutionStepValue(rVariable, i_order);
            un0 /= BDFBaseType::mBDF[0];
        } else if (itNode->HasDofFor(rDerivedVariable) && itNode->IsFixed(rDerivedVariable)) {
            un0 = dotun0;
            for (std::size_t i_order = 1; i_order < BDFBaseType::mOrder + 1; ++i_order)
                un0 -= BDFBaseType::mBDF[i_order] * itNode->FastGetSolutionStepValue(rVariable, i_order);
            un0 /= BDFBaseType::mBDF[0];
        } else if (!itNode->IsFixed(rVariable)) {
            un0 = un1 + DeltaTime * dotun1 + 0.5 * std::pow(DeltaTime, 2) * dot2un1;
        }
    }
 
    void CreateVariablesList(Parameters ThisParameters)
    {
        const std::size_t n_variables = ThisParameters["solution_variables"].size();
 
        // The current dimension
        mDomainSize = ThisParameters["domain_size"].GetInt();
 
        const auto variable_names = ThisParameters["solution_variables"].GetStringArray();
 
        for (std::size_t p_var = 0; p_var < n_variables; ++p_var){
            const std::string& variable_name = variable_names[p_var];
 
            if(KratosComponents<Variable<double>>::Has(variable_name)){
                const auto& r_var = KratosComponents<Variable<double>>::Get(variable_name);
                mDoubleVariable.push_back(&r_var);
                mFirstDoubleDerivatives.push_back(&(r_var.GetTimeDerivative()));
                mSecondDoubleDerivatives.push_back(&((r_var.GetTimeDerivative()).GetTimeDerivative()));
            } else if (KratosComponents< Variable< array_1d< double, 3> > >::Has(variable_name)) {
                // Components
                const auto& r_var_x = KratosComponents<Variable<double>>::Get(variable_name+"_X");
                const auto& r_var_y = KratosComponents<Variable<double>>::Get(variable_name+"_Y");
                mDoubleVariable.push_back(&r_var_x);
                mDoubleVariable.push_back(&r_var_y);
                mFirstDoubleDerivatives.push_back(&(r_var_x.GetTimeDerivative()));
                mFirstDoubleDerivatives.push_back(&(r_var_y.GetTimeDerivative()));
                mSecondDoubleDerivatives.push_back(&((r_var_x.GetTimeDerivative()).GetTimeDerivative()));
                mSecondDoubleDerivatives.push_back(&((r_var_y.GetTimeDerivative()).GetTimeDerivative()));
                if (mDomainSize == 3) {
                    const auto& r_var_z = KratosComponents<Variable<double>>::Get(variable_name+"_Z");
                    mDoubleVariable.push_back(&r_var_z);
                    mFirstDoubleDerivatives.push_back(&(r_var_z.GetTimeDerivative()));
                    mSecondDoubleDerivatives.push_back(&((r_var_z.GetTimeDerivative()).GetTimeDerivative()));
                }
            } else {
                KRATOS_ERROR << "Only double and vector variables are allowed in the variables list." ;
            }
        }
    }
 
 
 
}; /* Class ResidualBasedBDFCustomScheme */
}  /* namespace Kratos.*/
 
#endif /* KRATOS_RESIDUAL_BASED_BDF_CUSTOM_SCHEME defined */