From dc794f9dfb4e78aa3d84888ced34ffdca86e9a6f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 6 Mar 2024 13:34:11 +0100 Subject: [PATCH] example --- ...ixir_navierstokes_blast_wave_sutherland.jl | 111 ------------------ ...erstokes_taylor_green_vortex_sutherland.jl | 91 ++++++++++++++ test/test_parabolic_2d.jl | 18 +++ 3 files changed, 109 insertions(+), 111 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_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 deleted file mode 100644 index 370ec675b9b..00000000000 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_blast_wave_sutherland.jl +++ /dev/null @@ -1,111 +0,0 @@ - -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 - # Use 2.0E-3 instead of 1.0E-3 to avoid negative temperature - p = r > 0.5 ? 2.0E-3 : 1.245 - - 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/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl new file mode 100644 index 00000000000..a16858120f8 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 + +# Use Sutherland's law for a temperature-dependent viscosity +@inline function mu(u, equations) + T_ref = 291.15 + + R_specific_air = 287.052874 + T = R_specific_air * 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_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D) + +The classical viscous Taylor-Green vortex in 2D. +This forms the basis behind the 3D case found for instance in + - Jonathan R. Bull and Antony Jameson + Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes + [DOI: 10.2514/6.2014-3210](https://doi.org/10.2514/6.2014-3210) +""" +function initial_condition_taylor_green_vortex(x, t, + equations::CompressibleEulerEquations2D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) + v2 = -A * cos(x[1]) * sin(x[2]) + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0 / 4.0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0) .* pi +coordinates_max = (1.0, 1.0) .* pi +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-9 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 9f1382caa62..659033b8641 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -490,6 +490,24 @@ end tspan=(0.0, 0.7)) end +@trixi_testset "TreeMesh2D: elixir_navierstokes_taylor_green_vortex_sutherland.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_navierstokes_taylor_green_vortex_sutherland.jl"), + l2=[ + 0.001452856280034929, + 0.0007538775539989481, + 0.0007538775539988681, + 0.011035506549989587, + ], + linf=[ + 0.003291912841311362, + 0.002986462478096974, + 0.0029864624780958637, + 0.0231954665514138, + ], + tspan=(0.0, 1.0)) +end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"),