residual_criterion.hpp

Go to the documentation of this file.

//
//   Project Name:        KratosSolidMechanicsApplication $
//   Created by:          $Author:            JMCarbonell $
//   Last modified by:    $Co-Author:                     $
//   Date:                $Date:                 May 2018 $
//   Revision:            $Revision:                  0.0 $
//
//
 
#if !defined(KRATOS_RESIDUAL_CRITERION_H_INCLUDED )
#define  KRATOS_RESIDUAL_CRITERION_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "custom_solvers/convergence_criteria/convergence_criterion.hpp"
 
namespace Kratos
{
 
 
 
 
 
 
 
template<class TSparseSpace, class TDenseSpace>
class ResidualCriterion : public  ConvergenceCriterion< TSparseSpace, TDenseSpace >
{
 public:
 
 
  typedef ConvergenceCriterion< TSparseSpace, TDenseSpace >     BaseType;
  typedef typename BaseType::LocalFlagType                 LocalFlagType;
  typedef typename BaseType::DataType                           DataType;
  typedef typename BaseType::DofsArrayType                 DofsArrayType;
  typedef typename BaseType::SystemMatrixType           SystemMatrixType;
  typedef typename BaseType::SystemVectorType           SystemVectorType;
 
  typedef array_1d<double,3>                                  VectorType;
  typedef Variable<double>                           VectorComponentType;
  typedef Variable<VectorType>                        VariableVectorType;
  typedef Variable<double>                            VariableScalarType;
  typedef const VariableVectorType*                VariableVectorPointer;
  typedef const VariableScalarType*                VariableScalarPointer;
 
  typedef std::vector<Variable<double> >         ComponentVariableVector;
 
  KRATOS_CLASS_POINTER_DEFINITION(ResidualCriterion);
 
 
  ResidualCriterion(DataType RatioTolerance,
                    DataType AbsoluteTolerance)
      : BaseType(), mRatioTolerance(RatioTolerance), mAbsoluteTolerance(AbsoluteTolerance)
  {
    mpScalarVariable = nullptr;
    mpVectorVariable = nullptr;
    this->Set(LocalFlagType::SUPPLIED_DOF, false);
  }
 
  ResidualCriterion(const VariableScalarType& rScalarVariable,
                    DataType RatioTolerance,
                    DataType AbsoluteTolerance)
      : BaseType(), mRatioTolerance(RatioTolerance), mAbsoluteTolerance(AbsoluteTolerance)
  {
    mpScalarVariable = &rScalarVariable;
    mpVectorVariable = nullptr;
    this->Set(LocalFlagType::SUPPLIED_DOF, true);
  }
 
  ResidualCriterion(const VariableVectorType& rVectorVariable,
                    DataType RatioTolerance,
                    DataType AbsoluteTolerance)
      :BaseType(), mRatioTolerance(RatioTolerance), mAbsoluteTolerance(AbsoluteTolerance)
  {
    mpVectorVariable = &rVectorVariable;
    mpScalarVariable = nullptr;
    this->Set(LocalFlagType::SUPPLIED_DOF, true);
  }
 
  ResidualCriterion( ResidualCriterion const& rOther )
      :BaseType(rOther)
      ,mRatioTolerance(rOther.mRatioTolerance)
      ,mAbsoluteTolerance(rOther.mAbsoluteTolerance)
      ,mInitialResidualNorm(rOther.mInitialResidualNorm)
      ,mpScalarVariable(rOther.mpScalarVariable)
      ,mpVectorVariable(rOther.mpVectorVariable)
  {
  }
 
  ~ResidualCriterion() override {}
 
 
 
