residual_criteria.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
 
#pragma once
 
// System includes
 
// External includes
#include <limits>
 
// Project includes
#include "includes/model_part.h"
#include "utilities/constraint_utilities.h"
#include "utilities/parallel_utilities.h"
#include "utilities/reduction_utilities.h"
#include "solving_strategies/convergencecriterias/convergence_criteria.h"
 
namespace Kratos
{
 
 
 
 
 
template<class TSparseSpace,
         class TDenseSpace
         >
class ResidualCriteria
    : public  ConvergenceCriteria< TSparseSpace, TDenseSpace >
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( ResidualCriteria );
 
    using BaseType = ConvergenceCriteria< TSparseSpace, TDenseSpace >;
 
    using ClassType = ResidualCriteria< TSparseSpace, TDenseSpace >;
 
    using TDataType = typename BaseType::TDataType;
 
    using DofsArrayType = typename BaseType::DofsArrayType;
 
    using DofType = typename Node::DofType;
 
    using TSystemMatrixType = typename BaseType::TSystemMatrixType;
 
    using TSystemVectorType = typename BaseType::TSystemVectorType;
 
    using IndexType = std::size_t;
 
    using SizeType = std::size_t;
 
 
    explicit ResidualCriteria()
        : BaseType()
    {
    }
 
    explicit ResidualCriteria(Kratos::Parameters ThisParameters)
        : BaseType()
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
 
