Skip to content

Commit b13f5ff

Browse files
jbrillonjbrillon
authored
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
1 parent 41feecb commit b13f5ff

32 files changed

+752
-292
lines changed

‎src/dg/dg.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1736,7 +1736,7 @@ void DGBase<dim,real,MeshType>::output_face_results_vtk (const unsigned int cycl
17361736
data_out.set_flags(vtkflags);
17371737

17381738
const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
1739-
std::string filename = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1739+
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)+"-";
17401740
filename += dealii::Utilities::int_to_string(cycle, 4) + ".";
17411741
filename += dealii::Utilities::int_to_string(iproc, 4);
17421742
filename += ".vtu";
@@ -1747,13 +1747,13 @@ void DGBase<dim,real,MeshType>::output_face_results_vtk (const unsigned int cycl
17471747
if (iproc == 0) {
17481748
std::vector<std::string> filenames;
17491749
for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) {
1750-
std::string fn = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1750+
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)+"-";
17511751
fn += dealii::Utilities::int_to_string(cycle, 4) + ".";
17521752
fn += dealii::Utilities::int_to_string(iproc, 4);
17531753
fn += ".vtu";
17541754
filenames.push_back(fn);
17551755
}
1756-
std::string master_fn = "surface_solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1756+
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)+"-";
17571757
master_fn += dealii::Utilities::int_to_string(cycle, 4) + ".pvtu";
17581758
std::ofstream master_output(master_fn);
17591759
data_out.write_pvtu_record(master_output, filenames);
@@ -1853,7 +1853,7 @@ void DGBase<dim,real,MeshType>::output_results_vtk (const unsigned int cycle)//
18531853
data_out.set_flags(vtkflags);
18541854

18551855
const int iproc = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
1856-
std::string filename = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1856+
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)+"-";
18571857
filename += dealii::Utilities::int_to_string(cycle, 4) + ".";
18581858
filename += dealii::Utilities::int_to_string(iproc, 4);
18591859
filename += ".vtu";
@@ -1864,13 +1864,13 @@ void DGBase<dim,real,MeshType>::output_results_vtk (const unsigned int cycle)//
18641864
if (iproc == 0) {
18651865
std::vector<std::string> filenames;
18661866
for (unsigned int iproc = 0; iproc < dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); ++iproc) {
1867-
std::string fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1867+
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)+"-";
18681868
fn += dealii::Utilities::int_to_string(cycle, 4) + ".";
18691869
fn += dealii::Utilities::int_to_string(iproc, 4);
18701870
fn += ".vtu";
18711871
filenames.push_back(fn);
18721872
}
1873-
std::string master_fn = "solution-" + dealii::Utilities::int_to_string(dim, 1) +"D_maxpoly"+dealii::Utilities::int_to_string(max_degree, 2)+"-";
1873+
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)+"-";
18741874
master_fn += dealii::Utilities::int_to_string(cycle, 4) + ".pvtu";
18751875
std::ofstream master_output(master_fn);
18761876
data_out.write_pvtu_record(master_output, filenames);

‎src/parameters/all_parameters.cpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -180,6 +180,10 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
180180
"Dissipative numerical flux. "
181181
"Choices are <symm_internal_penalty | bassi_rebay_2>.");
182182

183+
prm.declare_entry("solution_vtk_files_directory_name", ".",
184+
dealii::Patterns::FileName(dealii::Patterns::FileName::FileType::input),
185+
"Name of directory for writing solution vtk files. Current directory by default.");
186+
183187
Parameters::LinearSolverParam::declare_parameters (prm);
184188
Parameters::ManufacturedConvergenceStudyParam::declare_parameters (prm);
185189
Parameters::ODESolverParam::declare_parameters (prm);
@@ -314,6 +318,7 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm)
314318
if (flux_reconstruction_aux_string == "kPlus") flux_reconstruction_aux_type = kPlus;
315319
if (flux_reconstruction_aux_string == "k10Thousand") flux_reconstruction_aux_type = k10Thousand;
316320

321+
solution_vtk_files_directory_name = prm.get("solution_vtk_files_directory_name");
317322

318323
pcout << "Parsing linear solver subsection..." << std::endl;
319324
linear_solver_param.parse_parameters (prm);

‎src/parameters/all_parameters.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -192,6 +192,9 @@ class AllParameters
192192
/// Store flux reconstruction type
193193
Flux_Reconstruction_Aux flux_reconstruction_aux_type;
194194

