Second-Order Finite Volume Integral in 1D

examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl

@@ -0,0 +1,57 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+surf_flux = flux_hllc
+solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                                     volume_flux_fv = surf_flux,
+                                                                     slope_limiter = monotonized_central))
+
+coordinates_min = 0.0
+coordinates_max = 2.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# 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 = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK22(),
+            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

src/Trixi.jl

@@ -235,13 +235,17 @@ export DG,
        FDSBP,
        VolumeIntegralWeakForm, VolumeIntegralStrongForm,
        VolumeIntegralFluxDifferencing,
-       VolumeIntegralPureLGLFiniteVolume,
+       VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
        MortarL2
 
+export reconstruction_small_stencil, reconstruction_constant,
+       minmod, monotonized_central, superbee, vanLeer_limiter,
+       central_slope
+
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
 

src/equations/linear_scalar_advection_1d.jl

@@ -216,6 +216,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation1D) = u

src/equations/linear_scalar_advection_2d.jl

@@ -277,6 +277,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation2D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation2D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u

src/equations/linear_scalar_advection_3d.jl

@@ -196,6 +196,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation3D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation3D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation3D) = u

src/solvers/dg.jl

@@ -184,6 +184,65 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
+                                        volume_flux_fv = flux_lax_friedrichs,
+                                        reconstruction_mode = reconstruction_small_stencil,
+                                        slope_limiter = minmod)
+
+This gives a formally O(2)-accurate finite volume scheme on an LGL-type subcell
+mesh (LGL = Legendre-Gauss-Lobatto).
+
+!!! 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)
+"""
+struct VolumeIntegralPureLGLFiniteVolumeO2{RealT, Basis, VolumeFluxFV, Reconstruction,
+                                           Limiter} <: AbstractVolumeIntegral
+    x_interfaces::Vector{RealT} # x-coordinates of the sub-cell element interfaces
+    volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
+    reconstruction_mode::Reconstruction # which type of FV reconstruction to use
+    slope_limiter::Limiter # which type of slope limiter function
+end
+
+function VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
+                                             volume_flux_fv = flux_lax_friedrichs,
+                                             reconstruction_mode = reconstruction_small_stencil,
+                                             slope_limiter = minmod) where {Basis}
+    # Suffices to store only the intermediate boundaries of the sub-cell elements                                             
+    x_interfaces = cumsum(basis.weights)[1:(end - 1)] .- 1
+    VolumeIntegralPureLGLFiniteVolumeO2{eltype(basis.weights),
+                                        typeof(basis),
+                                        typeof(volume_flux_fv),
+                                        typeof(reconstruction_mode),
+                                        typeof(slope_limiter)}(x_interfaces,
+                                                               volume_flux_fv,
+                                                               reconstruction_mode,
+                                                               slope_limiter)
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   integral::VolumeIntegralPureLGLFiniteVolumeO2)
+    @nospecialize integral # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        setup = [
+            "FV flux" => integral.volume_flux_fv,
+            "Reconstruction" => integral.reconstruction_mode,
+            "Slope limiter" => integral.slope_limiter,
+            "Subcell boundaries" => vcat([-1.0], integral.x_interfaces, [1.0]),
+        ]
+        summary_box(io, "VolumeIntegralPureLGLFiniteVolumeO2", setup)
+    end
+end
+
 """
     VolumeIntegralSubcellLimiting(limiter;
                                   volume_flux_dg, volume_flux_fv)

src/solvers/dgsem_tree/dg.jl

@@ -63,6 +63,9 @@ include("dg_parallel.jl")
 # Helper structs for parabolic AMR
 include("containers_viscous.jl")
 
+# Some standard functions for the canonical second-order FV (MUSCL) scheme
+include("finite_volume_O2.jl")
+
 # 1D DG implementation
 include("dg_1d.jl")
 include("dg_1d_parabolic.jl")

src/solvers/dgsem_tree/dg_1d.jl

@@ -61,8 +61,9 @@ end
 
 function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}, P4estMesh{1}},
                       equations,
-                      volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG,
-                      uEltype)
+                      volume_integral::Union{VolumeIntegralPureLGLFiniteVolume,
+                                             VolumeIntegralPureLGLFiniteVolumeO2},
+                      dg::DG, uEltype)
     A2dp1_x = Array{uEltype, 2}
     fstar1_L_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
                                 for _ in 1:Threads.nthreads()]
