root_finding_rrk_ode_solver.h

#ifndef __ENTROPY_RRK_ODESOLVER_H_
#define __ENTROPY_RRK_ODESOLVER_H_

#include "dg/dg_base.hpp"
#include "rrk_ode_solver_base.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 RootFindingRRKODESolver: public RRKODESolverBase<dim, real, MeshType>
{
public:
    explicit RootFindingRRKODESolver(
            std::shared_ptr<RKTableauBase<dim,real,MeshType>> rk_tableau_input);

protected:

    real compute_relaxation_parameter(const real dt,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg,
            const std::vector<dealii::LinearAlgebra::distributed::Vector<double>> &rk_stage,
            const dealii::LinearAlgebra::distributed::Vector<double> &solution_update
            ) override;

    real compute_root_function(
            const double gamma,
            const dealii::LinearAlgebra::distributed::Vector<double> &u_n,
            const dealii::LinearAlgebra::distributed::Vector<double> &d,
            const double eta_n,
            const double e,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg) const;

    real compute_numerical_entropy(
            const dealii::LinearAlgebra::distributed::Vector<double> &u,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg) const;
    
    real compute_integrated_numerical_entropy(
            const dealii::LinearAlgebra::distributed::Vector<double> &u,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg) const;
    
    real compute_entropy_change_estimate(const real dt,
            std::shared_ptr<DGBase<dim,real,MeshType>> dg,
            const std::vector<dealii::LinearAlgebra::distributed::Vector<double>> &rk_stage,
            const bool use_M_norm_for_entropy_change_est = true) const;
private:
    real FR_entropy_cumulative = 0;

};

} // ODE namespace
} // PHiLiP namespace

#endif