explicit_solving_strategy_bfecc.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Ruben Zorrilla, Eduard Gómez
//
//
 
#if !defined(KRATOS_EXPLICIT_SOLVING_STRATEGY_BFECC)
#define KRATOS_EXPLICIT_SOLVING_STRATEGY_BFECC
 
/* System includes */
 
 
/* External includes */
 
/* Project includes */
#include "includes/define.h"
#include "includes/model_part.h"
#include "factories/factory.h"
#include "solving_strategies/strategies/explicit_solving_strategy.h"
 
namespace Kratos
{
 
 
 
 
template <class TSparseSpace, class TDenseSpace>
class ExplicitSolvingStrategyBFECC : public ExplicitSolvingStrategy<TSparseSpace, TDenseSpace>
{
public:
 
    enum Substep { FORWARD, BACKWARD, FINAL };
 
    typedef ModelPart::SizeType SizeType;
 
    // The base solving strategy class definition
    typedef SolvingStrategy<TSparseSpace, TDenseSpace> SolvingStrategyType;
 
    // The base class definition
    typedef ExplicitSolvingStrategy<TSparseSpace, TDenseSpace> BaseType;
 
    typedef ExplicitSolvingStrategyBFECC<TSparseSpace, TDenseSpace> ClassType;
 
    // The explicit builder and solver definition
    typedef typename BaseType::ExplicitBuilderType ExplicitBuilderType;
 
    typedef typename BaseType::DofType DofType;
 
    typedef typename TDenseSpace::VectorType LocalSystemVectorType;
    typedef typename TDenseSpace::MatrixType LocalSystemMatrixType;
 
    KRATOS_CLASS_POINTER_DEFINITION(ExplicitSolvingStrategyBFECC);
 
    // Time-stepping settings
    struct SubstepData {
        enum Direction : int {BACK=-1, FORTH=1};
 
        SubstepData(const double Theta, const Direction Dir)
            : theta(Theta), direction(Dir)
        { }
 
        int TimeDirection() const noexcept { return static_cast<int>(direction); }
 
        double theta;
        Direction direction;
        LocalSystemVectorType u_fixed;
    };
 
 
    explicit ExplicitSolvingStrategyBFECC()
        : BaseType()
    {
    }
 
    explicit ExplicitSolvingStrategyBFECC(
        ModelPart &rModelPart,
        Parameters ThisParameters)
        : BaseType(rModelPart)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    explicit ExplicitSolvingStrategyBFECC(
        ModelPart &rModelPart,
        typename ExplicitBuilderType::Pointer pExplicitBuilder,
        bool MoveMeshFlag = false,
        int RebuildLevel = 0)
        : BaseType(rModelPart, pExplicitBuilder, MoveMeshFlag, RebuildLevel)
    {
    }
 
    explicit ExplicitSolvingStrategyBFECC(
        ModelPart &rModelPart,
        bool MoveMeshFlag = false,
        int RebuildLevel = 0)
        : BaseType(rModelPart, MoveMeshFlag, RebuildLevel)
    {
    }
 
    typename SolvingStrategyType::Pointer Create(
        ModelPart& rModelPart,
        Parameters ThisParameters
        ) const override
    {
        return Kratos::make_shared<ClassType>(rModelPart, ThisParameters);
    }
 
    ExplicitSolvingStrategyBFECC(const ExplicitSolvingStrategyBFECC &Other) = delete;
 
    ~ExplicitSolvingStrategyBFECC() override = default;
 
    virtual Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name" : "explicit_solving_strategy_bfecc"
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    virtual void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
    }
 
    static std::string Name()
    {
        std::stringstream s;
        s << "explicit_solving_strategy_bfecc";
        return s.str();
    }
 
 
 
 
    void Initialize() override
    {
        BaseType::Initialize();
        BaseType::GetModelPart().GetProcessInfo().SetValue(TIME_INTEGRATION_THETA, 0.0);
    }
 
 
    std::string Info() const override
    {
        std::stringstream ss;
        ss << "ExplicitSolvingStrategyBFECC";
        return ss.str();
    }
 
