convective_numerical_flux.cpp

#include <boost/preprocessor/seq/for_each.hpp>

#include "ADTypes.hpp"

#include "convective_numerical_flux.hpp"

namespace PHiLiP {
namespace NumericalFlux {

using AllParam = Parameters::AllParameters;

// Protyping low level functions
template<int nstate, typename real_tensor>
std::array<real_tensor, nstate> array_average(
    const std::array<real_tensor, nstate> &array1,
    const std::array<real_tensor, nstate> &array2)
{
    std::array<real_tensor,nstate> array_average;
    for (int s=0; s<nstate; s++) {
        array_average[s] = 0.5*(array1[s] + array2[s]);
    }
    return array_average;
}

template <int dim, int nspecies, int nstate, typename real>
NumericalFluxConvective<dim, nspecies, nstate, real>::NumericalFluxConvective(
    std::unique_ptr< BaselineNumericalFluxConvective<dim,nspecies,nstate,real> > baseline_input,
    std::unique_ptr< RiemannSolverDissipation<dim,nspecies,nstate,real> >   riemann_solver_dissipation_input)
    : baseline(std::move(baseline_input))
    , riemann_solver_dissipation(std::move(riemann_solver_dissipation_input))
{ }

template<int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> NumericalFluxConvective<dim,nspecies,nstate,real>
::evaluate_flux (
    const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const dealii::Tensor<1,dim,real> &normal_int) const
{
    // baseline flux (without upwind dissipation)
    const std::array<real, nstate> baseline_flux_dot_n 
        = this->baseline->evaluate_flux(soln_int, soln_ext, normal_int);

    // Riemann solver dissipation
    const std::array<real, nstate> riemann_solver_dissipation_dot_n 
        = this->riemann_solver_dissipation->evaluate_riemann_solver_dissipation(soln_int, soln_ext, normal_int);

    // convective numerical flux: sum of baseline and Riemann solver dissipation term
    std::array<real, nstate> numerical_flux_dot_n;
    for (int s=0; s<nstate; s++) {
        numerical_flux_dot_n[s] = baseline_flux_dot_n[s] + riemann_solver_dissipation_dot_n[s];
    }
    return numerical_flux_dot_n;
}

template <int dim, int nspecies, int nstate, typename real>
LaxFriedrichs<dim, nspecies, nstate, real>::LaxFriedrichs(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< CentralBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< LaxFriedrichsRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
RoePike<dim, nspecies, nstate, real>::RoePike(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< CentralBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< RoePikeRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
L2Roe<dim, nspecies, nstate, real>::L2Roe(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< CentralBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< L2RoeRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
Central<dim, nspecies, nstate, real>::Central(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< CentralBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< ZeroRiemannSolverDissipation<dim, nspecies, nstate, real> > ())
{}

template <int dim, int nspecies, int nstate, typename real>
EntropyConserving<dim, nspecies, nstate, real>::EntropyConserving(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< EntropyConservingBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< ZeroRiemannSolverDissipation<dim, nspecies, nstate, real> > ())
{}

template <int dim, int nspecies, int nstate, typename real>
EntropyConservingWithLaxFriedrichsDissipation<dim, nspecies, nstate, real>::EntropyConservingWithLaxFriedrichsDissipation(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< EntropyConservingBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< LaxFriedrichsRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
EntropyConservingWithRoeDissipation<dim, nspecies, nstate, real>::EntropyConservingWithRoeDissipation(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< EntropyConservingBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< RoePikeRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
EntropyConservingWithL2RoeDissipation<dim, nspecies, nstate, real>::EntropyConservingWithL2RoeDissipation(
    std::shared_ptr<Physics::PhysicsBase<dim, nspecies, nstate, real>> physics_input)
    : NumericalFluxConvective<dim,nspecies,nstate,real>(
        std::make_unique< EntropyConservingBaselineNumericalFluxConvective<dim, nspecies, nstate, real> > (physics_input), 
        std::make_unique< L2RoeRiemannSolverDissipation<dim, nspecies, nstate, real> > (physics_input))
{}

template <int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> CentralBaselineNumericalFluxConvective<dim,nspecies,nstate,real>::evaluate_flux(
 const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const dealii::Tensor<1,dim,real> &normal_int) const
{
    using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
    RealArrayVector conv_phys_flux_int;
    RealArrayVector conv_phys_flux_ext;

    conv_phys_flux_int = pde_physics->convective_flux (soln_int);
    conv_phys_flux_ext = pde_physics->convective_flux (soln_ext);
    
    RealArrayVector flux_avg;
    for (int s=0; s<nstate; s++) {
        flux_avg[s] = 0.0;
        for (int d=0; d<dim; ++d) {
            flux_avg[s][d] = 0.5*(conv_phys_flux_int[s][d] + conv_phys_flux_ext[s][d]);
        }
    }

    std::array<real, nstate> numerical_flux_dot_n;
    for (int s=0; s<nstate; s++) {
        real flux_dot_n = 0.0;
        for (int d=0; d<dim; ++d) {
            flux_dot_n += flux_avg[s][d]*normal_int[d];
        }
        numerical_flux_dot_n[s] = flux_dot_n;
    }
    return numerical_flux_dot_n;
}

template <int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> EntropyConservingBaselineNumericalFluxConvective<dim,nspecies,nstate,real>::evaluate_flux(
 const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const dealii::Tensor<1,dim,real> &normal_int) const
{
    using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
    RealArrayVector conv_phys_split_flux;

    conv_phys_split_flux = pde_physics->convective_numerical_split_flux (soln_int,soln_ext);

    // Scalar dissipation
    std::array<real, nstate> numerical_flux_dot_n;
    for (int s=0; s<nstate; s++) {
        real flux_dot_n = 0.0;
        for (int d=0; d<dim; ++d) {
            flux_dot_n += conv_phys_split_flux[s][d] * normal_int[d];
        }
        numerical_flux_dot_n[s] = flux_dot_n;
    }
    return numerical_flux_dot_n;
}

template<int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> ZeroRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_riemann_solver_dissipation (
    const std::array<real, nstate> &/*soln_int*/,
    const std::array<real, nstate> &/*soln_ext*/,
    const dealii::Tensor<1,dim,real> &/*normal_int*/) const
{
    // zero upwind dissipation
    std::array<real, nstate> numerical_flux_dot_n;
    numerical_flux_dot_n.fill(0.0);
    return numerical_flux_dot_n;
}

template<int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> LaxFriedrichsRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_riemann_solver_dissipation (
    const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const dealii::Tensor<1,dim,real> &normal_int) const
{
    const real conv_max_eig_int = pde_physics->max_convective_normal_eigenvalue(soln_int,normal_int);
    const real conv_max_eig_ext = pde_physics->max_convective_normal_eigenvalue(soln_ext,normal_int);
    // Replaced the std::max with an if-statement for the AD to work properly.
    //const real conv_max_eig = std::max(conv_max_eig_int, conv_max_eig_ext);
    real conv_max_eig;
    if (conv_max_eig_int > conv_max_eig_ext) {
        conv_max_eig = conv_max_eig_int;
    } else {
        conv_max_eig = conv_max_eig_ext;
    }
    //conv_max_eig = std::max(conv_max_eig_int, conv_max_eig_ext);
    // Scalar dissipation
    std::array<real, nstate> numerical_flux_dot_n;
    for (int s=0; s<nstate; s++) {
        numerical_flux_dot_n[s] = - 0.5 * conv_max_eig * (soln_ext[s]-soln_int[s]);
    }

    return numerical_flux_dot_n;
}

template <int dim, int nspecies, int nstate, typename real>
void RoePikeRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_entropy_fix (
    const std::array<real, 3> &eig_L,
    const std::array<real, 3> &eig_R,
    std::array<real, 3> &eig_RoeAvg,
    const real /*vel2_ravg*/,
    const real /*sound_ravg*/) const
{
    // Harten's entropy fix
    // -- See Blazek 2015, p.103-105
    for(int e=0;e<3;e++) {
        const real eps = std::max(abs(eig_RoeAvg[e] - eig_L[e]), abs(eig_R[e] - eig_RoeAvg[e]));
        if(eig_RoeAvg[e] < eps) {
            eig_RoeAvg[e] = 0.5*(eig_RoeAvg[e] * eig_RoeAvg[e]/eps + eps);
        }
    }
}

template <int dim, int nspecies, int nstate, typename real>
void RoePikeRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_additional_modifications (
    const std::array<real, nstate> &/*soln_int*/,
    const std::array<real, nstate> &/*soln_ext*/,
    const std::array<real, 3> &/*eig_L*/,
    const std::array<real, 3> &/*eig_R*/,
    real &/*dV_normal*/, 
    dealii::Tensor<1,dim,real> &/*dV_tangent*/
    ) const
{
    // No additional modifications for the Roe-Pike scheme
}

template <int dim, int nspecies, int nstate, typename real>
void L2RoeRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_shock_indicator (
    const std::array<real, 3> &eig_L,
    const std::array<real, 3> &eig_R,
    int &ssw_LEFT,
    int &ssw_RIGHT) const
{
    // Shock indicator of Wada & Liou (1994 Flux) -- Eq.(39)
    // -- See also p.74 of Osswald et al. (2016 L2Roe)
    
    ssw_LEFT=0; ssw_RIGHT=0; // initialize
    
    // ssw_L: i=L --> j=R
    if((eig_L[0]>0.0 && eig_R[0]<0.0) || (eig_L[2]>0.0 && eig_R[2]<0.0)) {
        ssw_LEFT = 1;
    }
    
    // ssw_R: i=R --> j=L
    if((eig_R[0]>0.0 && eig_L[0]<0.0) || (eig_R[2]>0.0 && eig_L[2]<0.0)) {
        ssw_RIGHT = 1;
    }
}

template <int dim, int nspecies, int nstate, typename real>
void L2RoeRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_entropy_fix (
    const std::array<real, 3> &eig_L,
    const std::array<real, 3> &eig_R,
    std::array<real, 3> &eig_RoeAvg,
    const real vel2_ravg,
    const real sound_ravg) const
{
    // Van Leer et al. (1989 Sonic) entropy fix for acoustic waves
    // -- p.74 of Osswald et al. (2016 L2Roe)
    for(int e=0;e<3;e++) {
        if(e!=1) {
            // const real deig = std::max((eig_R[e]-eig_L[e]), 0.0);
            const real deig = std::max(static_cast<real>(eig_R[e] - eig_L[e]), static_cast<real>(0.0));
            if(eig_RoeAvg[e] < 2.0*deig) {
                eig_RoeAvg[e] = 0.25*(eig_RoeAvg[e] * eig_RoeAvg[e]/deig) + deig;
            }
        }
    }
    
    // Entropy fix of Liou (2000 Mass)
    // -- p.74 of Osswald et al. (2016 L2Roe)
    int ssw_L, ssw_R;
    evaluate_shock_indicator(eig_L,eig_R,ssw_L,ssw_R);
    if(ssw_L!=0 || ssw_R!=0) {
        eig_RoeAvg[1] = std::max(sound_ravg, static_cast<real>(sqrt(vel2_ravg)));
    }
}

template <int dim, int nspecies, int nstate, typename real>
void L2RoeRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_additional_modifications  (
    const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const std::array<real, 3> &eig_L,
    const std::array<real, 3> &eig_R,
    real &dV_normal, 
    dealii::Tensor<1,dim,real> &dV_tangent) const
{
    const real mach_number_L = this->euler_physics->compute_mach_number(soln_int);
    const real mach_number_R = this->euler_physics->compute_mach_number(soln_ext);

    // Osswald's two modifications to Roe-Pike scheme --> L2Roe
    // - Blending factor (variable 'z' in reference)
    const real blending_factor = std::min(static_cast<real>(1.0), std::max(mach_number_L,mach_number_R));
    // - Scale jump in (1) normal and (2) tangential velocities
    int ssw_L, ssw_R;
    evaluate_shock_indicator(eig_L,eig_R,ssw_L,ssw_R);
    if(ssw_L==0 && ssw_R==0)
    {
        dV_normal *= blending_factor;
        for (int d=0;d<dim;d++)
        {
            dV_tangent[d] *= blending_factor;
        }
    }
}

template <int dim, int nspecies, int nstate, typename real>
std::array<real, nstate> RoeBaseRiemannSolverDissipation<dim,nspecies,nstate,real>
::evaluate_riemann_solver_dissipation (
    const std::array<real, nstate> &soln_int,
    const std::array<real, nstate> &soln_ext,
    const dealii::Tensor<1,dim,real> &normal_int) const
{
    // See Blazek 2015, p.103-105
    // -- Note: Modified calculation of alpha_{3,4} to use 
    //          dVt (jump in tangential velocities);
    //          expressions are equivalent
    
    // Blazek 2015
    // p. 103-105
    // Note: This is in fact the Roe-Pike method of Roe & Pike (1984 - Efficient)
    const std::array<real,nstate> prim_soln_int = euler_physics->convert_conservative_to_primitive(soln_int);
    const std::array<real,nstate> prim_soln_ext = euler_physics->convert_conservative_to_primitive(soln_ext);
    // Left cell
    const real density_L = prim_soln_int[0];
    const dealii::Tensor< 1,dim,real > velocities_L = euler_physics->extract_velocities_from_primitive(prim_soln_int);
    const real pressure_L = prim_soln_int[nstate-1];

    //const real normal_vel_L = velocities_L*normal_int;
    real normal_vel_L = 0.0;
    for (int d=0; d<dim; ++d) {
        normal_vel_L+= velocities_L[d]*normal_int[d];
    }
    const real specific_enthalpy_L = euler_physics->compute_specific_enthalpy(soln_int, pressure_L);

    // Right cell
    const real density_R = prim_soln_ext[0];
    const dealii::Tensor< 1,dim,real > velocities_R = euler_physics->extract_velocities_from_primitive(prim_soln_ext);
    const real pressure_R = prim_soln_ext[nstate-1];

    //const real normal_vel_R = velocities_R*normal_int;
    real normal_vel_R = 0.0;
    for (int d=0; d<dim; ++d) {
        normal_vel_R+= velocities_R[d]*normal_int[d];
    }
    const real specific_enthalpy_R = euler_physics->compute_specific_enthalpy(soln_ext, pressure_R);

    // Roe-averaged states
    const real r = sqrt(density_R/density_L);
    const real rp1 = r+1.0;

    const real density_ravg = r*density_L;
    //const dealii::Tensor< 1,dim,real > velocities_ravg = (r*velocities_R + velocities_L) / rp1;
    dealii::Tensor< 1,dim,real > velocities_ravg;
    for (int d=0; d<dim; ++d) {
        velocities_ravg[d] = (r*velocities_R[d] + velocities_L[d]) / rp1;
    }
    const real specific_total_enthalpy_ravg = (r*specific_enthalpy_R + specific_enthalpy_L) / rp1;

    const real vel2_ravg = euler_physics->compute_velocity_squared (velocities_ravg);
    //const real normal_vel_ravg = velocities_ravg*normal_int;
    real normal_vel_ravg = 0.0;
    for (int d=0; d<dim; ++d) {
        normal_vel_ravg += velocities_ravg[d]*normal_int[d];
    }

    const real sound2_ravg = euler_physics->gamm1*(specific_total_enthalpy_ravg-0.5*vel2_ravg);
    real sound_ravg = 1e10;
    if (sound2_ravg > 0.0) {
        sound_ravg = sqrt(sound2_ravg);
    }

    // Compute eigenvalues
    std::array<real, 3> eig_ravg;
    eig_ravg[0] = abs(normal_vel_ravg-sound_ravg);
    eig_ravg[1] = abs(normal_vel_ravg);
    eig_ravg[2] = abs(normal_vel_ravg+sound_ravg);

    const real sound_L = euler_physics->compute_sound(density_L, pressure_L);
    std::array<real, 3> eig_L;
    eig_L[0] = abs(normal_vel_L-sound_L);
    eig_L[1] = abs(normal_vel_L);
    eig_L[2] = abs(normal_vel_L+sound_L);

    const real sound_R = euler_physics->compute_sound(density_R, pressure_R);
    std::array<real, 3> eig_R;
    eig_R[0] = abs(normal_vel_R-sound_R);
    eig_R[1] = abs(normal_vel_R);
    eig_R[2] = abs(normal_vel_R+sound_R);

    // Jumps in pressure and density
    const real dp = pressure_R - pressure_L;
    const real drho = density_R - density_L;

    // Jump in normal velocity
    real dVn = normal_vel_R-normal_vel_L;

    // Jumps in tangential velocities
    dealii::Tensor<1,dim,real> dVt;
    for (int d=0;d<dim;d++) {
        dVt[d] = (velocities_R[d] - velocities_L[d]) - dVn*normal_int[d];
    }

    // Evaluate entropy fix on wave speeds
    evaluate_entropy_fix (eig_L, eig_R, eig_ravg, vel2_ravg, sound_ravg);

    // Evaluate additional modifications to the Roe-Pike scheme (if applicable)
    evaluate_additional_modifications (soln_int, soln_ext, eig_L, eig_R, dVn, dVt);

    // Product of eigenvalues and wave strengths
    real coeff[4];
    coeff[0] = eig_ravg[0]*(dp-density_ravg*sound_ravg*dVn)/(2.0*sound2_ravg);
    coeff[1] = eig_ravg[1]*(drho - dp/sound2_ravg);
    coeff[2] = eig_ravg[1]*density_ravg;
    coeff[3] = eig_ravg[2]*(dp+density_ravg*sound_ravg*dVn)/(2.0*sound2_ravg);

    // Evaluate |A_Roe| * (W_R - W_L)
    std::array<real,nstate> AdW;

    // Vn-c (i=1)
    AdW[0] = coeff[0] * 1.0;
    for (int d=0;d<dim;d++) {
        AdW[1+d] = coeff[0] * (velocities_ravg[d] - sound_ravg * normal_int[d]);
    }
    AdW[nstate-1] = coeff[0] * (specific_total_enthalpy_ravg - sound_ravg*normal_vel_ravg);

    // Vn (i=2)
    AdW[0] += coeff[1] * 1.0;
    for (int d=0;d<dim;d++) {
        AdW[1+d] += coeff[1] * velocities_ravg[d];
    }
    AdW[nstate-1] += coeff[1] * vel2_ravg * 0.5;

    // (i=3,4)
    AdW[0] += coeff[2] * 0.0;
    real dVt_dot_vel_ravg = 0.0;
    for (int d=0;d<dim;d++) {
        AdW[1+d] += coeff[2]*dVt[d];
        dVt_dot_vel_ravg += velocities_ravg[d]*dVt[d];
    }
    AdW[nstate-1] += coeff[2]*dVt_dot_vel_ravg;

    // Vn+c (i=5)
    AdW[0] += coeff[3] * 1.0;
    for (int d=0;d<dim;d++) {
        AdW[1+d] += coeff[3] * (velocities_ravg[d] + sound_ravg * normal_int[d]);
    }
    AdW[nstate-1] += coeff[3] * (specific_total_enthalpy_ravg + sound_ravg*normal_vel_ravg);

    std::array<real, nstate> numerical_flux_dot_n;
    for (int s=0; s<nstate; s++) {
        numerical_flux_dot_n[s] = - 0.5 * AdW[s];
    }

    return numerical_flux_dot_n;
}

#if PHILIP_SPECIES==1
    // Define a sequence of indices representing the range [1, 6]
    #define POSSIBLE_NSTATE (1)(2)(3)(4)(5)(6)

    // Define a macro to instantiate functions for a specific nstate
    #define INSTANTIATE_UPTO_NSTATE6(r, data, nstate) \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_UPTO_NSTATE6, _, POSSIBLE_NSTATE)

    // Removes nstate 6 from list, if you want to compile the function in INSTANTIATE_UPTO_NSTATE5 for nstate 6, then comment the two lines defining new nstate list
    #undef POSSIBLE_NSTATE
    #define POSSIBLE_NSTATE (1)(2)(3)(4)(5)

    #define INSTANTIATE_UPTO_NSTATE5(r, data, nstate) \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class EntropyConservingBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class EntropyConservingBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class EntropyConservingBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class EntropyConservingBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class EntropyConservingBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>; \
        \
        template class ZeroRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class ZeroRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class ZeroRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class ZeroRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class ZeroRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_UPTO_NSTATE5, _, POSSIBLE_NSTATE)

    // Instantiation for NSTATE = DIM + 2 with different types
    #define POSSIBLE_TYPE (double)(FadType)(RadType)(FadFadType)(RadFadType)

    #define INSTANTIATE_TYPES(r, data, type) \
        template class RoePike<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>; \
        template class L2Roe<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type >; \
        template class EntropyConservingWithRoeDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>; \
        template class EntropyConservingWithL2RoeDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>; \
        template class RoeBaseRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>; \
        template class RoePikeRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>; \
        template class L2RoeRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+2, type>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TYPES, _, POSSIBLE_TYPE)
#else
    #define POSSIBLE_TYPE (double)(FadType)(RadType)(FadFadType)(RadFadType)
    #define INSTANTIATE_TYPES(r, data, type) \
        template class NumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class Central<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class LaxFriedrichs<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class BaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class CentralBaselineNumericalFluxConvective<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class RiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class LaxFriedrichsRiemannSolverDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class EntropyConserving<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>; \
        template class EntropyConservingWithLaxFriedrichsDissipation<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TYPES, _, POSSIBLE_TYPE)
#endif
} // NumericalFlux namespace
} // PHiLiP namespace