positivity_preserving_limiter.h

#ifndef __POSITIVITY_PRESERVING_LIMITER__
#define __POSITIVITY_PRESERVING_LIMITER__

#include "bound_preserving_limiter.h"

namespace PHiLiP {
/**********************************
* Zhang, Xiangxiong, and Chi-Wang Shu. 
* "On positivity-preserving high order discontinuous Galerkin schemes 
* for compressible Euler equations on rectangular meshes." 
* Journal of Computational Physics 229.23 (2010): 8918-8934.
**********************************/
/**********************************
* Wang, Cheng, et al.
* "Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations."
* Journal of Computational Physics 231.2 (2012): 653-665.
**********************************/
template<int dim, int nstate, typename real>
class PositivityPreservingLimiter : public BoundPreservingLimiterState <dim, nstate, real>
{
public:
    explicit PositivityPreservingLimiter(
        const Parameters::AllParameters* const parameters_input);

    ~PositivityPreservingLimiter() = default;

    const Parameters::FlowSolverParam flow_solver_param; 

    std::shared_ptr<BoundPreservingLimiterState<dim, nstate, real>> tvbLimiter;

    std::shared_ptr < Physics::Euler<dim, nstate, double > > euler_physics;

    using BoundPreservingLimiterState<dim, nstate, real>::get_soln_cell_avg;

    void limit(
        dealii::LinearAlgebra::distributed::Vector<double>&     solution,
        const dealii::DoFHandler<dim>&                          dof_handler,
        const dealii::hp::FECollection<dim>&                    fe_collection,
        const dealii::hp::QCollection<dim>&                     volume_quadrature_collection,
        const unsigned int                                      grid_degree,
        const unsigned int                                      max_degree,
        const dealii::hp::FECollection<1>                       oneD_fe_collection_1state,
        const dealii::hp::QCollection<1>                        oneD_quadrature_collection,
        double                                                  dt) override;
protected:

    std::array<real, nstate> get_soln_cell_avg_PPL(
        const std::array<std::array<std::vector<real>, nstate>, dim>&        soln_at_q,
        const unsigned int                                                   n_quad_pts,
        const std::vector<real>&                                             quad_weights_GLL,
        const std::vector<real>&                                             quad_weights_GL,
        double&                                                              dt);

    std::vector<real> get_theta2_Zhang2010(
        const std::vector< real >&                      p_lim,
        const std::array<real, nstate>&                 soln_cell_avg,
        const std::array<std::vector<real>, nstate>&    soln_at_q,
        const unsigned int                              n_quad_pts,
        const double                                    eps,
        const double                                    gamma);

    real get_theta2_Wang2012(
        const std::array<std::vector<real>, nstate>&    soln_at_q,
        const unsigned int                              n_quad_pts,
        const double                                    p_avg);

    real get_density_scaling_value(
        const double    density_avg,
        const double    density_min,
        const double    pos_eps,
        const double    p_avg);

    void write_limited_solution(
        dealii::LinearAlgebra::distributed::Vector<double>&     solution,
        const std::array<std::vector<real>, nstate>&            soln_coeff,
        const unsigned int                                      n_shape_fns,
        const std::vector<dealii::types::global_dof_index>&     current_dofs_indices);

    // Values required to compute solution cell average in 2D/3D
    real dx; 
    real dy; 
    real dz; 
}; // End of PositivityPreservingLimiter Class
} // PHiLiP namespace

#endif