flow_solver.h

#ifndef __FLOW_SOLVER_H__
#define __FLOW_SOLVER_H__

// for FlowSolver class:
#include "flow_solver_cases/flow_solver_case_base.h"
#include "physics/physics.h"
#include "parameters/all_parameters.h"

// for generate_grid
#include <deal.II/base/function.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/vector_tools.h>

#include <stdlib.h>
#include <iostream>
#include <string>
#include <vector>

#include "dg/dg_base.hpp"
#include "dg/dg_factory.hpp"
#include "physics/physics_factory.h"
//#include "ode_solver/runge_kutta_ode_solver.h"
#include "ode_solver/runge_kutta_ode_solver.h"
#include "ode_solver/ode_solver_factory.h"

#include "reduced_order/pod_basis_online.h"

#include <deal.II/base/table_handler.h>

#include <deal.II/base/parameter_handler.h>

namespace PHiLiP {
namespace FlowSolver {

#if PHILIP_DIM==1
    using Triangulation = dealii::Triangulation<PHILIP_DIM>;
#else
    using Triangulation = dealii::parallel::distributed::Triangulation<PHILIP_DIM>;
#endif


class FlowSolverBase
{
public:
    virtual ~FlowSolverBase() = default;


    virtual int run() const = 0;
};


template <int dim, int nstate>
class FlowSolver : public FlowSolverBase
{
public:
    FlowSolver(
        const Parameters::AllParameters *const parameters_input, 
        std::shared_ptr<FlowSolverCaseBase<dim, nstate>> flow_solver_case_input,
        const dealii::ParameterHandler &parameter_handler_input);
    
    std::shared_ptr<FlowSolverCaseBase<dim, nstate>> flow_solver_case;

    const dealii::ParameterHandler &parameter_handler;

    int run () const override;

    void initialize_data_table_from_file(
        std::string data_table_filename_with_extension,
        const std::shared_ptr <dealii::TableHandler> data_table) const;

    std::string get_restart_filename_without_extension(const unsigned int restart_index_input) const;

protected:
    const MPI_Comm mpi_communicator; 
    const int mpi_rank; 
    const int n_mpi; 

    dealii::ConditionalOStream pcout;
    const Parameters::AllParameters all_param; 
    const Parameters::FlowSolverParam flow_solver_param; 
    const Parameters::ODESolverParam ode_param; 
    const unsigned int poly_degree; 
    const unsigned int grid_degree; 
    const double final_time; 

    const std::string input_parameters_file_reference_copy_filename;

    const bool do_output_solution_at_fixed_times; 
    const unsigned int number_of_fixed_times_to_output_solution; 
    const bool output_solution_at_exact_fixed_times;
    
public:
    std::shared_ptr<DGBase<dim, double>> dg;

    std::shared_ptr<ODE::ODESolverBase<dim, double>> ode_solver;

    std::shared_ptr<ProperOrthogonalDecomposition::OnlinePOD<dim>> time_pod;

private:
    std::vector<std::string> get_data_table_column_names(const std::string string_input) const;

    void write_restart_parameter_file(const unsigned int restart_index_input,
                                      const double constant_time_step_input) const;

    std::string double_to_string(const double value_input) const;

#if PHILIP_DIM>1
    void output_restart_files(
        const unsigned int current_restart_index,
        const double constant_time_step,
        const std::shared_ptr <dealii::TableHandler> unsteady_data_table) const;
#endif


    void perform_steady_state_mesh_adaptation() const;

    dealii::Table<1,double> output_solution_fixed_times;
};

} // FlowSolver namespace
} // PHiLiP namespace
#endif