195+
/// Name of directory for writing solution vtk files
196+
std::string solution_vtk_files_directory_name;
197+
195198
/// Declare parameters that can be set as inputs and set up the default options
196199
/** This subroutine should call the sub-parameter classes static declare_parameters()
197200
* such that each sub-parameter class is responsible to declare their own parameters.

‎src/parameters/parameters_flow_solver.cpp‎

Lines changed: 24 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -74,11 +74,24 @@ void FlowSolverParam::declare_parameters(dealii::ParameterHandler &prm)
7474
dealii::Patterns::FileName(dealii::Patterns::FileName::FileType::input),
7575
"Filename of the input mesh: input_mesh_filename.msh");
7676

77-
prm.enter_subsection("taylor_green_vortex_energy_check");
77+
prm.enter_subsection("taylor_green_vortex");
7878
{
7979
prm.declare_entry("expected_kinetic_energy_at_final_time", "1",
8080
dealii::Patterns::Double(0, dealii::Patterns::Double::max_double_value),
8181
"For integration test purposes, expected kinetic energy at final time.");
82+
83+
prm.declare_entry("expected_theoretical_dissipation_rate_at_final_time", "1",
84+
dealii::Patterns::Double(0, dealii::Patterns::Double::max_double_value),
85+
"For integration test purposes, expected theoretical kinetic energy dissipation rate at final time.");
86+
87+
prm.declare_entry("density_initial_condition_type", "uniform",
88+
dealii::Patterns::Selection(
89+
" uniform | "
90+
" isothermal "),
91+
"The type of density initialization. "
92+
"Choices are "
93+
" <uniform | "
94+
" isothermal>.");
8295
}
8396
prm.leave_subsection();
8497
}
@@ -90,10 +103,10 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm)
90103
prm.enter_subsection("flow_solver");
91104
{
92105
const std::string flow_case_type_string = prm.get("flow_case_type");
93-
if (flow_case_type_string == "taylor_green_vortex") {flow_case_type = taylor_green_vortex;}
94-
else if (flow_case_type_string == "burgers_viscous_snapshot") {flow_case_type = burgers_viscous_snapshot;}
95-
else if (flow_case_type_string == "burgers_rewienski_snapshot") {flow_case_type = burgers_rewienski_snapshot;}
96-
else if (flow_case_type_string == "naca0012") {flow_case_type = naca0012;}
106+
if (flow_case_type_string == "taylor_green_vortex") {flow_case_type = taylor_green_vortex;}
107+
else if (flow_case_type_string == "burgers_viscous_snapshot") {flow_case_type = burgers_viscous_snapshot;}
108+
else if (flow_case_type_string == "burgers_rewienski_snapshot") {flow_case_type = burgers_rewienski_snapshot;}
109+
else if (flow_case_type_string == "naca0012") {flow_case_type = naca0012;}
97110

98111
final_time = prm.get_double("final_time");
99112
courant_friedrich_lewy_number = prm.get_double("courant_friedrich_lewy_number");
@@ -108,9 +121,14 @@ void FlowSolverParam::parse_parameters(dealii::ParameterHandler &prm)
108121
output_restart_files_every_dt_time_intervals = prm.get_double("output_restart_files_every_dt_time_intervals");
109122
input_mesh_filename = prm.get("input_mesh_filename");
110123

111-
prm.enter_subsection("taylor_green_vortex_energy_check");
124+
prm.enter_subsection("taylor_green_vortex");
112125
{
113126
expected_kinetic_energy_at_final_time = prm.get_double("expected_kinetic_energy_at_final_time");
127+
expected_theoretical_dissipation_rate_at_final_time = prm.get_double("expected_theoretical_dissipation_rate_at_final_time");
128+
129+
const std::string density_initial_condition_type_string = prm.get("density_initial_condition_type");
130+
if (density_initial_condition_type_string == "uniform") {density_initial_condition_type = uniform;}
131+
else if (density_initial_condition_type_string == "isothermal") {density_initial_condition_type = isothermal;}
114132
}
115133
prm.leave_subsection();
116134
}

‎src/parameters/parameters_flow_solver.h‎

Lines changed: 16 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -40,16 +40,29 @@ class FlowSolverParam
4040
* will read file: input_mesh_filename.msh */
4141
std::string input_mesh_filename;
4242

