diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl new file mode 100644 index 0000000000..a1a9841eb0 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl @@ -0,0 +1,66 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations1D(1.4) + +initial_condition = initial_condition_convergence_test + +source_terms = source_terms_convergence_test + +# you can either use a single function to impose the BCs weakly in all +# 1*ndims == 2 directions or you can pass a tuple containing BCs for +# each direction +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = (x_neg = boundary_condition, + x_pos = boundary_condition) + +polydeg = 8 # Governs in this case only the number of subcells +basis = LobattoLegendreBasis(polydeg) +surf_flux = flux_hll +solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, + volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, + volume_flux_fv = surf_flux, + reconstruction_mode = reconstruction_small_stencil_inner, + slope_limiter = monotonized_central)) + +f1() = SVector(0.0) +f2() = SVector(2.0) +initial_refinement_level = 3 # required for convergence study +cells_per_dimension = 2^initial_refinement_level +mesh = StructuredMesh((cells_per_dimension,), (f1, f2), periodicity = false) + +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, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.3) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, ORK256(), + 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_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl index 7961bf5076..df8bdae970 100644 --- a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl +++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl @@ -9,12 +9,13 @@ equations = CompressibleEulerEquations1D(1.4) initial_condition = initial_condition_convergence_test -polydeg = 3 +polydeg = 3 # Governs in this case only the number of subcells basis = LobattoLegendreBasis(polydeg) surf_flux = flux_hllc solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux, volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = surf_flux, + reconstruction_mode = reconstruction_small_stencil_full, slope_limiter = monotonized_central)) coordinates_min = 0.0 @@ -37,7 +38,8 @@ summary_callback = SummaryCallback() analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, extra_analysis_errors = (:l2_error_primitive, - :linf_error_primitive)) + :linf_error_primitive, + :conservation_error)) alive_callback = AliveCallback(analysis_interval = analysis_interval) diff --git a/src/Trixi.jl b/src/Trixi.jl index c480473bab..5d53a58fd1 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -242,7 +242,8 @@ export DG, SurfaceIntegralUpwind, MortarL2 -export reconstruction_small_stencil, reconstruction_constant, +export reconstruction_small_stencil_inner, reconstruction_small_stencil_full, + reconstruction_constant, minmod, monotonized_central, superbee, vanLeer_limiter, central_slope diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index cf3dd08f25..07765eeb9a 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -187,20 +187,36 @@ end """ VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis; volume_flux_fv = flux_lax_friedrichs, - reconstruction_mode = reconstruction_small_stencil, + reconstruction_mode = reconstruction_small_stencil_inner, slope_limiter = minmod) -This gives a formally O(2)-accurate finite volume scheme on an LGL-type subcell +This gives an up to O(2)-accurate finite volume scheme on an LGL-type subcell mesh (LGL = Legendre-Gauss-Lobatto). +Depending on the `reconstruction_mode` and `slope_limiter`, experimental orders of convergence +between 1 and 2 can be expected in practice. +Currently, all reconstructions are purely cell-local, i.e., no neighboring elements are +queried at reconstruction stage. + +The non-boundary subcells are always reconstructed using the standard MUSCL-type reconstruction. +For the subcells at the boundaries, two options are available: + +1) The unlimited slope is used on these cells. This gives full O(2) accuracy, but may lead to overshoots between cells. + The `reconstruction_mode` corresponding to this is `reconstruction_small_stencil_full`. +2) On boundary subcells, the solution is represented using a constant value, thereby falling back to formally only O(1). + The `reconstruction_mode` corresponding to this is `reconstruction_small_stencil_inner`. + In the reference below, this is the recommended reconstruction mode and is thus used by default. !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. ## References -- Hennemann, 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) +See especially Sections 3.2 and 4 and Appendix D of the paper + +- Rueda-Ramírez, Hennemann, Hindenlang, Winters, & Gassner (2021). + "An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. + Part II: Subcell finite volume shock capturing" + [JCP: 2021.110580](https://doi.org/10.1016/j.jcp.2021.110580) """ struct VolumeIntegralPureLGLFiniteVolumeO2{RealT, Basis, VolumeFluxFV, Reconstruction, Limiter} <: AbstractVolumeIntegral diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 681b396eed..2b2e6cb5b1 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -461,7 +461,7 @@ end # Cell index: # 1 2 3 4 - # Obtain unlimited values in primal variables + ## Obtain unlimited values in primal variables ## # Note: If i - 2 = 0 we do not go to neighbor element, as one would do in a finite volume scheme. # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out. @@ -474,12 +474,12 @@ end u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1), element), equations) - # Reconstruct values at interfaces with limiting + ## Reconstruct values at interfaces with limiting ## u_ll, u_rr = reconstruction_mode(u_mm, u_ll, u_rr, u_pp, x_interfaces, i, slope_limiter, dg) - # Convert primal variables back to conservative variables + ## Convert primal variables back to conservative variables ## flux = volume_flux_fv(prim2cons(u_ll, equations), prim2cons(u_rr, equations), 1, equations) # orientation 1: x direction diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl index 119f86c378..0ee8a01b5a 100644 --- a/src/solvers/dgsem_tree/finite_volume_O2.jl +++ b/src/solvers/dgsem_tree/finite_volume_O2.jl @@ -1,12 +1,3 @@ -@inline function linear_reconstruction(u_ll, u_rr, s_l, s_r, - x_ll, x_rr, x_interfaces, node_index) - # Linear reconstruction at the interface - u_ll = u_ll + s_l * (x_interfaces[node_index - 1] - x_ll) - u_rr = u_rr + s_r * (x_interfaces[node_index - 1] - x_rr) - - return u_ll, u_rr -end - """ reconstruction_constant(u_mm, u_ll, u_rr, u_pp, x_interfaces, @@ -22,6 +13,15 @@ Formally O(1) accurate. return u_ll, u_rr end +@inline function linear_reconstruction(u_ll, u_rr, s_l, s_r, + x_ll, x_rr, x_interfaces, node_index) + # Linear reconstruction at the interface + u_ll = u_ll + s_l * (x_interfaces[node_index - 1] - x_ll) + u_rr = u_rr + s_r * (x_interfaces[node_index - 1] - x_rr) + + return u_ll, u_rr +end + # Reference element: # -1 -----------------0----------------- 1 -> x # Gauss-Lobatto-Legendre nodes (schematic for k = 3): @@ -37,17 +37,18 @@ end # 1 2 3 4 """ - reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp, - x_interfaces, node_index, limiter, dg) + reconstruction_small_stencil_full(u_mm, u_ll, u_rr, u_pp, + x_interfaces, node_index, limiter, dg) Computes limited (linear) slopes on the subcells for a DGSEM element. The supplied `limiter` governs the choice of slopes given the nodal values `u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes. -Formally O(2) accurate. +Formally O(2) accurate when used without a limiter, i.e., the `central_slope`. +For the subcell boundaries the reconstructed slopes are not limited. """ -@inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp, - x_interfaces, node_index, - limiter, dg) +@inline function reconstruction_small_stencil_full(u_mm, u_ll, u_rr, u_pp, + x_interfaces, node_index, + limiter, dg) @unpack nodes = dg.basis x_ll = nodes[node_index - 1] x_rr = nodes[node_index] @@ -74,6 +75,47 @@ Formally O(2) accurate. linear_reconstruction(u_ll, u_rr, s_l, s_r, x_ll, x_rr, x_interfaces, node_index) end +""" + reconstruction_small_stencil_inner(u_mm, u_ll, u_rr, u_pp, + x_interfaces, node_index, limiter, dg) + +Computes limited (linear) slopes on the *inner* subcells for a DGSEM element. +The supplied `limiter` governs the choice of slopes given the nodal values +`u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes. +For the outer, i.e., boundary subcells, constant values are used, i.e, no reconstruction. +This reduces the order of the scheme below 2. +""" +@inline function reconstruction_small_stencil_inner(u_mm, u_ll, u_rr, u_pp, + x_interfaces, node_index, + limiter, dg) + @unpack nodes = dg.basis + x_ll = nodes[node_index - 1] + x_rr = nodes[node_index] + + # Middle element slope + s_mm = (u_rr - u_ll) / (x_rr - x_ll) + + if node_index == 2 # Catch case mm == ll + # Do not reconstruct at the boundary + s_l = zero(s_mm) + else + # Left element slope + s_ll = (u_ll - u_mm) / (x_ll - nodes[node_index - 2]) + s_l = limiter.(s_ll, s_mm) + end + + if node_index == nnodes(dg) # Catch case rr == pp + # Do not reconstruct at the boundary + s_r = zero(s_mm) + else + # Right element slope + s_rr = (u_pp - u_rr) / (nodes[node_index + 1] - x_rr) + s_r = limiter.(s_mm, s_rr) + end + + linear_reconstruction(u_ll, u_rr, s_l, s_r, x_ll, x_rr, x_interfaces, node_index) +end + """ central_slope(sl, sr)