periodic_turbulence.h

#ifndef __PERIODIC_TURBULENCE_H__
#define __PERIODIC_TURBULENCE_H__

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

#include "dg/dg_base.hpp"
#include "periodic_cube_flow.h"
#include "physics/navier_stokes.h"

namespace PHiLiP {
namespace FlowSolver {

template <int dim, int nstate>
class PeriodicTurbulence : public PeriodicCubeFlow<dim,nstate>
{
    static const int NUMBER_OF_INTEGRATED_QUANTITIES = 5;
public:
    explicit PeriodicTurbulence(const Parameters::AllParameters *const parameters_input);

    void compute_and_update_integrated_quantities(DGBase<dim, double> &dg);

    double get_integrated_kinetic_energy() const;

    double get_integrated_enstrophy() const;

    double get_vorticity_based_dissipation_rate() const;

    double get_pressure_dilatation_based_dissipation_rate () const;

    double get_deviatoric_strain_rate_tensor_based_dissipation_rate() const;

    double get_strain_rate_tensor_based_dissipation_rate() const;

    void output_velocity_field(
            std::shared_ptr<DGBase<dim,double>> dg,
            const unsigned int output_file_index,
            const double current_time) const;

    double compute_current_integrated_numerical_entropy(
            const std::shared_ptr <DGBase<dim, double>> dg) const;

    void update_numerical_entropy(
            const double FR_entropy_contribution_RRK_solver,
            const unsigned int current_iteration,
            const std::shared_ptr <DGBase<dim, double>> dg);
    
    double get_numerical_entropy(const std::shared_ptr <DGBase<dim, double>> /*dg*/) const;

protected:
    const std::string unsteady_data_table_filename_with_extension;

    const unsigned int number_of_times_to_output_velocity_field;

    const bool output_velocity_field_at_fixed_times;

    const bool output_vorticity_magnitude_field_in_addition_to_velocity;

    const std::string output_flow_field_files_directory_name;

    const bool output_solution_at_exact_fixed_times;

    std::shared_ptr< Physics::NavierStokes<dim,dim+2,double> > navier_stokes_physics;

    bool is_taylor_green_vortex = false; 
    bool is_decaying_homogeneous_isotropic_turbulence = false; 
    bool is_viscous_flow = true; 
    bool do_calculate_numerical_entropy = false; 

    void display_additional_flow_case_specific_parameters() const override;

    double get_constant_time_step(std::shared_ptr<DGBase<dim,double>> dg) const override;

    using CubeFlow_UniformGrid<dim, nstate>::get_adaptive_time_step;

    using CubeFlow_UniformGrid<dim, nstate>::get_adaptive_time_step_initial;

    using CubeFlow_UniformGrid<dim, nstate>::update_maximum_local_wave_speed;

    using FlowSolverCaseBase<dim,nstate>::compute_unsteady_data_and_write_to_table;
    void compute_unsteady_data_and_write_to_table(
            const std::shared_ptr<ODE::ODESolverBase<dim, double>> ode_solver, 
            const std::shared_ptr <DGBase<dim, double>> dg,
            const std::shared_ptr<dealii::TableHandler> unsteady_data_table) override;

    enum IntegratedQuantitiesEnum {
        kinetic_energy,
        enstrophy,
        pressure_dilatation,
        deviatoric_strain_rate_tensor_magnitude_sqr,
        strain_rate_tensor_magnitude_sqr
    };
    std::array<double,NUMBER_OF_INTEGRATED_QUANTITIES> integrated_quantities;

    double maximum_local_wave_speed;

    double previous_numerical_entropy = 0;

    double cumulative_numerical_entropy_change_FRcorrected = 0;

    double initial_numerical_entropy_abs = 0;

    dealii::Table<1,double> output_velocity_field_times;

    unsigned int index_of_current_desired_time_to_output_velocity_field;

    std::string flow_field_quantity_filename_prefix;

    std::shared_ptr<dealii::TableHandler> exact_output_times_of_velocity_field_files_table;
};

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