displacement_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
 
// Project includes
#include "includes/model_part.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 DisplacementCriteria
    : public ConvergenceCriteria< TSparseSpace, TDenseSpace >
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION( DisplacementCriteria );
 
    using BaseType = ConvergenceCriteria< TSparseSpace, TDenseSpace >;
 
    using ClassType = DisplacementCriteria< TSparseSpace, TDenseSpace >;
 
    using SparseSpaceType = TSparseSpace;
 
    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 DisplacementCriteria()
        : BaseType()
    {
    }
 
    explicit DisplacementCriteria(Kratos::Parameters ThisParameters)
        : BaseType()
    {
        // Validate and assign defaults
        ThisParameters = this->ValidateAndAssignParameters(ThisParameters, this->GetDefaultParameters());
        this->AssignSettings(ThisParameters);
    }
 
    explicit DisplacementCriteria(
        TDataType NewRatioTolerance,
        TDataType AlwaysConvergedNorm)
        : BaseType(),
          mRatioTolerance(NewRatioTolerance),
          mAlwaysConvergedNorm(AlwaysConvergedNorm)
    {
    }
 
    explicit DisplacementCriteria( DisplacementCriteria const& rOther )
      :BaseType(rOther)
      ,mRatioTolerance(rOther.mRatioTolerance)
      ,mAlwaysConvergedNorm(rOther.mAlwaysConvergedNorm)
      ,mReferenceDispNorm(rOther.mReferenceDispNorm)
    {
    }
 
    ~DisplacementCriteria() override {}
 
 
    DisplacementCriteria& operator=(DisplacementCriteria 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 (SparseSpaceType::Size(rDx) != 0) { //if we are solving for something
            // Some values
            const int rank = rModelPart.GetCommunicator().GetDataCommunicator().Rank();
            const TDataType approx_zero_tolerance = std::numeric_limits<TDataType>::epsilon();
 
            // The final correction norm calculation
            SizeType size_solution = 0;
            const TDataType final_correction_norm = CalculateFinalCorrectionNorm(size_solution, rDofSet, rDx, rModelPart);
 
            TDataType ratio = 0.0;
 
            CalculateReferenceNorm(rDofSet, rModelPart);
            if (mReferenceDispNorm < approx_zero_tolerance) {
                KRATOS_WARNING("DisplacementCriteria") << "NaN norm is detected. Setting reference to convergence criteria" << std::endl;
                mReferenceDispNorm = final_correction_norm;
            }
 
            if(final_correction_norm < approx_zero_tolerance) {
                ratio = 0.0;
            } else {
                ratio = final_correction_norm/mReferenceDispNorm;
            }
 
            const TDataType float_size_solution = static_cast<TDataType>(size_solution);
 
            const TDataType absolute_norm = (final_correction_norm/std::sqrt(float_size_solution));
 
            KRATOS_INFO_IF("DISPLACEMENT 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; //  || (final_correction_norm/x.size())<=1e-7)
            KRATOS_INFO_IF("DISPLACEMENT CRITERION", has_achieved_convergence && this->GetEchoLevel() > 0 && rank == 0) << "Convergence is achieved" << std::endl;
            return has_achieved_convergence;
        } else { //in this case all the displacements are imposed!
            return true;
        }
    }
 
    void Initialize(
        ModelPart& rModelPart
        ) override
    {
        BaseType::mConvergenceCriteriaIsInitialized = true;
    }
 
    void InitializeSolutionStep(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        const TSystemMatrixType& rA,
        const TSystemVectorType& rDx,
        const TSystemVectorType& rb
        ) override
    {
        BaseType::InitializeSolutionStep(rModelPart, rDofSet, rA, rDx, rb);
    }
 
    void FinalizeSolutionStep(
        ModelPart& rModelPart,
        DofsArrayType& rDofSet,
        const TSystemMatrixType& rA,
        const TSystemVectorType& rDx,
        const TSystemVectorType& rb
        ) override
    {
        BaseType::FinalizeSolutionStep(rModelPart, rDofSet, rA, rDx, rb);
    }
 
    static std::string Name()
    {
        return "displacement_criteria";
    }
 
    Parameters GetDefaultParameters() const override
    {
        Parameters default_parameters = Parameters(R"(
        {
            "name"                            : "displacement_criteria",
            "displacement_relative_tolerance" : 1.0e-4,
            "displacement_absolute_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;
    }
 
 
 
 
    std::string Info() const override
    {
        return "DisplacementCriteria";
    }
 
    void PrintInfo(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
    void PrintData(std::ostream& rOStream) const override
    {
        rOStream << Info();
    }
 
 
protected:
 
 
    TDataType mRatioTolerance = 0.0;      
 
    TDataType mAlwaysConvergedNorm = 0.0; 
 
    TDataType mReferenceDispNorm = 0.0;   
 
 
 
    void AssignSettings(const Parameters ThisParameters) override
    {
        BaseType::AssignSettings(ThisParameters);
        mAlwaysConvergedNorm = ThisParameters["displacement_absolute_tolerance"].GetDouble();
        mRatioTolerance = ThisParameters["displacement_relative_tolerance"].GetDouble();
    }
 
 
 
 
private:
 
 
 
 
    bool IsFreeAndLocalDof(
        const DofType& rDof,
        const int Rank
        )
    {
        if constexpr (!TSparseSpace::IsDistributedSpace()) {
            return rDof.IsFree();
        } else {
            return (rDof.IsFree() && (rDof.GetSolutionStepValue(PARTITION_INDEX) == Rank));
        }
    }
 
    void CalculateReferenceNorm(
        DofsArrayType& rDofSet,
        ModelPart& rModelPart
        )
    {
        // Retrieve the data communicator
        const auto& r_data_communicator = rModelPart.GetCommunicator().GetDataCommunicator();
 
        // The current MPI rank
        const int rank = r_data_communicator.Rank();
 
        // Auxiliary struct
        struct TLS {TDataType dof_value{};};
 
        const TDataType reference_disp_norm = block_for_each<SumReduction<TDataType>>(rDofSet, TLS(), [this, &rank](auto& rDof, TLS& rTLS) {
            if (ClassType::IsFreeAndLocalDof(rDof, rank)) {
                rTLS.dof_value = rDof.GetSolutionStepValue();
                return std::pow(rTLS.dof_value, 2);
            } else {
                return TDataType();
            }
        });
        mReferenceDispNorm = std::sqrt(r_data_communicator.SumAll(reference_disp_norm));
    }
 
    TDataType CalculateFinalCorrectionNorm(
        SizeType& rDofNum,
        DofsArrayType& rDofSet,
        const TSystemVectorType& rDx,
        ModelPart& rModelPart
        )
    {
        // 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 final_correction_norm = TDataType();
        unsigned int dof_num = 0;
 
        // Custom reduction
        using CustomReduction = CombinedReduction<SumReduction<TDataType>,SumReduction<unsigned int>>;
 
        // Auxiliary struct
        struct TLS {TDataType dof_value{}; TDataType variation_dof_value{};};
 
        // Loop over Dofs
        std::tie(final_correction_norm, dof_num) = block_for_each<CustomReduction>(rDofSet, TLS(), [this, &rDx, &rank](auto& rDof, TLS& rTLS) {
            if (this->IsFreeAndLocalDof(rDof, rank)) {
                rTLS.variation_dof_value = SparseSpaceType::GetValue(rDx, rDof.EquationId());
                return std::make_tuple(std::pow(rTLS.variation_dof_value, 2), 1);
            } else {
                return std::make_tuple(TDataType(), 0);
            }
        });
 
        rDofNum = static_cast<SizeType>(r_data_communicator.SumAll(dof_num));
        return std::sqrt(r_data_communicator.SumAll(final_correction_norm));
    }
 
}; /* Class DisplacementCriteria */
 
 
 
}  /* namespace Kratos.*/