burgers.h

#ifndef __BURGERS__
#define __BURGERS__

#include <deal.II/base/tensor.h>

#include "parameters/all_parameters.h"
#include "parameters/parameters_manufactured_solution.h"

#include "physics.h"

namespace PHiLiP {
namespace Physics {

template <int dim, int nspecies, int nstate, typename real>
class Burgers : 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;
protected:
    double diffusion_scaling_coeff;
public:

    const bool hasConvection;
    const bool hasDiffusion;
    const Parameters::AllParameters::TestType test_type; 

    Burgers(
        const Parameters::AllParameters *const                    parameters_input,
        const double                                              diffusion_coefficient,
        const bool                                                convection = true, 
        const bool                                                diffusion = true,
        const dealii::Tensor<2,3,double>                          input_diffusion_tensor = Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(),
        std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function = nullptr,
        const Parameters::AllParameters::TestType                 parameters_test = Parameters::AllParameters::TestType::run_control,
        const bool                                                has_nonzero_physical_source = false);

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (const std::array<real,nstate> &solution) 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> 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;

    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;

    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> convective_eigenvalues (
        const std::array<real,nstate> &/*solution*/,
        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;

    //  /// Diffusion matrix is identity
    //  std::array<dealii::Tensor<1,dim,real>,nstate> apply_diffusion_matrix (
    //      const std::array<real,nstate> &solution,
    //      const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_grad) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> dissipative_flux (
        const std::array<real,nstate> &solution,
        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> &solution,
        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> &solution,
        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> &solution,
        const real current_time) 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;

protected:
    real diffusion_coefficient () const;
};


} // Physics namespace
} // PHiLiP namespace

#endif