linear_strategy.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:               March 2018 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_LINEAR_STRATEGY_H_INCLUDED)
#define KRATOS_LINEAR_STRATEGY_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_solvers/solution_strategies/solution_strategy.hpp"
 
//default builder and solver
#include "custom_solvers/solution_builders_and_solvers/block_builder_and_solver.hpp"
 
namespace Kratos
{
 
 
 
 
 
 
template <class TSparseSpace,
          class TDenseSpace,  // = DenseSpace<double>,
          class TLinearSolver // = LinearSolver<TSparseSpace,TDenseSpace>
          >
class LinearStrategy : public SolutionStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
{
 public:
 
  // Counted pointer of ClassName
  KRATOS_CLASS_POINTER_DEFINITION(LinearStrategy);
 
  typedef SolutionStrategy<TSparseSpace, TDenseSpace, TLinearSolver>            BaseType;
 
  typedef typename BaseType::LocalFlagType                                 LocalFlagType;
 
  typedef typename BaseType::BuilderAndSolverType                   BuilderAndSolverType;
 
  typedef typename BaseType::SchemeType                                       SchemeType;
 
  typedef TLinearSolver                                                 LinearSolverType;
 
  typedef TSparseSpace                                                   SparseSpaceType;
 
  typedef typename BaseType::DofsArrayType                                 DofsArrayType;
 
  typedef typename BaseType::SystemMatrixType                           SystemMatrixType;
 
  typedef typename BaseType::SystemVectorType                           SystemVectorType;
 
  typedef typename BaseType::SystemMatrixPointerType             SystemMatrixPointerType;
 
  typedef typename BaseType::SystemVectorPointerType             SystemVectorPointerType;
 
 
 
