From e1e680ca8574acd10daa2e5bc5e1f49e1ce008f9 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 28 Jul 2023 04:35:33 +0200 Subject: [PATCH] Enstrophy for 2D Navier-Stokes (#1591) * Doubly periodic shear layer * test if prject toml shows up in git diff * remove chnages * Enstrophy for 2D Navier-Stokes --- Project.toml | 2 +- .../elixir_navierstokes_shear_layer.jl | 71 +++++++++++++++++++ src/callbacks_step/analysis_dg2d.jl | 18 +++++ .../compressible_navier_stokes_2d.jl | 15 ++++ test/test_parabolic_2d.jl | 4 ++ 5 files changed, 109 insertions(+), 1 deletion(-) create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl diff --git a/Project.toml b/Project.toml index db410317851..b3ca99be9ec 100644 --- a/Project.toml +++ b/Project.toml @@ -60,8 +60,8 @@ HDF5 = "0.14, 0.15, 0.16" IfElse = "0.1" LinearMaps = "2.7, 3.0" LoopVectorization = "0.12.118" -Makie = "0.19" MPI = "0.20" +Makie = "0.19" MuladdMacro = "0.2.2" Octavian = "0.3.5" OffsetArrays = "1.3" diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl new file mode 100644 index 00000000000..a7cb2fc89f1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl @@ -0,0 +1,71 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), + Prandtl=prandtl_number()) + +function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) + k = 80 + delta = 0.05 + u0 = 1.0 + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = x[2] <= 0.5 ? u0*tanh(k*(x[2]*0.5 - 0.25)) : tanh(k*(0.75 -x[2]*0.5)) + v2 = u0*delta * sin(2*pi*(x[1]*0.5 + 0.25)) + p = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_shear_layer + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +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, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval,) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary \ No newline at end of file diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 4e456f79872..aecabf0e4b7 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -213,6 +213,24 @@ function integrate(func::Func, u, end end +function integrate(func::Func, u, + mesh::Union{TreeMesh{2}, P4estMesh{2}}, + equations, equations_parabolic, + dg::DGSEM, + cache, cache_parabolic; normalize = true) where {Func} + gradients_x, gradients_y = cache_parabolic.gradients + integrate_via_indices(u, mesh, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, element) + gradients_1_local = get_node_vars(gradients_x, equations_parabolic, dg, i, j, + element) + gradients_2_local = get_node_vars(gradients_y, equations_parabolic, dg, i, j, + element) + return func(u_local, (gradients_1_local, gradients_2_local), + equations_parabolic) + end +end + function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 9b06e0b5abf..a1f11717e69 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -300,6 +300,21 @@ end return T end +@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D) + # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v + + omega = vorticity(u, gradients, equations) + return 0.5 * u[1] * omega^2 +end + +@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D) + # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function. + _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations) + _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations) + + return dv2dx - dv1dy +end + # TODO: can we generalize this to MHD? """ struct BoundaryConditionNavierStokesWall diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 471b976e990..57f296b55fe 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -136,6 +136,10 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), initial_refinement_level = 2, tspan=(0.0, 0.1), + analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_integrals=(energy_kinetic, + energy_internal, + enstrophy)), l2 = [0.002111672530658797, 0.0034322351490857846, 0.0038742528195910416, 0.012469246082568561], linf = [0.012006418939223495, 0.035520871209746126, 0.024512747492231427, 0.11191122588756564] )