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 nspecies, int nstate>
class PeriodicTurbulence : public PeriodicCubeFlow<dim,nspecies,nstate>
{
    static const int NUMBER_OF_INTEGRATED_QUANTITIES = 9;

public:
    explicit PeriodicTurbulence(const Parameters::AllParameters *const parameters_input);

    void output_velocity_field_if_current_time_is_output_time(
        const double current_time,
        const std::shared_ptr <DGBase<dim, nspecies, double>> dg);
    
    void compute_and_update_integrated_quantities(DGBase<dim, nspecies, double> &dg);

    double get_integrated_kinetic_energy() const;

    double get_integrated_enstrophy() const;

    double get_integrated_incompressible_kinetic_energy() const;

    double get_integrated_incompressible_enstrophy() const;

    double get_integrated_incompressible_palinstrophy() const;

    double get_integrated_angular_momentum() 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;

    virtual unsigned int get_number_of_degrees_of_freedom_per_state_from_poly_degree(const unsigned int poly_degree_input) const;

    void output_velocity_field(
            std::shared_ptr<DGBase<dim,nspecies,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, nspecies, 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, nspecies, double>> dg);
    
    double get_numerical_entropy(const std::shared_ptr <DGBase<dim, nspecies, 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 bool output_density_field_in_addition_to_velocity;

    const bool output_viscosity_field_in_addition_to_velocity;

    const std::string output_flow_field_files_directory_name;

    const bool output_solution_at_exact_fixed_times;

    const unsigned int output_velocity_number_of_subvisions; 

    const bool do_compute_angular_momentum;

    std::shared_ptr< Physics::NavierStokes<dim,nspecies,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; 

    virtual void display_additional_flow_case_specific_parameters() const override;

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

protected:
    double get_adaptive_time_step(std::shared_ptr<DGBase<dim,nspecies,double>> dg) const override;

    void update_maximum_local_wave_speed(DGBase<dim, nspecies, double> &dg) override;

    double compute_angular_momentum(const dealii::Point<dim> position, const dealii::Tensor<1,3,double> vorticity) const;

public:
    virtual void compute_unsteady_data_and_write_to_table(
            const std::shared_ptr<ODE::ODESolverBase<dim, nspecies, double>> ode_solver, 
            const std::shared_ptr <DGBase<dim, nspecies, double>> dg,
            const std::shared_ptr<dealii::TableHandler> unsteady_data_table,
            const bool do_write_unsteady_data_table_file) override;

protected:
    enum IntegratedQuantitiesEnum {
        kinetic_energy,
        enstrophy,
        pressure_dilatation,
        viscosity_times_deviatoric_strain_rate_tensor_magnitude_sqr,
        viscosity_times_strain_rate_tensor_magnitude_sqr,
        incompressible_kinetic_energy,
        incompressible_enstrophy,
        incompressible_palinstrophy,
        angular_momentum
    };
    std::array<double,NUMBER_OF_INTEGRATED_QUANTITIES> integrated_quantities;

    double integrated_kinetic_energy_at_previous_time_step;

    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