  LinearStrategy(ModelPart& rModelPart,
                 typename SchemeType::Pointer pScheme,
                 typename BuilderAndSolverType::Pointer pBuilderAndSolver,
                 Flags& rOptions)
      : SolutionStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, rOptions)
  {
    KRATOS_TRY
 
    // Saving the scheme
    mpScheme = pScheme;
 
    // Setting up the default builder and solver
    mpBuilderAndSolver = pBuilderAndSolver;
 
    // Set Options
    this->mOptions = rOptions;
 
    // Tells to the builder and solver if the reactions have to be Calculated or not
    mpBuilderAndSolver->GetOptions().Set(LocalFlagType::COMPUTE_REACTIONS, this->mOptions.Is(LocalFlagType::COMPUTE_REACTIONS));
 
    // Tells to the Builder And Solver if the system matrix and vectors need to be reshaped at each step or not
    mpBuilderAndSolver->GetOptions().Set(LocalFlagType::REFORM_DOFS, this->mOptions.Is(LocalFlagType::REFORM_DOFS));
 
    mpBuilderAndSolver->SetEchoLevel(this->mEchoLevel);
 
    mpA  = TSparseSpace::CreateEmptyMatrixPointer();
    mpDx = TSparseSpace::CreateEmptyVectorPointer();
    mpb  = TSparseSpace::CreateEmptyVectorPointer();
 
    KRATOS_CATCH("")
  }
 
 
  LinearStrategy(ModelPart& rModelPart,
                 typename SchemeType::Pointer pScheme,
                 typename LinearSolverType::Pointer pLinearSolver,
                 Flags& rOptions)
      : LinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, Kratos::make_shared<BlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver> >(pLinearSolver), rOptions)
  {}
 
 
  ~LinearStrategy() override
  {
    Clear();
  }
 
 
 
 
  void InitializeSolutionStep() override
  {
    KRATOS_TRY
 
    //initialize system elements, conditions..
    if(this->IsNot(LocalFlagType::INITIALIZED))
      this->Initialize();
 
    //prints informations about the current time
    //KRATOS_INFO("") << "  [STEP:" << this->GetModelPart().GetProcessInfo()[STEP] << "  TIME: "<< this->GetModelPart().GetProcessInfo()[TIME]<< "]\n" << LoggerMessage::Category::STATUS;
 
    //set up the system
    if(this->mOptions.Is(LocalFlagType::REFORM_DOFS)){
      this->Clear();
      this->SetSystemDofs();
    }
 
    //setting up the vectors involved to the correct size
    // // double begin_time = OpenMPUtils::GetCurrentTime();
    // mpBuilderAndSolver->SetUpSystemMatrices(mpScheme, this->GetModelPart(), mpA, mpDx, mpb);
    // double end_time = OpenMPUtils::GetCurrentTime();
 
    // if (this->mEchoLevel >= 2)
    //   KRATOS_INFO("system_resize_time") << ": system_resize_time : " << end_time - begin_time << "\n" << LoggerMessage::Category::STATISTICS;
 
    //initial operations ... things that are constant over the Solution Step
    mpBuilderAndSolver->InitializeSolutionStep(mpScheme, this->GetModelPart(), mpA, mpDx, mpb);
 
    //initial operations ... things that are constant over the Solution Step
    mpScheme->InitializeSolutionStep(this->GetModelPart());
 
    this->Predict();
 
    KRATOS_CATCH("")
  }
 
  void FinalizeSolutionStep() override
  {
    KRATOS_TRY
 
    //Finalization of the solution step, operations to be done after achieving convergence
 
    //calculate reactions if required
    if(this->mOptions.Is(LocalFlagType::COMPUTE_REACTIONS))
      mpBuilderAndSolver->CalculateReactions(mpScheme, this->GetModelPart(), (*mpA), (*mpDx), (*mpb));
 
    //finalize scheme anb builder and solver
    mpScheme->FinalizeSolutionStep(this->GetModelPart());
    mpBuilderAndSolver->FinalizeSolutionStep(mpScheme, this->GetModelPart(), mpA, mpDx, mpb);
 
    KRATOS_CATCH("")
  }
 
 
  bool SolveSolutionStep() override
  {
    KRATOS_TRY
 
    return this->SolveIteration();
 
    KRATOS_CATCH("")
  }
 
 
  bool SolveIteration() override
  {
    KRATOS_TRY
 
    // Warning info
    if(!SparseSpaceType::Size((*mpDx)))
      KRATOS_WARNING("DOFS") << "solution has zero size, no free DOFs" << std::endl;
 
    // Initialize Iteration
    mpScheme->InitializeNonLinearIteration(this->GetModelPart());
 
    //function to perform the building and the solving phase.
    if(this->mOptions.IsNot(LocalFlagType::CONSTANT_SYSTEM_MATRIX)){
 
      TSparseSpace::SetToZero((*mpA));
      TSparseSpace::SetToZero((*mpDx));
      TSparseSpace::SetToZero((*mpb));
 
      mpBuilderAndSolver->BuildAndSolve(mpScheme, this->GetModelPart(), (*mpA), (*mpDx), (*mpb));
    }
    else{
 
      TSparseSpace::SetToZero((*mpDx));
      TSparseSpace::SetToZero((*mpb));
 
      mpBuilderAndSolver->BuildRHSAndSolve(mpScheme, this->GetModelPart(), (*mpA), (*mpDx), (*mpb));
    }
 
    // EchoInfo
    mpBuilderAndSolver->EchoInfo(this->GetModelPart(), (*mpA), (*mpDx), (*mpb));
 
    // Updating the results
    this->Update();
 
    // Finalize Iteration
    mpScheme->FinalizeNonLinearIteration(this->GetModelPart());
 
    return true;
 
    KRATOS_CATCH("")
  }
 
 
  void Clear() override
  {
    KRATOS_TRY
 
    // if the preconditioner is saved between solves, it should be cleared here.
    mpBuilderAndSolver->GetLinearSystemSolver()->Clear();
 
    //deallocate the systemvectors
    if(mpA != nullptr)
      SparseSpaceType::Clear(mpA);
    if(mpDx != nullptr)
      SparseSpaceType::Clear(mpDx);
    if(mpb != nullptr)
      SparseSpaceType::Clear(mpb);
 
    mpBuilderAndSolver->Clear();
    mpScheme->Clear();
 
    KRATOS_CATCH("")
  }
 
  int Check() override
  {
    KRATOS_TRY
 
    //check the model part
    BaseType::Check();
 
    //check the scheme
    mpScheme->Check(this->GetModelPart());
 
    //check the builder and solver
    mpBuilderAndSolver->Check(this->GetModelPart());
 
    return 0;
 
    KRATOS_CATCH("")
  }
 
 
 
  void SetEchoLevel(const int Level) override
  {
    BaseType::SetEchoLevel(Level);
    mpBuilderAndSolver->SetEchoLevel(Level);
  }
 
 
  void SetScheme(typename SchemeType::Pointer pScheme)
  {
    mpScheme = pScheme;
  };
 
  typename SchemeType::Pointer GetScheme()
  {
    return mpScheme;
  };
 
  void SetBuilderAndSolver(typename BuilderAndSolverType::Pointer pBuilderAndSolver)
  {
    mpBuilderAndSolver = pBuilderAndSolver;
  };
 
  typename BuilderAndSolverType::Pointer GetBuilderAndSolver()
  {
    return mpBuilderAndSolver;
  };
 
 
  SystemMatrixType &GetSystemMatrix()
  {
    SystemMatrixType &mA = *mpA;
 
    return mA;
  }
 
  void GetDirectSystemMatrix(SystemMatrixType& A)
  {
    A = *mpA;
  }
 
 
 
 
 
 protected:
 
 
  typename SchemeType::Pointer mpScheme; 
  typename BuilderAndSolverType::Pointer mpBuilderAndSolver; 
 
  SystemVectorPointerType mpDx; 
  SystemVectorPointerType mpb; 
  SystemMatrixPointerType mpA; 
 
 
 
  void Initialize() override
  {
    KRATOS_TRY
 
    //Initialize The Scheme - OPERATIONS TO BE DONE ONCE
    if( mpScheme->IsNot(LocalFlagType::INITIALIZED) )
      mpScheme->Initialize(this->GetModelPart());
 
    //set up the system
    this->SetSystemDofs();
 
    this->Set(LocalFlagType::INITIALIZED,true);
 
    KRATOS_CATCH("")
  }
 
 
  void Predict() override
  {
    KRATOS_TRY
 
    mpScheme->Predict(this->GetModelPart(), mpBuilderAndSolver->GetDofSet(), (*mpDx));
 
    KRATOS_CATCH("")
  }
 
 
  void Update() override
  {
    KRATOS_TRY
 
    mpScheme->Update(this->GetModelPart(), mpBuilderAndSolver->GetDofSet(), (*mpDx));
 
    KRATOS_CATCH("")
  }
 
 
  void SetSystemDofs()
  {
    KRATOS_TRY
 
   if (this->mEchoLevel >= 2)
      KRATOS_INFO(" Reform Dofs ") << " Flag = " <<this->mOptions.Is(LocalFlagType::REFORM_DOFS) << std::endl;
 
    //set up the system, operation performed just once unless it is required to reform the dof set at each iteration
 
    //setting up the list of the DOFs to be solved
    // double begin_time = OpenMPUtils::GetCurrentTime();
    mpBuilderAndSolver->SetUpDofSet(mpScheme, this->GetModelPart());
    // double end_time = OpenMPUtils::GetCurrentTime();
    // if (this->mEchoLevel >= 2)
    //   KRATOS_INFO("setup_dofs_time") << "setup_dofs_time : " << end_time - begin_time << "\n" << LoggerMessage::Category::STATISTICS;
 
    //shaping correctly the system
    // begin_time = OpenMPUtils::GetCurrentTime();
    mpBuilderAndSolver->SetUpSystem();
    // end_time = OpenMPUtils::GetCurrentTime();
    // if (this->mEchoLevel >= 2)
    //   KRATOS_INFO("setup_system_time") << ": setup_system_time : " << end_time - begin_time << "\n" << LoggerMessage::Category::STATISTICS;    //set up the system, operation performed just once unless it is required to reform the dof set at each iteration
 
    KRATOS_CATCH("")
  }
 
 
 
 private:
 
  LinearStrategy(const LinearStrategy &Other){};
 
 
}; 
 
 
 
 
}  // namespace Kratos.
#endif // KRATOS_LINEAR_STRATEGY_H_INCLUDED defined