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

@@ -1736,7 +1736,7 @@ void DGBase<dim,real,MeshType>::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<dim,real,MeshType>::output_face_results_vtk (const unsigned int cycl
     if (iproc == 0) {
         std::vector<std::string> 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<dim,real,MeshType>::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<dim,real,MeshType>::output_results_vtk (const unsigned int cycle)//
     if (iproc == 0) {
         std::vector<std::string> 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);

src/parameters/all_parameters.cpp

@@ -180,6 +180,10 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
                       "Dissipative numerical flux. "
                       "Choices are <symm_internal_penalty | bassi_rebay_2>.");
 
+    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);

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.

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 "
+                              " <uniform | "
+                              " isothermal>.");
         }
         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();
     }

src/parameters/parameters_flow_solver.h

@@ -40,16 +40,29 @@ 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
     int restart_file_index; ///< Index of desired restart file for restarting the computation from
     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);
 

src/physics/euler.cpp

@@ -204,13 +204,30 @@ template <int dim, int nstate, typename real>
 inline real Euler<dim,nstate,real>
 ::compute_total_energy ( const std::array<real,nstate> &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 <int dim, int nstate, typename real>
+inline real Euler<dim,nstate,real>
+::compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const
+{
+    const real density = primitive_soln[0];
     const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive<real>(primitive_soln);
     const real vel2 = compute_velocity_squared<real>(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 <int dim, int nstate, typename real>
+inline real Euler<dim,nstate,real>
+::compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &conservative_soln ) const
+{
+    const std::array<real,nstate> primitive_soln = convert_conservative_to_primitive<real>(conservative_soln);
+    const real kinetic_energy = compute_kinetic_energy_from_primitive_solution(primitive_soln);
+    return kinetic_energy;
 }
 
 template <int dim, int nstate, typename real>

src/physics/euler.h

@@ -207,6 +207,12 @@ class Euler : public PhysicsBase <dim, nstate, real>
      */
     real compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const;
 
+    /// Given primitive variables, returns kinetic energy
+    real compute_kinetic_energy_from_primitive_solution ( const std::array<real,nstate> &primitive_soln ) const;
+
+    /// Given conservative variables, returns kinetic energy
+    real compute_kinetic_energy_from_conservative_solution ( const std::array<real,nstate> &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

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 <int dim, int nstate, typename real>
 InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
@@ -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 <int dim, int nstate, typename real>
 real InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
@@ -29,7 +27,7 @@ ::primitive_value(const dealii::Point<dim,real> &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<dim,real> &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<dim,real> &point, const unsigned int istate) const
     return value;
 }
 
+template <int dim, int nstate, typename real>
+real InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>
+::density(const dealii::Point<dim,real> &/*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 <int dim, int nstate, typename real>
+InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real>
+::InitialConditionFunction_TaylorGreenVortex_Isothermal (
+    const double       gamma_gas,
+    const double       mach_inf)
+    : InitialConditionFunction_TaylorGreenVortex<dim,nstate,real>(gamma_gas,mach_inf)
+{}
+
+template <int dim, int nstate, typename real>
+real InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real>
+::density(const dealii::Point<dim,real> &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<dim,nstate, real>::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<InitialConditionFunction_TaylorGreenVortex<dim,nstate,real> >(
+        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<InitialConditionFunction_TaylorGreenVortex<dim,nstate,real> >(
                     param->euler_param.gamma_gas,
                     param->euler_param.mach_inf);
+            } else if (density_initial_condition_type == DensityInitialConditionEnum::isothermal) {
+                return std::make_shared<InitialConditionFunction_TaylorGreenVortex_Isothermal<dim,nstate,real> >(
+                    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<InitialConditionFunction_BurgersRewienski<dim,nstate,real> > ();
     } else if (flow_type == FlowCaseEnum::burgers_viscous_snapshot) {
@@ -185,6 +227,7 @@ template class InitialConditionFactory <PHILIP_DIM, PHILIP_DIM, double>;
 template class InitialConditionFunction_BurgersViscous<PHILIP_DIM, PHILIP_DIM, double>;
 template class InitialConditionFunction_BurgersRewienski<PHILIP_DIM, PHILIP_DIM, double>;
 template class InitialConditionFunction_TaylorGreenVortex <PHILIP_DIM,PHILIP_DIM+2,double>;
+template class InitialConditionFunction_TaylorGreenVortex_Isothermal <PHILIP_DIM,PHILIP_DIM+2,double>;
 template class InitialConditionFunction_Zero <PHILIP_DIM,1, double>;
 template class InitialConditionFunction_Zero <PHILIP_DIM,2, double>;
 template class InitialConditionFunction_Zero <PHILIP_DIM,3, double>;