From b13f5ff03c9ee050d21d1716516908e0b6513072 Mon Sep 17 00:00:00 2001 From: Julien Brillon <43619715+jbrillon@users.noreply.github.com> Date: Tue, 10 May 2022 10:22:05 -0400 Subject: [PATCH] Computation of Theoretical Dissipation Rates for Turbulent Flows (#137) * adding compute_vorticity function to NavierStokes * using physics object to compute kinetic energy in periodic cube flow * adding enstrophy and related kinetic energy dissipation rate computation; moving initial_condition_function to flow_solver_case_base; computing solution gradient at quad point inside integrate_over_domain * correcting nondimensionalization of integrated quantities * adding TGV isothermal initialization in addition to the existing uniform density initialization * bug fix in soln_grad_at_q computation (did not initialize as zero); adding parameters to output enstrophy and vorticity_based_dissipation_rate * adding parameter solution_vtk_files_directory_name * adding pressure dilatation and deviatoric strain-rate tensor based dissipation rates * computing all quantities simultaneously for computational efficiency; removing L2 error initial condition calculation from periodic cube flow * cleaning up parameters; adding viscous flow flag * resolving additional merge conflicts * addressing PR comments; moving PeriodicTurbulence to separate files * using dealii::Tensor<2,dim,real> in NavierStokes to signify tensor * adding periodic_turbulence.cpp/.h files -- forgot to commit earlier --- src/dg/dg.cpp | 12 +- src/parameters/all_parameters.cpp | 5 + src/parameters/all_parameters.h | 3 + src/parameters/parameters_flow_solver.cpp | 30 ++- src/parameters/parameters_flow_solver.h | 19 +- src/physics/euler.cpp | 23 +- src/physics/euler.h | 6 + .../initial_conditions/initial_condition.cpp | 59 ++++- .../initial_conditions/initial_condition.h | 34 ++- src/physics/navier_stokes.cpp | 161 +++++++++++-- src/physics/navier_stokes.h | 59 ++++- src/testing/CMakeLists.txt | 1 + src/testing/flow_solver.cpp | 5 +- src/testing/flow_solver.h | 2 +- .../1D_burgers_rewienski_snapshot.cpp | 4 +- .../1D_burgers_rewienski_snapshot.h | 4 +- .../1d_burgers_viscous_snapshot.cpp | 4 +- .../1d_burgers_viscous_snapshot.h | 4 +- .../flow_solver_case_base.cpp | 19 +- .../flow_solver_cases/flow_solver_case_base.h | 14 +- src/testing/flow_solver_cases/naca0012.cpp | 4 +- src/testing/flow_solver_cases/naca0012.h | 2 +- .../flow_solver_cases/periodic_cube_flow.cpp | 162 +------------ .../flow_solver_cases/periodic_cube_flow.h | 61 +---- .../flow_solver_cases/periodic_turbulence.cpp | 228 ++++++++++++++++++ .../flow_solver_cases/periodic_turbulence.h | 101 ++++++++ .../taylor_green_vortex_energy_check.cpp | 5 +- .../taylor_green_vortex_restart_check.cpp | 5 +- ..._taylor_green_vortex_energy_check_long.prm | 2 +- ...taylor_green_vortex_energy_check_quick.prm | 2 +- ...cous_taylor_green_vortex_restart_check.prm | 2 +- ...rtex_restart_from_parameter_file_check.prm | 2 +- 32 files changed, 752 insertions(+), 292 deletions(-) create mode 100644 src/testing/flow_solver_cases/periodic_turbulence.cpp create mode 100644 src/testing/flow_solver_cases/periodic_turbulence.h diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp index 13b2f53d2..55cd56c6d 100644 --- a/src/dg/dg.cpp +++ b/src/dg/dg.cpp @@ -1736,7 +1736,7 @@ void DGBase::output_face_results_vtk (const unsigned int cycl data_out.set_flags(vtkflags); const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); - std::string filename = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string filename = this->all_parameters->solution_vtk_files_directory_name + "/" + "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; filename += dealii::Utilities::int_to_string(cycle, 4) + "."; filename += dealii::Utilities::int_to_string(iproc, 4); filename += ".vtu"; @@ -1747,13 +1747,13 @@ void DGBase::output_face_results_vtk (const unsigned int cycl if (iproc == 0) { std::vector filenames; for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) { - std::string fn = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string fn = this->all_parameters->solution_vtk_files_directory_name + "/" + "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; fn += dealii::Utilities::int_to_string(cycle, 4) + "."; fn += dealii::Utilities::int_to_string(iproc, 4); fn += ".vtu"; filenames.push_back(fn); } - std::string master_fn = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string master_fn = this->all_parameters->solution_vtk_files_directory_name + "/" + "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; master_fn += dealii::Utilities::int_to_string(cycle, 4) + ".pvtu"; std::ofstream master_output(master_fn); data_out.write_pvtu_record(master_output, filenames); @@ -1853,7 +1853,7 @@ void DGBase::output_results_vtk (const unsigned int cycle)// data_out.set_flags(vtkflags); const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); - std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string filename = this->all_parameters->solution_vtk_files_directory_name + "/" + "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; filename += dealii::Utilities::int_to_string(cycle, 4) + "."; filename += dealii::Utilities::int_to_string(iproc, 4); filename += ".vtu"; @@ -1864,13 +1864,13 @@ void DGBase::output_results_vtk (const unsigned int cycle)// if (iproc == 0) { std::vector filenames; for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) { - std::string fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string fn = this->all_parameters->solution_vtk_files_directory_name + "/" + "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; fn += dealii::Utilities::int_to_string(cycle, 4) + "."; fn += dealii::Utilities::int_to_string(iproc, 4); fn += ".vtu"; filenames.push_back(fn); } - std::string master_fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; + std::string master_fn = this->all_parameters->solution_vtk_files_directory_name + "/" + "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-"; master_fn += dealii::Utilities::int_to_string(cycle, 4) + ".pvtu"; std::ofstream master_output(master_fn); data_out.write_pvtu_record(master_output, filenames); diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp index 986a3e183..2f57a2027 100644 --- a/src/parameters/all_parameters.cpp +++ b/src/parameters/all_parameters.cpp @@ -180,6 +180,10 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm) "Dissipative numerical flux. " "Choices are ."); + prm.declare_entry("solution_vtk_files_directory_name", ".", + dealii::Patterns::FileName(dealii::Patterns::FileName::FileType::input), + "Name of directory for writing solution vtk files. Current directory by default."); + Parameters::LinearSolverParam::declare_parameters (prm); Parameters::ManufacturedConvergenceStudyParam::declare_parameters (prm); Parameters::ODESolverParam::declare_parameters (prm); @@ -314,6 +318,7 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm) if (flux_reconstruction_aux_string == "kPlus") flux_reconstruction_aux_type = kPlus; if (flux_reconstruction_aux_string == "k10Thousand") flux_reconstruction_aux_type = k10Thousand; + solution_vtk_files_directory_name = prm.get("solution_vtk_files_directory_name"); pcout << "Parsing linear solver subsection..." << std::endl; linear_solver_param.parse_parameters (prm); diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h index 44ac0a6af..27ab7b8b2 100644 --- a/src/parameters/all_parameters.h +++ b/src/parameters/all_parameters.h @@ -192,6 +192,9 @@ class AllParameters /// Store flux reconstruction type Flux_Reconstruction_Aux flux_reconstruction_aux_type; + /// Name of directory for writing solution vtk files + std::string solution_vtk_files_directory_name; + /// Declare parameters that can be set as inputs and set up the default options /** This subroutine should call the sub-parameter classes static declare_parameters() * such that each sub-parameter class is responsible to declare their own parameters. diff --git a/src/parameters/parameters_flow_solver.cpp b/src/parameters/parameters_flow_solver.cpp index 30768603b..3d973a0bb 100644 --- a/src/parameters/parameters_flow_solver.cpp +++ b/src/parameters/parameters_flow_solver.cpp @@ -74,11 +74,24 @@ void FlowSolverParam::declare_parameters(dealii::ParameterHandler &prm) dealii::Patterns::FileName(dealii::Patterns::FileName::FileType::input), "Filename of the input mesh: input_mesh_filename.msh"); - prm.enter_subsection("taylor_green_vortex_energy_check"); + prm.enter_subsection("taylor_green_vortex"); { prm.declare_entry("expected_kinetic_energy_at_final_time", "1", dealii::Patterns::Double(0, dealii::Patterns::Double::max_double_value), "For integration test purposes, expected kinetic energy at final time."); + + prm.declare_entry("expected_theoretical_dissipation_rate_at_final_time", "1", + dealii::Patterns::Double(0, dealii::Patterns::Double::max_double_value), + "For integration test purposes, expected theoretical kinetic energy dissipation rate at final time."); + + prm.declare_entry("density_initial_condition_type", "uniform", + dealii::Patterns::Selection( + " uniform | " + " isothermal "), + "The type of density initialization. " + "Choices are " + " ."); } prm.leave_subsection(); } @@ -90,10 +103,10 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm) prm.enter_subsection("flow_solver"); { const std::string flow_case_type_string = prm.get("flow_case_type"); - if (flow_case_type_string == "taylor_green_vortex") {flow_case_type = taylor_green_vortex;} - else if (flow_case_type_string == "burgers_viscous_snapshot") {flow_case_type = burgers_viscous_snapshot;} - else if (flow_case_type_string == "burgers_rewienski_snapshot") {flow_case_type = burgers_rewienski_snapshot;} - else if (flow_case_type_string == "naca0012") {flow_case_type = naca0012;} + if (flow_case_type_string == "taylor_green_vortex") {flow_case_type = taylor_green_vortex;} + else if (flow_case_type_string == "burgers_viscous_snapshot") {flow_case_type = burgers_viscous_snapshot;} + else if (flow_case_type_string == "burgers_rewienski_snapshot") {flow_case_type = burgers_rewienski_snapshot;} + else if (flow_case_type_string == "naca0012") {flow_case_type = naca0012;} final_time = prm.get_double("final_time"); courant_friedrich_lewy_number = prm.get_double("courant_friedrich_lewy_number"); @@ -108,9 +121,14 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm) output_restart_files_every_dt_time_intervals = prm.get_double("output_restart_files_every_dt_time_intervals"); input_mesh_filename = prm.get("input_mesh_filename"); - prm.enter_subsection("taylor_green_vortex_energy_check"); + prm.enter_subsection("taylor_green_vortex"); { expected_kinetic_energy_at_final_time = prm.get_double("expected_kinetic_energy_at_final_time"); + expected_theoretical_dissipation_rate_at_final_time = prm.get_double("expected_theoretical_dissipation_rate_at_final_time"); + + const std::string density_initial_condition_type_string = prm.get("density_initial_condition_type"); + if (density_initial_condition_type_string == "uniform") {density_initial_condition_type = uniform;} + else if (density_initial_condition_type_string == "isothermal") {density_initial_condition_type = isothermal;} } prm.leave_subsection(); } diff --git a/src/parameters/parameters_flow_solver.h b/src/parameters/parameters_flow_solver.h index bdf3c0856..ff37548a2 100644 --- a/src/parameters/parameters_flow_solver.h +++ b/src/parameters/parameters_flow_solver.h @@ -40,9 +40,6 @@ class FlowSolverParam * will read file: input_mesh_filename.msh */ std::string input_mesh_filename; - /** For integration test purposes, expected kinetic energy at final time. */ - double expected_kinetic_energy_at_final_time; - bool restart_computation_from_file; ///< Restart computation from restart file bool output_restart_files; ///< Output the restart files std::string restart_files_directory_name; ///< Name of directory for writing and reading restart files @@ -50,6 +47,22 @@ class FlowSolverParam int output_restart_files_every_x_steps; ///< Outputs the restart files every x steps double output_restart_files_every_dt_time_intervals; ///< Outputs the restart files at time intervals of dt + /** For taylor green vortex integration tests, expected kinetic energy at final time. */ + double expected_kinetic_energy_at_final_time; + + /** For taylor green vortex integration tests, + * expected theoretical kinetic energy dissipation + * rate at final time. */ + double expected_theoretical_dissipation_rate_at_final_time; + + /// For taylor green vortex, selects the type of density initialization + enum DensityInitialConditionType{ + uniform, + isothermal, + }; + /// Selected DensityInitialConditionType from the input file + DensityInitialConditionType density_initial_condition_type; + /// Declares the possible variables and sets the defaults. static void declare_parameters (dealii::ParameterHandler &prm); diff --git a/src/physics/euler.cpp b/src/physics/euler.cpp index 71421d090..fd31d0930 100644 --- a/src/physics/euler.cpp +++ b/src/physics/euler.cpp @@ -204,13 +204,30 @@ template inline real Euler ::compute_total_energy ( const std::array &primitive_soln ) const { - const real density = primitive_soln[0]; const real pressure = primitive_soln[nstate-1]; + const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln); + const real tot_energy = pressure / this->gamm1 + kinetic_energy; + return tot_energy; +} + +template +inline real Euler +::compute_kinetic_energy_from_primitive_solution ( const std::array &primitive_soln ) const +{ + const real density = primitive_soln[0]; const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive(primitive_soln); const real vel2 = compute_velocity_squared(velocities); + const real kinetic_energy = 0.5*density*vel2; + return kinetic_energy; +} - const real tot_energy = pressure / gamm1 + 0.5*density*vel2; - return tot_energy; +template +inline real Euler +::compute_kinetic_energy_from_conservative_solution ( const std::array &conservative_soln ) const +{ + const std::array primitive_soln = convert_conservative_to_primitive(conservative_soln); + const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln); + return kinetic_energy; } template diff --git a/src/physics/euler.h b/src/physics/euler.h index 540069194..0f7fbb130 100644 --- a/src/physics/euler.h +++ b/src/physics/euler.h @@ -207,6 +207,12 @@ class Euler : public PhysicsBase */ real compute_total_energy ( const std::array &primitive_soln ) const; + /// Given primitive variables, returns kinetic energy + real compute_kinetic_energy_from_primitive_solution ( const std::array &primitive_soln ) const; + + /// Given conservative variables, returns kinetic energy + real compute_kinetic_energy_from_conservative_solution ( const std::array &conservative_soln ) const; + /// Evaluate entropy from conservative variables /** Note that it is not the actual entropy since it's missing some constants. * Used to check entropy convergence diff --git a/src/physics/initial_conditions/initial_condition.cpp b/src/physics/initial_conditions/initial_condition.cpp index 72ba8aa0a..3e0a81202 100644 --- a/src/physics/initial_conditions/initial_condition.cpp +++ b/src/physics/initial_conditions/initial_condition.cpp @@ -3,7 +3,7 @@ namespace PHiLiP { // ======================================================== -// TAYLOR GREEN VORTEX -- Initial Condition +// TAYLOR GREEN VORTEX -- Initial Condition (Uniform density) // ======================================================== template InitialConditionFunction_TaylorGreenVortex @@ -14,9 +14,7 @@ ::InitialConditionFunction_TaylorGreenVortex ( , gamma_gas(gamma_gas) , mach_inf(mach_inf) , mach_inf_sqr(mach_inf*mach_inf) -{ - // Nothing to do here yet -} +{} template real InitialConditionFunction_TaylorGreenVortex @@ -29,7 +27,7 @@ ::primitive_value(const dealii::Point &point, const unsigned int istat if(istate==0) { // density - value = 1.0; + value = this->density(point); } if(istate==1) { // x-velocity @@ -45,7 +43,7 @@ ::primitive_value(const dealii::Point &point, const unsigned int istat } if(istate==4) { // pressure - value = 1.0/(gamma_gas*mach_inf_sqr) + (1.0/16.0)*(cos(2.0*x)+cos(2.0*y))*(cos(2.0*z)+2.0); + value = 1.0/(this->gamma_gas*this->mach_inf_sqr) + (1.0/16.0)*(cos(2.0*x)+cos(2.0*y))*(cos(2.0*z)+2.0); } } return value; @@ -69,7 +67,7 @@ ::convert_primitive_to_conversative_value( if(istate==1) value = rho*u; // x-momentum if(istate==2) value = rho*v; // y-momentum if(istate==3) value = rho*w; // z-momentum - if(istate==4) value = p/(1.4-1.0) + 0.5*rho*(u*u + v*v + w*w); // total energy + if(istate==4) value = p/(this->gamma_gas-1.0) + 0.5*rho*(u*u + v*v + w*w); // total energy } return value; @@ -84,6 +82,40 @@ ::value(const dealii::Point &point, const unsigned int istate) const return value; } +template +real InitialConditionFunction_TaylorGreenVortex +::density(const dealii::Point &/*point*/) const +{ + // Note: This is in non-dimensional form (free-stream values as reference) + real value = 0.; + // density + value = 1.0; + return value; +} + +// ======================================================== +// TAYLOR GREEN VORTEX -- Initial Condition (Isothermal density) +// ======================================================== +template +InitialConditionFunction_TaylorGreenVortex_Isothermal +::InitialConditionFunction_TaylorGreenVortex_Isothermal ( + const double gamma_gas, + const double mach_inf) + : InitialConditionFunction_TaylorGreenVortex(gamma_gas,mach_inf) +{} + +template +real InitialConditionFunction_TaylorGreenVortex_Isothermal +::density(const dealii::Point &point) const +{ + // Note: This is in non-dimensional form (free-stream values as reference) + real value = 0.; + // density + value = this->primitive_value(point, 4); // get pressure + value *= this->gamma_gas*this->mach_inf_sqr; + return value; +} + // ======================================================== // 1D BURGERS REWIENSKI -- Initial Condition // ======================================================== @@ -145,9 +177,19 @@ InitialConditionFactory::create_InitialConditionFunction( // Get the flow case type const FlowCaseEnum flow_type = param->flow_solver_param.flow_case_type; if (flow_type == FlowCaseEnum::taylor_green_vortex) { - if constexpr (dim==3 && nstate==dim+2) return std::make_shared >( + if constexpr (dim==3 && nstate==dim+2){ + // Get the density initial condition type + const DensityInitialConditionEnum density_initial_condition_type = param->flow_solver_param.density_initial_condition_type; + if(density_initial_condition_type == DensityInitialConditionEnum::uniform) { + return std::make_shared >( param->euler_param.gamma_gas, param->euler_param.mach_inf); + } else if (density_initial_condition_type == DensityInitialConditionEnum::isothermal) { + return std::make_shared >( + param->euler_param.gamma_gas, + param->euler_param.mach_inf); + } + } } else if (flow_type == FlowCaseEnum::burgers_rewienski_snapshot) { if constexpr (dim==1 && nstate==dim) return std::make_shared > (); } else if (flow_type == FlowCaseEnum::burgers_viscous_snapshot) { @@ -185,6 +227,7 @@ template class InitialConditionFactory ; template class InitialConditionFunction_BurgersViscous; template class InitialConditionFunction_BurgersRewienski; template class InitialConditionFunction_TaylorGreenVortex ; +template class InitialConditionFunction_TaylorGreenVortex_Isothermal ; template class InitialConditionFunction_Zero ; template class InitialConditionFunction_Zero ; template class InitialConditionFunction_Zero ; diff --git a/src/physics/initial_conditions/initial_condition.h b/src/physics/initial_conditions/initial_condition.h index 05b15e12e..788de099f 100644 --- a/src/physics/initial_conditions/initial_condition.h +++ b/src/physics/initial_conditions/initial_condition.h @@ -57,7 +57,7 @@ class FreeStreamInitialConditions : public InitialConditionFunction class InitialConditionFunction_TaylorGreenVortex: public InitialConditionFunction { @@ -65,7 +65,7 @@ class InitialConditionFunction_TaylorGreenVortex: public InitialConditionFunctio using dealii::Function::value; ///< dealii::Function we are templating on public: - /// Constructor for TaylorGreenVortex_InitialCondition + /// Constructor for TaylorGreenVortex_InitialCondition with uniform density /** Calls the Function(const unsigned int n_components) constructor in deal.II * This sets the public attribute n_components = nstate, which can then be accessed * by all the other functions @@ -89,6 +89,34 @@ class InitialConditionFunction_TaylorGreenVortex: public InitialConditionFunctio /// Converts value from: primitive to conservative real convert_primitive_to_conversative_value(const dealii::Point &point, const unsigned int istate = 0) const; + + /// Value of initial condition for density + virtual real density(const dealii::Point &point) const; +}; + +/// Initial Condition Function: Taylor Green Vortex (isothermal density) +template +class InitialConditionFunction_TaylorGreenVortex_Isothermal + : public InitialConditionFunction_TaylorGreenVortex +{ +public: + /// Constructor for TaylorGreenVortex_InitialCondition with isothermal density + /** Calls the Function(const unsigned int n_components) constructor in deal.II + * This sets the public attribute n_components = nstate, which can then be accessed + * by all the other functions + * -- Reference: (1) Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * (2) Brian Vermeire 2014 Thesis + * These initial conditions are given in nondimensional form (free-stream as reference) + */ + InitialConditionFunction_TaylorGreenVortex_Isothermal ( + const double gamma_gas = 1.4, + const double mach_inf = 0.1); + +protected: + /// Value of initial condition for density + real density(const dealii::Point &point) const override; }; /// Initial Condition Function: 1D Burgers Rewienski @@ -130,6 +158,8 @@ class InitialConditionFactory protected: /// Enumeration of all flow solver initial conditions types defined in the Parameters class using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType; + /// Enumeration of all taylor green vortex initial condition sub-types defined in the Parameters class + using DensityInitialConditionEnum = Parameters::FlowSolverParam::DensityInitialConditionType; public: /// Construct InitialConditionFunction object from global parameter file diff --git a/src/physics/navier_stokes.cpp b/src/physics/navier_stokes.cpp index e0547ab46..61bb123a7 100644 --- a/src/physics/navier_stokes.cpp +++ b/src/physics/navier_stokes.cpp @@ -125,9 +125,7 @@ ::compute_scaled_viscosity_coefficient (const std::array &primitiv */ const real2 viscosity_coefficient = compute_viscosity_coefficient(primitive_soln); const real2 scaled_viscosity_coefficient = viscosity_coefficient/reynolds_number_inf; - // print the value for Re - // std::cout << "\n Reynolds number inside compute_scaled_viscosity_coefficient(): " << reynolds_number_inf << "\n" << std::endl; - + return scaled_viscosity_coefficient; } @@ -166,13 +164,148 @@ ::compute_heat_flux ( return heat_flux; } +template +dealii::Tensor<1,3,real> NavierStokes +::compute_vorticity ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const +{ + // Compute the vorticity + dealii::Tensor<1,3,real> vorticity; + for(int d=0; d<3; ++d) { + vorticity[d] = 0.0; + } + if constexpr(dim>1) { + // Get velocity gradient + const std::array,nstate> primitive_soln_gradient = convert_conservative_gradient_to_primitive_gradient(conservative_soln, conservative_soln_gradient); + const dealii::Tensor<2,dim,real> velocities_gradient = extract_velocities_gradient_from_primitive_solution_gradient(primitive_soln_gradient); + if constexpr(dim==2) { + // vorticity exists only in z-component + vorticity[2] = velocities_gradient[1][0] - velocities_gradient[0][1]; // z-component + } + if constexpr(dim==3) { + vorticity[0] = velocities_gradient[2][1] - velocities_gradient[1][2]; // x-component + vorticity[1] = velocities_gradient[0][2] - velocities_gradient[2][0]; // y-component + vorticity[2] = velocities_gradient[1][0] - velocities_gradient[0][1]; // z-component + } + } + return vorticity; +} + +template +real NavierStokes +::compute_enstrophy ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const +{ + // Compute the vorticity + dealii::Tensor<1,3,real> vorticity = compute_vorticity(conservative_soln, conservative_soln_gradient); + // Compute enstrophy + real enstrophy = 0.0; + for(int d=0; d<3; ++d) { + enstrophy += vorticity[d]*vorticity[d]; + } + const real density = conservative_soln[0]; + enstrophy *= 0.5*density; + return enstrophy; +} + +template +real NavierStokes +::compute_vorticity_based_dissipation_rate_from_integrated_enstrophy ( + const real integrated_enstrophy) const +{ + real dissipation_rate = 2.0*integrated_enstrophy/(this->reynolds_number_inf); + return dissipation_rate; +} + +template +real NavierStokes +::compute_pressure_dilatation ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const +{ + // Get pressure + const real pressure = this->template compute_pressure(conservative_soln); + + // Get velocity gradient + const std::array,nstate> primitive_soln_gradient = convert_conservative_gradient_to_primitive_gradient(conservative_soln, conservative_soln_gradient); + const dealii::Tensor<2,dim,real> velocities_gradient = extract_velocities_gradient_from_primitive_solution_gradient(primitive_soln_gradient); + + // Compute the pressure dilatation + real pressure_dilatation = 0.0; + for(int d=0; d +dealii::Tensor<2,dim,real> NavierStokes +::compute_deviatoric_strain_rate_tensor ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const +{ + // Get velocity gradient + const std::array,nstate> primitive_soln_gradient = convert_conservative_gradient_to_primitive_gradient(conservative_soln, conservative_soln_gradient); + const dealii::Tensor<2,dim,real> velocities_gradient = extract_velocities_gradient_from_primitive_solution_gradient(primitive_soln_gradient); + + // Compute the deviatoric strain rate tensor + dealii::Tensor<2,dim,real> deviatoric_strain_rate_tensor; + for(int d1=0; d1 +real NavierStokes +::get_tensor_magnitude_sqr ( + const dealii::Tensor<2,dim,real> &tensor) const +{ + real tensor_magnitude_sqr = 0.0; + for (int i=0; i +real NavierStokes +::compute_deviatoric_strain_rate_tensor_magnitude_sqr ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const +{ + // Compute the deviatoric strain rate tensor + const dealii::Tensor<2,dim,real> deviatoric_strain_rate_tensor = compute_deviatoric_strain_rate_tensor(conservative_soln,conservative_soln_gradient); + // Get magnitude squared + real deviatoric_strain_rate_tensor_magnitude_sqr = get_tensor_magnitude_sqr(deviatoric_strain_rate_tensor); + + return deviatoric_strain_rate_tensor_magnitude_sqr; +} + +template +real NavierStokes +::compute_deviatoric_strain_rate_tensor_based_dissipation_rate_from_integrated_deviatoric_strain_rate_tensor_magnitude_sqr ( + const real integrated_deviatoric_strain_rate_tensor_magnitude_sqr) const +{ + real dissipation_rate = 2.0*integrated_deviatoric_strain_rate_tensor_magnitude_sqr/(this->reynolds_number_inf); + return dissipation_rate; +} + template template -std::array,dim> NavierStokes +dealii::Tensor<2,dim,real2> NavierStokes ::extract_velocities_gradient_from_primitive_solution_gradient ( const std::array,nstate> &primitive_soln_gradient) const { - std::array,dim> velocities_gradient; + dealii::Tensor<2,dim,real2> velocities_gradient; for (int d1=0; d1 template -std::array,dim> NavierStokes +dealii::Tensor<2,dim,real2> NavierStokes ::compute_viscous_stress_tensor ( const std::array &primitive_soln, const std::array,nstate> &primitive_soln_gradient) const @@ -191,7 +324,7 @@ ::compute_viscous_stress_tensor ( /* Nondimensionalized viscous stress tensor, $\bm{\tau}^{*}$ * Reference: Masatsuka 2018 "I do like CFD", p.148, eq.(4.14.12) */ - const std::array,dim> vel_gradient = extract_velocities_gradient_from_primitive_solution_gradient(primitive_soln_gradient); + const dealii::Tensor<2,dim,real2> vel_gradient = extract_velocities_gradient_from_primitive_solution_gradient(primitive_soln_gradient); const real2 scaled_viscosity_coefficient = compute_scaled_viscosity_coefficient(primitive_soln); // Divergence of velocity @@ -204,7 +337,7 @@ ::compute_viscous_stress_tensor ( } // Viscous stress tensor, \tau_{i,j} - std::array,dim> viscous_stress_tensor; + dealii::Tensor<2,dim,real2> viscous_stress_tensor; for (int d1=0; d1,nstate> primitive_soln_gradient = convert_conservative_gradient_to_primitive_gradient(conservative_soln, solution_gradient); // Step 3: Viscous stress tensor, Velocities, Heat flux - const std::array,dim> viscous_stress_tensor = compute_viscous_stress_tensor(primitive_soln, primitive_soln_gradient); + const dealii::Tensor<2,dim,real2> viscous_stress_tensor = compute_viscous_stress_tensor(primitive_soln, primitive_soln_gradient); const dealii::Tensor<1,dim,real2> vel = this->template extract_velocities_from_primitive(primitive_soln); // from Euler const dealii::Tensor<1,dim,real2> heat_flux = compute_heat_flux(primitive_soln, primitive_soln_gradient); @@ -670,10 +803,10 @@ template dealii::Tensor<1,PHILIP_DIM,RadType> NavierStokes < PHILIP_DIM, PHILIP_ template dealii::Tensor<1,PHILIP_DIM,FadFadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, FadFadType>::compute_heat_flux(const std::array &primitive_soln, const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; template dealii::Tensor<1,PHILIP_DIM,RadFadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, RadFadType>::compute_heat_flux(const std::array &primitive_soln, const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; -template std::array,PHILIP_DIM> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, double>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; -template std::array,PHILIP_DIM> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, FadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; -template std::array,PHILIP_DIM> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, RadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; -template std::array,PHILIP_DIM> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, FadFadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; -template std::array,PHILIP_DIM> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, RadFadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; +template dealii::Tensor<2,PHILIP_DIM,double> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, double>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; +template dealii::Tensor<2,PHILIP_DIM,FadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, FadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; +template dealii::Tensor<2,PHILIP_DIM,RadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, RadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; +template dealii::Tensor<2,PHILIP_DIM,FadFadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, FadFadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; +template dealii::Tensor<2,PHILIP_DIM,RadFadType> NavierStokes < PHILIP_DIM, PHILIP_DIM+2, RadFadType>::extract_velocities_gradient_from_primitive_solution_gradient(const std::array,PHILIP_DIM+2> &primitive_soln_gradient) const; } // Physics namespace } // PHiLiP namespace diff --git a/src/physics/navier_stokes.h b/src/physics/navier_stokes.h index d8ec7d842..d022bd457 100644 --- a/src/physics/navier_stokes.h +++ b/src/physics/navier_stokes.h @@ -81,9 +81,60 @@ class NavierStokes : public Euler const std::array &primitive_soln, const std::array,nstate> &primitive_soln_gradient) const; + /// Evaluate vorticity from conservative variables and gradient of conservative variables + dealii::Tensor<1,3,real> compute_vorticity ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const; + + /// Evaluate enstrophy from conservative variables and gradient of conservative variables + real compute_enstrophy ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const; + + /** Evaluate non-dimensional theoretical vorticity-based dissipation rate integrated enstrophy. + * Note: For incompressible flows or when dilatation effects are negligible + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * -- Equation (56) with free-stream nondimensionalization applied + * */ + real compute_vorticity_based_dissipation_rate_from_integrated_enstrophy ( + const real integrated_enstrophy) const; + + /** Evaluate pressure dilatation from conservative variables and gradient of conservative variables + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + real compute_pressure_dilatation ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const; + + /** Evaluate the deviatoric strain-rate tensor from conservative variables and gradient of conservative variables + * -- Reference: Pope, Stephen B., "Turbulent Flows", Cambridge University Press (2000). Eq.(2.70) + * */ + dealii::Tensor<2,dim,real> compute_deviatoric_strain_rate_tensor ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const; + + /// Evaluate the square of the deviatoric strain-rate tensor magnitude (i.e. double dot product) from conservative variables and gradient of conservative variables + real compute_deviatoric_strain_rate_tensor_magnitude_sqr ( + const std::array &conservative_soln, + const std::array,nstate> &conservative_soln_gradient) const; + + /** Evaluate non-dimensional theoretical deviatoric strain-rate tensor based dissipation rate from integrated + * deviatoric strain-rate tensor magnitude squared. + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * -- Equation (57a) with free-stream nondimensionalization applied + * */ + real compute_deviatoric_strain_rate_tensor_based_dissipation_rate_from_integrated_deviatoric_strain_rate_tensor_magnitude_sqr ( + const real integrated_deviatoric_strain_rate_tensor_magnitude_sqr) const; + /** Extract gradient of velocities */ template - std::array,dim> + dealii::Tensor<2,dim,real2> extract_velocities_gradient_from_primitive_solution_gradient ( const std::array,nstate> &primitive_soln_gradient) const; @@ -91,7 +142,7 @@ class NavierStokes : public Euler * Reference: Masatsuka 2018 "I do like CFD", p.148, eq.(4.14.12) */ template - std::array,dim> + dealii::Tensor<2,dim,real2> compute_viscous_stress_tensor ( const std::array &primitive_soln, const std::array,nstate> &primitive_soln_gradient) const; @@ -190,6 +241,10 @@ class NavierStokes : public Euler std::array &soln_bc, std::array,nstate> &soln_grad_bc) const override; +private: + /// Returns the square of the magnitude of the tensor (i.e. the double dot product of a tensor with itself) + real get_tensor_magnitude_sqr (const dealii::Tensor<2,dim,real> &tensor) const; + }; } // Physics namespace diff --git a/src/testing/CMakeLists.txt b/src/testing/CMakeLists.txt index e32d7c510..9d6c47cef 100644 --- a/src/testing/CMakeLists.txt +++ b/src/testing/CMakeLists.txt @@ -24,6 +24,7 @@ set(TESTS_SOURCE fd_state_sensitivity_wrt_parameter.cpp flow_solver_cases/flow_solver_case_base.cpp flow_solver_cases/periodic_cube_flow.cpp + flow_solver_cases/periodic_turbulence.cpp flow_solver_cases/1D_burgers_rewienski_snapshot.cpp flow_solver_cases/1d_burgers_viscous_snapshot.cpp flow_solver.cpp diff --git a/src/testing/flow_solver.cpp b/src/testing/flow_solver.cpp index 73c203bf1..74b410de3 100644 --- a/src/testing/flow_solver.cpp +++ b/src/testing/flow_solver.cpp @@ -21,7 +21,6 @@ FlowSolver::FlowSolver( : TestsBase::TestsBase(parameters_input) , flow_solver_case(flow_solver_case) , parameter_handler(parameter_handler_input) -, initial_condition_function(InitialConditionFactory::create_InitialConditionFunction(parameters_input)) , all_param(*parameters_input) , flow_solver_param(all_param.flow_solver_param) , ode_param(all_param.ode_solver_param) @@ -33,7 +32,7 @@ FlowSolver::FlowSolver( { flow_solver_case->set_higher_order_grid(dg); dg->allocate_system(); - flow_solver_case->display_flow_solver_setup(initial_condition_function); + flow_solver_case->display_flow_solver_setup(); dealii::LinearAlgebra::distributed::Vector solution_no_ghost; solution_no_ghost.reinit(dg->locally_owned_dofs, MPI_COMM_WORLD); @@ -63,7 +62,7 @@ FlowSolver::FlowSolver( } else { // Initialize solution from initial_condition_function pcout << "Initializing solution with initial condition function... " << std::flush; - dealii::VectorTools::interpolate(dg->dof_handler, *initial_condition_function, solution_no_ghost); + dealii::VectorTools::interpolate(dg->dof_handler, *(flow_solver_case->initial_condition_function), solution_no_ghost); } dg->solution = solution_no_ghost; //< assignment dg->solution.update_ghost_values(); diff --git a/src/testing/flow_solver.h b/src/testing/flow_solver.h index ddad12e63..95537dd62 100644 --- a/src/testing/flow_solver.h +++ b/src/testing/flow_solver.h @@ -25,6 +25,7 @@ #include "ode_solver/explicit_ode_solver.h" #include "ode_solver/ode_solver_factory.h" #include "flow_solver_cases/periodic_cube_flow.h" +#include "flow_solver_cases/periodic_turbulence.h" #include "flow_solver_cases/1D_burgers_rewienski_snapshot.h" #include "flow_solver_cases/1d_burgers_viscous_snapshot.h" #include "flow_solver_cases/naca0012.h" @@ -74,7 +75,6 @@ class FlowSolver : public TestsBase std::string get_restart_filename_without_extension(const int restart_index_input) const; protected: - std::shared_ptr> initial_condition_function; ///< Initial condition function const Parameters::AllParameters all_param; ///< All parameters const Parameters::FlowSolverParam flow_solver_param; ///< Flow solver parameters const Parameters::ODESolverParam ode_param; ///< ODE solver parameters diff --git a/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.cpp b/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.cpp index ab0d5bbbf..345c9d05e 100644 --- a/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.cpp +++ b/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.cpp @@ -42,7 +42,7 @@ std::shared_ptr BurgersRewienskiSnapshot::generate_gr } template -void BurgersRewienskiSnapshot::display_additional_flow_case_specific_parameters(std::shared_ptr> /*initial_condition*/) const +void BurgersRewienskiSnapshot::display_additional_flow_case_specific_parameters() const { // Display the information about the grid this->pcout << "\n- GRID INFORMATION:" << std::endl; @@ -57,7 +57,7 @@ void BurgersRewienskiSnapshot::compute_unsteady_data_and_write_to_t const unsigned int current_iteration, const double current_time, const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const + const std::shared_ptr unsteady_data_table) { if (this->all_param.ode_solver_param.output_solution_vector_modulo > 0) { if (current_iteration % this->all_param.ode_solver_param.output_solution_vector_modulo == 0) { diff --git a/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.h b/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.h index e8d4022cb..e3dd14ddd 100644 --- a/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.h +++ b/src/testing/flow_solver_cases/1D_burgers_rewienski_snapshot.h @@ -38,7 +38,7 @@ class BurgersRewienskiSnapshot: public FlowSolverCaseBase const unsigned int current_iteration, const double current_time, const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const override; + const std::shared_ptr unsteady_data_table) override; /// Function for postprocessing when solving for steady state void steady_state_postprocessing(std::shared_ptr > dg) const override; @@ -49,7 +49,7 @@ class BurgersRewienskiSnapshot: public FlowSolverCaseBase const double domain_right; ///< Domain right-boundary value for generating the grid /// Display additional more specific flow case parameters - void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const override; + void display_additional_flow_case_specific_parameters() const override; }; } // Tests namespace diff --git a/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.cpp b/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.cpp index fe1d9bc6d..2d570ee4a 100644 --- a/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.cpp +++ b/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.cpp @@ -41,7 +41,7 @@ std::shared_ptr BurgersViscousSnapshot::generate_grid } template -void BurgersViscousSnapshot::display_additional_flow_case_specific_parameters(std::shared_ptr> /*initial_condition*/) const +void BurgersViscousSnapshot::display_additional_flow_case_specific_parameters() const { // Display the information about the grid this->pcout << "\n- GRID INFORMATION:" << std::endl; @@ -56,7 +56,7 @@ void BurgersViscousSnapshot::compute_unsteady_data_and_write_to_tab const unsigned int current_iteration, const double current_time, const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const + const std::shared_ptr unsteady_data_table) { if (this->all_param.ode_solver_param.output_solution_vector_modulo > 0) { if (current_iteration % this->all_param.ode_solver_param.output_solution_vector_modulo == 0) { diff --git a/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.h b/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.h index 48a200c44..14918df1d 100644 --- a/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.h +++ b/src/testing/flow_solver_cases/1d_burgers_viscous_snapshot.h @@ -38,7 +38,7 @@ class BurgersViscousSnapshot: public FlowSolverCaseBase const unsigned int current_iteration, const double current_time, const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const override; + const std::shared_ptr unsteady_data_table) override; protected: const int number_of_refinements; ///< Number of cells per direction for the grid @@ -46,7 +46,7 @@ class BurgersViscousSnapshot: public FlowSolverCaseBase const double domain_right; ///< Domain right-boundary value for generating the grid /// Display additional more specific flow case parameters - void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const override; + void display_additional_flow_case_specific_parameters() const override; }; } // Tests namespace diff --git a/src/testing/flow_solver_cases/flow_solver_case_base.cpp b/src/testing/flow_solver_cases/flow_solver_case_base.cpp index 5ecf85dbc..0043ad9a3 100644 --- a/src/testing/flow_solver_cases/flow_solver_case_base.cpp +++ b/src/testing/flow_solver_cases/flow_solver_case_base.cpp @@ -5,7 +5,8 @@ namespace Tests { template FlowSolverCaseBase::FlowSolverCaseBase(const PHiLiP::Parameters::AllParameters *const parameters_input) - : all_param(*parameters_input) + : initial_condition_function(InitialConditionFactory::create_InitialConditionFunction(parameters_input)) + , all_param(*parameters_input) , mpi_communicator(MPI_COMM_WORLD) , mpi_rank(dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) , n_mpi(dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) @@ -50,7 +51,7 @@ std::string FlowSolverCaseBase::get_flow_case_string() const } template -void FlowSolverCaseBase::display_flow_solver_setup(std::shared_ptr> initial_condition) const +void FlowSolverCaseBase::display_flow_solver_setup() const { const std::string pde_string = this->get_pde_string(); pcout << "- PDE Type: " << pde_string << " " << "(dim=" << dim << ", nstate=" << nstate << ")" << std::endl; @@ -64,7 +65,7 @@ void FlowSolverCaseBase::display_flow_solver_setup(std::shared_ptrall_param.flow_solver_param.final_time << std::endl; } - this->display_additional_flow_case_specific_parameters(initial_condition); + this->display_additional_flow_case_specific_parameters(); } template @@ -91,11 +92,21 @@ void FlowSolverCaseBase::compute_unsteady_data_and_write_to_table( const unsigned int /*current_iteration*/, const double /*current_time*/, const std::shared_ptr > /*dg*/, - const std::shared_ptr /*unsteady_data_table*/) const + const std::shared_ptr /*unsteady_data_table*/) { // do nothing by default } +template +void FlowSolverCaseBase::add_value_to_data_table( + const double value, + const std::string value_string, + const std::shared_ptr data_table) const +{ + data_table->add_value(value_string, value); + data_table->set_precision(value_string, 16); + data_table->set_scientific(value_string, true); +} #if PHILIP_DIM==1 template class FlowSolverCaseBase; diff --git a/src/testing/flow_solver_cases/flow_solver_case_base.h b/src/testing/flow_solver_cases/flow_solver_case_base.h index 366ad5546..25e29bd37 100644 --- a/src/testing/flow_solver_cases/flow_solver_case_base.h +++ b/src/testing/flow_solver_cases/flow_solver_case_base.h @@ -25,11 +25,13 @@ class FlowSolverCaseBase ///Constructor FlowSolverCaseBase(const Parameters::AllParameters *const parameters_input); + std::shared_ptr> initial_condition_function; ///< Initial condition function + /// Destructor ~FlowSolverCaseBase() {}; /// Displays the flow setup parameters - void display_flow_solver_setup(std::shared_ptr> initial_condition) const; + void display_flow_solver_setup() const; /// Pure Virtual function to generate the grid virtual std::shared_ptr generate_grid() const = 0; @@ -42,7 +44,7 @@ class FlowSolverCaseBase const unsigned int current_iteration, const double current_time, const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const; + const std::shared_ptr unsteady_data_table); /// Virtual function to compute the constant time step virtual double get_constant_time_step(std::shared_ptr > dg) const; @@ -61,8 +63,14 @@ class FlowSolverCaseBase */ dealii::ConditionalOStream pcout; + /// Add a value to a given data table with scientific format + void add_value_to_data_table( + const double value, + const std::string value_string, + const std::shared_ptr data_table) const; + /// Display additional more specific flow case parameters - virtual void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const = 0; + virtual void display_additional_flow_case_specific_parameters() const = 0; private: /// Returns the pde type string from the all_param class member diff --git a/src/testing/flow_solver_cases/naca0012.cpp b/src/testing/flow_solver_cases/naca0012.cpp index 8ababeb68..11e5d807f 100644 --- a/src/testing/flow_solver_cases/naca0012.cpp +++ b/src/testing/flow_solver_cases/naca0012.cpp @@ -28,7 +28,7 @@ NACA0012::NACA0012(const PHiLiP::Parameters::AllParameters *const p } template -void NACA0012::display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const +void NACA0012::display_additional_flow_case_specific_parameters() const { using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation; const PDE_enum pde_type = this->all_param.pde_type; @@ -43,7 +43,7 @@ void NACA0012::display_additional_flow_case_specific_parameters(std: this->pcout << "- - Farfield conditions: " << std::endl; const dealii::Point dummy_point; for (int s=0;spcout << "- - - State " << s << "; Value: " << initial_condition->value(dummy_point, s) << std::endl; + this->pcout << "- - - State " << s << "; Value: " << this->initial_condition_function->value(dummy_point, s) << std::endl; } } diff --git a/src/testing/flow_solver_cases/naca0012.h b/src/testing/flow_solver_cases/naca0012.h index a7db420c2..a21d3b1af 100644 --- a/src/testing/flow_solver_cases/naca0012.h +++ b/src/testing/flow_solver_cases/naca0012.h @@ -40,7 +40,7 @@ class NACA0012 : public FlowSolverCaseBase protected: /// Display additional more specific flow case parameters - void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const override; + void display_additional_flow_case_specific_parameters() const override; }; } // Tests namespace diff --git a/src/testing/flow_solver_cases/periodic_cube_flow.cpp b/src/testing/flow_solver_cases/periodic_cube_flow.cpp index 1492d328b..44ea35355 100644 --- a/src/testing/flow_solver_cases/periodic_cube_flow.cpp +++ b/src/testing/flow_solver_cases/periodic_cube_flow.cpp @@ -1,15 +1,8 @@ #include "periodic_cube_flow.h" -#include + #include #include -#include -#include -#include -#include -#include "physics/physics_factory.h" -#include "dg/dg.h" #include "mesh/grids/straight_periodic_cube.hpp" -#include namespace PHiLiP { @@ -55,164 +48,13 @@ void PeriodicCubeFlow::display_grid_parameters() const } template -void PeriodicCubeFlow::display_additional_flow_case_specific_parameters(std::shared_ptr> /*initial_condition*/) const -{ - - this->display_grid_parameters(); -} - -//========================================================= -// TURBULENCE IN PERIODIC CUBE DOMAIN -//========================================================= -template -PeriodicTurbulence::PeriodicTurbulence(const PHiLiP::Parameters::AllParameters *const parameters_input) - : PeriodicCubeFlow(parameters_input) - , unsteady_data_table_filename_with_extension(this->all_param.flow_solver_param.unsteady_data_table_filename+".txt") -{ - // Get the flow case type - using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType; - const FlowCaseEnum flow_type = parameters_input->flow_solver_param.flow_case_type; - - // Flow case identifiers - this->is_taylor_green_vortex = (flow_type == FlowCaseEnum::taylor_green_vortex); -} - -template -void PeriodicTurbulence::display_additional_flow_case_specific_parameters(std::shared_ptr> /*initial_condition*/) const +void PeriodicCubeFlow::display_additional_flow_case_specific_parameters() const { - this->pcout << "- - Courant-Friedrich-Lewy number: " << this->all_param.flow_solver_param.courant_friedrich_lewy_number << std::endl; - std::string flow_type_string; - if(this->is_taylor_green_vortex) { - this->pcout << "- - Freestream Reynolds number: " << this->all_param.navier_stokes_param.reynolds_number_inf << std::endl; - this->pcout << "- - Freestream Mach number: " << this->all_param.euler_param.mach_inf << std::endl; - } this->display_grid_parameters(); } -template -double PeriodicTurbulence::get_constant_time_step(std::shared_ptr> dg) const -{ - const unsigned int number_of_degrees_of_freedom = dg->dof_handler.n_dofs(); - const double approximate_grid_spacing = (this->domain_right-this->domain_left)/pow(number_of_degrees_of_freedom,(1.0/dim)); - const double constant_time_step = this->all_param.flow_solver_param.courant_friedrich_lewy_number * approximate_grid_spacing; - return constant_time_step; -} - -template -void PeriodicTurbulence::compute_unsteady_data_and_write_to_table( - const unsigned int current_iteration, - const double current_time, - const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const -{ - // Compute kinetic energy - const double kinetic_energy = this->compute_kinetic_energy(*dg); - - if(this->mpi_rank==0) { - // Add time to table - std::string time_string = "time"; - unsteady_data_table->add_value(time_string, current_time); - unsteady_data_table->set_precision(time_string, 16); - unsteady_data_table->set_scientific(time_string, true); - // Add kinetic energy to table - std::string kinetic_energy_string = "kinetic_energy"; - unsteady_data_table->add_value(kinetic_energy_string, kinetic_energy); - unsteady_data_table->set_precision(kinetic_energy_string, 16); - unsteady_data_table->set_scientific(kinetic_energy_string, true); - // Write to file - std::ofstream unsteady_data_table_file(this->unsteady_data_table_filename_with_extension); - unsteady_data_table->write_text(unsteady_data_table_file); - } - // Print to console - this->pcout << " Iter: " << current_iteration - << " Time: " << current_time - << " Energy: " << kinetic_energy - << std::endl; - - // Abort if energy is nan - if(std::isnan(kinetic_energy)) { - this->pcout << " ERROR: Kinetic energy at time " << current_time << " is nan." << std::endl; - this->pcout << " Consider decreasing the time step / CFL number." << std::endl; - std::abort(); - } -} - -template -double PeriodicTurbulence::integrand_kinetic_energy(const std::array &soln_at_q) const -{ - // Description: Returns nondimensional kinetic energy - const double nondimensional_density = soln_at_q[0]; - double dot_product_of_nondimensional_momentum = 0.0; - for (int d=0; ddomain_size; -} - -template -double PeriodicTurbulence::integrand_l2_error_initial_condition(const std::array &soln_at_q, const dealii::Point qpoint) const -{ - // Description: Returns l2 error with the initial condition function - // Purpose: For checking the initialization - std::shared_ptr> initial_condition_function = InitialConditionFactory::create_InitialConditionFunction(&this->all_param); - double integrand_value = 0.0; - for (int istate=0; istatevalue(qpoint, istate); - integrand_value += pow(soln_at_q[istate] - exact_soln_at_q, 2.0); - } - return integrand_value; -} - -template -double PeriodicTurbulence::integrate_over_domain(DGBase &dg,const IntegratedQuantitiesEnum integrated_quantity) const -{ - double integral_value = 0.0; - - // Overintegrate the error to make sure there is not integration error in the error estimate - int overintegrate = 10; - dealii::QGauss quad_extra(dg.max_degree+1+overintegrate); - dealii::FEValues fe_values_extra(*(dg.high_order_grid->mapping_fe_field), dg.fe_collection[dg.max_degree], quad_extra, - dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points); - - const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; - std::array soln_at_q; - - std::vector dofs_indices (fe_values_extra.dofs_per_cell); - for (auto cell : dg.dof_handler.active_cell_iterators()) { - if (!cell->is_locally_owned()) continue; - fe_values_extra.reinit (cell); - cell->get_dof_indices (dofs_indices); - - for (unsigned int iquad=0; iquad qpoint = (fe_values_extra.quadrature_point(iquad)); - - double integrand_value = 0.0; - if(integrated_quantity == IntegratedQuantitiesEnum::kinetic_energy) {integrand_value = integrand_kinetic_energy(soln_at_q);} - if(integrated_quantity == IntegratedQuantitiesEnum::l2_error_initial_condition) {integrand_value = integrand_l2_error_initial_condition(soln_at_q,qpoint);} - - integral_value += integrand_value * fe_values_extra.JxW(iquad); - } - } - const double integral_value_mpi_sum = dealii::Utilities::MPI::sum(integral_value, this->mpi_communicator); - return integral_value_mpi_sum; -} - -template -double PeriodicTurbulence::compute_kinetic_energy(DGBase &dg) const -{ - return integrate_over_domain(dg, IntegratedQuantitiesEnum::kinetic_energy); -} - #if PHILIP_DIM==3 template class PeriodicCubeFlow ; -template class PeriodicTurbulence ; #endif } // Tests namespace diff --git a/src/testing/flow_solver_cases/periodic_cube_flow.h b/src/testing/flow_solver_cases/periodic_cube_flow.h index 7a27d41e2..3ab0fdbcf 100644 --- a/src/testing/flow_solver_cases/periodic_cube_flow.h +++ b/src/testing/flow_solver_cases/periodic_cube_flow.h @@ -1,23 +1,15 @@ #ifndef __PERIODIC_CUBE_FLOW_H__ #define __PERIODIC_CUBE_FLOW_H__ -// for FlowSolver class: -#include "physics/initial_conditions/initial_condition.h" -#include "dg/dg.h" -#include "physics/physics.h" -#include "parameters/all_parameters.h" -#include -#include -#include #include "flow_solver_case_base.h" namespace PHiLiP { namespace Tests { #if PHILIP_DIM==1 - using Triangulation = dealii::Triangulation; +using Triangulation = dealii::Triangulation; #else - using Triangulation = dealii::parallel::distributed::Triangulation; +using Triangulation = dealii::parallel::distributed::Triangulation; #endif template @@ -40,59 +32,12 @@ class PeriodicCubeFlow : public FlowSolverCaseBase const double domain_size; ///< Domain size (length in 1D, area in 2D, and volume in 3D) /// Display additional more specific flow case parameters - virtual void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const override; + virtual void display_additional_flow_case_specific_parameters() const override; /// Display grid parameters void display_grid_parameters() const; }; -template -class PeriodicTurbulence : public PeriodicCubeFlow -{ -public: - /// Constructor. - PeriodicTurbulence(const Parameters::AllParameters *const parameters_input); - - /// Destructor - ~PeriodicTurbulence() {}; - - /// Computes the kinetic energy - double compute_kinetic_energy(DGBase &dg) const; - -protected: - /// Filename (with extension) for the unsteady data table - const std::string unsteady_data_table_filename_with_extension; - - /// Display additional more specific flow case parameters - void display_additional_flow_case_specific_parameters(std::shared_ptr> initial_condition) const override; - - bool is_taylor_green_vortex = false; ///< Identifies if taylor green vortex case; initialized as false. - - /// Function to compute the constant time step - double get_constant_time_step(std::shared_ptr> dg) const override; - - /// Compute the desired unsteady data and write it to a table - void compute_unsteady_data_and_write_to_table( - const unsigned int current_iteration, - const double current_time, - const std::shared_ptr > dg, - const std::shared_ptr unsteady_data_table) const override; - - /// List of possible integrated quantities over the domain - enum IntegratedQuantitiesEnum { - kinetic_energy, - l2_error_initial_condition - }; - - /// Kinetic energy integrand used for integrating over the entire domain - double integrand_kinetic_energy(const std::array &soln_at_q) const; - /// Integrand for computing the L2-error of the initialization with the initial condition - double integrand_l2_error_initial_condition(const std::array &soln_at_q, const dealii::Point qpoint) const; - - /// Integrates over the entire domain - double integrate_over_domain(DGBase &dg,const IntegratedQuantitiesEnum integrated_quantity) const; -}; - } // Tests namespace } // PHiLiP namespace #endif \ No newline at end of file diff --git a/src/testing/flow_solver_cases/periodic_turbulence.cpp b/src/testing/flow_solver_cases/periodic_turbulence.cpp new file mode 100644 index 000000000..fd11db6ab --- /dev/null +++ b/src/testing/flow_solver_cases/periodic_turbulence.cpp @@ -0,0 +1,228 @@ +#include "periodic_turbulence.h" + +#include +#include +#include +#include +#include +#include +#include +#include "physics/physics_factory.h" +#include +#include +#include "math.h" + +namespace PHiLiP { + +namespace Tests { + +//========================================================= +// TURBULENCE IN PERIODIC CUBE DOMAIN +//========================================================= +template +PeriodicTurbulence::PeriodicTurbulence(const PHiLiP::Parameters::AllParameters *const parameters_input) + : PeriodicCubeFlow(parameters_input) + , unsteady_data_table_filename_with_extension(this->all_param.flow_solver_param.unsteady_data_table_filename+".txt") +{ + // Get the flow case type + using FlowCaseEnum = Parameters::FlowSolverParam::FlowCaseType; + const FlowCaseEnum flow_type = this->all_param.flow_solver_param.flow_case_type; + + // Flow case identifiers + this->is_taylor_green_vortex = (flow_type == FlowCaseEnum::taylor_green_vortex); + this->is_viscous_flow = (this->all_param.pde_type == Parameters::AllParameters::PartialDifferentialEquation::navier_stokes); + + // Navier-Stokes object; create using dynamic_pointer_cast and the create_Physics factory + PHiLiP::Parameters::AllParameters parameters_navier_stokes = this->all_param; + parameters_navier_stokes.pde_type = Parameters::AllParameters::PartialDifferentialEquation::navier_stokes; + this->navier_stokes_physics = std::dynamic_pointer_cast>( + Physics::PhysicsFactory::create_Physics(¶meters_navier_stokes)); + + /* Initialize integrated quantities as NAN; + done as a precaution in the case compute_integrated_quantities() is not called + before a member function of kind get_integrated_quantity() is called + */ + std::fill(this->integrated_quantities.begin(), this->integrated_quantities.end(), NAN); +} + +template +void PeriodicTurbulence::display_additional_flow_case_specific_parameters() const +{ + this->pcout << "- - Courant-Friedrich-Lewy number: " << this->all_param.flow_solver_param.courant_friedrich_lewy_number << std::endl; + std::string flow_type_string; + if(this->is_taylor_green_vortex) { + this->pcout << "- - Freestream Reynolds number: " << this->all_param.navier_stokes_param.reynolds_number_inf << std::endl; + this->pcout << "- - Freestream Mach number: " << this->all_param.euler_param.mach_inf << std::endl; + } + this->display_grid_parameters(); +} + +template +double PeriodicTurbulence::get_constant_time_step(std::shared_ptr> dg) const +{ + const unsigned int number_of_degrees_of_freedom = dg->dof_handler.n_dofs(); + const double approximate_grid_spacing = (this->domain_right-this->domain_left)/pow(number_of_degrees_of_freedom,(1.0/dim)); + const double constant_time_step = this->all_param.flow_solver_param.courant_friedrich_lewy_number * approximate_grid_spacing; + return constant_time_step; +} + +template +void PeriodicTurbulence::compute_and_update_integrated_quantities(DGBase &dg) +{ + std::array integral_values; + std::fill(integral_values.begin(), integral_values.end(), 0.0); + + // Overintegrate the error to make sure there is not integration error in the error estimate + int overintegrate = 10; + dealii::QGauss quad_extra(dg.max_degree+1+overintegrate); + dealii::FEValues fe_values_extra(*(dg.high_order_grid->mapping_fe_field), dg.fe_collection[dg.max_degree], quad_extra, + dealii::update_values | dealii::update_gradients | dealii::update_JxW_values | dealii::update_quadrature_points); + + const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points; + std::array soln_at_q; + std::array,nstate> soln_grad_at_q; + + std::vector dofs_indices (fe_values_extra.dofs_per_cell); + for (auto cell : dg.dof_handler.active_cell_iterators()) { + if (!cell->is_locally_owned()) continue; + fe_values_extra.reinit (cell); + cell->get_dof_indices (dofs_indices); + + for (unsigned int iquad=0; iquad qpoint = (fe_values_extra.quadrature_point(iquad)); + + std::array integrand_values; + std::fill(integrand_values.begin(), integrand_values.end(), 0.0); + integrand_values[IntegratedQuantitiesEnum::kinetic_energy] = this->navier_stokes_physics->compute_kinetic_energy_from_conservative_solution(soln_at_q); + integrand_values[IntegratedQuantitiesEnum::enstrophy] = this->navier_stokes_physics->compute_enstrophy(soln_at_q,soln_grad_at_q); + integrand_values[IntegratedQuantitiesEnum::pressure_dilatation] = this->navier_stokes_physics->compute_pressure_dilatation(soln_at_q,soln_grad_at_q); + integrand_values[IntegratedQuantitiesEnum::deviatoric_strain_rate_tensor_magnitude_sqr] = this->navier_stokes_physics->compute_deviatoric_strain_rate_tensor_magnitude_sqr(soln_at_q,soln_grad_at_q); + + for(int i_quantity=0; i_quantityintegrated_quantities[i_quantity] = dealii::Utilities::MPI::sum(integral_values[i_quantity], this->mpi_communicator); + this->integrated_quantities[i_quantity] /= this->domain_size; // divide by total domain volume + } +} + +template +double PeriodicTurbulence::get_integrated_kinetic_energy() const +{ + const double integrated_kinetic_energy = this->integrated_quantities[IntegratedQuantitiesEnum::kinetic_energy]; + // // Abort if energy is nan + // if(std::isnan(integrated_kinetic_energy)) { + // this->pcout << " ERROR: Kinetic energy at time " << current_time << " is nan." << std::endl; + // this->pcout << " Consider decreasing the time step / CFL number." << std::endl; + // std::abort(); + // } + return integrated_kinetic_energy; +} + +template +double PeriodicTurbulence::get_integrated_enstrophy() const +{ + return this->integrated_quantities[IntegratedQuantitiesEnum::enstrophy]; +} + +template +double PeriodicTurbulence::get_vorticity_based_dissipation_rate() const +{ + const double integrated_enstrophy = this->integrated_quantities[IntegratedQuantitiesEnum::enstrophy]; + double vorticity_based_dissipation_rate = 0.0; + if (is_viscous_flow){ + vorticity_based_dissipation_rate = this->navier_stokes_physics->compute_vorticity_based_dissipation_rate_from_integrated_enstrophy(integrated_enstrophy); + } + return vorticity_based_dissipation_rate; +} + +template +double PeriodicTurbulence::get_pressure_dilatation_based_dissipation_rate() const +{ + const double integrated_pressure_dilatation = this->integrated_quantities[IntegratedQuantitiesEnum::pressure_dilatation]; + return (-1.0*integrated_pressure_dilatation); // See reference (listed in header file), equation (57b) +} + +template +double PeriodicTurbulence::get_deviatoric_strain_rate_tensor_based_dissipation_rate() const +{ + const double integrated_deviatoric_strain_rate_tensor_magnitude_sqr = this->integrated_quantities[IntegratedQuantitiesEnum::deviatoric_strain_rate_tensor_magnitude_sqr]; + double deviatoric_strain_rate_tensor_based_dissipation_rate = 0.0; + if (is_viscous_flow){ + deviatoric_strain_rate_tensor_based_dissipation_rate = + this->navier_stokes_physics->compute_deviatoric_strain_rate_tensor_based_dissipation_rate_from_integrated_deviatoric_strain_rate_tensor_magnitude_sqr(integrated_deviatoric_strain_rate_tensor_magnitude_sqr); + } + return deviatoric_strain_rate_tensor_based_dissipation_rate; +} + +template +void PeriodicTurbulence::compute_unsteady_data_and_write_to_table( + const unsigned int current_iteration, + const double current_time, + const std::shared_ptr > dg, + const std::shared_ptr unsteady_data_table) +{ + // Compute and update integrated quantities + this->compute_and_update_integrated_quantities(*dg); + // Get computed quantities + const double integrated_kinetic_energy = this->get_integrated_kinetic_energy(); + const double integrated_enstrophy = this->get_integrated_enstrophy(); + const double vorticity_based_dissipation_rate = this->get_vorticity_based_dissipation_rate(); + const double pressure_dilatation_based_dissipation_rate = this->get_pressure_dilatation_based_dissipation_rate(); + const double deviatoric_strain_rate_tensor_based_dissipation_rate = this->get_deviatoric_strain_rate_tensor_based_dissipation_rate(); + + if(this->mpi_rank==0) { + // Add values to data table + this->add_value_to_data_table(current_time,"time",unsteady_data_table); + this->add_value_to_data_table(integrated_kinetic_energy,"kinetic_energy",unsteady_data_table); + this->add_value_to_data_table(integrated_enstrophy,"enstrophy",unsteady_data_table); + if(is_viscous_flow) this->add_value_to_data_table(vorticity_based_dissipation_rate,"eps_vorticity",unsteady_data_table); + this->add_value_to_data_table(pressure_dilatation_based_dissipation_rate,"eps_pressure",unsteady_data_table); + if(is_viscous_flow) this->add_value_to_data_table(deviatoric_strain_rate_tensor_based_dissipation_rate,"eps_strain",unsteady_data_table); + // Write to file + std::ofstream unsteady_data_table_file(this->unsteady_data_table_filename_with_extension); + unsteady_data_table->write_text(unsteady_data_table_file); + } + // Print to console + this->pcout << " Iter: " << current_iteration + << " Time: " << current_time + << " Energy: " << integrated_kinetic_energy + << " Enstrophy: " << integrated_enstrophy; + if(is_viscous_flow) { + this->pcout << " eps_vorticity: " << vorticity_based_dissipation_rate + << " eps_p+eps_strain: " << (pressure_dilatation_based_dissipation_rate + deviatoric_strain_rate_tensor_based_dissipation_rate); + } + this->pcout << std::endl; + + // Abort if energy is nan + if(std::isnan(integrated_kinetic_energy)) { + this->pcout << " ERROR: Kinetic energy at time " << current_time << " is nan." << std::endl; + this->pcout << " Consider decreasing the time step / CFL number." << std::endl; + std::abort(); + } +} + +#if PHILIP_DIM==3 +template class PeriodicTurbulence ; +#endif + +} // Tests namespace +} // PHiLiP namespace + diff --git a/src/testing/flow_solver_cases/periodic_turbulence.h b/src/testing/flow_solver_cases/periodic_turbulence.h new file mode 100644 index 000000000..b7422e7e3 --- /dev/null +++ b/src/testing/flow_solver_cases/periodic_turbulence.h @@ -0,0 +1,101 @@ +#ifndef __PERIODIC_TURBULENCE_H__ +#define __PERIODIC_TURBULENCE_H__ + +#include "periodic_cube_flow.h" +#include "dg/dg.h" +#include "physics/navier_stokes.h" + +namespace PHiLiP { +namespace Tests { + +template +class PeriodicTurbulence : public PeriodicCubeFlow +{ + /** Number of different computed quantities + * Corresponds to the number of items in IntegratedQuantitiesEnum + * */ + static const int NUMBER_OF_INTEGRATED_QUANTITIES = 4; + +public: + /// Constructor. + PeriodicTurbulence(const Parameters::AllParameters *const parameters_input); + + /// Destructor + ~PeriodicTurbulence() {}; + + /// Computes the integrated quantities over the domain simultaneously and updates the array storing them + void compute_and_update_integrated_quantities(DGBase &dg); + + /** Gets the nondimensional integrated kinetic energy given a DG object from dg->solution + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + double get_integrated_kinetic_energy() const; + + /** Gets the nondimensional integrated enstrophy given a DG object from dg->solution + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + double get_integrated_enstrophy() const; + + /** Gets non-dimensional theoretical vorticity tensor based dissipation rate + * Note: For incompressible flows or when dilatation effects are negligible + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + double get_vorticity_based_dissipation_rate() const; + + /** Evaluate non-dimensional theoretical pressure-dilatation dissipation rate + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + double get_pressure_dilatation_based_dissipation_rate () const; + + /** Gets non-dimensional theoretical deviatoric strain-rate tensor based dissipation rate + * -- Reference: Cox, Christopher, et al. "Accuracy, stability, and performance comparison + * between the spectral difference and flux reconstruction schemes." + * Computers & Fluids 221 (2021): 104922. + * */ + double get_deviatoric_strain_rate_tensor_based_dissipation_rate() const; + +protected: + /// Filename (with extension) for the unsteady data table + const std::string unsteady_data_table_filename_with_extension; + + /// Pointer to Navier-Stokes physics object for computing things on the fly + std::shared_ptr< Physics::NavierStokes > navier_stokes_physics; + + bool is_taylor_green_vortex = false; ///< Identifies if taylor green vortex case; initialized as false. + bool is_viscous_flow = true; ///< Identifies if viscous flow; initialized as true. + + /// Display additional more specific flow case parameters + void display_additional_flow_case_specific_parameters() const override; + + /// Function to compute the constant time step + double get_constant_time_step(std::shared_ptr> dg) const override; + + /// Compute the desired unsteady data and write it to a table + void compute_unsteady_data_and_write_to_table( + const unsigned int current_iteration, + const double current_time, + const std::shared_ptr > dg, + const std::shared_ptr unsteady_data_table) override; + + /// List of possible integrated quantities over the domain + enum IntegratedQuantitiesEnum { + kinetic_energy, + enstrophy, + pressure_dilatation, + deviatoric_strain_rate_tensor_magnitude_sqr + }; + /// Array for storing the integrated quantities; done for computational efficiency + std::array integrated_quantities; +}; + +} // Tests namespace +} // PHiLiP namespace +#endif \ No newline at end of file diff --git a/src/testing/taylor_green_vortex_energy_check.cpp b/src/testing/taylor_green_vortex_energy_check.cpp index 9ee22c73a..a83943a14 100644 --- a/src/testing/taylor_green_vortex_energy_check.cpp +++ b/src/testing/taylor_green_vortex_energy_check.cpp @@ -1,6 +1,6 @@ #include "taylor_green_vortex_energy_check.h" #include "flow_solver.h" -#include "flow_solver_cases/periodic_cube_flow.h" +#include "flow_solver_cases/periodic_turbulence.h" namespace PHiLiP { namespace Tests { @@ -23,7 +23,8 @@ int TaylorGreenVortexEnergyCheck::run_test() const // Compute kinetic energy std::unique_ptr> flow_solver_case = std::make_unique>(this->all_parameters); - const double kinetic_energy_computed = flow_solver_case->compute_kinetic_energy(*(flow_solver->dg)); + flow_solver_case->compute_and_update_integrated_quantities(*(flow_solver->dg)); + const double kinetic_energy_computed = flow_solver_case->get_integrated_kinetic_energy(); const double relative_error = abs(kinetic_energy_computed - kinetic_energy_expected)/kinetic_energy_expected; if (relative_error > 1.0e-10) { diff --git a/src/testing/taylor_green_vortex_restart_check.cpp b/src/testing/taylor_green_vortex_restart_check.cpp index d20b34ea0..b941d7485 100644 --- a/src/testing/taylor_green_vortex_restart_check.cpp +++ b/src/testing/taylor_green_vortex_restart_check.cpp @@ -1,6 +1,6 @@ #include "taylor_green_vortex_restart_check.h" #include "flow_solver.h" -#include "flow_solver_cases/periodic_cube_flow.h" +#include "flow_solver_cases/periodic_turbulence.h" #include #include #include @@ -111,7 +111,8 @@ int TaylorGreenVortexRestartCheck::run_test() const // Compute kinetic energy at final time achieved by restarting the computation std::unique_ptr> flow_solver_case = std::make_unique>(this->all_parameters); - const double kinetic_energy_computed = flow_solver_case->compute_kinetic_energy(*(flow_solver_restart_to_complete_run->dg)); + flow_solver_case->compute_and_update_integrated_quantities(*(flow_solver_restart_to_complete_run->dg)); + const double kinetic_energy_computed = flow_solver_case->get_integrated_kinetic_energy(); const double relative_error = abs(kinetic_energy_computed - kinetic_energy_expected)/kinetic_energy_expected; if (relative_error > 1.0e-10) { diff --git a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_long.prm b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_long.prm index acc4e3d11..60759bc72 100644 --- a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_long.prm +++ b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_long.prm @@ -53,7 +53,7 @@ subsection flow_solver set final_time = 4.0005346370569432e+00 set courant_friedrich_lewy_number = 0.005 set unsteady_data_table_filename = tgv_kinetic_energy_vs_time_table_for_energy_check - subsection taylor_green_vortex_energy_check + subsection taylor_green_vortex set expected_kinetic_energy_at_final_time = 1.0848940550514036e-01 end end diff --git a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_quick.prm b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_quick.prm index 90ce06d6f..5624227cf 100644 --- a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_quick.prm +++ b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_energy_check_quick.prm @@ -53,7 +53,7 @@ subsection flow_solver set final_time = 1.2248096860487504e-02 set courant_friedrich_lewy_number = 0.005 set unsteady_data_table_filename = tgv_kinetic_energy_vs_time_table_for_energy_check - subsection taylor_green_vortex_energy_check + subsection taylor_green_vortex set expected_kinetic_energy_at_final_time = 1.2074066465958100e-01 end end diff --git a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_check.prm b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_check.prm index b50edf00c..b5867036b 100644 --- a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_check.prm +++ b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_check.prm @@ -53,7 +53,7 @@ subsection flow_solver set final_time = 1.2248096860487504e-02 set courant_friedrich_lewy_number = 0.005 set unsteady_data_table_filename = tgv_kinetic_energy_vs_time_table_for_restart_check - subsection taylor_green_vortex_energy_check + subsection taylor_green_vortex set expected_kinetic_energy_at_final_time = 1.2074066465958100e-01 end end diff --git a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_from_parameter_file_check.prm b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_from_parameter_file_check.prm index e778f953f..843551405 100644 --- a/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_from_parameter_file_check.prm +++ b/tests/integration_tests_control_files/taylor_green_vortex_integration/viscous_taylor_green_vortex_restart_from_parameter_file_check.prm @@ -231,7 +231,7 @@ subsection flow_solver set unsteady_data_table_filename = tgv_kinetic_energy_vs_time_table_for_restart_check # default: unsteady_data_table - subsection taylor_green_vortex_energy_check + subsection taylor_green_vortex # For integration test purposes, expected kinetic energy at final time. set expected_kinetic_energy_at_final_time = 1.2074066465958100e-01 # default: 1 end