residualbased_linear_strategy.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
//
 
#if !defined(KRATOS_RESIDUALBASED_LINEAR_STRATEGY )
#define  KRATOS_RESIDUALBASED_LINEAR_STRATEGY
 
// System includes
 
// External includes
 
// Project includes
#include "includes/define.h"
#include "solving_strategies/strategies/implicit_solving_strategy.h"
#include "utilities/builtin_timer.h"
 
//default builder and solver
#include "solving_strategies/builder_and_solvers/builder_and_solver.h"
#include "solving_strategies/builder_and_solvers/residualbased_block_builder_and_solver.h"
 
namespace Kratos
{
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace, //= DenseSpace<double>,
         class TLinearSolver //= LinearSolver<TSparseSpace,TDenseSpace>
         >
class ResidualBasedLinearStrategy
    : public ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(ResidualBasedLinearStrategy);
 
    typedef SolvingStrategy<TSparseSpace, TDenseSpace> SolvingStrategyType;
 
    typedef ImplicitSolvingStrategy<TSparseSpace, TDenseSpace, TLinearSolver> BaseType;
 
    typedef ResidualBasedLinearStrategy<TSparseSpace,TDenseSpace,TLinearSolver> ClassType;
 
    typedef typename BaseType::TDataType TDataType;
 
    typedef TSparseSpace SparseSpaceType;
 
    typedef typename BaseType::TSchemeType TSchemeType;
 
    typedef typename BaseType::TBuilderAndSolverType TBuilderAndSolverType;
 
    typedef typename BaseType::DofsArrayType DofsArrayType;
 
    typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
 
    typedef typename BaseType::TSystemVectorType TSystemVectorType;
 
    typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
 
    typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
 
    typedef typename BaseType::TSystemMatrixPointerType TSystemMatrixPointerType;
 
    typedef typename BaseType::TSystemVectorPointerType TSystemVectorPointerType;
 
 
 
    explicit ResidualBasedLinearStrategy() : BaseType()
    {
    }
 
    explicit ResidualBasedLinearStrategy(ModelPart& rModelPart, Parameters ThisParameters)
        : BaseType(rModelPart)
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
 
        // Set flags to start the calculations correctly
        mSolutionStepIsInitialized = false;
        mInitializeWasPerformed = false;
 
        // Tells to the builder and solver if the reactions have to be Calculated or not
        GetBuilderAndSolver()->SetCalculateReactionsFlag(mCalculateReactionsFlag);
 
