Second-Order Finite Volume Integral in 1D

examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl

@@ -0,0 +1,65 @@
+
+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)
+cells_per_dimension = (8,)
+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

@@ -0,0 +1,59 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+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
+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,
+                                                              :conservation_error))
+
+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,18 @@ export DG,
        FDSBP,
        VolumeIntegralWeakForm, VolumeIntegralStrongForm,
        VolumeIntegralFluxDifferencing,
-       VolumeIntegralPureLGLFiniteVolume,
+       VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
        MortarL2
 
+export reconstruction_small_stencil_inner, reconstruction_small_stencil_full,
+       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,81 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
+                                        volume_flux_fv = flux_lax_friedrichs,
+                                        reconstruction_mode = reconstruction_small_stencil_inner,
+                                        slope_limiter = minmod)
+
+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
+
+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
+    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

@@ -67,6 +67,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")