Viscous terms from entropy gradients (#1202)

* first draft of entropy stable gradients

* update TODOs

* added elixir with discontinuous solution for ES testing

* forgot one TODO marker

* adding GradientVariables**** type parameter

* fix bug (missing type param in constructor)

* add elixir with periodic CNS manufactured solution

* fix gradient conversion

* fix entropy variables w/BCs

* updating elixirs with `gradient_variables=GradientVariablesPrimitive()`

* Update src/equations/compressible_navier_stokes_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* removing unused routines

* removing specialized entropy2cons and cons2entropy routines

* adding more tests

* removing unused elixirs

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* adding docs/warning for gradient variable type

* extracting w_inner[4] into a more clearly named variable

Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Jesse Chan <[email protected]>
Co-authored-by: Jesse Chan <[email protected]>
Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
Co-authored-by: Jesse Chan <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>

examples/dgmulti_2d/elixir_navier_stokes_convergence.jl

@@ -11,7 +11,7 @@ equations = CompressibleEulerEquations2D(1.4)
 # Note: If you change the Navier-Stokes parameters here, also change them in the initial condition
 # I really do not like this structure but it should work for now
 equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(),
-                                                          Mach_freestream=0.5)
+                                                          Mach_freestream=0.5, gradient_variables=GradientVariablesPrimitive())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),

examples/tree_2d_dgsem/elixir_navier_stokes_convergence.jl

@@ -9,7 +9,7 @@ prandtl_number() = 0.72
 
 equations = CompressibleEulerEquations2D(1.4)
 equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(),
-                                                          Mach_freestream=0.5)
+                                                          Mach_freestream=0.5, gradient_variables=GradientVariablesPrimitive())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,

src/Trixi.jl

@@ -137,6 +137,8 @@ export AcousticPerturbationEquations2D,
 export LaplaceDiffusion2D,
        CompressibleNavierStokesDiffusion2D
 
+export GradientVariablesPrimitive, GradientVariablesEntropy
+
 export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov,
        flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner,
        flux_nonconservative_powell,

src/equations/compressible_navier_stokes_2d.jl