        // Tells to the Builder And Solver if the system matrix and vectors need to
        // be reshaped at each step or not
        GetBuilderAndSolver()->SetReshapeMatrixFlag(mReformDofSetAtEachStep);
    }
 
    explicit ResidualBasedLinearStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TLinearSolver::Pointer pNewLinearSolver,
        bool CalculateReactionFlag = false,
        bool ReformDofSetAtEachStep = false,
        bool CalculateNormDxFlag = false,
        bool MoveMeshFlag = false
        ) : BaseType(rModelPart, MoveMeshFlag),
            mpScheme(pScheme),
            mReformDofSetAtEachStep(ReformDofSetAtEachStep),
            mCalculateNormDxFlag(CalculateNormDxFlag),
            mCalculateReactionsFlag(CalculateReactionFlag)
    {
        KRATOS_TRY
 
        // Setting up the default builder and solver
        mpBuilderAndSolver = typename TBuilderAndSolverType::Pointer
                             (
                                 new ResidualBasedBlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver > (pNewLinearSolver)
                             );
 
        // Set flag to start the calculations correctly
        mSolutionStepIsInitialized = false;
        mInitializeWasPerformed = false;
 
        // Tells to the builder and solver if the reactions have to be Calculated or not
        GetBuilderAndSolver()->SetCalculateReactionsFlag(mCalculateReactionsFlag);
 
        // Tells to the Builder And Solver if the system matrix and vectors need to
        //be reshaped at each step or not
        GetBuilderAndSolver()->SetReshapeMatrixFlag(mReformDofSetAtEachStep);
 
        // Set EchoLevel to the default value (only time is displayed)
        this->SetEchoLevel(1);
 
        // By default the matrices are rebuilt at each solution step
        BaseType::SetRebuildLevel(1);
 
        KRATOS_CATCH("")
    }
 
    explicit ResidualBasedLinearStrategy(
        ModelPart& rModelPart,
        typename TSchemeType::Pointer pScheme,
        typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver,
        bool CalculateReactionFlag = false,
        bool ReformDofSetAtEachStep = false,
        bool CalculateNormDxFlag = false,
        bool MoveMeshFlag = false
        ) : BaseType(rModelPart, MoveMeshFlag),
            mpScheme(pScheme),
            mpBuilderAndSolver(pNewBuilderAndSolver),
            mReformDofSetAtEachStep(ReformDofSetAtEachStep),
            mCalculateNormDxFlag(CalculateNormDxFlag),
            mCalculateReactionsFlag(CalculateReactionFlag)
    {
        KRATOS_TRY
 
        // Set flag to start the calculations correctly
        mSolutionStepIsInitialized = false;
        mInitializeWasPerformed = false;
 
        // Tells to the builder and solver if the reactions have to be Calculated or not
        GetBuilderAndSolver()->SetCalculateReactionsFlag(mCalculateReactionsFlag);
 
        // Tells to the Builder And Solver if the system matrix and vectors need to
        //be reshaped at each step or not
        GetBuilderAndSolver()->SetReshapeMatrixFlag(mReformDofSetAtEachStep);
 
        //set EchoLevel to the default value (only time is displayed)
        this->SetEchoLevel(1);
 
        // By default the matrices are rebuilt at each solution step
        BaseType::SetRebuildLevel(1);
 
        KRATOS_CATCH("")
    }
 
    ~ResidualBasedLinearStrategy() override
    {
        // If the linear solver has not been deallocated, clean it before
        // deallocating mpA. This prevents a memory error with the the ML
        // solver (which holds a reference to it).
        // NOTE: The linear solver is hold by the B&S
        auto p_builder_and_solver = this->GetBuilderAndSolver();
        if (p_builder_and_solver != nullptr) {
            p_builder_and_solver->Clear();
        }
 
        // Deallocating system vectors to avoid errors in MPI. Clear calls
        // TrilinosSpace::Clear for the vectors, which preserves the Map of
        // current vectors, performing MPI calls in the process. Due to the
        // way Python garbage collection works, this may happen after
        // MPI_Finalize has already been called and is an error. Resetting
        // the pointers here prevents Clear from operating with the
        // (now deallocated) vectors.
        mpA.reset();
        mpDx.reset();
        mpb.reset();
 
        this->Clear();
    }
 
    void SetScheme(typename TSchemeType::Pointer pScheme)
    {
        mpScheme = pScheme;
    };
 
    typename TSchemeType::Pointer GetScheme()
    {
        return mpScheme;
    };
 
    void SetBuilderAndSolver(typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver)
    {
        mpBuilderAndSolver = pNewBuilderAndSolver;
    };
 
    typename TBuilderAndSolverType::Pointer GetBuilderAndSolver()
    {
        return mpBuilderAndSolver;
    };
 
    void SetCalculateReactionsFlag(bool CalculateReactionsFlag)
    {
        mCalculateReactionsFlag = CalculateReactionsFlag;
        GetBuilderAndSolver()->SetCalculateReactionsFlag(mCalculateReactionsFlag);
    }
 
    bool GetCalculateReactionsFlag()
    {
        return mCalculateReactionsFlag;
    }
 
    void SetReformDofSetAtEachStepFlag(bool Flag)
    {
        mReformDofSetAtEachStep = Flag;
        GetBuilderAndSolver()->SetReshapeMatrixFlag(mReformDofSetAtEachStep);
    }
 
    bool GetReformDofSetAtEachStepFlag()
    {
        return mReformDofSetAtEachStep;
    }
 
    void SetEchoLevel(int Level) override
    {
        BaseType::SetEchoLevel(Level);
        GetBuilderAndSolver()->SetEchoLevel(Level);
    }
 
    //*********************************************************************************
    typename SolvingStrategyType::Pointer Create(
        ModelPart& rModelPart,
        Parameters ThisParameters
        ) const override
    {
        return Kratos::make_shared<ClassType>(rModelPart, ThisParameters);
    }
 
    void Predict() override
    {
        KRATOS_TRY
        const DataCommunicator &r_comm = BaseType::GetModelPart().GetCommunicator().GetDataCommunicator();
        //OPERATIONS THAT SHOULD BE DONE ONCE - internal check to avoid repetitions
        //if the operations needed were already performed this does nothing
        if(mInitializeWasPerformed == false)
            Initialize();
 
        //initialize solution step
        if (mSolutionStepIsInitialized == false)
            InitializeSolutionStep();
 
        TSystemMatrixType& rA  = *mpA;
        TSystemVectorType& rDx = *mpDx;
        TSystemVectorType& rb  = *mpb;
 
        DofsArrayType& r_dof_set = GetBuilderAndSolver()->GetDofSet();
 
        this->GetScheme()->Predict(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);
 
        // Applying constraints if needed
        auto& r_constraints_array = BaseType::GetModelPart().MasterSlaveConstraints();
        const int local_number_of_constraints = r_constraints_array.size();
        const int global_number_of_constraints = r_comm.SumAll(local_number_of_constraints);
        if(global_number_of_constraints != 0) {
            const auto& r_process_info = BaseType::GetModelPart().GetProcessInfo();
 
            block_for_each(r_constraints_array, [&r_process_info](MasterSlaveConstraint& rConstraint){
                rConstraint.ResetSlaveDofs(r_process_info);
            });
            block_for_each(r_constraints_array, [&r_process_info](MasterSlaveConstraint& rConstraint){
                rConstraint.Apply(r_process_info);
            });
 
            //the following is needed since we need to eventually compute time derivatives after applying
            //Master slave relations
            TSparseSpace::SetToZero(rDx);
            this->GetScheme()->Update(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);
        }
 
        if (BaseType::MoveMeshFlag() == true) BaseType::MoveMesh();
 
        KRATOS_CATCH("")
    }
 
    void Initialize() override
    {
        KRATOS_TRY
 
        if (mInitializeWasPerformed == false)
        {
            //pointers needed in the solution
            typename TSchemeType::Pointer p_scheme = GetScheme();
 
            //Initialize The Scheme - OPERATIONS TO BE DONE ONCE
            if (p_scheme->SchemeIsInitialized() == false)
                p_scheme->Initialize(BaseType::GetModelPart());
 
            //Initialize The Elements - OPERATIONS TO BE DONE ONCE
            if (p_scheme->ElementsAreInitialized() == false)
                p_scheme->InitializeElements(BaseType::GetModelPart());
 
            //Initialize The Conditions - OPERATIONS TO BE DONE ONCE
            if (p_scheme->ConditionsAreInitialized() == false)
                p_scheme->InitializeConditions(BaseType::GetModelPart());
 
            mInitializeWasPerformed = true;
        }
 
        KRATOS_CATCH("")
    }
 
    double Solve() override
    {
        BaseType::Solve();
 
        //calculate if needed the norm of Dx
        double norm_dx = 0.00;
        if (mCalculateNormDxFlag == true)
            norm_dx = TSparseSpace::TwoNorm(*mpDx);
 
        return norm_dx;
    }
 
    void Clear() override
    {
        KRATOS_TRY;
 
        // Setting to zero the internal flag to ensure that the dof sets are recalculated. Also clear the linear solver stored in the B&S
        auto p_builder_and_solver = GetBuilderAndSolver();
        if (p_builder_and_solver != nullptr) {
            p_builder_and_solver->SetDofSetIsInitializedFlag(false);
            p_builder_and_solver->Clear();
        }
 
        // Clearing the system of equations
        if (mpA != nullptr)
            SparseSpaceType::Clear(mpA);
        if (mpDx != nullptr)
            SparseSpaceType::Clear(mpDx);
        if (mpb != nullptr)
            SparseSpaceType::Clear(mpb);
 
        // Clearing scheme
        auto p_scheme = GetScheme();
        if (p_scheme != nullptr) {
            GetScheme()->Clear();
        }
 
        mInitializeWasPerformed = false;
        mSolutionStepIsInitialized = false;
 
        KRATOS_CATCH("");
    }
 
    void CalculateOutputData() override
    {
        TSystemMatrixType& rA  = *mpA;
        TSystemVectorType& rDx = *mpDx;
        TSystemVectorType& rb  = *mpb;
 
        GetScheme()->CalculateOutputData(BaseType::GetModelPart(),
                                         GetBuilderAndSolver()->GetDofSet(),
                                         rA, rDx, rb);
    }
 
    void InitializeSolutionStep() override
    {
        KRATOS_TRY
 
        if (mSolutionStepIsInitialized == false)
        {
            //pointers needed in the solution
            typename TSchemeType::Pointer p_scheme = GetScheme();
            typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver();
 
            //set up the system, operation performed just once unless it is required
            //to reform the dof set at each iteration
            BuiltinTimer system_construction_time;
            if (p_builder_and_solver->GetDofSetIsInitializedFlag() == false ||
                    mReformDofSetAtEachStep == true)
            {
                //setting up the list of the DOFs to be solved
                BuiltinTimer setup_dofs_time;
                p_builder_and_solver->SetUpDofSet(p_scheme, BaseType::GetModelPart());
                KRATOS_INFO_IF("ResidualBasedLinearStrategy", BaseType::GetEchoLevel() > 0) << "Setup Dofs Time: " << setup_dofs_time << std::endl;
 
                //shaping correctly the system
                BuiltinTimer setup_system_time;
                p_builder_and_solver->SetUpSystem(BaseType::GetModelPart());
                KRATOS_INFO_IF("ResidualBasedLinearStrategy", BaseType::GetEchoLevel() > 0) << "Setup System Time: " << setup_system_time << std::endl;
 
                //setting up the Vectors involved to the correct size
                BuiltinTimer system_matrix_resize_time;
                p_builder_and_solver->ResizeAndInitializeVectors(p_scheme, mpA, mpDx, mpb,
                                                                 BaseType::GetModelPart());
                KRATOS_INFO_IF("ResidualBasedLinearStrategy", BaseType::GetEchoLevel() > 0) << "System Matrix Resize Time: " << system_matrix_resize_time << std::endl;
            }
 
            KRATOS_INFO_IF("ResidualBasedLinearStrategy", BaseType::GetEchoLevel() > 0) << "System Construction Time: " << system_construction_time << std::endl;
 
            TSystemMatrixType& rA  = *mpA;
            TSystemVectorType& rDx = *mpDx;
            TSystemVectorType& rb  = *mpb;
 
            //initial operations ... things that are constant over the Solution Step
            p_builder_and_solver->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb);
 
            //initial operations ... things that are constant over the Solution Step
            p_scheme->InitializeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb);
 
            mSolutionStepIsInitialized = true;
        }
 
        KRATOS_CATCH("")
    }
 
    void FinalizeSolutionStep() override
    {
        KRATOS_TRY;
        typename TSchemeType::Pointer p_scheme = GetScheme();
        typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver();
 
        TSystemMatrixType &rA  = *mpA;
        TSystemVectorType &rDx = *mpDx;
        TSystemVectorType &rb  = *mpb;
 
        //Finalisation of the solution step,
        //operations to be done after achieving convergence, for example the
        //Final Residual Vector (mb) has to be saved in there
        //to avoid error accumulation
 
        p_scheme->FinalizeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb);
        p_builder_and_solver->FinalizeSolutionStep(BaseType::GetModelPart(), rA, rDx, rb);
 
        //Cleaning memory after the solution
        p_scheme->Clean();
 
        //reset flags for next step
        mSolutionStepIsInitialized = false;
 
        //deallocate the systemvectors if needed
        if (mReformDofSetAtEachStep == true)
        {
            SparseSpaceType::Clear(mpA);
            SparseSpaceType::Clear(mpDx);
            SparseSpaceType::Clear(mpb);
 
            this->Clear();
        }
 
        KRATOS_CATCH("");
    }
 
    bool SolveSolutionStep() override
    {
        //pointers needed in the solution
        typename TSchemeType::Pointer p_scheme = GetScheme();
        typename TBuilderAndSolverType::Pointer p_builder_and_solver = GetBuilderAndSolver();
 
        TSystemMatrixType& rA  = *mpA;
        TSystemVectorType& rDx = *mpDx;
        TSystemVectorType& rb  = *mpb;
 
        p_scheme->InitializeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb);
 
        if (BaseType::mRebuildLevel > 0 || BaseType::mStiffnessMatrixIsBuilt == false)
        {
            TSparseSpace::SetToZero(rA);
            TSparseSpace::SetToZero(rDx);
            TSparseSpace::SetToZero(rb);
            // passing smart pointers instead of references here
            // to prevent dangling pointer to system matrix when
            // reusing ml preconditioners in the trilinos tpl
            p_builder_and_solver->BuildAndSolve(p_scheme, BaseType::GetModelPart(), rA, rDx, rb);
            BaseType::mStiffnessMatrixIsBuilt = true;
        }
        else
        {
            TSparseSpace::SetToZero(rDx);
            TSparseSpace::SetToZero(rb);
            p_builder_and_solver->BuildRHSAndSolve(p_scheme, BaseType::GetModelPart(), rA, rDx, rb);
        }
 
        // Debugging info
        EchoInfo();
 
        //update results
        DofsArrayType& r_dof_set = p_builder_and_solver->GetDofSet();
        p_scheme->Update(BaseType::GetModelPart(), r_dof_set, rA, rDx, rb);
 
        //move the mesh if needed
        if (BaseType::MoveMeshFlag() == true) BaseType::MoveMesh();
 
        p_scheme->FinalizeNonLinIteration(BaseType::GetModelPart(), rA, rDx, rb);
 
        // Calculate reactions if required
        if (mCalculateReactionsFlag == true)
            p_builder_and_solver->CalculateReactions(p_scheme,
                                                     BaseType::GetModelPart(),
                                                     rA, rDx, rb);
 
        return true;
    }
 
    TSystemMatrixType& GetSystemMatrix() override
    {
        TSystemMatrixType& mA = *mpA;
 
        return mA;
    }
 
    TSystemVectorType& GetSystemVector() override
    {
        TSystemVectorType& mb = *mpb;
 
        return mb;
    }
 
    TSystemVectorType& GetSolutionVector() override
    {
        TSystemVectorType& mDx = *mpDx;
 
        return mDx;
    }
 
    double GetResidualNorm() override
    {
        if (TSparseSpace::Size(*mpb) != 0)
            return TSparseSpace::TwoNorm(*mpb);
        else
            return 0.0;
    }
 
    int Check() override
    {
        KRATOS_TRY
 
        BaseType::Check();
 
        GetBuilderAndSolver()->Check(BaseType::GetModelPart());
 
        GetScheme()->Check(BaseType::GetModelPart());
 
        return 0;
 
        KRATOS_CATCH("")
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"                         : "linear_strategy",
            "compute_norm_dx"              : false,
            "reform_dofs_at_each_step"     : false,
            "compute_reactions"            : false,
            "builder_and_solver_settings"  : {},
            "linear_solver_settings"       : {},
            "scheme_settings"              : {}
        })");
 
        // Getting base class default parameters
        const Parameters base_default_parameters = BaseType::GetDefaultParameters();
        default_parameters.RecursivelyAddMissingParameters(base_default_parameters);
        return default_parameters;
    }
 
    static std::string Name()
    {
        return "linear_strategy";
    }
 
 
 
 
 
 
 
 
    std::string Info() const override
    {
        return "ResidualBasedLinearStrategy";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
 
 
protected:
 
 
 
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
        mCalculateNormDxFlag = ThisParameters["compute_norm_dx"].GetBool();
        mReformDofSetAtEachStep = ThisParameters["reform_dofs_at_each_step"].GetBool();
        mCalculateReactionsFlag = ThisParameters["compute_reactions"].GetBool();
 
        // Saving the scheme
        if (ThisParameters["scheme_settings"].Has("name")) {
            KRATOS_ERROR << "IMPLEMENTATION PENDING IN CONSTRUCTOR WITH PARAMETERS" << std::endl;
        }
 
        // Setting up the default builder and solver
        if (ThisParameters["builder_and_solver_settings"].Has("name")) {
            KRATOS_ERROR << "IMPLEMENTATION PENDING IN CONSTRUCTOR WITH PARAMETERS" << std::endl;
        }
    }
 
 
 
 
 