@@ -309,6 +310,24 @@ function calc_volume_integral!(du, u,
     return nothing
 end
 
+function calc_volume_integral!(du, u,
+                               mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+                               nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralPureLGLFiniteVolumeO2,
+                               dg::DGSEM, cache)
+    @unpack x_interfaces, volume_flux_fv, reconstruction_mode, slope_limiter = volume_integral
+
+    # Calculate LGL second-order FV volume integral
+    @threaded for element in eachelement(dg, cache)
+        fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
+                   dg, cache, element,
+                   x_interfaces, reconstruction_mode, slope_limiter,
+                   true)
+    end
+
+    return nothing
+end
+
 @inline function fv_kernel!(du, u,
                             mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                             nonconservative_terms, equations,
@@ -335,6 +354,35 @@ end
     return nothing
 end
 
+@inline function fv_kernel!(du, u,
+                            mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+                            nonconservative_terms, equations,
+                            volume_flux_fv, dg::DGSEM, cache, element,
+                            x_interfaces, reconstruction_mode, slope_limiter,
+                            alpha = true)
+    @unpack fstar1_L_threaded, fstar1_R_threaded = cache
+    @unpack inverse_weights = dg.basis
+
+    # Calculate FV two-point fluxes
+    fstar1_L = fstar1_L_threaded[Threads.threadid()]
+    fstar1_R = fstar1_R_threaded[Threads.threadid()]
+    calcflux_fv!(fstar1_L, fstar1_R, u, mesh, nonconservative_terms, equations,
+                 volume_flux_fv,
+                 dg, element, cache,
+                 x_interfaces, reconstruction_mode, slope_limiter)
+
+    # Calculate FV volume integral contribution
+    for i in eachnode(dg)
+        for v in eachvariable(equations)
+            du[v, i, element] += (alpha *
+                                  (inverse_weights[i] *
+                                   (fstar1_L[v, i + 1] - fstar1_R[v, i])))
+        end
+    end
+
+    return nothing
+end
+
 @inline function calcflux_fv!(fstar1_L, fstar1_R, u::AbstractArray{<:Any, 3},
                               mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                               nonconservative_terms::False,
@@ -388,6 +436,60 @@ end
     return nothing
 end
 
+@inline function calcflux_fv!(fstar1_L, fstar1_R, u::AbstractArray{<:Any, 3},
+                              mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+                              nonconservative_terms::False,
+                              equations, volume_flux_fv, dg::DGSEM, element, cache,
+                              x_interfaces, reconstruction_mode, slope_limiter)
+    fstar1_L[:, 1] .= zero(eltype(fstar1_L))
+    fstar1_L[:, nnodes(dg) + 1] .= zero(eltype(fstar1_L))
+    fstar1_R[:, 1] .= zero(eltype(fstar1_R))
+    fstar1_R[:, nnodes(dg) + 1] .= zero(eltype(fstar1_R))
+
+    for i in 2:nnodes(dg)
+        #             Reference element:             
+        #  -1 -----------------0----------------- 1 -> x
+        # Gauss Lobatto Legendre nodes (schematic for k = 3):
+        #   .         .                 .         .
+        #   ^         ^                 ^         ^
+        # i - 2,    i - 1,              i,      i + 1
+        # mm        ll                  rr      pp
+        # Cell boundaries (schematic for k = 3): 
+        # (note that only the inner three boundaries are stored)
+        #  -1 -----------------0----------------- 1 -> x
+        #   |     |            |             |    |
+        # Cell index:
+        #      1         2              3       4
+
+        # 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.
+        u_mm = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), element),
+                         equations)
+        u_ll = cons2prim(get_node_vars(u, equations, dg, i - 1, element), equations)
+        u_rr = cons2prim(get_node_vars(u, equations, dg, i, element), equations)
+        # Note: If i + 1 > nnodes(dg) 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.
+        u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1),
+                                       element), equations)
+
+        # 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
+        flux = volume_flux_fv(prim2cons(u_ll, equations), prim2cons(u_rr, equations),
+                              1, equations) # orientation 1: x direction
+
+        set_node_vars!(fstar1_L, flux, equations, dg, i)
+        set_node_vars!(fstar1_R, flux, equations, dg, i)
+    end
+
+    return nothing
+end
+
 # We pass the `surface_integral` argument solely for dispatch
 function prolong2interfaces!(cache, u,
                              mesh::TreeMesh{1}, equations, surface_integral, dg::DG)