43-
/** For integration test purposes, expected kinetic energy at final time. */
44-
double expected_kinetic_energy_at_final_time;
45-
4643
bool restart_computation_from_file; ///< Restart computation from restart file
4744
bool output_restart_files; ///< Output the restart files
4845
std::string restart_files_directory_name; ///< Name of directory for writing and reading restart files
4946
int restart_file_index; ///< Index of desired restart file for restarting the computation from
5047
int output_restart_files_every_x_steps; ///< Outputs the restart files every x steps
5148
double output_restart_files_every_dt_time_intervals; ///< Outputs the restart files at time intervals of dt
5249

50+
/** For taylor green vortex integration tests, expected kinetic energy at final time. */
51+
double expected_kinetic_energy_at_final_time;
52+
53+
/** For taylor green vortex integration tests,
54+
* expected theoretical kinetic energy dissipation
55+
* rate at final time. */
56+
double expected_theoretical_dissipation_rate_at_final_time;
57+
58+
/// For taylor green vortex, selects the type of density initialization
59+
enum DensityInitialConditionType{
60+
uniform,
61+
isothermal,
62+
};
63+
/// Selected DensityInitialConditionType from the input file
64+
DensityInitialConditionType density_initial_condition_type;
65+
5366
/// Declares the possible variables and sets the defaults.
5467
static void declare_parameters (dealii::ParameterHandler &prm);
5568

‎src/physics/euler.cpp‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -204,13 +204,30 @@ template <int dim, int nstate, typename real>
204204
inline real Euler<dim,nstate,real>
205205
::compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const
206206
{
207-
const real density = primitive_soln[0];
208207
const real pressure = primitive_soln[nstate-1];
208+
const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln);
209+
const real tot_energy = pressure / this->gamm1 + kinetic_energy;
210+
return tot_energy;
211+
}
212+
213+
template <int dim, int nstate, typename real>
214+
inline real Euler<dim,nstate,real>
215+
::compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const
216+
{
217+
const real density = primitive_soln[0];
209218
const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive<real>(primitive_soln);
210219
const real vel2 = compute_velocity_squared<real>(velocities);
220+
const real kinetic_energy = 0.5*density*vel2;
221+
return kinetic_energy;
222+
}
211223

212-
const real tot_energy = pressure / gamm1 + 0.5*density*vel2;
213-
return tot_energy;
224+
template <int dim, int nstate, typename real>
225+
inline real Euler<dim,nstate,real>
226+
::compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const
227+
{
228+
const std::array<real,nstate> primitive_soln = convert_conservative_to_primitive<real>(conservative_soln);
229+
const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln);
230+
return kinetic_energy;
214231
}
215232

216233
template <int dim, int nstate, typename real>

‎src/physics/euler.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,12 @@ class Euler : public PhysicsBase <dim, nstate, real>
207207
*/
208208
real compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const;
209209

