physics.h

#ifndef __PHYSICS__
#define __PHYSICS__

#include <deal.II/base/tensor.h>
#include <deal.II/numerics/data_component_interpretation.h>
#include <deal.II/fe/fe_update_flags.h>
#include <deal.II/base/types.h>
#include <deal.II/base/conditional_ostream.h>

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


namespace PHiLiP {
namespace Physics {


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

    using NonPhysicalBehaviorEnum = Parameters::AllParameters::NonPhysicalBehaviorEnum;

    PhysicsBase(
        const Parameters::AllParameters *const                    parameters_input,
        const bool                                                has_nonzero_diffusion_input,
        const bool                                                has_nonzero_physical_source_input,
        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_input = nullptr);

    PhysicsBase(
        const Parameters::AllParameters *const                    parameters_input,
        const bool                                                has_nonzero_diffusion_input,
        const bool                                                has_nonzero_physical_source_input,
        std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function_input = nullptr);

    virtual ~PhysicsBase() = default;

    const bool has_nonzero_diffusion;

    const bool has_nonzero_physical_source;

    const Parameters::AllParameters *const all_parameters;
    
    const NonPhysicalBehaviorEnum non_physical_behavior_type;    

    std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function;

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

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

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

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

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

    virtual 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 = 0;

    virtual 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 = 0;

    virtual std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
        const std::array<real,nstate> &solution) const = 0;

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

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

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


    virtual std::array<real,nstate> convective_eigenvalues (
        const std::array<real,nstate> &/*solution*/,
        const dealii::Tensor<1,dim,real> &/*normal*/) const = 0;

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

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

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

    // /// Evaluate the diffusion matrix \f$ A \f$ such that \f$F_v = A \nabla u\f$.
    // virtual 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 = 0;

    virtual std::array<real,nstate> dissipative_flux_dot_normal (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const bool on_boundary,
        const dealii::types::global_dof_index cell_index,
        const dealii::Tensor<1,dim,real> &normal,
        const int boundary_type);

    virtual 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> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const dealii::types::global_dof_index cell_index);

    virtual 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 = 0;


/*    virtual std::array<dealii::Tensor<1,dim,real>,nstate> artificial_dissipative_flux (
        const real viscosity_coefficient,
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) 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 dealii::types::global_dof_index cell_index) const = 0;

    virtual std::array<real,nstate> physical_source_term (
        const dealii::Point<dim,real> &pos,
        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;

    virtual std::array<real,nstate> artificial_source_term (
        const real viscosity_coefficient,
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &solution) const;

    virtual 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*/,
        const std::array<real,nstate> &/*filtered_soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_soln_grad_int*/,
        std::array<real,nstate> &/*soln_bc*/,
        std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;

    virtual void boundary_face_values_viscous_flux (
        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*/,
        const std::array<real,nstate> &/*filtered_soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_soln_grad_int*/,
        std::array<real,nstate> &/*soln_bc*/,
        std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;

    virtual 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 = 0;


    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 dealii::Vector<double> post_compute_derived_quantities_scalar (
        const double                &uh,
        const dealii::Tensor<1,dim> &/*duh*/,
        const 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;


    template<typename real2> 
    real2 handle_non_physical_result (const std::string message = "") const;

public:


    const double BIG_NUMBER = 1e100;
    
protected:

    dealii::ConditionalOStream pcout;


    dealii::Tensor<2,dim,double> diffusion_tensor;
};
} // Physics namespace
} // PHiLiP namespace

#endif