estimate_dt_utilities.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Jordi Cotela
//                   Ruben Zorrilla
//
//
 
#ifndef KRATOS_ESTIMATE_DT_UTILITIES_H
#define KRATOS_ESTIMATE_DT_UTILITIES_H
 
// System includes
 
 
// External includes
 
 
// Project includes
#include "includes/define.h"
#include "includes/model_part.h"
#include "includes/kratos_parameters.h"
#include "includes/serializer.h"
 
// Application includes
 
 
namespace Kratos
{
 
 
 
class KRATOS_API(FLUID_DYNAMICS_APPLICATION) EstimateDtUtility
{
public:
 
 
    KRATOS_CLASS_POINTER_DEFINITION(EstimateDtUtility);
 
    typedef std::function<double(const Geometry<Node>&)> ElementSizeFunctionType;
 
    // Function type for the CFL calculator function
    typedef std::function<double(const Element &, const ElementSizeFunctionType &, const double)> CFLCalculatorType;
 
 
 
    EstimateDtUtility(
        ModelPart &ModelPart,
        const double CFL,
        const double DtMin,
        const double DtMax)
        : mrModelPart(ModelPart)
    {
        mCFL = CFL;
        mDtMin = DtMin;
        mDtMax = DtMax;
        mViscousFourier = 0.0;
        mThermalFourier = 0.0;
        mConsiderArtificialDiffusion = false;
        mConsiderCompressibilityInCFL = false;
 
        SetDtEstimationMagnitudesFlag();
    }
 
 
    EstimateDtUtility(
        ModelPart &ModelPart,
        const double CFL,
        const double ViscousFourier,
        const double ThermalFourier,
        const bool ConsiderArtificialDiffusion,
        const bool NodalDensityFormulation,
        const double DtMin,
        const double DtMax,
        const bool ConsiderCompressibilityInCFL = false)
        : mrModelPart(ModelPart)
    {
        mCFL = CFL;
        mDtMin = DtMin;
        mDtMax = DtMax;
        mViscousFourier = ViscousFourier;
        mThermalFourier = ThermalFourier;
        mConsiderArtificialDiffusion = ConsiderArtificialDiffusion;
        mNodalDensityFormulation = NodalDensityFormulation;
        mConsiderCompressibilityInCFL = ConsiderCompressibilityInCFL;
 
        SetDtEstimationMagnitudesFlag();
    }
 
 
    EstimateDtUtility(
        ModelPart& ModelPart,
        Parameters& rParameters)
        : mrModelPart(ModelPart)
    {
        Parameters defaultParameters(R"({
            "automatic_time_step"             : true,
            "CFL_number"                      : 1.0,
            "Viscous_Fourier_number"          : 0.0,
            "Thermal_Fourier_number"          : 0.0,
            "consider_artificial_diffusion"   : false,
            "nodal_density_formulation"       : false,
            "minimum_delta_time"              : 1e-4,
            "maximum_delta_time"              : 0.1,
            "consider_compressibility_in_CFL" : false
        })");
 
        rParameters.ValidateAndAssignDefaults(defaultParameters);
 
        mCFL = rParameters["CFL_number"].GetDouble();
        mViscousFourier = rParameters["Viscous_Fourier_number"].GetDouble();
        mThermalFourier = rParameters["Thermal_Fourier_number"].GetDouble();
        mConsiderArtificialDiffusion = rParameters["consider_artificial_diffusion"].GetBool();
        mNodalDensityFormulation = rParameters["nodal_density_formulation"].GetBool();
        mDtMin = rParameters["minimum_delta_time"].GetDouble();
        mDtMax = rParameters["maximum_delta_time"].GetDouble();
        mConsiderCompressibilityInCFL = rParameters["consider_compressibility_in_CFL"].GetBool();
 
        SetDtEstimationMagnitudesFlag();
    }
 
    ~EstimateDtUtility()
    {}
 
 
    void SetCFL(const double CFL);
 
    void SetViscousFourier(const double ViscousFourier);
 
    void SetThermalFourier(const double ThermalFourier);
 
    void SetDtMin(const double DtMin);
    
    void SetDtMax(const double DtMax);
 
    double EstimateDt() const;
 
 
private:
 
 
    KRATOS_DEFINE_LOCAL_FLAG(CFL_ESTIMATION);
    KRATOS_DEFINE_LOCAL_FLAG(VISCOUS_FOURIER_ESTIMATION);
    KRATOS_DEFINE_LOCAL_FLAG(THERMAL_FOURIER_ESTIMATION);
 
 
    double    mCFL;                         // User-defined CFL number
    double    mViscousFourier;              // User-defined viscous Fourier number 
    double    mThermalFourier;              // User-defined thermal Fourier number
    bool      mConsiderArtificialDiffusion; // Speficies if the artificial diffusion values are considered in the Peclet numbers
    bool      mNodalDensityFormulation;     // Specifies if the density is nodally stored (only required for the Peclet number)
    double    mDtMax;                       // User-defined maximum time increment allowed
    double    mDtMin;                       // User-defined minimum time increment allowed
    bool      mConsiderCompressibilityInCFL;// User-defined formulation. CFL number depends on this parameter.
    Flags     mDtEstimationMagnitudesFlags; // Flags indicating the reference magnitudes used in the Dt estimation
    ModelPart &mrModelPart;                 // The problem's model part
 
 
 
 
    void SetDtEstimationMagnitudesFlag();
 
    template<const bool ConsiderCFL, const bool ConsiderViscousFourier, const bool ConsiderThermalFourier>
    double InternalEstimateDt() const;
 
    template <class... CharacteristicNumbersPairsType>
    double CalculateNewDeltaTime(
        const double CurrentDeltaTime,
        const CharacteristicNumbersPairsType&... rCharacteristicNumbersPairs) const
    {
        KRATOS_TRY
 
        // Calculate the corresponding new time increments from the provided pairs
        const double zero_tol = std::min(0.1 * mDtMin, 1e-10);
        double new_dt_list[sizeof...(CharacteristicNumbersPairsType)] = {(
            (std::get<0>(rCharacteristicNumbersPairs) > zero_tol) ? CurrentDeltaTime * std::get<1>(rCharacteristicNumbersPairs) / std::get<0>(rCharacteristicNumbersPairs) : mDtMin
        )...};
 
        // Get the minimum one among all the obtained ones and check the user-defined bounds
        double new_dt = *(std::min_element(new_dt_list, new_dt_list + sizeof...(CharacteristicNumbersPairsType)));
        LimitNewDeltaTime(new_dt);
 
        // Perform MPI sync if needed
        new_dt = mrModelPart.GetCommunicator().GetDataCommunicator().MinAll(new_dt);
 
        return new_dt;
 
        KRATOS_CATCH("");
    }
 
    void LimitNewDeltaTime(double& rNewDeltaTime) const;
 
    
    const CFLCalculatorType GetCFLCalculatorUtility() const;
 
};
 
 
 
} // namespace Kratos.
 
 
#endif  /* KRATOS_ESTIMATE_DT_UTILITIES_H */