  //Criterion that needs to be called after getting the solution
  bool PostCriteria(ModelPart& rModelPart,
                    DofsArrayType& rDofSet,
                    const SystemMatrixType& rA,
                    const SystemVectorType& rDx,
                    const SystemVectorType& rb
                    ) override
  {
    if (TSparseSpace::Size(rb) != 0) //if we are solving for something
    {
 
      std::size_t size = 0;
 
      if( this->IsNot(LocalFlagType::INITIALIZED) )
      {
        size = CalculateResidualNorm(rDofSet,rb,mInitialResidualNorm);
        this->Set(LocalFlagType::INITIALIZED,true);
      }
 
      DataType ratio;
      DataType CurrentResidualNorm;
      size = CalculateResidualNorm(rDofSet,rb,CurrentResidualNorm);
 
      if( size == 0 ){
        KRATOS_WARNING("") << GetDofName() << " Residual Dofs vector has size: " << size << std::endl;
        return true;
      }
 
      if(mInitialResidualNorm == 0.00)
      {
        ratio = 0.00;
      }
      else
      {
        ratio = CurrentResidualNorm/mInitialResidualNorm;
      }
 
      const DataType absolute_norm = (CurrentResidualNorm/static_cast<DataType>(size));
 
      if (rModelPart.GetCommunicator().MyPID() == 0)
      {
        if (this->GetEchoLevel() >= 1)
        {
          std::cout << "RESIDUAL (" << GetDofName() << ") ["<<rModelPart.GetProcessInfo()[NL_ITERATION_NUMBER]<<"] :: Ratio = "<< ratio  << "; Norm = " << absolute_norm <<std::endl;
        }
      }
 
      rModelPart.GetProcessInfo()[CONVERGENCE_RATIO] = ratio;
      rModelPart.GetProcessInfo()[RESIDUAL_NORM] = absolute_norm;
 
      if (ratio <= mRatioTolerance || absolute_norm < mAbsoluteTolerance)
      {
        if (rModelPart.GetCommunicator().MyPID() == 0)
        {
          if (this->GetEchoLevel() >= 1)
          {
            std::cout << "Convergence is achieved" << std::endl;
          }
        }
        return true;
      }
      else
      {
        return false;
      }
    }
    else
    {
      return true;
    }
  }
 
 
  void InitializeSolutionStep(ModelPart& rModelPart,
                              DofsArrayType& rDofSet,
                              const SystemMatrixType& rA,
                              const SystemVectorType& rDx,
                              const SystemVectorType& rb) override
  {
    this->Set(LocalFlagType::INITIALIZED,false);
  }
 
 
 
 
 
 
 
 protected:
 
 
 
 
 
 
 
 
 
 private:
 
 
  DataType mRatioTolerance;
 
  DataType mAbsoluteTolerance;
 
  DataType mInitialResidualNorm;
 
  VariableScalarPointer  mpScalarVariable;
 
  VariableVectorPointer  mpVectorVariable;
 
 
 
  std::size_t CalculateResidualNorm(DofsArrayType& rDofSet, const SystemVectorType& rb, DataType& rResidualNorm)
  {
    KRATOS_TRY
 
    std::size_t size = 0;
    rResidualNorm = 0.0;
 
    if( this->Is(CriterionLocalFlags::SUPPLIED_DOF) ){
 
      if( mpVectorVariable != nullptr ){
 
        DataType temp;
        for(typename DofsArrayType::iterator i_dof = rDofSet.begin() ; i_dof != rDofSet.end() ; ++i_dof)
        {
          if(i_dof->IsFree())
          {
            if(CheckVectorDof(i_dof))
            {
              temp = rb[i_dof->EquationId()];
              rResidualNorm +=  temp*temp;
              ++size;
            }
          }
        }
      }
      else if( mpScalarVariable != nullptr ){
 
        DataType temp;
 
        for(typename DofsArrayType::iterator i_dof = rDofSet.begin() ; i_dof != rDofSet.end() ; ++i_dof)
        {
          if(i_dof->IsFree())
          {
            if(i_dof->GetVariable() == *mpScalarVariable)
            {
              temp = rb[i_dof->EquationId()];
              rResidualNorm += temp*temp;
              ++size;
            }
          }
        }
 
      }
      else{
        KRATOS_ERROR << " No variable dof supplied to the convergence criterion " << std::endl;
      }
 
      rResidualNorm = std::sqrt(rResidualNorm);
    }
    else{
      rResidualNorm = TSparseSpace::TwoNorm(rb);
      size = TSparseSpace::Size(rb);
    }
 
    return size;
 
    KRATOS_CATCH("")
  }
 
  bool CheckVectorDof(typename DofsArrayType::iterator& rDofIter)
  {
    KRATOS_TRY
 
    const std::string& variable_name = mpVectorVariable->Name();
    const Variable<double>& var_x = KratosComponents<Variable<double> >::Get(variable_name+"_X");
    const Variable<double>& var_y = KratosComponents<Variable<double> >::Get(variable_name+"_Y");
    const Variable<double>& var_z = KratosComponents<Variable<double> >::Get(variable_name+"_Z");
 
    if( rDofIter->GetVariable() == var_x || rDofIter->GetVariable() == var_y || rDofIter->GetVariable() == var_z )
      return true;
    else
      return false;
 
    KRATOS_CATCH("")
  }
 
  std::string GetDofName()
  {
    std::string name = "DOFS";
 
    if( this->Is(CriterionLocalFlags::SUPPLIED_DOF) ){
 
      if( mpVectorVariable != nullptr ){
        name =  mpVectorVariable->Name();
      }
      else if( mpScalarVariable != nullptr ){
        name = mpScalarVariable->Name();
      }
    }
 
    return name;
  }
 
 
 
 
 
}; // Class ResidualCriterion
 
 
 
 
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_RESIDUAL_CRITERION_H_INCLUDED defined