Viscous terms from entropy gradients

examples/tree_2d_dgsem/elixir_navier_stokes_es.jl

@@ -0,0 +1,75 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+reynolds_number() = 1600
+prandtl_number() = 0.72
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, Reynolds=reynolds_number(), Prandtl=prandtl_number(),
+                                                          Mach_freestream=0.5)
+
+# TODO: For entropy stability testing
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0, -2.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 2.0,  2.0) # maximum coordinates (max(x), max(y))
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+function initial_condition_blast_wave(x, t, equations)
+  # Set up polar coordinates
+  inicenter = SVector(0.0, 0.0)
+  x_norm = x[1] - inicenter[1]
+  y_norm = x[2] - inicenter[2]
+  r = sqrt(x_norm^2 + y_norm^2)
+  phi = atan(y_norm, x_norm)
+  sin_phi, cos_phi = sincos(phi)
+
+  # Calculate primitive variables
+  rho = r > 0.5 ? 1.0 : 2.0
+  v1  = r > 0.5 ? 0.0 : 0.2 * cos_phi
+  v2  = r > 0.5 ? 0.0 : 0.2 * sin_phi
+  p   = r > 0.5 ? 1.0 : 3.0
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_blast_wave
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span `tspan`
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary
+

src/equations/compressible_navier_stokes_2d.jl

@@ -132,7 +132,11 @@ 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) = cons2prim
+
+# TODO: parabolic; entropy stable viscous terms
+gradient_variable_transformation(::CompressibleNavierStokesDiffusion2D) = 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
@@ -204,6 +208,39 @@ end
   return SVector(rho, v1, v2, T)
 end
 
+
+# Convert conservative variables to entropy
+# Note, only w_2, w_3, w_4 are needed for the viscous fluxes so we avoid computing
+# w_1 and simply copy over rho.
+# TODO: parabolic; entropy stable viscous terms
+@inline function cons2entropy(u, equations::CompressibleNavierStokesDiffusion2D)
+  rho, rho_v1, rho_v2, rho_e = u
+
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho)
+
+  return SVector(rho, rho_v1/p, rho_v2/p, -rho/p)
+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_gradient_variables(u, gradient_entropy_vars, equations::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.

src/solvers/dgsem_tree/dg_2d_parabolic.jl

@@ -26,6 +26,11 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::Abstra
     gradients, u_transformed, t, mesh, equations_parabolic, boundary_conditions_parabolic, dg,
     cache, cache_parabolic)
 
+  # Convert the gradients to a form more suitable for viscous flux calculations
+  # TODO: parabolic; entropy stable viscous terms
+  @trixi_timeit timer() "transform gradients" transform_gradients!(
+    gradients, u, mesh, equations_parabolic, dg, parabolic_scheme, cache, cache_parabolic)
+
   # Compute and store the viscous fluxes
   @trixi_timeit timer() "calculate viscous fluxes" calc_viscous_fluxes!(
     flux_viscous, gradients, u_transformed, mesh, equations_parabolic, dg, cache, cache_parabolic)
@@ -99,6 +104,34 @@ function transform_variables!(u_transformed, u, mesh::TreeMesh{2},
 end
 
 
+# Transform gradients from the `u_transformed` variable set prior to computing the viscous fluxes
+# (e.g. entropy to primitive variables).
+# TODO: parabolic; entropy stable viscous terms
+function transform_gradients!(gradients, u, mesh::TreeMesh{2},
+                              equations_parabolic::AbstractEquationsParabolic,
+                              dg::DG, parabolic_scheme, cache, cache_parabolic)
+  gradients_x, gradients_y = gradients
+
+  @threaded for element in eachelement(dg, cache)
+    for j in eachnode(dg), i in eachnode(dg)
+      # Get the solution and gradients
+      u_node = get_node_vars(u, equations_parabolic, dg, i, j, element)
+      entropy_gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, element)
+      entropy_gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, element)
+
+      # Convert entropy gradients into their primitive variable form
+      primitive_gradients_1_node = convert_gradient_variables(u_node, entropy_gradients_1_node, equations_parabolic)
+      primitive_gradients_2_node = convert_gradient_variables(u_node, entropy_gradients_2_node, equations_parabolic)
+
+      # Save them back in same arrays to avoid extra storage
+      set_node_vars!(gradients_x, primitive_gradients_1_node, equations_parabolic, dg, i, j, element)
+      set_node_vars!(gradients_y, primitive_gradients_2_node, equations_parabolic, dg, i, j, element)
+    end
+  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,