example using sutherland

examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl

@@ -0,0 +1,110 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+
+# Compute dynamic viscosity according to Sutherland's law for air
+@inline function mu(u, equations)
+    T_ref = 291.15
+    T = Trixi.temperature(u, equations)
+    C_air = 120.0
+    mu0_air = 18.27e-6
+
+    return mu0_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5
+end
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
+                                                          Prandtl = prandtl_number())
+
+"""
+    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+A medium blast wave taken from
+- Sebastian Hennemann, Gregor J. Gassner (2020)
+  A provably entropy stable subcell shock capturing approach for high order split form DG
+  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+"""
+function initial_condition_blast_wave(x, t, equations)
+    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
+    # 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 : 1.1691
+    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
+    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
+    p = r > 0.5 ? 2.0E-3 : 1.245 # Use 2.0E-3 instead of 1.0E-3 to avoid negative temperature
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Unstructured mesh with 48 cells of the square domain [-1, 1]^n
+mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp",
+                           joinpath(@__DIR__, "square_unstructured_1.inp"))
+
+mesh = P4estMesh{2}(mesh_file, polydeg = 3, initial_refinement_level = 3)
+
+boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver,
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK54(),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/equations/compressible_navier_stokes_3d.jl

@@ -80,7 +80,7 @@ where
 w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p}
 ```
 """
-struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real,
+struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu,
                                            E <: AbstractCompressibleEulerEquations{3}} <:
        AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables}
     # TODO: parabolic
@@ -111,8 +111,9 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio
     kappa = gamma * inv_gamma_minus_one / Prandtl
 
     CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma),
+                                        typeof(mu),
                                         typeof(equations)}(gamma, inv_gamma_minus_one,
-                                                           μ, Prandtl, kappa,
+                                                           mu, Prandtl, kappa,
                                                            equations,
                                                            gradient_variables)
 end