From bf88ef8c8785a12300e3cc5c58d4605add0af748 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Mar 2024 11:42:00 +0100 Subject: [PATCH] example using sutherland --- ...ixir_navierstokes_blast_wave_sutherland.jl | 110 ++++++++++++++++++ .../compressible_navier_stokes_3d.jl | 5 +- 2 files changed, 113 insertions(+), 2 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl new file mode 100644 index 00000000000..1a78796338c --- /dev/null +++ b/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 diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 556d883f5ee..42b66fadb81 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/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