diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded_fvO2.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded_fvO2.jl new file mode 100644 index 00000000000..20dd6650362 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_extruded_fvO2.jl @@ -0,0 +1,88 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)) + +polydeg = 3 # governs in this case only the number of subcells +basis = LobattoLegendreBasis(polydeg) +surface_flux = flux_lax_friedrichs + +volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, + volume_flux_fv = surface_flux, + reconstruction_mode = reconstruction_O2_full, + slope_limiter = monotonized_central) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D. +# This particular mesh is unstructured in the yz-plane, but extruded in x-direction. +# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded +# in x-direction to ensure free stream preservation on a non-conforming mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +function mapping(xi, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + z = zeta + + 1 / 6 * (cos(1.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + y = eta + 1 / 6 * (cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(2 * pi * (2 * z - 3) / 3)) + + return SVector(xi, y, z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + joinpath(@__DIR__, "cube_unstructured_2.inp")) + +mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, + mapping = mapping) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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_restart = SaveRestartCallback(interval = 100, + save_final_restart = true) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_restart, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., callback = callbacks); diff --git a/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl new file mode 100644 index 00000000000..32a9905ae13 --- /dev/null +++ b/examples/structured_3d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl @@ -0,0 +1,60 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test +source_terms = source_terms_convergence_test + +polydeg = 2 # governs in this case only the number of subcells +basis = LobattoLegendreBasis(polydeg) +surface_flux = flux_hll + +volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, + volume_flux_fv = surface_flux, + reconstruction_mode = reconstruction_O2_full, + slope_limiter = monotonized_central) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (2.0, 2.0, 2.0) +cells_per_dimension = (2, 2, 2) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = false) + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = boundary_condition_default(mesh, boundary_condition) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) +############################################################################### +# run the simulation + +sol = solve(ode, ParsaniKetchesonDeconinck3S82(); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., callback = callbacks); diff --git a/src/solvers/dgsem/calc_volume_integral.jl b/src/solvers/dgsem/calc_volume_integral.jl index 7271a502a92..a337d8b6eeb 100644 --- a/src/solvers/dgsem/calc_volume_integral.jl +++ b/src/solvers/dgsem/calc_volume_integral.jl @@ -80,11 +80,7 @@ function calc_volume_integral!(du, u, mesh, return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}, - TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}, - UnstructuredMesh2D, T8codeMesh{2}, - TreeMesh{3}}, +function calc_volume_integral!(du, u, mesh, have_nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingRRG, dg::DGSEM, cache) @@ -143,11 +139,7 @@ function calc_volume_integral!(du, u, mesh, return nothing end -function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}, - TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}, - UnstructuredMesh2D, T8codeMesh{2}, - TreeMesh{3}}, +function calc_volume_integral!(du, u, mesh, have_nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolumeO2, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index 7de87ad75e6..960a082d25a 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -528,6 +528,92 @@ end return nothing end +@inline function calcflux_fvO2!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, + fstar3_L, fstar3_R, u, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, + have_nonconservative_terms::False, + equations, + volume_flux_fv, dg::DGSEM, element, cache, + sc_interface_coords, reconstruction_mode, slope_limiter) + @unpack normal_vectors_1, normal_vectors_2, normal_vectors_3 = cache.normal_vectors + + for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg) + u_ll = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), j, k, element), + equations) + u_lr = cons2prim(get_node_vars(u, equations, dg, i - 1, j, k, element), + equations) + u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element), + equations) + + u_rr = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1), j, k, + element), equations) + + u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr, + sc_interface_coords, i, + slope_limiter, dg) + + normal_direction = get_normal_vector(normal_vectors_1, i - 1, j, k, element) + + contravariant_flux = volume_flux_fv(prim2cons(u_l, equations), + prim2cons(u_r, equations), + normal_direction, equations) + + set_node_vars!(fstar1_L, contravariant_flux, equations, dg, i, j, k) + set_node_vars!(fstar1_R, contravariant_flux, equations, dg, i, j, k) + end + + for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg) + u_ll = cons2prim(get_node_vars(u, equations, dg, i, max(1, j - 2), k, element), + equations) + u_lr = cons2prim(get_node_vars(u, equations, dg, i, j - 1, k, element), + equations) + u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element), + equations) + u_rr = cons2prim(get_node_vars(u, equations, dg, i, min(nnodes(dg), j + 1), k, + element), equations) + + u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr, + sc_interface_coords, j, + slope_limiter, dg) + + normal_direction = get_normal_vector(normal_vectors_2, i, j - 1, k, element) + + contravariant_flux = volume_flux_fv(prim2cons(u_l, equations), + prim2cons(u_r, equations), + normal_direction, equations) + + set_node_vars!(fstar2_L, contravariant_flux, equations, dg, i, j, k) + set_node_vars!(fstar2_R, contravariant_flux, equations, dg, i, j, k) + end + + for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg) + u_ll = cons2prim(get_node_vars(u, equations, dg, i, j, max(1, k - 2), element), + equations) + u_lr = cons2prim(get_node_vars(u, equations, dg, i, j, k - 1, element), + equations) + u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, k, element), + equations) + u_rr = cons2prim(get_node_vars(u, equations, dg, i, j, min(nnodes(dg), k + 1), + element), equations) + + u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr, + sc_interface_coords, k, + slope_limiter, dg) + + normal_direction = get_normal_vector(normal_vectors_3, i, j, k - 1, element) + + contravariant_flux = volume_flux_fv(prim2cons(u_l, equations), + prim2cons(u_r, equations), + normal_direction, equations) + + set_node_vars!(fstar3_L, contravariant_flux, equations, dg, i, j, k) + set_node_vars!(fstar3_R, contravariant_flux, equations, dg, i, j, k) + end + + return nothing +end + function calc_interface_flux!(cache, u, mesh::StructuredMesh{3}, have_nonconservative_terms, # can be True/False equations, surface_integral, dg::DG) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 10101dd34d3..12de969bf2f 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -216,6 +216,29 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_free_stream_extruded_fvO2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_extruded_fvO2.jl"), + l2=[ + 1.5286679406524262e-16, + 4.950311927136664e-16, + 4.69913693630953e-16, + 4.012046698403924e-16, + 2.0419659904946413e-15 + ], + linf=[ + 2.220446049250313e-15, + 3.3029134982598407e-15, + 7.716050021144838e-15, + 1.27675647831893e-14, + 3.019806626980426e-14 + ], + tspan=(0.0, 0.1)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_euler_free_stream_boundaries.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_boundaries.jl"), diff --git a/test/test_structured_3d.jl b/test/test_structured_3d.jl index 4d0464921e6..1ae0c13fb16 100644 --- a/test/test_structured_3d.jl +++ b/test/test_structured_3d.jl @@ -147,6 +147,28 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_source_terms_nonperiodic_fvO2.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic_fvO2.jl"), + l2=[ + 0.046860892952192236, + 0.04269641872975366, + 0.04269641872975368, + 0.042696418729753556, + 0.1443421265167028 + ], + linf=[ + 0.42250001347847244, + 0.2975151811754566, + 0.29751518117545483, + 0.29751518117545617, + 0.3982472383589144 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_euler_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), l2=[