newton_raphson_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_NEWTON_RAPHSON_STRATEGY_H_INCLUDED)
#define KRATOS_NEWTON_RAPHSON_STRATEGY_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_solvers/solution_strategies/linear_strategy.hpp"
 
//default convergence criterion
#include "custom_solvers/convergence_criteria/convergence_criterion.hpp"
 
#include "solid_mechanics_application_variables.h"
 
namespace Kratos
{
 
 
 
 
 
 
template <class TSparseSpace,
          class TDenseSpace,  // = DenseSpace<double>,
          class TLinearSolver // = LinearSolver<TSparseSpace,TDenseSpace>
          >
class NewtonRaphsonStrategy : public LinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>
{
 public:
 
  // Counted pointer of ClassName
  KRATOS_CLASS_POINTER_DEFINITION(NewtonRaphsonStrategy);
 
  typedef LinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>              BaseType;
 
  typedef typename BaseType::LocalFlagType                                 LocalFlagType;
 
  typedef ConvergenceCriterion<TSparseSpace, TDenseSpace>       ConvergenceCriterionType;
 
  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;
 
 
 
  NewtonRaphsonStrategy(ModelPart& rModelPart,
                        typename SchemeType::Pointer pScheme,
                        typename BuilderAndSolverType::Pointer pBuilderAndSolver,
                        typename ConvergenceCriterionType::Pointer pConvergenceCriterion,
                        Flags& rOptions,
                        unsigned int MaxIterations = 30)
      : LinearStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, pBuilderAndSolver, rOptions)
  {
    KRATOS_TRY
 
    // Saving the convergence criteria to be used
    mpConvergenceCriteria = pConvergenceCriterion;
 
    // Maximum iterations allowed
    mMaxIterationNumber = MaxIterations;
 
    KRATOS_CATCH("")
  }
 
 
  NewtonRaphsonStrategy(ModelPart& rModelPart,
                        typename SchemeType::Pointer pScheme,
                        typename LinearSolverType::Pointer pLinearSolver,
                        typename ConvergenceCriterionType::Pointer pConvergenceCriterion,
                        Flags& rOptions,
                        unsigned int MaxIterations = 30)
      : NewtonRaphsonStrategy<TSparseSpace, TDenseSpace, TLinearSolver>(rModelPart, pScheme, Kratos::make_shared<BlockBuilderAndSolver<TSparseSpace, TDenseSpace, TLinearSolver> >(pLinearSolver), pConvergenceCriterion, rOptions, MaxIterations)
  {}
 
 
  ~NewtonRaphsonStrategy() override
  {
    Clear();
  }
 
 
 
 
  void InitializeSolutionStep() override
  {
    KRATOS_TRY
 
    //set implex
    if(this->mOptions.Is(LocalFlagType::IMPLEX))
      this->GetModelPart().GetProcessInfo().SetValue(IMPLEX, true);
 
    BaseType::InitializeSolutionStep();
 
    KRATOS_CATCH("")
  }
 
  void FinalizeSolutionStep() override
  {
    KRATOS_TRY
 
    //Finalization of the solution step, operations to be done after achieving convergence
 
    //set implex calculation
    if(this->mOptions.Is(LocalFlagType::IMPLEX)){
      this->GetModelPart().GetProcessInfo().SetValue(IMPLEX, false);
      if(this->mOptions.IsNot(LocalFlagType::COMPUTE_REACTIONS))
        this->mpBuilderAndSolver->BuildRHS(this->mpScheme, this->GetModelPart(), (*this->mpb));
    }
 
    BaseType::FinalizeSolutionStep();
 
    KRATOS_CATCH("")
  }
 
 
  bool SolveSolutionStep() override
  {
    KRATOS_TRY
 
    //initializing the parameters of the Newton-Raphson cicle
    unsigned int iteration_number = 1;
 
    //setting the iteration number
    this->GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] = iteration_number;
 
    this->Set(LocalFlagType::CONVERGED, this->SolveIteration());
 
    //iteration cycle... performed only for NonLinearProblems
    while( this->IsNot(LocalFlagType::CONVERGED) && iteration_number++ < mMaxIterationNumber)
    {
      //setting the iteration number
      this->GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] = iteration_number;
 
