mhd.h

#ifndef __MHD__
#define __MHD__

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

namespace PHiLiP {
namespace Physics {


template <int dim, int nstate, typename real>
class MHD : public PhysicsBase <dim, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using PhysicsBase<dim,nstate,real>::dissipative_flux;
    using PhysicsBase<dim,nstate,real>::source_term;
public:
    MHD(
        const Parameters::AllParameters *const                    parameters_input,
        const double                                              gamma_gas, 
        const dealii::Tensor<2,3,double>                          input_diffusion_tensor = Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(),
        std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function = nullptr,
        const bool                                                has_nonzero_diffusion = false,
        const bool                                                has_nonzero_physical_source = false)
    : PhysicsBase<dim,nstate,real>(parameters_input, has_nonzero_diffusion, has_nonzero_physical_source, input_diffusion_tensor, manufactured_solution_function)
    , gam(gamma_gas)
    , gamm1(gam-1.0)
    {
        static_assert(nstate==8, "Physics::MHD() should be created with nstate=8");

    };

    const double gam;
    const double gamm1;

   // double mach_inf_sqr = 1;

    //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;

    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;

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

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

    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;

    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;

    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> 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;

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

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

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

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

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

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

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


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

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

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


    real compute_temperature ( const std::array<real,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;
    
    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;


    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_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;

//    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;
//    virtual std::vector<std::string> post_get_names () const;
//    virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const;
//    virtual dealii::UpdateFlags post_get_needed_update_flags () const;
protected:


};

} // Physics namespace
} // PHiLiP namespace

#endif