convection_diffusion.h

#ifndef __CONVECTION_DIFFUSION__
#define __CONVECTION_DIFFUSION__

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

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

namespace PHiLiP {
namespace Physics {

template <int dim, int nspecies, int nstate, typename real>
class ConvectionDiffusion : 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 linear_advection_velocity[3] = { 1.1, -atan(1)*4.0 / exp(1), exp(1)/(atan(1)*4.0) };
    double diffusion_scaling_coeff = 0.1*atan(1)*4.0/exp(1);
public:
    const bool hasConvection; 

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

    ConvectionDiffusion (
        const Parameters::AllParameters *const                    parameters_input,
        const bool                                                convection = true, 
        const bool                                                diffusion = true,
        const dealii::Tensor<2,3,double>                          input_diffusion_tensor = Parameters::ManufacturedSolutionParam::get_default_diffusion_tensor(),
        const dealii::Tensor<1,3,double>                          input_advection_vector = Parameters::ManufacturedSolutionParam::get_default_advection_vector(),
        const double                                              input_diffusion_coefficient = Parameters::ManufacturedSolutionParam::get_default_diffusion_coefficient(),
        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) : 
            PhysicsBase<dim,nspecies,nstate,real>(parameters_input, diffusion, has_nonzero_physical_source, input_diffusion_tensor, manufactured_solution_function), 
            linear_advection_velocity{input_advection_vector[0], input_advection_vector[1], input_advection_vector[2]},
            diffusion_scaling_coeff(input_diffusion_coefficient),
            hasConvection(convection), 
            hasDiffusion(diffusion),
            test_type(parameters_test)
    {};

    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> &soln1,
        const std::array<real,nstate> &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;

    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:
    dealii::Tensor<1,dim,real> advection_speed () const;
    real diffusion_coefficient () const;
};


} // Physics namespace
} // PHiLiP namespace

#endif