        this->mActualizeRHSIsNeeded = true;
    }
 
    explicit ResidualCriteria(
        TDataType NewRatioTolerance,
        TDataType AlwaysConvergedNorm)
        : BaseType(),
          mRatioTolerance(NewRatioTolerance),
          mAlwaysConvergedNorm(AlwaysConvergedNorm)
    {
        this->mActualizeRHSIsNeeded = true;
    }
 
    explicit ResidualCriteria( ResidualCriteria const& rOther )
      :BaseType(rOther)
      ,mRatioTolerance(rOther.mRatioTolerance)
      ,mInitialResidualNorm(rOther.mInitialResidualNorm)
      ,mCurrentResidualNorm(rOther.mCurrentResidualNorm)
      ,mAlwaysConvergedNorm(rOther.mAlwaysConvergedNorm)
      ,mReferenceDispNorm(rOther.mReferenceDispNorm)
      ,mActiveDofs(rOther.mActiveDofs)
      ,mInitialDoFId(rOther.mInitialDoFId)
    {
        this->mActualizeRHSIsNeeded = true;
    }
 
    ~ResidualCriteria() override {}
 
 
    ResidualCriteria& operator=(ResidualCriteria const& rOther) = delete;
 
 
    typename BaseType::Pointer Create(Parameters ThisParameters) const override
    {
        return Kratos::make_shared<ClassType>(ThisParameters);
    }
 
    bool PostCriteria(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        const TSystemMatrixType& rA,
        const TSystemVectorType& rDx,
        const TSystemVectorType& rb
        ) override
    {
        if (TSparseSpace::Size(rb) != 0) { //if we are solving for something
            // Some values
            const int rank = rModelPart.GetCommunicator().GetDataCommunicator().Rank();
 
            // Calculate the residual norm
            SizeType size_residual;
            CalculateResidualNorm(rModelPart, mCurrentResidualNorm, size_residual, rDofSet, rb);
 
            TDataType ratio{};
            if(mInitialResidualNorm < std::numeric_limits<TDataType>::epsilon()) {
                ratio = 0.0;
            } else {
                ratio = mCurrentResidualNorm/mInitialResidualNorm;
            }
 
            const TDataType float_size_residual = static_cast<TDataType>(size_residual);
            const TDataType absolute_norm = (mCurrentResidualNorm/float_size_residual);
 
            KRATOS_INFO_IF("RESIDUAL CRITERION", this->GetEchoLevel() > 1 && rank == 0) << " :: [ Initial residual norm = " << mInitialResidualNorm << "; Current residual norm =  " << mCurrentResidualNorm << "]" << std::endl;
            KRATOS_INFO_IF("RESIDUAL CRITERION", this->GetEchoLevel() > 0 && rank == 0) << " :: [ Obtained ratio = " << ratio << "; Expected ratio = " << mRatioTolerance << "; Absolute norm = " << absolute_norm << "; Expected norm =  " << mAlwaysConvergedNorm << "]" << std::endl;
 
            rModelPart.GetProcessInfo()[CONVERGENCE_RATIO] = ratio;
            rModelPart.GetProcessInfo()[RESIDUAL_NORM] = absolute_norm;
 
            const bool has_achieved_convergence = ratio <= mRatioTolerance || absolute_norm < mAlwaysConvergedNorm;
            KRATOS_INFO_IF("RESIDUAL CRITERION", has_achieved_convergence && this->GetEchoLevel() > 0 && rank == 0) << "Convergence is achieved" << std::endl;
            return has_achieved_convergence;
        } else {
            return true;
        }
    }
 
    void InitializeSolutionStep(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        const TSystemMatrixType& rA,
        const TSystemVectorType& rDx,
        const TSystemVectorType& rb
        ) override
    {
        BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb);
 
        // Filling mActiveDofs when MPC exist
        if (rModelPart.NumberOfMasterSlaveConstraints() > 0) {
            ComputeActiveDofs(rModelPart, rDofSet);
        }
 
        SizeType size_residual;
        CalculateResidualNorm(rModelPart, mInitialResidualNorm, size_residual, rDofSet, rb);
    }
 
    void FinalizeSolutionStep(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        const TSystemMatrixType& rA,
        const TSystemVectorType& rDx,
        const TSystemVectorType& rb
        ) override
    {
        BaseType::FinalizeSolutionStep(rModelPart, rDofSet, rA, rDx, rb);
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"                        : "residual_criteria",
            "residual_absolute_tolerance" : 1.0e-4,
            "residual_relative_tolerance" : 1.0e-9
        })");
 
        // 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 "residual_criteria";
    }
 
 
 
 
    std::string Info() const override
    {
        return "ResidualCriteria";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
    TDataType mRatioTolerance{};                                     
 
    TDataType mInitialResidualNorm{};                                
 
    TDataType mCurrentResidualNorm{};                                
 
    TDataType mAlwaysConvergedNorm{};                                
 
    TDataType mReferenceDispNorm{};                                  
 
    std::vector<int> mActiveDofs;                                    
 
    IndexType mInitialDoFId = std::numeric_limits<IndexType>::max(); 
 
 
 
    virtual void ComputeActiveDofs(
        ModelPart& rModelPart,
        const ModelPart::DofsArrayType& rDofSet
        )
    {
        // Filling mActiveDofs when MPC exist
        ConstraintUtilities::ComputeActiveDofs(rModelPart, mActiveDofs, rDofSet);
    }
 
    bool IsActiveAndLocalDof(
        const DofType& rDof,
        const int Rank
        )
    {
        const IndexType dof_id = rDof.EquationId();
        if constexpr (!TSparseSpace::IsDistributedSpace()) {
            return mActiveDofs[dof_id] == 1;
        } else {
            KRATOS_DEBUG_ERROR_IF((dof_id - mInitialDoFId) >= mActiveDofs.size() && (rDof.GetSolutionStepValue(PARTITION_INDEX) == Rank)) << "DofId is greater than the size of the active Dofs vector. DofId: " << dof_id << "\tInitialDoFId: " << mInitialDoFId << "\tActiveDofs size: " << mActiveDofs.size() << std::endl;
            return (mActiveDofs[dof_id - mInitialDoFId] == 1 && (rDof.GetSolutionStepValue(PARTITION_INDEX) == Rank));
        }
    }
 
    bool IsFreeAndLocalDof(
        const DofType& rDof,
        const int Rank
        )
    {
        if constexpr (!TSparseSpace::IsDistributedSpace()) {
            return rDof.IsFree();
        } else {
            return (rDof.IsFree() && (rDof.GetSolutionStepValue(PARTITION_INDEX) == Rank));
        }
    }
 
    virtual void CalculateResidualNorm(
        ModelPart& rModelPart,
        TDataType& rResidualSolutionNorm,
        SizeType& rDofNum,
        DofsArrayType& rDofSet,
        const TSystemVectorType& rb
        )
    {
        // Retrieve the data communicator
        const auto& r_data_communicator = rModelPart.GetCommunicator().GetDataCommunicator();
 
        // The current MPI rank
        const int rank = r_data_communicator.Rank();
 
        // Initialize
        TDataType residual_solution_norm = TDataType();
        unsigned int dof_num = 0;
 
        // Custom reduction
        using CustomReduction = CombinedReduction<SumReduction<TDataType>,SumReduction<unsigned int>>;
 
        // Auxiliary struct
        struct TLS {TDataType residual_dof_value{};};
 
        // Loop over Dofs
        if (rModelPart.NumberOfMasterSlaveConstraints() > 0) {
            std::tie(residual_solution_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [this, &rb, &rank](auto& rDof, TLS& rTLS) {
                if (this->IsActiveAndLocalDof(rDof, rank)) {
                    rTLS.residual_dof_value = TSparseSpace::GetValue(rb, rDof.EquationId());
                    return std::make_tuple(std::pow(rTLS.residual_dof_value, 2), 1);
                } else {
                    return std::make_tuple(TDataType(), 0);
                }
            });
        } else {
            std::tie(residual_solution_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [this, &rb, &rank](auto& rDof, TLS& rTLS) {
                if (this->IsFreeAndLocalDof(rDof, rank)) {
                    rTLS.residual_dof_value = TSparseSpace::GetValue(rb, rDof.EquationId());
                    return std::make_tuple(std::pow(rTLS.residual_dof_value, 2), 1);
                } else {
                    return std::make_tuple(TDataType(), 0);
                }
            });
        }
 
        rDofNum = static_cast<SizeType>(r_data_communicator.SumAll(dof_num));
        rResidualSolutionNorm = std::sqrt(r_data_communicator.SumAll(residual_solution_norm));
    }
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
        mAlwaysConvergedNorm = ThisParameters["residual_absolute_tolerance"].GetDouble();
        mRatioTolerance = ThisParameters["residual_relative_tolerance"].GetDouble();
    }
 
 
}; // Class ClassName
 
 
 
 
}  // namespace Kratos.