initial_condition_function.h

#ifndef __INITIAL_CONDITION_FUNCTION_H__
#define __INITIAL_CONDITION_FUNCTION_H__

// for the initial condition function:
#include <deal.II/lac/vector.h>
#include <deal.II/base/function.h>
#include "parameters/all_parameters.h"
#include "../euler.h" // for FreeStreamInitialConditions

namespace PHiLiP {

template <int dim, int nstate, typename real>
class InitialConditionFunction : public dealii::Function<dim,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:
    InitialConditionFunction();

    virtual real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const = 0;

};

template <int dim, int nstate, typename real>
class FreeStreamInitialConditions : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:
    std::array<double,nstate> farfield_conservative;


    explicit FreeStreamInitialConditions (const Physics::Euler<dim,nstate,double> euler_physics)
            : InitialConditionFunction<dim,nstate,real>()
    {
        //const double density_bc = 2.33333*euler_physics.density_inf;
        const double density_bc = euler_physics.density_inf;
        const double pressure_bc = 1.0/(euler_physics.gam*euler_physics.mach_inf_sqr);
        std::array<double,nstate> primitive_boundary_values;
        primitive_boundary_values[0] = density_bc;
        for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = euler_physics.velocities_inf[d]; }
        primitive_boundary_values[nstate-1] = pressure_bc;
        farfield_conservative = euler_physics.convert_primitive_to_conservative(primitive_boundary_values);
    }

    double value (const dealii::Point<dim> &/*point*/, const unsigned int istate) const
    {
        return farfield_conservative[istate];
    }
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_EulerBase : public InitialConditionFunction<dim, nstate, real>
{
protected:
    using dealii::Function<dim, real>::value; 

public:
    explicit InitialConditionFunction_EulerBase(
        Parameters::AllParameters const* const param);

    real value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

protected:
    virtual real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const = 0;

    real convert_primitive_to_conversative_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const;

private:
    // Euler physics pointer. Used to convert primitive to conservative.
    std::shared_ptr < Physics::Euler<dim, nstate, double > > euler_physics;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_TaylorGreenVortex : public InitialConditionFunction_EulerBase<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    explicit InitialConditionFunction_TaylorGreenVortex (
            Parameters::AllParameters const *const param);

    const double gamma_gas; 
    const double mach_inf; 
    const double mach_inf_sqr; 

protected:
    real primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;

    virtual real density(const dealii::Point<dim,real> &point) const;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_TaylorGreenVortex_Isothermal
    : public InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
{
public:

    explicit InitialConditionFunction_TaylorGreenVortex_Isothermal (
            Parameters::AllParameters const *const param);

protected:
    real density(const dealii::Point<dim,real> &point) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_BurgersRewienski: public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_BurgersRewienski ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_BurgersViscous: public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_BurgersViscous ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_BurgersInviscid: public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_BurgersInviscid ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_BurgersInviscidEnergy
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_BurgersInviscidEnergy ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_Advection
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_Advection ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_AdvectionEnergy
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_AdvectionEnergy ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_ConvDiff
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_ConvDiff ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_ConvDiffEnergy
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    InitialConditionFunction_ConvDiffEnergy ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_1DSine
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:
    InitialConditionFunction_1DSine ();

    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_IsentropicVortex
        : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:

    explicit InitialConditionFunction_IsentropicVortex (
            Parameters::AllParameters const *const param);

    real value (const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

protected:

    // Euler physics pointer. Used to convert primitive to conservative.
    std::shared_ptr < Physics::Euler<dim, nstate, double > > euler_physics;

};




template <int dim, int nstate, typename real>
class InitialConditionFunction_KHI : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 
    
public:
    explicit InitialConditionFunction_KHI(
            Parameters::AllParameters const *const param);

    real value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;

protected:
    
    const real atwood_number;

    // Euler physics pointer. Used to convert primitive to conservative.
    std::shared_ptr < Physics::Euler<dim, nstate, double > > euler_physics;

};


template <int dim, int nstate, typename real>
class InitialConditionFunction_SodShockTube: public InitialConditionFunction_EulerBase<dim,nstate,real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

    explicit InitialConditionFunction_SodShockTube (
            Parameters::AllParameters const* const param);
};


template <int dim, int nstate, typename real>
class InitialConditionFunction_LowDensity2D: public InitialConditionFunction_EulerBase<dim,nstate,real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

    explicit InitialConditionFunction_LowDensity2D(
        Parameters::AllParameters const* const param);
};


template <int dim, int nstate, typename real>
class InitialConditionFunction_LeblancShockTube : public InitialConditionFunction_EulerBase<dim, nstate, real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

    explicit InitialConditionFunction_LeblancShockTube(
        Parameters::AllParameters const* const param);
};


template <int dim, int nstate, typename real>
class InitialConditionFunction_ShuOsherProblem : public InitialConditionFunction_EulerBase<dim, nstate, real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

    explicit InitialConditionFunction_ShuOsherProblem(
        Parameters::AllParameters const* const param);
};

template <int dim, int nstate, typename real>
class InitialConditionFunction_Zero : public InitialConditionFunction<dim,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 
    
public:
    InitialConditionFunction_Zero();

    real value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const override;
};

template <int dim, int nstate, typename real>
class InitialConditionFactory
{
protected:    
    using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType;
    using DensityInitialConditionEnum = Parameters::FlowSolverParam::DensityInitialConditionType;

public:
    static std::shared_ptr<InitialConditionFunction<dim,nstate,real>>
    create_InitialConditionFunction(
        Parameters::AllParameters const *const param);
};

} // PHiLiP namespace
#endif