LO & HO flux comparison

src/Trixi.jl

@@ -127,7 +127,8 @@ export TreeMesh
 export DG,
        DGSEM, LobattoLegendreBasis,
        VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing,
-       VolumeIntegralLocalComparison,
+       VolumeIntegralLocalComparison, VolumeIntegralFluxComparison,
+       VolumeIntegralLocalFluxComparison,
        VolumeIntegralPureLGLFiniteVolume,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
        MortarL2

src/solvers/dg/dg.jl

@@ -91,6 +91,52 @@ 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_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)
+  @nospecialize integral # reduce precompilation time
+
+  if get(io, :compact, false)
+    show(io, integral)
+  else
+    setup = [
+            "volume flux a" => integral.volume_flux_a,
+            "volume flux b" => integral.volume_flux_b
+            ]
+    summary_box(io, "VolumeIntegralFluxComparison", setup)
+  end
+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{VolumeFluxA, VolumeFluxB} <: AbstractVolumeIntegral
+  volume_flux_a::VolumeFluxA
+  volume_flux_b::VolumeFluxB
+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 a" => integral.volume_flux_a,
+            "volume flux b" => integral.volume_flux_b
+            ]
+    summary_box(io, "VolumeIntegralLocalFluxComparison", setup)
+  end
+end
+
+
 """
     VolumeIntegralShockCapturingHG
 

src/solvers/dg/dg_1d.jl

@@ -60,6 +60,30 @@ 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_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{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[]
@@ -268,6 +292,116 @@ 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_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()]
+
+    # 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_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_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_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_b - flux_a)
+      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
+      δ = (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, 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, element)
+    end
+    add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[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 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
+  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_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 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_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_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_a + δ * (flux_b - flux_a))
+      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,