euler.h

#ifndef __EULER__
#define __EULER__

#include <deal.II/base/tensor.h>
#include "physics.h"
#include "parameters/all_parameters.h"
#include "parameters/parameters_manufactured_solution.h"

namespace PHiLiP {
namespace Physics {


template <int dim, int nspecies, int nstate, typename real>
class Euler : public PhysicsBase <dim, nspecies, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using PhysicsBase<dim,nspecies,nstate,real>::dissipative_flux;
    using PhysicsBase<dim,nspecies,nstate,real>::source_term;
    using PhysicsBase<dim,nspecies,nstate,real>::boundary_face_values;
public:
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    Euler ( 
        const Parameters::AllParameters *const                    parameters_input,
        const double                                              ref_length,
        const double                                              gamma_gas,
        const double                                              mach_inf,
        const double                                              angle_of_attack,
        const double                                              side_slip_angle,
        std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function = nullptr,
        const two_point_num_flux_enum                             two_point_num_flux_type = two_point_num_flux_enum::KG,
        const bool                                                has_nonzero_diffusion = false,
        const bool                                                has_nonzero_physical_source = false);

    const double ref_length; 
    const double gam; 
    const double gamm1; 

    const double density_inf;

    const double mach_inf; 
    const double mach_inf_sqr; 

    const double angle_of_attack;

    const double side_slip_angle;


    const double sound_inf; 
    const double pressure_inf; 
    const double entropy_inf; 
    const two_point_num_flux_enum two_point_num_flux_type; 
    double temperature_inf; 
    double dynamic_pressure_inf; 

    //const double internal_energy_inf;

    dealii::Tensor<1,dim,double> velocities_inf; // should be const


    // dealii::Tensor<1,dim,double> compute_velocities_inf() const;

    // std::array<real,nstate> manufactured_solution (const dealii::Point<dim,double> &pos) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
        const std::array<real,nstate> &conservative_soln) const override;

    std::array<real,nstate> convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const;

    dealii::Tensor<2,nstate,real> convective_flux_directional_jacobian (
        const std::array<real,nstate> &conservative_soln,
        const dealii::Tensor<1,dim,real> &normal) const;

    std::array<real,nstate> convective_eigenvalues (
        const std::array<real,nstate> &/*conservative_soln*/,
        const dealii::Tensor<1,dim,real> &/*normal*/) const override;

    real max_convective_eigenvalue (const std::array<real,nstate> &soln) const override;


    real max_convective_normal_eigenvalue (
        const std::array<real,nstate> &soln,
        const dealii::Tensor<1,dim,real> &normal) const override;

    real max_viscous_eigenvalue (const std::array<real,nstate> &soln) const override;

    std::array<dealii::Tensor<1,dim,real>,nstate> dissipative_flux (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const dealii::types::global_dof_index cell_index) const;

    virtual std::array<dealii::Tensor<1,dim,real>,nstate> dissipative_flux (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const;

    std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_soln,
        const real current_time,
        const dealii::types::global_dof_index cell_index) const;

    virtual std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_soln,
        const real current_time) const;

    std::array<real,nstate> convective_source_term (
        const dealii::Point<dim,real> &pos) const;

protected:

    template<typename real2>
    bool check_positive_quantity(real2 &quantity, const std::string qty_name) const;

