ode_solver_base.h

#ifndef __ODE_SOLVER_BASE__
#define __ODE_SOLVER_BASE__

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/table_handler.h>
#include <deal.II/lac/vector.h>
#include "parameters/all_parameters.h"
#include "limiter/bound_preserving_limiter_factory.hpp"
#include <iostream>
#include <stdexcept>

#include "dg/dg_base.hpp"
#include "parameters/all_parameters.h"
#include "reduced_order/pod_basis_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 ODESolverBase
{
public:
    explicit ODESolverBase(std::shared_ptr< DGBase<dim, real, MeshType> > dg_input,
                           std::shared_ptr< ProperOrthogonalDecomposition::PODBase<dim>> pod); 

    ODESolverBase(std::shared_ptr< DGBase<dim, real, MeshType> >  dg_input);

    virtual ~ODESolverBase() = default; 

    dealii::TableHandler solutions_table;

    void write_ode_solver_steady_state_convergence_data_to_table(
        const unsigned int current_iteration,
        const double current_residual,
        const std::shared_ptr <dealii::TableHandler> data_table) const;

    virtual int steady_state ();



    void initialize_steady_polynomial_ramping (const unsigned int global_final_poly_degree);



    void valid_initial_conditions () const;

    int advance_solution_time (double time_advance);

    virtual void step_in_time(real dt, const bool pseudotime) = 0;

    virtual double get_automatic_error_adaptive_step_size (real dt, const bool pseudotime);

    virtual double get_automatic_initial_step_size (real dt, const bool pseudotime);

    virtual void allocate_ode_system () = 0;

    double residual_norm; 
    double residual_norm_decrease; 

protected:

    double CFL_factor;

    double update_norm; 
    double initial_residual_norm; 

    dealii::LinearAlgebra::distributed::Vector<double> solution_update;

public:
    std::shared_ptr<DGBase<dim,real,MeshType>> dg;

    std::shared_ptr<ProperOrthogonalDecomposition::PODBase<dim>> pod;

    std::unique_ptr<BoundPreservingLimiter<dim,real>> limiter;

protected:
    const Parameters::AllParameters *const all_parameters;

    const Parameters::ODESolverParam ode_param;

public:

    double current_time;


    unsigned int current_iteration;


    double current_desired_time_for_output_solution_every_dt_time_intervals;

    double get_original_time_step() const;

    double get_modified_time_step() const;

protected:
    double original_time_step;
    double modified_time_step;
public:
    

    double FR_entropy_contribution_RRK_solver = 0;
    

    double relaxation_parameter_RRK_solver=1;

protected:
    const MPI_Comm mpi_communicator; 
    const int mpi_rank; 
    dealii::ConditionalOStream pcout; 
};
} // ODE namespace
} // PHiLiP namespace

#endif