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
#include "../real_gas.h" // for RealGasBase
#include "../navier_stokes.h" // for InitialConditionFunction_TurbulentChannelFlow

namespace PHiLiP {

template <int dim, int nspecies, 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 nspecies, int nstate, typename real>
class FreeStreamInitialConditions : public InitialConditionFunction<dim,nspecies,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

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


    explicit FreeStreamInitialConditions (const Physics::Euler<dim,nspecies,nstate,double> euler_physics)
            : InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_TurbulentChannelFlow : public InitialConditionFunction<dim,nspecies,nstate,real>
{
protected:
    using dealii::Function<dim,real>::value; 

public:
    InitialConditionFunction_TurbulentChannelFlow (
        const Physics::NavierStokes<dim,nspecies,nstate,double> navier_stokes_physics_,
        const double channel_friction_velocity_reynolds_number_,
        const double domain_length_x_,
        const double domain_length_y_,
        const double domain_length_z_);

    const Physics::NavierStokes<dim,nspecies,nstate,double> navier_stokes_physics; 
    const double channel_friction_velocity_reynolds_number; 
    const double domain_length_x; 
    const double domain_length_y; 
    const double domain_length_z; 
    const double channel_height; 
    const double half_channel_height; 

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

protected:
    real get_distance_from_wall(const dealii::Point<dim,real> &point) const;
    virtual real x_velocity (const dealii::Point<dim,real> &point, const real density, const real temperature) const;
    virtual real y_velocity (const dealii::Point<dim,real> &point) const;
};

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_TurbulentChannelFlow_Turbulent : public InitialConditionFunction_TurbulentChannelFlow<dim,nspecies,nstate,real>
{
public:

    InitialConditionFunction_TurbulentChannelFlow_Turbulent (
        const Physics::NavierStokes<dim,nspecies,nstate,double> navier_stokes_physics_,
        const double channel_friction_velocity_reynolds_number_,
        const double domain_length_x_,
        const double domain_length_y_,
        const double domain_length_z_);

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

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_TurbulentChannelFlow_Manufactured : public InitialConditionFunction_TurbulentChannelFlow_Turbulent<dim,nspecies,nstate,real>
{
public:

    InitialConditionFunction_TurbulentChannelFlow_Manufactured (
        const Physics::NavierStokes<dim,nspecies,nstate,double> navier_stokes_physics_,
        const double channel_friction_velocity_reynolds_number_,
        const double domain_length_x_,
        const double domain_length_y_,
        const double domain_length_z_);

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

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

public:
    InitialConditionFunction_NavierStokesBase (
            Parameters::AllParameters const *const param);

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

    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;

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

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

public:
    explicit InitialConditionFunction_RealGasBase(
        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;

    // Real Gas physics pointer. Used to convert primitive to conservative.
    std::shared_ptr < Physics::RealGas<dim, nspecies, nstate, double > > real_gas_physics;
};

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

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

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

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

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_TaylorGreenVortex_Isothermal
    : public InitialConditionFunction_TaylorGreenVortex<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_DipoleWallCollision : public InitialConditionFunction_NavierStokesBase<dim,nspecies,nstate,real>
{
public:

    InitialConditionFunction_DipoleWallCollision (
            Parameters::AllParameters const *const param,
            const real extremum_vorticity_value_,
            const real dipole_radius,
            const real dipole_axis_angle_wrt_x_axis_in_degrees);

    const real extremum_vorticity_value; // Extremum vorticity value
    const real r0; // dipole radius
    const real x1; // x-coordinate of dipole 1
    const real y1; // y-coordinate of dipole 1
    const real x2; // x-coordinate of dipole 2
    const real y2; // y-coordinate of dipole 2

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

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_DipoleWallCollision_Normal : public InitialConditionFunction_DipoleWallCollision<dim,nspecies,nstate,real>
{
public:

    InitialConditionFunction_DipoleWallCollision_Normal (
            Parameters::AllParameters const *const param);
};

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_DipoleWallCollision_Oblique : public InitialConditionFunction_DipoleWallCollision<dim,nspecies,nstate,real>
{
public:

    InitialConditionFunction_DipoleWallCollision_Oblique (
            Parameters::AllParameters const *const param);
};

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_BurgersRewienski: public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_BurgersViscous: public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_BurgersInviscid: public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_BurgersInviscidEnergy
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_Advection
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_AdvectionEnergy
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_ConvDiff
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_ConvDiffEnergy
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_1DSine
        : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_IsentropicVortex
        : public InitialConditionFunction<dim,nspecies,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, nspecies, nstate, double > > euler_physics;

};




template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_KHI : public InitialConditionFunction<dim,nspecies,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, nspecies, nstate, double > > euler_physics;

};


template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_SodShockTube: public InitialConditionFunction_NavierStokesBase<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFunction_LeblancShockTube : public InitialConditionFunction_NavierStokesBase<dim, nspecies, 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 nspecies, int nstate, typename real>
class InitialConditionFunction_LowDensity: public InitialConditionFunction_NavierStokesBase<dim,nspecies,nstate,real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

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


template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_ShuOsherProblem : public InitialConditionFunction_NavierStokesBase<dim, nspecies, 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 nspecies, int nstate, typename real>
class InitialConditionFunction_DoubleMachReflection : public InitialConditionFunction_NavierStokesBase<dim, nspecies, nstate, real>
{
protected:
    real primitive_value(const dealii::Point<dim, real>& point, const unsigned int istate = 0) const override;

public:

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


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

public:

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



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

public:

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



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

public:

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


template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_Multispecies_VortexAdvection: public InitialConditionFunction_RealGasBase<dim,nspecies,nstate,real>
{
public:
    InitialConditionFunction_Multispecies_VortexAdvection (
            Parameters::AllParameters const *const param,
            bool high_temperature);
protected:
    real primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
    bool use_high_temp_ic;
};


template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_Multispecies_SodShockTube: public InitialConditionFunction_RealGasBase<dim,nspecies,nstate,real>
{
public:
    InitialConditionFunction_Multispecies_SodShockTube (
            Parameters::AllParameters const *const param);
protected:
    real primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
};


template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_Multispecies_IsentropicVortex: public InitialConditionFunction_RealGasBase<dim,nspecies,nstate,real>
{
public:
    InitialConditionFunction_Multispecies_IsentropicVortex (
            Parameters::AllParameters const *const param);
protected:
    real primitive_value(const dealii::Point<dim,real> &point, const unsigned int istate = 0) const;
};

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

public:
    explicit InitialConditionFunction_Multispecies_TaylorGreenVortex (
            Parameters::AllParameters const *const param, const bool use_smooth_interface);

    const double gamma_gas; 
    const double mach_inf; 
    const double mach_inf_sqr; 
    const bool smooth_interface; 

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

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

template <int dim, int nspecies, int nstate, typename real>
class InitialConditionFunction_Zero : public InitialConditionFunction<dim,nspecies,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 nspecies, int nstate, typename real>
class InitialConditionFactory
{
protected:    
    using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType;
    using DensityInitialConditionEnum = Parameters::FlowSolverParam::DensityInitialConditionType;
    using XVelocityInitialConditionEnum = Parameters::FlowSolverParam::XVelocityInitialConditionType;

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

} // PHiLiP namespace
#endif