Second order finite volume integral in 2D

examples/tree_2d_dgsem/elixir_euler_convergence_pure_fvO2.jl

@@ -0,0 +1,58 @@
+
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+polydeg = 3 # Governs in this case only the number of subcells
+basis = LobattoLegendreBasis(polydeg)
+surface_flux = flux_hllc
+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)
+coordinates_max = (2.0, 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 = 1.1)
+
+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);

src/solvers/dgsem_tree/dg_2d.jl

@@ -62,7 +62,7 @@ end
 
 function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
                                   P4estMesh{2}, T8codeMesh{2}}, equations,
-                      volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG,
+                      volume_integral::AbstractVolumeIntegralPureLGLFiniteVolume, dg::DG,
                       uEltype)
     A3d = Array{uEltype, 3}
 
@@ -295,6 +295,154 @@ end
     end
 end
 
+function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+                               have_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)
+        fvO2_kernel!(du, u, mesh,
+                     have_nonconservative_terms, equations,
+                     volume_flux_fv, dg, cache, element,
+                     x_interfaces, reconstruction_mode, slope_limiter, true)
+    end
+
+    return nothing
+end
+
+
+@inline function fvO2_kernel!(du, u,
+                              mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+                              have_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, fstar2_L_threaded, fstar2_R_threaded = cache
+    @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes
+
+    # Calculate FV two-point fluxes
+    fstar1_L = fstar1_L_threaded[Threads.threadid()]
+    fstar2_L = fstar2_L_threaded[Threads.threadid()]
+    fstar1_R = fstar1_R_threaded[Threads.threadid()]
+    fstar2_R = fstar2_R_threaded[Threads.threadid()]
+    calcflux_fvO2!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, 
+                   have_nonconservative_terms, equations,
+                   volume_flux_fv, dg, element, cache,
+                   x_interfaces, reconstruction_mode, slope_limiter)
+
+    # Calculate FV volume integral contribution
+    for j in eachnode(dg), i in eachnode(dg)
+        for v in eachvariable(equations)
+            du[v, i, j, element] += (alpha *
+                                     (inverse_weights[i] *
+                                      (fstar1_L[v, i + 1, j] - fstar1_R[v, i, j]) +
+                                      inverse_weights[j] *
+                                      (fstar2_L[v, i, j + 1] - fstar2_R[v, i, j])))
+        end
+    end
+
+    return nothing
+end
+
+@inline function calcflux_fvO2!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, 
+                                u::AbstractArray{<:Any, 4},
+                                mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+                                have_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 j in eachnode(dg), i in 2:nnodes(dg) # We compute FV02 fluxes at the (nnodes(dg) - 1) subcell boundaries
+        #             Reference element:
+        #  -1 ------------------0------------------ 1 -> x
+        # Gauss-Lobatto-Legendre nodes (schematic for k = 3):
+        #   .          .                  .         .
+        #   ^          ^                  ^         ^
+        # Node indices:
+        #   1          2                  3         4
+        # The inner subcell boundaries are governed by the
+        # cumulative sum of the quadrature weights - 1 .
+        #  -1 ------------------0------------------ 1 -> x
+        #        w1-1      (w1+w2)-1   (w1+w2+w3)-1
+        #   |     |             |             |     |
+        # Note that only the inner boundaries are stored.
+        # Subcell interface indices, loop only over 2 -> nnodes(dg) = 4
+        #   1     2             3             4     5
+        #
+        # In general a four-point stencil is required, since we reconstruct the
+        # piecewise linear solution in both subcells next to the subcell interface.
+        # Since these subcell boundaries are not aligned with the DG nodes,
+        # on each neighboring subcell two linear solutions are reconstructed => 4 point stencil.
+        # For the outer interfaces the stencil shrinks since we do not consider values
+        # outside the element (this is a volume integral).
+        #
+        # The left subcell node values are labelled `_ll` (left-left) and `_lr` (left-right), while
+        # the right subcell node values are labelled `_rl` (right-left) and `_rr` (right-right).
+
+        ## Obtain unlimited values in primitive 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_ll = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), j, element),
+                         equations)
+        u_lr = cons2prim(get_node_vars(u, equations, dg, i - 1, j, element),
+                         equations)
+        u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, 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_rr = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1), j,
+                                       element), equations)
+
+        ## Reconstruct values at interfaces with limiting ##
+        u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr,
+                                       x_interfaces, i,
+                                       slope_limiter, dg)
+
+        ## Convert primitive variables back to conservative variables ##
+        flux = volume_flux_fv(prim2cons(u_l, equations), prim2cons(u_r, equations),
+                              1, equations) # orientation 1: x direction
+
+        set_node_vars!(fstar1_L, flux, equations, dg, i, j)
+        set_node_vars!(fstar1_R, flux, equations, dg, i, j)
+    end
+
+    fstar2_L[:, :, 1] .= zero(eltype(fstar2_L))
+    fstar2_L[:, :, nnodes(dg) + 1] .= zero(eltype(fstar2_L))
+    fstar2_R[:, :, 1] .= zero(eltype(fstar2_R))
+    fstar2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fstar2_R))
+
+    for j in 2:nnodes(dg), i in eachnode(dg) # We compute FV02 fluxes at the (nnodes(dg) - 1) subcell boundaries
+        u_ll = cons2prim(get_node_vars(u, equations, dg, i, max(1, j - 2), element),
+                         equations)
+        u_lr = cons2prim(get_node_vars(u, equations, dg, i, j - 1, element),
+                         equations)
+        u_rl = cons2prim(get_node_vars(u, equations, dg, i, j, element),
+                         equations)
+        u_rr = cons2prim(get_node_vars(u, equations, dg, i, min(nnodes(dg), j + 1),
+                                       element), equations)
+
+        ## Reconstruct values at interfaces with limiting ##
+        u_l, u_r = reconstruction_mode(u_ll, u_lr, u_rl, u_rr,
+                                       x_interfaces, j,
+                                       slope_limiter, dg)
+
+        ## Convert primitive variables back to conservative variables ##
+        flux = volume_flux_fv(prim2cons(u_l, equations), prim2cons(u_r, equations),
+                              2, equations) # orientation 1: x direction
+
+        set_node_vars!(fstar2_L, flux, equations, dg, i, j)
+        set_node_vars!(fstar2_R, flux, equations, dg, i, j)
+    end
+
+    return nothing
+end
+
 @inline function fv_kernel!(du, u,
                             mesh::Union{TreeMesh{2}, StructuredMesh{2},
                                         UnstructuredMesh2D, P4estMesh{2},