    void PrintInfo(std::ostream &rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream &rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
 
 
 
 
    void SolveWithLumpedMassMatrix() override
    {
        KRATOS_TRY
 
        // Validate time step size
        const double dt = BaseType::GetDeltaTime();
        KRATOS_ERROR_IF(dt < 1.0e-12) << "ProcessInfo DELTA_TIME is close to zero." << std::endl;
 
        // Stash the prev step solution
        const auto original_starting_values = ExtractSolutionStepData(1);
 
        PerformSubstep(FORWARD);
        PerformSubstep(BACKWARD);
 
        CorrectErrorAfterForwardsAndBackwards(original_starting_values);
        PerformSubstep(FINAL);
 
        KRATOS_CATCH("");
    }
 
    LocalSystemVectorType CopySolutionStepData(const SubstepData SData)
    {
        KRATOS_TRY
 
        const SizeType source      = SData.direction == SubstepData::FORTH ? 1 : 0;
        const SizeType destination = SData.direction == SubstepData::FORTH ? 0 : 1;
 
        // Get the required data from the explicit builder and solver
        auto& r_explicit_bs = BaseType::GetExplicitBuilder();
        auto& r_dof_set = r_explicit_bs.GetDofSet();
        const SizeType dof_size = r_explicit_bs.GetEquationSystemSize();
 
        LocalSystemVectorType u_fixed(dof_size);
        IndexPartition<SizeType>(r_dof_set.size()).for_each(
            [&](const SizeType i_dof){
                auto r_dof = *(r_dof_set.begin() + i_dof);
                double& r_u_0 = r_dof.GetSolutionStepValue(destination);
                if (r_dof.IsFixed()) {
                    u_fixed(i_dof) = r_u_0;
                }
                r_u_0 = r_dof.GetSolutionStepValue(source);
            }
        );
 
        return u_fixed;
 
        KRATOS_CATCH("")
    }
 
    LocalSystemVectorType ExtractSolutionStepData(const SizeType BufferPosition) const
    {
        // Get the required data from the explicit builder and solver
        auto& r_explicit_bs = BaseType::GetExplicitBuilder();
        auto& r_dof_set = r_explicit_bs.GetDofSet();
        const SizeType dof_size = r_explicit_bs.GetEquationSystemSize();
 
        LocalSystemVectorType u_copy(dof_size);
        IndexPartition<SizeType>(r_dof_set.size()).for_each(
            [&](const SizeType i_dof){
                const auto it_dof = r_dof_set.cbegin() + i_dof;
                u_copy[i_dof] = it_dof->GetSolutionStepValue(BufferPosition);
            }
        );
 
        return u_copy;
    }
 
    virtual void InitializeBFECCForwardSubstep() {};
 
    virtual void FinalizeBFECCForwardSubstep() {};
 
    virtual void InitializeBFECCBackwardSubstep() {};
 
    virtual void FinalizeBFECCBackwardSubstep() {};
 
    virtual void InitializeBFECCFinalSubstep() {};
 
    virtual void FinalizeBFECCFinalSubstep() {};
 
    virtual void PerformSubstep(const Substep SubstepType)
    {
        KRATOS_TRY
 
        // Get the required data from the explicit builder and solver
        auto& r_explicit_bs = BaseType::GetExplicitBuilder();
        auto& r_dof_set = r_explicit_bs.GetDofSet();
        const auto& r_lumped_mass_vector = r_explicit_bs.GetLumpedMassMatrixVector();
 
        // Get model part and information
        const double dt = BaseType::GetDeltaTime();
        KRATOS_ERROR_IF(dt < 1.0e-12) << "ProcessInfo DELTA_TIME is close to zero." << std::endl;
        auto& r_model_part = BaseType::GetModelPart();
        auto& r_process_info = r_model_part.GetProcessInfo();
 
        // Clone the previous step and initialize values
        const auto substep_settings = InitializeSubstep(SubstepType);
 
        r_process_info.GetValue(TIME_INTEGRATION_THETA) = substep_settings.theta;
        r_explicit_bs.BuildRHS(r_model_part);
 
        IndexPartition<SizeType>(r_dof_set.size()).for_each(
            [&](SizeType i_dof){
                auto it_dof = r_dof_set.begin() + i_dof;
 
                // Save current value in the corresponding vector
                const double residual = it_dof->GetSolutionStepReactionValue();
 
                // Do the DOF update
                double& r_u = it_dof->GetSolutionStepValue(0);
                double& r_u_prev_step = it_dof->GetSolutionStepValue(1);
 
                if (!it_dof->IsFixed()) {
                    const double mass = r_lumped_mass_vector[i_dof];
                    r_u += substep_settings.TimeDirection() * (dt / mass) * residual;
                } else {
                    r_u = substep_settings.theta*substep_settings.u_fixed[i_dof] + (1 - substep_settings.theta)*r_u_prev_step;
                }
            }
        );
 
        FinalizeSubstep(SubstepType);
 
        KRATOS_CATCH("Substep type: " + [](Substep substep) -> std::string {
            switch(substep) {
            case FORWARD:  return "FORWARD";
            case BACKWARD: return "BACKWARD";
            case FINAL:    return "FINAL";
            default: return std::to_string(static_cast<int>(substep));
            }
        }(SubstepType));
    }
 
    void CorrectErrorAfterForwardsAndBackwards(const LocalSystemVectorType& rPrevStepSolution)
    {
        auto& r_explicit_bs = BaseType::GetExplicitBuilder();
        auto& r_dof_set = r_explicit_bs.GetDofSet();
 
        IndexPartition<SizeType>(r_dof_set.size()).for_each(
            [&](const SizeType dof_index)
            {
                const double prev_step_solution = rPrevStepSolution[dof_index];
                auto& r_dof = *(r_dof_set.begin() + dof_index);
 
                const double error = (r_dof.GetSolutionStepValue(0) - prev_step_solution)/2.0;
 
                r_dof.GetSolutionStepValue(1) = prev_step_solution - error;
                // Will be copied to buffer position (0) during the last substep initialize
 
            }
        );
    }
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
    SubstepData InitializeSubstep(const Substep SubstepType)
    {
        switch(SubstepType)
        {
            case FORWARD:
            {
                SubstepData s = {0.0, SubstepData::FORTH};
                s.u_fixed = CopySolutionStepData(s);
                InitializeBFECCForwardSubstep();
                return s;
            }
            case BACKWARD:
            {
                SubstepData s = {1.0, SubstepData::BACK};
                s.u_fixed = CopySolutionStepData(s);
                InitializeBFECCBackwardSubstep();
                return s;
            }
            case FINAL:
            {
                SubstepData s = {0.0, SubstepData::FORTH};
                s.u_fixed = CopySolutionStepData(s);
                InitializeBFECCFinalSubstep();
                return s;
            }
            default:
                KRATOS_ERROR << "Invalid value for Substep" << std::endl;
        }
    }
 
    void  FinalizeSubstep(const Substep SubstepType)
    {
        switch(SubstepType) {
            case FORWARD:  FinalizeBFECCForwardSubstep();   break;
            case BACKWARD: FinalizeBFECCBackwardSubstep();  break;
            case FINAL:    FinalizeBFECCFinalSubstep();     break;
            default: KRATOS_ERROR << "Invalid value for Substep" << std::endl;
        }
    }
 
 
 
 
 
 
 
};
 
 
 
 
} /* namespace Kratos.*/
 
#endif /* KRATOS_EXPLICIT_SOLVING_STRATEGY_BFECC  defined */