@@ -81,7 +81,7 @@ Other normalization strategies exist, see the reference below for details.
   [CERFACS Technical report](https://www.cerfacs.fr/~montagna/TR-CFD-13-77.pdf)
 The scaling used herein is Section 4.5 of the reference.
 """
-struct CompressibleNavierStokesDiffusion2D{RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4}
+struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4}
   # TODO: parabolic
   # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations
   # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function
@@ -97,9 +97,31 @@ struct CompressibleNavierStokesDiffusion2D{RealT <: Real, E <: AbstractCompressi
   R::RealT                   # gas constant (depends on nondimensional scaling!)
 
   equations_hyperbolic::E    # CompressibleEulerEquations2D
+  gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy
 end
 
-function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D; Reynolds, Prandtl, Mach_freestream)
+"""
+#!!! warning "Experimental code"
+#    This code is experimental and may be changed or removed in any future release.
+
+`GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters
+for `CompressibleNavierStokesDiffusion2D`. By default, the gradient variables are set to be
+`GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable
+formulation from
+- Hughes, Mallet, Franca (1986)
+  A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the
+  compressible Euler and Navier-Stokes equations and the second law of thermodynamics.
+  [https://doi.org/10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
+
+Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably entropy stable.
+"""
+struct GradientVariablesPrimitive end
+struct GradientVariablesEntropy end
+
+# default to primitive gradient variables
+function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquations2D;
+                                             Reynolds, Prandtl, Mach_freestream,
+                                             gradient_variables = GradientVariablesPrimitive())
   gamma = equations.gamma
   inv_gamma_minus_one = equations.inv_gamma_minus_one
   Re, Pr, Ma = promote(Reynolds, Prandtl, Mach_freestream)
@@ -115,13 +137,12 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio
   p_inf = 1 / gamma
   u_inf = Mach_freestream
   R     = 1 / gamma
-  CompressibleNavierStokesDiffusion2D{typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one,
-                                                                        Re, Pr, Ma, kappa,
-                                                                        p_inf, u_inf, R,
-                                                                        equations)
+  CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), typeof(equations)}(gamma, inv_gamma_minus_one,
+                                                                                                    Re, Pr, Ma, kappa,
+                                                                                                    p_inf, u_inf, R,
+                                                                                                    equations, gradient_variables)
 end
 
-
 # TODO: parabolic
 # This is the flexibility a user should have to select the different gradient variable types
 # varnames(::typeof(cons2prim)   , ::CompressibleNavierStokesDiffusion2D) = ("v1", "v2", "T")
@@ -132,7 +153,8 @@ varnames(variable_mapping, equations_parabolic::CompressibleNavierStokesDiffusio
 
 # we specialize this function to compute gradients of primitive variables instead of
 # conservative variables.
-gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D) = cons2prim
+gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive}) = cons2prim
+gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy}) = cons2entropy
 
 # Explicit formulas for the diffussive Navier-Stokes fluxes are available, e.g. in Section 2
 # of the paper by Svärd, Carpenter and Nordström
@@ -144,14 +166,13 @@ gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D) = cons2p
 # Note, could be generalized to use Sutherland's law to get the molecular and thermal
 # diffusivity
 function flux(u, gradients, orientation::Integer, equations::CompressibleNavierStokesDiffusion2D)
+  # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
+  rho, v1, v2, _ = convert_transformed_to_primitive(u, equations)
   # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, T)
   # either computed directly or reverse engineered from the gradient of the entropy vairables
-  # by way of the `convert_gradient_variables` function
-  rho, v1, v2, _ = u
-
-  # gradients contains derivatives of each hyperbolic variable
-  _, dv1dx, dv2dx, dTdx = gradients[1]
-  _, dv1dy, dv2dy, dTdy = gradients[2]
+  # by way of the `convert_gradient_variables` function.
+  _, dv1dx, dv2dx, dTdx = convert_derivative_to_primitive(u, gradients[1], equations)
+  _, dv1dy, dv2dy, dTdy = convert_derivative_to_primitive(u, gradients[2], equations)
 
   # Components of viscous stress tensor
 
@@ -204,6 +225,57 @@ end
   return SVector(rho, v1, v2, T)
 end
 
+# Convert conservative variables to entropy
+# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms
+# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion2D`,
+# but this may be confusing to new users.
+cons2entropy(u, equations::CompressibleNavierStokesDiffusion2D) = cons2entropy(u, equations.equations_hyperbolic)
+entropy2cons(w, equations::CompressibleNavierStokesDiffusion2D) = entropy2cons(w, equations.equations_hyperbolic)
+
+# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables.
+# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed
+# variables into primitive variables.
+@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
+  return u_transformed
+end
+
+# TODO: parabolic. Make this more efficient!
+@inline function convert_transformed_to_primitive(u_transformed, equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+  # note: this uses CompressibleNavierStokesDiffusion2D versions of cons2prim and entropy2cons
+  return cons2prim(entropy2cons(u_transformed, equations), equations)
+end
+
+
+# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4) and
+# reverse engineers the gradients to be terms of the primitive variables (v1, v2, T).
+# Helpful because then the diffusive fluxes have the same form as on paper.
+# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused.
+# TODO: parabolic; entropy stable viscous terms
+@inline function convert_derivative_to_primitive(u, gradient, ::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
+  return gradient
+end
+
+# the first argument is always the "transformed" variables.
+@inline function convert_derivative_to_primitive(w, gradient_entropy_vars,
+                                                 equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+
+  # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back.
+  # We can fix this if we directly compute v1, v2, T from the entropy variables
+  u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion2D
+  rho, rho_v1, rho_v2, _ = u
+
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  T  = temperature(u, equations)
+
+  return SVector(gradient_entropy_vars[1],
+                 equations.R * T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[4]), # grad(u) = R*T*(grad(w_2)+v1*grad(w_4))
+                 equations.R * T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[4]), # grad(v) = R*T*(grad(w_3)+v2*grad(w_4))
+                 equations.R * T * T * gradient_entropy_vars[4]                                # grad(T) = R*T^2*grad(w_4))
+                )
+end
+
+
 # This routine is required because `prim2cons` is called in `initial_condition`, which
 # is called with `equations::CompressibleEulerEquations2D`. This means it is inconsistent
 # with `cons2prim(..., ::CompressibleNavierStokesDiffusion2D)` as defined above.
@@ -275,14 +347,14 @@ end
 
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector,
                                                                                            x, t, operator_type::Gradient,
-                                                                                           equations::CompressibleNavierStokesDiffusion2D)
+                                                                                           equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
   v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
   return SVector(u_inner[1], v1, v2, u_inner[4])
 end
 
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, u_inner, normal::AbstractVector,
                                                                                            x, t, operator_type::Divergence,
-                                                                                           equations::CompressibleNavierStokesDiffusion2D)
+                                                                                           equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
   # rho, v1, v2, _ = u_inner
   normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations)
   v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
@@ -294,15 +366,58 @@ end
 
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector,
                                                                                             x, t, operator_type::Gradient,
-                                                                                            equations::CompressibleNavierStokesDiffusion2D)
+                                                                                            equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
   v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
   T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations)
   return SVector(u_inner[1], v1, v2, T)
 end
 
 @inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, u_inner, normal::AbstractVector,
                                                                                             x, t, operator_type::Divergence,
-                                                                                            equations::CompressibleNavierStokesDiffusion2D)
+                                                                                            equations::CompressibleNavierStokesDiffusion2D{GradientVariablesPrimitive})
   return flux_inner
 end
 
+# specialized BC impositions for GradientVariablesEntropy.
+
+# This should return a SVector containing the boundary values of entropy variables.
+# Here, `w_inner` are the transformed variables (e.g., entropy variables).
+#
+# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions
+#             for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022.
+# DOI: 10.1016/j.jcp.2021.110723
+@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector,
+                                                                                                x, t, operator_type::Gradient,
+                                                                                                equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+  v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
+  negative_rho_inv_p = w_inner[4] # w_4 = -rho / p
+  return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p, negative_rho_inv_p)
+end
+
+# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness.
+@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Adiabatic})(flux_inner, w_inner, normal::AbstractVector,
+                                                                                                x, t, operator_type::Divergence,
+                                                                                                equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+  normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x, t, equations)
+  v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
+  _, tau_1n, tau_2n, _ = flux_inner # extract fluxes for 2nd and 3rd equations
+  normal_energy_flux = v1 * tau_1n + v2 * tau_2n + normal_heat_flux
+  return SVector(flux_inner[1], flux_inner[2], flux_inner[3], normal_energy_flux)
+end
+
+@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector,
+                                                                                                 x, t, operator_type::Gradient,
+                                                                                                 equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+  v1, v2 = boundary_condition.boundary_condition_velocity.boundary_value_function(x, t, equations)
+  T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t, equations)
+
+  # the entropy variables w2 = rho * v1 / p = v1 / (equations.R * T) = -v1 * w4. Similarly for w3
+  w4 = -1 / (equations.R * T)
+  return SVector(w_inner[1], -v1 * w4, -v2 * w4, w4)
+end
+
+@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip, <:Isothermal})(flux_inner, w_inner, normal::AbstractVector,
+                                                                                           x, t, operator_type::Divergence,
+                                                                                           equations::CompressibleNavierStokesDiffusion2D{GradientVariablesEntropy})
+  return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4])
+end

src/solvers/dgsem_tree/dg_2d_parabolic.jl

@@ -98,7 +98,6 @@ function transform_variables!(u_transformed, u, mesh::TreeMesh{2},
   end
 end
 
-
 # This is the version used when calculating the divergence of the viscous fluxes
 function calc_volume_integral!(du, flux_viscous,
                                mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic,

test/test_parabolic_2d.jl

@@ -134,6 +134,32 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (isothermal walls)" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"),
+      initial_refinement_level = 2, tspan=(0.0, 0.1),
+      heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)),
+      l2 = [0.002103629650384378, 0.0034358439333976123, 0.0038673598780978413, 0.012670355349347209],
+      linf = [0.012006261793021222, 0.035502125190110666, 0.025107947320650532, 0.11647078036915026]
+    )
+  end
+
+  @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables)" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"),
+      initial_refinement_level=2, tspan=(0,.1), gradient_variables=GradientVariablesEntropy(),
+      l2 = [0.002140374251726679, 0.0034258287094981717, 0.0038915122887464865, 0.012506862342821999],
+      linf = [0.012244412004805971, 0.035075591861236655, 0.02458089234452718, 0.11425600757951138]
+    )
+  end
+
+  @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (Entropy gradient variables, isothermal walls)" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"),
+      initial_refinement_level=2, tspan=(0,.1), gradient_variables=GradientVariablesEntropy(),
+      heat_bc_top_bottom=Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), equations)),
+      l2 = [0.002134973734788134, 0.0034301388278191753, 0.0038928324474145994, 0.012693611436279086],
+      linf = [0.012244236275815057, 0.035054066314196344, 0.02509959850525358, 0.1179561632485715]
+    )
+  end
+
   @trixi_testset "TreeMesh2D: elixir_navier_stokes_convergence.jl (flux differencing)" begin
     @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navier_stokes_convergence.jl"),
       initial_refinement_level = 2, tspan=(0.0, 0.1),