shock_capturing_entropy_viscosity_process.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main author:     Eduard Gómez
//
 
#if !defined(KRATOS_SHOCK_CAPTURING_ENTROPY_VISCOSITY_H_INCLUDED)
#define  KRATOS_SHOCK_CAPTURING_ENTROPY_VISCOSITY_H_INCLUDED
 
// System includes
 
// External includes
 
// Project includes
#include "containers/model.h"
#include "processes/process.h"
 
// Application includes
 
 
namespace Kratos
{
 
 
 
 
 
 
class KRATOS_API(FLUID_DYNAMICS_APPLICATION) ShockCapturingEntropyViscosityProcess : public Process
{
public:
 
 
    typedef ModelPart::NodeType NodeType;
 
    struct InfNormData
    {
        double EntropyResidual;
        double Density;
        double TotalVelocity;
    };
 
 
    struct TotalDerivativeUtil
    {
        Vector Value;
        Vector TimeDerivative;
        Matrix Flux;
 
        TotalDerivativeUtil()
            : TotalDerivativeUtil(0, 0)
        {}
 
        TotalDerivativeUtil(
            const std::size_t NumberOfDimensions,
            const std::size_t NumberOfNodes)
            : Value          {NumberOfNodes, 0.0},
              TimeDerivative {NumberOfNodes, 0.0},
              Flux           {NumberOfNodes, NumberOfDimensions, 0.0}
        {
        }
 
        void LoadNodalValues(
            const Variable<double>& rVariable,
            const NodeType& rNode,
            const std::size_t NodeIndex,
            const array_1d<double, 3>& Velocity,
            const double DeltaTime)
        {
            KRATOS_TRY
 
            const auto old_value = rNode.FastGetSolutionStepValue(rVariable, 1);
            Value[NodeIndex] = rNode.FastGetSolutionStepValue(rVariable);
            TimeDerivative[NodeIndex] = (Value[NodeIndex] - old_value) / DeltaTime;
 
            Flux(NodeIndex, 0) = Value[NodeIndex] * Velocity[0];
            Flux(NodeIndex, 1) = Value[NodeIndex] * Velocity[1];
            if(Flux.size2()==3){ // (if 3D)
                Flux(NodeIndex, 2) = Value[NodeIndex] * Velocity[2];
            }
 
            KRATOS_CATCH("")
        }
 
        double ComputeAtGaussPoint(
            const MatrixRow<const Matrix>& rShapeFunctions,
            const Matrix& rShapeFunctionsGradents) const
        {
            KRATOS_TRY
 
            const double divergence = Divergence(rShapeFunctionsGradents, Flux);
            const double time_derivative = inner_prod(TimeDerivative, rShapeFunctions);
            return time_derivative + divergence;
 
            KRATOS_CATCH("") // Catching possible error in divergence computation
        }
 
    private:
        static double Divergence(const Matrix& rShapeFunGradients, const Matrix& rNodalValues);
    };
 
 
    KRATOS_CLASS_POINTER_DEFINITION(ShockCapturingEntropyViscosityProcess);
 
 
    ShockCapturingEntropyViscosityProcess(
        Model& rModel,
        Parameters rParameters)
        : ShockCapturingEntropyViscosityProcess(rModel.GetModelPart(rParameters["model_part_name"].GetString()), rParameters) {};
 
    ShockCapturingEntropyViscosityProcess(
        ModelPart& rModelPart,
        Parameters rParameters);
 
 
 
 
    void ExecuteInitializeSolutionStep() override;
 
    int Check() override;
 
    const Parameters GetDefaultParameters() const override;
 
 
 
 
 
 
    std::string Info() const override
    {
        std::stringstream buffer;
        buffer << "ShockCapturingEntropyViscosityProcess";
        return buffer.str();
    }
 
    void PrintInfo(std::ostream& rOStream) const override { rOStream << "ShockCapturingEntropyViscosityProcess"; }
 
    void PrintData(std::ostream&) const override {}
 
 
 
private:
 
 
 
    ModelPart& mrModelPart;
    bool mIsInitialized = false;
    bool mComputeAreasEveryStep = false;
    double mEntropyConstant = 0.0;
    double mEnergyConstant = 0.0;
    double mArtificialMassDiffusivityPrandtl = 0.0;
    double mArtificialConductivityPrandtl = 0.0;
 
 
 
 
    void ValidateAndAssignParameters(Parameters &rParameters);
 
    void UpdateNodalAreaProcess();
 
    void ComputeNodalEntropies(const unsigned int WriteBufferIndex = 0);
 
    void ComputeArtificialMagnitudes();
 
    void DistributeVariablesToNodes(
        Element& rElement,
        const double ArtificialDynamicViscosity,
        const double ArtificialMassDiffusivity,
        const double ArtificialConductivity,
        const std::function<double(Geometry<Node>*)>& rGeometrySize) const;
 
    static double ComputeEntropy(
        const double Density, 
        const double Pressure, 
        const double Gamma)
    {
        return Density / (Gamma - 1.0) * std::log(Pressure / std::pow(Density, Gamma));
    }
 
    template <unsigned int TDim>
    static double MinimumEdgeLengthSquared(const Element &rElement);
 
    
    static InfNormData ComputeElementalInfNormData(
        const Element& rElement,
        const double DeltaTime,
        const double HeatCapacityRatio,
        const double SpecificHeatCV);
 
    static std::tuple<TotalDerivativeUtil, TotalDerivativeUtil, Vector> BuildTotalDerivativeUtils(
        const Element& rElement,
        const double DeltaTime,
        const double HeatCapacityRatio,
        const double SpecificHeatCV);
 
    static InfNormData ComputeInfNorms(
        const Geometry<NodeType>& rGeometry,
        const TotalDerivativeUtil& rEntropyTotalDerivative,
        const TotalDerivativeUtil& rDensityTotalDerivative,
        const Vector& rTotalVelocities);
 
 
 
 
 
 
 
 
    ShockCapturingEntropyViscosityProcess& operator=(ShockCapturingEntropyViscosityProcess const& rOther);
 
    ShockCapturingEntropyViscosityProcess(ShockCapturingEntropyViscosityProcess const& rOther);
 
 
}; // Class ShockCapturingEntropyViscosityProcess
 
 
 
 
inline std::ostream& operator << (
    std::ostream& rOStream,
    const ShockCapturingEntropyViscosityProcess& rThis)
{
    rThis.PrintInfo(rOStream);
    return rOStream;
}
 
 
 
}  // namespace Kratos.
 
#endif // KRATOS_SHOCK_CAPTURING_ENTROPY_VISCOSITY_H_INCLUDED  defined