210+
/// Given primitive variables, returns kinetic energy
211+
real compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const;
212+
213+
/// Given conservative variables, returns kinetic energy
214+
real compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const;
215+
210216
/// Evaluate entropy from conservative variables
211217
/** Note that it is not the actual entropy since it's missing some constants.
212218
* Used to check entropy convergence

‎src/physics/initial_conditions/initial_condition.cpp‎

Lines changed: 51 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33

44
namespace PHiLiP {
55
// ========================================================
6-
// TAYLOR GREEN VORTEX -- Initial Condition
6+
// TAYLOR GREEN VORTEX -- Initial Condition (Uniform density)
77
// ========================================================
88
template <int dim, int nstate, typename real>
99
InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
@@ -14,9 +14,7 @@ ::InitialConditionFunction_TaylorGreenVortex (
1414
, gamma_gas(gamma_gas)
1515
, mach_inf(mach_inf)
1616
, mach_inf_sqr(mach_inf*mach_inf)
17-
{
18-
// Nothing to do here yet
19-
}
17+
{}
2018

2119
template <int dim, int nstate, typename real>
2220
real InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
@@ -29,7 +27,7 @@ ::primitive_value(const dealii::Point<dim,real> &point, const unsigned int istat
2927

3028
if(istate==0) {
3129
// density
32-
value = 1.0;
30+
value = this->density(point);
3331
}
3432
if(istate==1) {
3533
// x-velocity
@@ -45,7 +43,7 @@ ::primitive_value(const dealii::Point<dim,real> &point, const unsigned int istat
4543
}
4644
if(istate==4) {
4745
// pressure
48-
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);
46+
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);
4947
}
5048
}
5149
return value;
@@ -69,7 +67,7 @@ ::convert_primitive_to_conversative_value(
6967
if(istate==1) value = rho*u; // x-momentum
7068
if(istate==2) value = rho*v; // y-momentum
7169
if(istate==3) value = rho*w; // z-momentum
72-
if(istate==4) value = p/(1.4-1.0) + 0.5*rho*(u*u + v*v + w*w); // total energy
70+
if(istate==4) value = p/(this->gamma_gas-1.0) + 0.5*rho*(u*u + v*v + w*w); // total energy
7371
}
7472

7573
return value;
@@ -84,6 +82,40 @@ ::value(const dealii::Point<dim,real> &point, const unsigned int istate) const
8482
return value;
8583
}
8684

85+
template <int dim, int nstate, typename real>
86+
real InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
87+
::density(const dealii::Point<dim,real> &/*point*/) const
88+
{
89+
// Note: This is in non-dimensional form (free-stream values as reference)
90+
real value = 0.;
91+
// density
92+
value = 1.0;
93+
return value;
94+
}
95+
96+
// ========================================================
97+
// TAYLOR GREEN VORTEX -- Initial Condition (Isothermal density)
98+
// ========================================================
99+
template <int dim, int nstate, typename real>
100+
InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real>
101+
::InitialConditionFunction_TaylorGreenVortex_Isothermal (
102+
const double gamma_gas,
103+
const double mach_inf)
104+
: InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>(gamma_gas,mach_inf)
105+
{}
106+
107+
template <int dim, int nstate, typename real>
108+
real InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real>
109+
::density(const dealii::Point<dim,real> &point) const
110+
{
111+
// Note: This is in non-dimensional form (free-stream values as reference)
112+
real value = 0.;
113+
// density
114+
value = this->primitive_value(point, 4); // get pressure
115+
value *= this->gamma_gas*this->mach_inf_sqr;
116+
return value;
117+
}
118+
87119
// ========================================================
88120
// 1D BURGERS REWIENSKI -- Initial Condition
89121
// ========================================================
@@ -145,9 +177,19 @@ InitialConditionFactory<dim,nstate, real>::create_InitialConditionFunction(
145177
// Get the flow case type
146178
const FlowCaseEnum flow_type = param->flow_solver_param.flow_case_type;
147179
if (flow_type == FlowCaseEnum::taylor_green_vortex) {
148-
if constexpr (dim==3 && nstate==dim+2) return std::make_shared<InitialConditionFunction_TaylorGreenVortex<dim,nstate,real> >(
180+
if constexpr (dim==3 && nstate==dim+2){
181+
// Get the density initial condition type
182+
const DensityInitialConditionEnum density_initial_condition_type = param->flow_solver_param.density_initial_condition_type;
183+
if(density_initial_condition_type == DensityInitialConditionEnum::uniform) {
184+
return std::make_shared<InitialConditionFunction_TaylorGreenVortex<dim,nstate,real> >(
149185
param->euler_param.gamma_gas,
150186
param->euler_param.mach_inf);
187+
} else if (density_initial_condition_type == DensityInitialConditionEnum::isothermal) {
188+
return std::make_shared<InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real> >(
189+
param->euler_param.gamma_gas,
190+
param->euler_param.mach_inf);
191+
}
192+
}
151193
} else if (flow_type == FlowCaseEnum::burgers_rewienski_snapshot) {
152194
if constexpr (dim==1 && nstate==dim) return std::make_shared<InitialConditionFunction_BurgersRewienski<dim,nstate,real> > ();
153195
} else if (flow_type == FlowCaseEnum::burgers_viscous_snapshot) {
@@ -185,6 +227,7 @@ template class InitialConditionFactory <PHILIP_DIM, PHILIP_DIM, double>;
185227
template class InitialConditionFunction_BurgersViscous<PHILIP_DIM, PHILIP_DIM, double>;
186228
template class InitialConditionFunction_BurgersRewienski<PHILIP_DIM, PHILIP_DIM, double>;
187229
template class InitialConditionFunction_TaylorGreenVortex <PHILIP_DIM,PHILIP_DIM+2,double>;
230+
template class InitialConditionFunction_TaylorGreenVortex_Isothermal <PHILIP_DIM,PHILIP_DIM+2,double>;
188231
template class InitialConditionFunction_Zero <PHILIP_DIM,1, double>;
189232
template class InitialConditionFunction_Zero <PHILIP_DIM,2, double>;
190233
template class InitialConditionFunction_Zero <PHILIP_DIM,3, double>;

0 commit comments

Comments
 (0)