      this->Set(LocalFlagType::CONVERGED, this->SolveIteration());
    }
 
    //plots a warning if the maximum number of iterations is exceeded
    if(iteration_number >= mMaxIterationNumber)
    {
      if( this->GetEchoLevel() >= 0 )
        KRATOS_INFO("  [Iterative loop interrupted] ") << "[" << iteration_number << " iterations performed] \n";
    }
 
    return (this->Is(LocalFlagType::CONVERGED));
 
    KRATOS_CATCH("")
 
  }
 
 
  bool SolveIteration() override
  {
    KRATOS_TRY
 
    bool is_converged = false;
 
    // Warning info
    if(!SparseSpaceType::Size((*this->mpDx)))
      KRATOS_WARNING("DOFS") << "solution has zero size, no free DOFs" << std::endl;
 
    // Initialize Iteration
    this->mpScheme->InitializeNonLinearIteration(this->GetModelPart());
 
    is_converged = mpConvergenceCriteria->PreCriteria(this->GetModelPart(), this->mpBuilderAndSolver->GetDofSet(), (*this->mpA), (*this->mpDx), (*this->mpb));
 
    // Function to perform the building and the solving phase.
    if(this->mOptions.IsNot(LocalFlagType::CONSTANT_SYSTEM_MATRIX)){
 
      TSparseSpace::SetToZero((*this->mpA));
      TSparseSpace::SetToZero((*this->mpDx));
      TSparseSpace::SetToZero((*this->mpb));
 
      this->mpBuilderAndSolver->BuildAndSolve(this->mpScheme, this->GetModelPart(), (*this->mpA), (*this->mpDx), (*this->mpb));
    }
    else{
 
      TSparseSpace::SetToZero((*this->mpDx));
      TSparseSpace::SetToZero((*this->mpb));
 
      this->mpBuilderAndSolver->BuildRHSAndSolve(this->mpScheme, this->GetModelPart(), (*this->mpA), (*this->mpDx), (*this->mpb));
    }
 
    // EchoInfo
    this->mpBuilderAndSolver->EchoInfo(this->GetModelPart(), (*this->mpA), (*this->mpDx), (*this->mpb));
 
    // Updating the results
    this->Update();
 
    // Finalize Iteration
    this->mpScheme->FinalizeNonLinearIteration(this->GetModelPart());
 
    if(is_converged == true)
    {
      //initialisation of the convergence criteria (after first calculation only)
      if( this->GetModelPart().GetProcessInfo()[NL_ITERATION_NUMBER] == 1 ){
        mpConvergenceCriteria->InitializeSolutionStep(this->GetModelPart(), this->mpBuilderAndSolver->GetDofSet(), (*this->mpA), (*this->mpDx), (*this->mpb));
      }
 
      if(mpConvergenceCriteria->Is(CriterionLocalFlags::UPDATE_RHS))
      {
        TSparseSpace::SetToZero((*this->mpb));
 
        this->mpBuilderAndSolver->BuildRHS(this->mpScheme, this->GetModelPart(), (*this->mpb));
      }
 
      is_converged = mpConvergenceCriteria->PostCriteria(this->GetModelPart(), this->mpBuilderAndSolver->GetDofSet(), (*this->mpA), (*this->mpDx), (*this->mpb));
    }
 
    return is_converged;
 
    KRATOS_CATCH("")
  }
 
 
  void Clear() override
  {
    KRATOS_TRY
 
    BaseType::Clear();
 
    KRATOS_CATCH("")
  }
 
  int Check() override
  {
    KRATOS_TRY
 
    //check linear strategy
    BaseType::Check();
 
    //check the convergence criterion
    mpConvergenceCriteria->Check(this->GetModelPart());
 
    return 0;
 
    KRATOS_CATCH("")
  }
 
 
 
  void SetEchoLevel(const int Level) override
  {
    BaseType::SetEchoLevel(Level);
  }
 
  void SetMaxIterationNumber(unsigned int MaxIterationNumber)
  {
    mMaxIterationNumber = MaxIterationNumber;
  }
 
  unsigned int GetMaxIterationNumber() override
  {
    return mMaxIterationNumber;
  }
 
 
 
 
 protected:
 
 
  typename ConvergenceCriterionType::Pointer mpConvergenceCriteria; 
 
  unsigned int mMaxIterationNumber; 
 
 
  void Initialize() override
  {
    KRATOS_TRY
 
    BaseType::Initialize();
 
    KRATOS_CATCH("")
  }
 
 
 private:
 
  NewtonRaphsonStrategy(const NewtonRaphsonStrategy &Other){};
 
 
}; 
 
 
 
 
}  // namespace Kratos.
#endif // KRATOS_NEWTON_RAPHSON_STRATEGY_H_INCLUDED defined