runge_kutta_store_entropy.h

#ifndef __RK_NUM_ENTROPY_H__
#define __RK_NUM_ENTROPY_H__

#include "dg/dg_base.hpp"
#include "ode_solver/runge_kutta_ode_solver.h"

namespace PHiLiP {
namespace ODE {

#if PHILIP_DIM==1
template <int dim, typename real, typename MeshType = dealii::Triangulation<dim>>
#else
template <int dim, typename real, typename MeshType = dealii::parallel::distributed::Triangulation<dim>>
#endif
class RKNumEntropy: public EmptyRRKBase <dim, real, MeshType>
{
public:
    explicit RKNumEntropy(
            std::shared_ptr<RKTableauBase<dim,real,MeshType>> rk_tableau_input);


    real compute_FR_entropy_contribution(const real dt,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg,
            const std::vector<dealii::LinearAlgebra::distributed::Vector<double>> &rk_stage,
            const bool compute_K_norm) const override;
    
    // "using" keyword to prevent compiler complaining
    using EmptyRRKBase<dim, real, MeshType>::compute_FR_entropy_contribution;
    using EmptyRRKBase<dim, real, MeshType>::update_relaxation_parameter;
    using EmptyRRKBase<dim, real, MeshType>::store_stage_solutions;
    

protected:

    std::shared_ptr<RKTableauBase<dim,real,MeshType>> butcher_tableau;

    const int n_rk_stages;

    const MPI_Comm mpi_communicator; 
    const int mpi_rank; 
    dealii::ConditionalOStream pcout; 


    std::vector<dealii::LinearAlgebra::distributed::Vector<double>> rk_stage_solution;


    void store_stage_solutions(const int istage,
            const dealii::LinearAlgebra::distributed::Vector<double> rk_stage_i) override;
    
    dealii::LinearAlgebra::distributed::Vector<double> compute_entropy_vars(
            const dealii::LinearAlgebra::distributed::Vector<double> &u,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg) const;

};

} // ODE namespace
} // PHiLiP namespace

#endif