private:
 
 
 
    typename TSchemeType::Pointer mpScheme = nullptr; 
    typename TBuilderAndSolverType::Pointer mpBuilderAndSolver = nullptr; 
 
    TSystemVectorPointerType mpDx; 
    TSystemVectorPointerType mpb; 
    TSystemMatrixPointerType mpA; 
 
    bool mReformDofSetAtEachStep;
 
    bool mCalculateNormDxFlag; 
 
    bool mCalculateReactionsFlag;
 
    bool mSolutionStepIsInitialized; 
 
    bool mInitializeWasPerformed; 
 
 
    virtual void EchoInfo()
    {
        TSystemMatrixType& rA  = *mpA;
        TSystemVectorType& rDx = *mpDx;
        TSystemVectorType& rb  = *mpb;
 
        if (BaseType::GetEchoLevel() == 3) //if it is needed to print the debug info
        {
            KRATOS_INFO("LHS") << "SystemMatrix = " << rA << std::endl;
            KRATOS_INFO("Dx")  << "Solution obtained = " << rDx << std::endl;
            KRATOS_INFO("RHS") << "RHS  = " << rb << std::endl;
        }
        if (this->GetEchoLevel() == 4) //print to matrix market file
        {
            std::stringstream matrix_market_name;
            matrix_market_name << "A_" << BaseType::GetModelPart().GetProcessInfo()[TIME] <<  ".mm";
            TSparseSpace::WriteMatrixMarketMatrix((char*) (matrix_market_name.str()).c_str(), rA, false);
 
            std::stringstream matrix_market_vectname;
            matrix_market_vectname << "b_" << BaseType::GetModelPart().GetProcessInfo()[TIME] << ".mm.rhs";
            TSparseSpace::WriteMatrixMarketVector((char*) (matrix_market_vectname.str()).c_str(), rb);
        }
    }
 
 
 
 
 
 
 
 
    ResidualBasedLinearStrategy(const ResidualBasedLinearStrategy& Other);
 
 
}; /* Class ResidualBasedLinearStrategy */
 
 
 
 
 
} /* namespace Kratos.*/
 
#endif /* KRATOS_RESIDUALBASED_LINEAR_STRATEGY  defined */