From db01ebec0702da190ce261947965e4b451c92c6d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 15 Mar 2021 15:49:43 +0100 Subject: [PATCH 1/8] WIP: compare entropy production of local HO fluxes in VolumeIntegralFluxComparison --- src/Trixi.jl | 2 +- src/solvers/dg/dg.jl | 23 +++++++++++++ src/solvers/dg/dg_1d.jl | 75 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 99 insertions(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index fb307712da8..49b93ae62be 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -127,7 +127,7 @@ export TreeMesh export DG, DGSEM, LobattoLegendreBasis, VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing, - VolumeIntegralLocalComparison, + VolumeIntegralLocalComparison, VolumeIntegralFluxComparison, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, MortarL2 diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index 684a6fb9834..fa9a79059b8 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -91,6 +91,29 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalComp end +# TODO: It would be really nice if we could convert other volume integrals to +# the flux form required here easily... +# TODO: This needs a better name... +struct VolumeIntegralFluxComparison{VolumeFlux일, VolumeFlux이} <: AbstractVolumeIntegral + volume_flux_일::VolumeFlux일 + volume_flux_이::VolumeFlux이 +end + +function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxComparison) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "volume flux 1" => integral.volume_flux_일, + "volume flux 2" => integral.volume_flux_이 + ] + summary_box(io, "VolumeIntegralFluxComparison", setup) + end +end + + """ VolumeIntegralShockCapturingHG diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index dd2040d0791..1e4637ca87c 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -60,6 +60,19 @@ function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_일_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_이_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_일_threaded, fluxes_이_threaded) +end + + function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -268,6 +281,68 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_일, volume_flux_이 = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_일_threaded, fluxes_이_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_일 = fluxes_일_threaded[Threads.threadid()] + fluxes_이 = fluxes_이_threaded[Threads.threadid()] + + # compute entropy variables + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + w[i] = cons2entropy(u_node, equations) + end + + # compute local high-order fluxes + local_fluxes_1!(fluxes_일, u, nonconservative_terms, equations, volume_flux_일, dg, element) + local_fluxes_1!(fluxes_이, u, nonconservative_terms, equations, volume_flux_이, dg, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_일, fluxes_이) + flux일 = fluxes_일[i] + flux이 = fluxes_이[i] + if dot(w[i+1] - w[i], flux일 - flux이) <= 0 + fluxes_일[i] = flux일 # flux일 is more dissipative + else + fluxes_일[i] = flux이 # flux이 is more dissipative + end + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_일[1], equations, dg, 1, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_일[i] - fluxes_일[i-1]), equations, dg, i, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_일[end], equations, dg, nnodes(dg), element) + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,3}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for k in i+1:nnodes(dg) + uk = get_node_vars(u, equations, dg, k, element) + for l in 1:i + ul = get_node_vars(u, equations, dg, l, element) + flux1 += weights[l] * derivative_split[l,k] * volume_flux(ul, uk, 1, equations) + end + end + fluxes[i] = flux1 + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, From a438fcdd54ed637387cbddfaabc9f423eed5ebc9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 15 Mar 2021 17:01:46 +0100 Subject: [PATCH 2/8] use some smoothening for flux comparison --- src/solvers/dg/dg_1d.jl | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index 1e4637ca87c..0129a05d619 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -308,11 +308,20 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, for i in eachindex(fluxes_일, fluxes_이) flux일 = fluxes_일[i] flux이 = fluxes_이[i] - if dot(w[i+1] - w[i], flux일 - flux이) <= 0 - fluxes_일[i] = flux일 # flux일 is more dissipative - else - fluxes_일[i] = flux이 # flux이 is more dissipative - end + # if dot(w[i+1] - w[i], flux일 - flux이) <= 0 + # fluxes_일[i] = flux일 # flux일 is more dissipative + # # fluxes_일[i] = 2 * flux일 - flux이 # flux일 is more dissipative + # else + # fluxes_일[i] = flux이 # flux이 is more dissipative + # # fluxes_일[i] = 2 * flux이 - flux일 # flux이 is more dissipative + # end + b = dot(w[i+1] - w[i], flux이 - flux일) + c = 1.0e-12 + # hyp = sqrt(b^2 + c^2) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_일[i] = flux일 + δ * (flux이 - flux일) end # update volume contribution in locally conservative form From 175e284a7697d437b0cf8127d32d808eb04264c0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 15 Mar 2021 18:23:24 +0100 Subject: [PATCH 3/8] WIP:VolumeIntegralLocalFluxComparison --- src/Trixi.jl | 1 + src/solvers/dg/dg.jl | 23 +++++++++++++++++++ src/solvers/dg/dg_1d.jl | 50 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 74 insertions(+) diff --git a/src/Trixi.jl b/src/Trixi.jl index 49b93ae62be..78ae8ad121f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -128,6 +128,7 @@ export DG, DGSEM, LobattoLegendreBasis, VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing, VolumeIntegralLocalComparison, VolumeIntegralFluxComparison, + VolumeIntegralLocalFluxComparison, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, MortarL2 diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index fa9a79059b8..d77540d9c44 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -114,6 +114,29 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxCompa end +# TODO: It would be really nice if we could convert other volume integrals to +# the flux form required here easily... +# TODO: This needs a better name... +struct VolumeIntegralLocalFluxComparison{VolumeFlux일, VolumeFlux이} <: AbstractVolumeIntegral + volume_flux_일::VolumeFlux일 + volume_flux_이::VolumeFlux이 +end + +function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalFluxComparison) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "volume flux 1" => integral.volume_flux_일, + "volume flux 2" => integral.volume_flux_이 + ] + summary_box(io, "VolumeIntegralLocalFluxComparison", setup) + end +end + + """ VolumeIntegralShockCapturingHG diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index 0129a05d619..ef7421d7b5f 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -73,6 +73,17 @@ function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -352,6 +363,45 @@ end end +function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_일, volume_flux_이 = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + w[i] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, element) + du_i = zero(u_i) + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, element) + flux일 = volume_flux_일(u_i, u_ii, 1, equations) + flux이 = volume_flux_이(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i] - w[ii], flux이 - flux일) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux일 + δ * (flux이 - flux일)) + end + add_to_node_vars!(du, du_i, equations, dg, i, element) + end + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, From 6e3425924935217fd71db6737599a2621e08b1a6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Mar 2021 06:19:38 +0100 Subject: [PATCH 4/8] volume fluxes 'numbered' by alphabet --- src/solvers/dg/dg.jl | 20 +++++++-------- src/solvers/dg/dg_1d.jl | 54 ++++++++++++++++++++--------------------- 2 files changed, 37 insertions(+), 37 deletions(-) diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index d77540d9c44..988419d443a 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -94,9 +94,9 @@ end # TODO: It would be really nice if we could convert other volume integrals to # the flux form required here easily... # TODO: This needs a better name... -struct VolumeIntegralFluxComparison{VolumeFlux일, VolumeFlux이} <: AbstractVolumeIntegral - volume_flux_일::VolumeFlux일 - volume_flux_이::VolumeFlux이 +struct VolumeIntegralFluxComparison{VolumeFlux_a, VolumeFlux_b} <: AbstractVolumeIntegral + volume_flux_a::VolumeFlux_a + volume_flux_b::VolumeFlux_b end function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxComparison) @@ -106,8 +106,8 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxCompa show(io, integral) else setup = [ - "volume flux 1" => integral.volume_flux_일, - "volume flux 2" => integral.volume_flux_이 + "volume flux 1" => integral.volume_flux_a, + "volume flux 2" => integral.volume_flux_b ] summary_box(io, "VolumeIntegralFluxComparison", setup) end @@ -117,9 +117,9 @@ end # TODO: It would be really nice if we could convert other volume integrals to # the flux form required here easily... # TODO: This needs a better name... -struct VolumeIntegralLocalFluxComparison{VolumeFlux일, VolumeFlux이} <: AbstractVolumeIntegral - volume_flux_일::VolumeFlux일 - volume_flux_이::VolumeFlux이 +struct VolumeIntegralLocalFluxComparison{VolumeFluxA, VolumeFluxB} <: AbstractVolumeIntegral + volume_flux_a::VolumeFluxA + volume_flux_b::VolumeFluxB end function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalFluxComparison) @@ -129,8 +129,8 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalFlux show(io, integral) else setup = [ - "volume flux 1" => integral.volume_flux_일, - "volume flux 2" => integral.volume_flux_이 + "volume flux a" => integral.volume_flux_a, + "volume flux b" => integral.volume_flux_b ] summary_box(io, "VolumeIntegralLocalFluxComparison", setup) end diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index ef7421d7b5f..e3d0e4c67de 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -66,10 +66,10 @@ function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeInteg FluxType = SVector{nvariables(equations), uEltype} w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()] - fluxes_일_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] - fluxes_이_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] - return (; w_threaded, fluxes_일_threaded, fluxes_이_threaded) + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) end @@ -296,14 +296,14 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms::Val{false}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DGSEM, cache) - @unpack volume_flux_일, volume_flux_이 = volume_integral + @unpack volume_flux_a, volume_flux_b = volume_integral @unpack inverse_weights = dg.basis - @unpack w_threaded, fluxes_일_threaded, fluxes_이_threaded = cache + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache @threaded for element in eachelement(dg, cache) w = w_threaded[Threads.threadid()] - fluxes_일 = fluxes_일_threaded[Threads.threadid()] - fluxes_이 = fluxes_이_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] # compute entropy variables for i in eachnode(dg) @@ -312,35 +312,35 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, end # compute local high-order fluxes - local_fluxes_1!(fluxes_일, u, nonconservative_terms, equations, volume_flux_일, dg, element) - local_fluxes_1!(fluxes_이, u, nonconservative_terms, equations, volume_flux_이, dg, element) + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, element) # compare entropy production of both fluxes and choose the more dissipative one - for i in eachindex(fluxes_일, fluxes_이) - flux일 = fluxes_일[i] - flux이 = fluxes_이[i] - # if dot(w[i+1] - w[i], flux일 - flux이) <= 0 - # fluxes_일[i] = flux일 # flux일 is more dissipative - # # fluxes_일[i] = 2 * flux일 - flux이 # flux일 is more dissipative + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + # if dot(w[i+1] - w[i], flux_a - flux_b) <= 0 + # fluxes_a[i] = flux_a # flux_a is more dissipative + # # fluxes_a[i] = 2 * flux_a - flux_b # flux_a is more dissipative # else - # fluxes_일[i] = flux이 # flux이 is more dissipative - # # fluxes_일[i] = 2 * flux이 - flux일 # flux이 is more dissipative + # fluxes_a[i] = flux_b # flux_b is more dissipative + # # fluxes_a[i] = 2 * flux_b - flux_a # flux_b is more dissipative # end - b = dot(w[i+1] - w[i], flux이 - flux일) + b = dot(w[i+1] - w[i], flux_b - flux_a) c = 1.0e-12 # hyp = sqrt(b^2 + c^2) hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way # δ = (hyp - b) / hyp # add anti-dissipation as dissipation δ = (hyp - b) / 2hyp # just use the more dissipative flux - fluxes_일[i] = flux일 + δ * (flux이 - flux일) + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) end # update volume contribution in locally conservative form - add_to_node_vars!(du, inverse_weights[1] * fluxes_일[1], equations, dg, 1, element) + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, element) for i in 2:nnodes(dg)-1 - add_to_node_vars!(du, inverse_weights[i] * (fluxes_일[i] - fluxes_일[i-1]), equations, dg, i, element) + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, element) end - add_to_node_vars!(du, -inverse_weights[end] * fluxes_일[end], equations, dg, nnodes(dg), element) + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), element) end end @@ -367,7 +367,7 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms::Val{false}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DGSEM, cache) - @unpack volume_flux_일, volume_flux_이 = volume_integral + @unpack volume_flux_a, volume_flux_b = volume_integral @unpack derivative_split = dg.basis @unpack w_threaded = cache @@ -386,15 +386,15 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, du_i = zero(u_i) for ii in eachnode(dg) u_ii = get_node_vars(u, equations, dg, ii, element) - flux일 = volume_flux_일(u_i, u_ii, 1, equations) - flux이 = volume_flux_이(u_i, u_ii, 1, equations) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) d_i_ii = derivative_split[i, ii] - b = -sign(d_i_ii) * dot(w[i] - w[ii], flux이 - flux일) + b = -sign(d_i_ii) * dot(w[i] - w[ii], flux_b - flux_a) c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way # δ = (hyp - b) / hyp # add anti-dissipation as dissipation δ = (hyp - b) / 2hyp # just use the more dissipative flux - du_i += d_i_ii * (flux일 + δ * (flux이 - flux일)) + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) end add_to_node_vars!(du, du_i, equations, dg, i, element) end From d06da8e1ac706b75eaf8b78be29e216164290354 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Mar 2021 06:26:19 +0100 Subject: [PATCH 5/8] adapt switch constant c --- src/solvers/dg/dg_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index e3d0e4c67de..d41192b96f6 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -327,7 +327,7 @@ function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, # # fluxes_a[i] = 2 * flux_b - flux_a # flux_b is more dissipative # end b = dot(w[i+1] - w[i], flux_b - flux_a) - c = 1.0e-12 + c = 1.0e-7 # hyp = sqrt(b^2 + c^2) hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way # δ = (hyp - b) / hyp # add anti-dissipation as dissipation From 5d7e500083f08b7f264551f6ce941f76f0493869 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Mar 2021 06:37:02 +0100 Subject: [PATCH 6/8] =?UTF-8?q?MR.=20Consistency=20=F0=9F=98=8E?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Michael Schlottke-Lakemper --- src/solvers/dg/dg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index 988419d443a..8e0d7b248a1 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -106,8 +106,8 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxCompa show(io, integral) else setup = [ - "volume flux 1" => integral.volume_flux_a, - "volume flux 2" => integral.volume_flux_b + "volume flux a" => integral.volume_flux_a, + "volume flux b" => integral.volume_flux_b ] summary_box(io, "VolumeIntegralFluxComparison", setup) end From 7d3018f8306591e6ee5c085ebb9bb2d243703862 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Mar 2021 07:01:28 +0100 Subject: [PATCH 7/8] VolumeIntegralLocalFluxComparison in 2D and 3D --- src/solvers/dg/dg_2d.jl | 67 ++++++++++++++++++++++++++++++++++ src/solvers/dg/dg_3d.jl | 80 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 147 insertions(+) diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index f7da5f8bc2a..c342a882f69 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -75,6 +75,17 @@ function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -456,6 +467,62 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + w[i,j] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for j in eachnode(dg), i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, j, element) + du_i = zero(u_i) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + + # x + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, j, element) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i,j] - w[ii,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) + end + + # y + for jj in eachnode(dg) + u_jj = get_node_vars(u, equations, dg, i, jj, element) + flux_a = volume_flux_a(u_i, u_jj, 2, equations) + flux_b = volume_flux_b(u_i, u_jj, 2, equations) + d_j_jj = derivative_split[j, jj] + b = -sign(d_j_jj) * dot(w[i,j] - w[i,jj], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a)) + end + + add_to_node_vars!(du, du_i, equations, dg, i, j, element) + end + + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, diff --git a/src/solvers/dg/dg_3d.jl b/src/solvers/dg/dg_3d.jl index 03c5b6cc973..6ef6a4b0ef3 100644 --- a/src/solvers/dg/dg_3d.jl +++ b/src/solvers/dg/dg_3d.jl @@ -94,6 +94,17 @@ function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -476,6 +487,75 @@ function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + w[i,j,k] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, j, k, element) + du_i = zero(u_i) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + + # x + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, j, k, element) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i,j,k] - w[ii,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) + end + + # y + for jj in eachnode(dg) + u_jj = get_node_vars(u, equations, dg, i, jj, k, element) + flux_a = volume_flux_a(u_i, u_jj, 2, equations) + flux_b = volume_flux_b(u_i, u_jj, 2, equations) + d_j_jj = derivative_split[j, jj] + b = -sign(d_j_jj) * dot(w[i,j,k] - w[i,jj,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a)) + end + + # z + for kk in eachnode(dg) + u_kk = get_node_vars(u, equations, dg, i, j, kk, element) + flux_a = volume_flux_a(u_i, u_kk, 3, equations) + flux_b = volume_flux_b(u_i, u_kk, 3, equations) + d_k_kk = derivative_split[k, kk] + b = -sign(d_k_kk) * dot(w[i,j,k] - w[i,j,kk], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_k_kk * (flux_a + δ * (flux_b - flux_a)) + end + + add_to_node_vars!(du, du_i, equations, dg, i, j, k, element) + end + + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, From e89bdf25c1a505fd119b6ee8a227b1051734d286 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 16 Mar 2021 07:20:10 +0100 Subject: [PATCH 8/8] VolumeIntegralFluxComparison in 2D, 3D --- src/solvers/dg/dg_1d.jl | 10 +-- src/solvers/dg/dg_2d.jl | 123 +++++++++++++++++++++++++++++ src/solvers/dg/dg_3d.jl | 166 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 294 insertions(+), 5 deletions(-) diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index d41192b96f6..c6d9879c26a 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -351,11 +351,11 @@ end for i in eachindex(fluxes) flux1 = zero(eltype(fluxes)) - for k in i+1:nnodes(dg) - uk = get_node_vars(u, equations, dg, k, element) - for l in 1:i - ul = get_node_vars(u, equations, dg, l, element) - flux1 += weights[l] * derivative_split[l,k] * volume_flux(ul, uk, 1, equations) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) end end fluxes[i] = flux1 diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index c342a882f69..9d6b6ce6c76 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -75,6 +75,19 @@ function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) +end + + function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") @@ -467,6 +480,116 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] + c = 1.0e-7 + + # compute entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + w[i,j] = cons2entropy(u_node, equations) + end + + # x + for j in eachnode(dg) + # compute local high-order fluxes + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + b = dot(w[i+1,j] - w[i,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, element) + end + + # y + for i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, element) + local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for j in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[j] + flux_b = fluxes_b[j] + b = dot(w[i,j+1] - w[i,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[j] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, element) + for j in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), element) + end + + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, j, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, j, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, j, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) + end + end + fluxes[i] = flux1 + end +end + +@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, element) + @unpack weights, derivative_split = dg.basis + + for j in eachindex(fluxes) + flux2 = zero(eltype(fluxes)) + for jjp in j+1:nnodes(dg) + ujjp = get_node_vars(u, equations, dg, i, jjp, element) + for jjm in 1:j + ujjm = get_node_vars(u, equations, dg, i, jjm, element) + flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations) + end + end + fluxes[j] = flux2 + end +end + + function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms::Val{false}, equations, volume_integral::VolumeIntegralLocalFluxComparison, diff --git a/src/solvers/dg/dg_3d.jl b/src/solvers/dg/dg_3d.jl index 6ef6a4b0ef3..fb1c6143c79 100644 --- a/src/solvers/dg/dg_3d.jl +++ b/src/solvers/dg/dg_3d.jl @@ -94,6 +94,19 @@ function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeInteg end +function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) +end + + function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") @@ -487,6 +500,159 @@ function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, end +function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] + c = 1.0e-7 + + # compute entropy variables + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + w[i,j,k] = cons2entropy(u_node, equations) + end + + # x + for k in eachnode(dg), j in eachnode(dg) + # compute local high-order fluxes + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, k, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, k, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + b = dot(w[i+1,j,k] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, k, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, k, element) + end + + # y + for k in eachnode(dg), i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, k, element) + local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, k, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for j in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[j] + flux_b = fluxes_b[j] + b = dot(w[i,j+1,k] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[j] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, k, element) + for j in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), k, element) + end + + # z + for j in eachnode(dg), i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_3!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, j, element) + local_fluxes_3!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, j, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for k in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[k] + flux_b = fluxes_b[k] + b = dot(w[i,j,k+1] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[k] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, j, 1, element) + for k in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[k] * (fluxes_a[k] - fluxes_a[k-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, j, nnodes(dg), element) + end + + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, j, k, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, j, k, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, j, k, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) + end + end + fluxes[i] = flux1 + end +end + +@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, k, element) + @unpack weights, derivative_split = dg.basis + + for j in eachindex(fluxes) + flux2 = zero(eltype(fluxes)) + for jjp in j+1:nnodes(dg) + ujjp = get_node_vars(u, equations, dg, i, jjp, k, element) + for jjm in 1:j + ujjm = get_node_vars(u, equations, dg, i, jjm, k, element) + flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations) + end + end + fluxes[j] = flux2 + end +end + +@inline function local_fluxes_3!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, j, element) + @unpack weights, derivative_split = dg.basis + + for k in eachindex(fluxes) + flux3 = zero(eltype(fluxes)) + for kkp in k+1:nnodes(dg) + ukkp = get_node_vars(u, equations, dg, i, j, kkp, element) + for kkm in 1:k + ukkm = get_node_vars(u, equations, dg, i, j, kkm, element) + flux3 += weights[kkm] * derivative_split[kkm,kkp] * volume_flux(ukkm, ukkp, 3, equations) + end + end + fluxes[k] = flux3 + end +end + + function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, nonconservative_terms::Val{false}, equations, volume_integral::VolumeIntegralLocalFluxComparison,