add different recontruction mode

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

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)
 

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
 

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

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
 

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)