public:
    template<typename real2>
    std::array<real2,nstate> convert_conservative_to_primitive_templated ( const std::array<real2,nstate> &conservative_soln ) const;

    std::array<real,nstate> convert_conservative_to_primitive ( const std::array<real,nstate> &conservative_soln ) const;

    std::array<real,nstate> convert_primitive_to_conservative ( const std::array<real,nstate> &primitive_soln ) const;

    template<typename real2>
    std::array<dealii::Tensor<1,dim,real2>,nstate> 
    convert_conservative_gradient_to_primitive_gradient_templated (
        const std::array<real2,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &conservative_soln_gradient) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> 
    convert_conservative_gradient_to_primitive_gradient (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> 
    convert_primitive_gradient_to_conservative_gradient (
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient) const;

    virtual real compute_gamma ( const std::array<real,nstate> &conservative_soln ) const;
    
    template<typename real2>
    real2 compute_pressure_templated ( const std::array<real2,nstate> &conservative_soln ) const;

    real compute_pressure ( const std::array<real,nstate> &conservative_soln ) const;

    template<typename real2>
    real2 compute_entropy (const real2 density, const real2 pressure) const;

    real compute_specific_enthalpy ( const std::array<real,nstate> &conservative_soln, const real pressure) const;

    template<typename real2>
    real2 compute_entropy_templated (const std::array<real2,nstate> &conservative_soln) const;

    real compute_entropy ( const std::array<real,nstate> &conservative_soln ) const;

    real compute_numerical_entropy_function(const std::array<real,nstate> &conservative_soln) const;

    real compute_sound ( const std::array<real,nstate> &conservative_soln ) const;
    real compute_sound ( const real density, const real pressure ) const;

    template<typename real2>
    dealii::Tensor<1,dim,real2> compute_velocities ( const std::array<real2,nstate> &conservative_soln ) const;
    template<typename real2>
    real2 compute_velocity_squared ( const dealii::Tensor<1,dim,real2> &velocities ) const;

    template<typename real2>
    dealii::Tensor<1,dim,real2> extract_velocities_from_primitive ( const std::array<real2,nstate> &primitive_soln ) const;

    real compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const;

    real compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const;

    real compute_incompressible_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const;

    real compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const;

    real compute_incompressible_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const;


    real compute_entropy_measure ( const std::array<real,nstate> &conservative_soln ) const;

    real compute_entropy_measure ( const real density, const real pressure ) const;

    real compute_mach_number ( const std::array<real,nstate> &conservative_soln ) const;


    template<typename real2>
    real2 compute_temperature ( const std::array<real2,nstate> &primitive_soln ) const;


    real compute_density_from_pressure_temperature ( const real pressure, const real temperature ) const;


    real compute_temperature_from_density_pressure ( const real density, const real pressure ) const;


    real compute_pressure_from_density_temperature ( const real density, const real temperature ) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const override;

    std::array<real,nstate> compute_entropy_variables (
                const std::array<real,nstate> &conservative_soln) const;

    std::array<real,nstate> compute_conservative_variables_from_entropy_variables (
                const std::array<real,nstate> &entropy_var) const;

    std::array<real,nstate> compute_kinetic_energy_variables (
                const std::array<real,nstate> &conservative_soln) const;


    real compute_mean_density(
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &convervative_soln2) const;


    real compute_mean_pressure(
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &convervative_soln2) const;


    dealii::Tensor<1,dim,real> compute_mean_velocities(
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &convervative_soln2) const;


    real compute_mean_specific_total_energy(
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &convervative_soln2) const;

    void boundary_face_values (
        const int /*boundary_type*/,
        const dealii::Point<dim, real> &/*pos*/,
        const dealii::Tensor<1,dim,real> &/*normal*/,
        const std::array<real,nstate> &/*soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
        std::array<real,nstate> &/*soln_bc*/,
        std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;

    virtual dealii::Vector<double> post_compute_derived_quantities_vector (
        const dealii::Vector<double>              &uh,
        const std::vector<dealii::Tensor<1,dim> > &duh,
        const std::vector<dealii::Tensor<2,dim> > &dduh,
        const dealii::Tensor<1,dim>               &normals,
        const dealii::Point<dim>                  &evaluation_points) const override;
    
    virtual std::vector<std::string> post_get_names () const override;
    
    virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const override;
    
    virtual dealii::UpdateFlags post_get_needed_update_flags () const override;

protected:
    void boundary_slip_wall (
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;

    void boundary_wall (
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;

    virtual void boundary_manufactured_solution (
        const dealii::Point<dim, real> &pos,
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;

    void boundary_pressure_outflow (
        const real total_inlet_pressure,
        const real back_pressure,
        const std::array<real,nstate> &soln_int,
        std::array<real,nstate> &soln_bc) const;

    void boundary_inflow (
        const real total_inlet_pressure,
        const real total_inlet_temperature,
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        std::array<real,nstate> &soln_bc) const;

    void boundary_riemann (
       const dealii::Tensor<1,dim,real> &normal_int,
       const std::array<real,nstate> &soln_int,
       std::array<real,nstate> &soln_bc) const;

    void boundary_farfield (
        std::array<real,nstate> &soln_bc) const;

    void boundary_p0_extrapolation (
        const std::array<real,nstate> &soln_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const;

    void boundary_custom (
        std::array<real,nstate> &soln_bc) const;

    void boundary_astrophysical_inflow (
        std::array<real,nstate> &soln_bc) const;

    std::array<real,nstate> get_manufactured_solution_value(
        const dealii::Point<dim,real> &pos) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> get_manufactured_solution_gradient(
        const dealii::Point<dim,real> &pos) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux_kennedy_gruber (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const;

    std::array<real,nstate> compute_ismail_roe_parameter_vector_from_primitive(
        const std::array<real,nstate> &primitive_soln) const;

    real compute_ismail_roe_logarithmic_mean(const real val1, const real val2) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux_ismail_roe (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux_chandrashekar (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux_ranocha (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const;
};

} // Physics namespace
} // PHiLiP namespace

#endif