LO & HO flux comparison

[ranocha]

High-order flux comparison

This basically implement a version comparing the local entropy production based on local high-order fluxes as in the work of Fisher and Carpenter. Right now, it is restricted to 1D.

The density wave works in 1D.

julia> trixi_include("examples/1d/elixir_euler_density_wave.jl", 
           solver=DGSEM(5, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ranocha)), 
           initial_refinement_level=4, 
           tspan=(0.0, 100.0))

EOC does not seem to be too bad (but not necessarily with a similarly clear superconvergence as central DG).

julia> convergence_test("examples/1d/elixir_euler_source_terms.jl", 4, 
           solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ranocha)))
...
####################################################################################################
l2
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
6.23e-06  -         3.65e-06  -         1.06e-05  -         
4.42e-07  3.82      2.83e-07  3.69      6.11e-07  4.11      
2.03e-08  4.45      1.69e-08  4.07      2.58e-08  4.57      
5.29e-10  5.26      4.32e-10  5.29      1.13e-09  4.51      

mean      4.51      mean      4.35      mean      4.39      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
1.73e-05  -         1.23e-05  -         3.77e-05  -         
1.91e-06  3.18      1.13e-06  3.44      3.41e-06  3.47      
9.14e-08  4.38      8.09e-08  3.81      1.31e-07  4.71      
2.05e-09  5.48      1.82e-09  5.47      6.05e-09  4.43      

mean      4.35      mean      4.24      mean      4.20      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("rho", "rho_v1", "rho_e")
  :l2        => [4.50854, 4.34757, 4.39473]
  :linf      => [4.34864, 4.24124, 4.20194]

julia> convergence_test("examples/1d/elixir_euler_source_terms.jl", 4, 
           solver=DGSEM(4, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ranocha)))
...
####################################################################################################
l2
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
1.07e-07  -         9.58e-08  -         1.25e-07  -         
3.87e-09  4.79      3.19e-09  4.91      4.76e-09  4.72      
1.14e-10  5.09      2.42e-11  7.04      2.02e-10  4.56      
6.30e-12  4.17      1.01e-12  4.58      1.08e-11  4.22      

mean      4.68      mean      5.51      mean      4.50      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_e               
error     EOC       error     EOC       error     EOC       
2.88e-07  -         2.30e-07  -         7.54e-07  -         
1.67e-08  4.11      8.13e-09  4.82      3.52e-08  4.42      
5.83e-10  4.84      1.37e-10  5.90      1.36e-09  4.69      
3.27e-11  4.16      5.20e-12  4.71      6.71e-11  4.35      

mean      4.37      mean      5.14      mean      4.49      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("rho", "rho_v1", "rho_e")
  :l2        => [4.68368, 5.51074, 4.49975]
  :linf      => [4.36869, 5.14326, 4.48567]

julia> convergence_test("examples/1d/elixir_burgers_basic.jl", 4, 
           solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ec)))
...
####################################################################################################
l2
scalar              
error     EOC       
1.17e-04  -         
1.24e-05  3.24      
1.09e-06  3.51      
8.08e-08  3.75      

mean      3.50      
----------------------------------------------------------------------------------------------------
linf
scalar              
error     EOC       
4.71e-04  -         
7.03e-05  2.74      
8.26e-06  3.09      
8.10e-07  3.35      

mean      3.06      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("scalar",)
  :l2        => [3.49955]
  :linf      => [3.06089]

julia> convergence_test("examples/1d/elixir_burgers_basic.jl", 4, 
           solver=DGSEM(4, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ec)))
####################################################################################################
l2
scalar              
error     EOC       
3.83e-06  -         
1.40e-07  4.77      
4.91e-09  4.84      
3.79e-11  7.02      

mean      5.54      
----------------------------------------------------------------------------------------------------
linf
scalar              
error     EOC       
1.63e-05  -         
7.48e-07  4.45      
3.49e-08  4.42      
2.68e-10  7.02      

mean      5.30      
----------------------------------------------------------------------------------------------------
Dict{Symbol, Any} with 3 entries:
  :variables => ("scalar",)
  :l2        => [5.54138]
  :linf      => [5.29776]

Interestingly, we get linear stability for high wave numbers for the Burgers' example (where the element-based comparison results in eigenvalues with positive real parts)

julia> trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ec)), k=20); J = jacobian_ad_forward(semi); λ = eigvals(J); maximum(real, λ)
...
-6.14732499940009e-15

julia> trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(flux_ec))), k=20); J = jacobian_ad_forward(semi); λ = eigvals(J); maximum(real, λ)
...
1.1973486561303046

However, we get eigenvalues with positive real parts for low wave numbers (where the element-based comparison was linearly stable).

julia> trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ec)), k=1); J = jacobian_ad_forward(semi); λ = eigvals(J); maximum(real, λ)
...
3.676973569222142e-5

julia> trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(flux_ec))), k=1); J = jacobian_ad_forward(semi); λ = eigvals(J); maximum(real, λ)
...
4.349791949228003e-15

Local flux comparison

I also implemented the version where each two-point flux is compared using the entropy dissipation approximation in 1D. It seems to be linearly stable for Burgers 🥳

julia> max_real_parts = Float64[]; for k in 1:20
           trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralLocalFluxComparison(flux_central, flux_ec)), k=k)
           J = jacobian_ad_forward(semi)
           λ = eigvals(J)
           push!(max_real_parts, maximum(real, λ))
       end; max_real_parts
...
20-element Vector{Float64}:
  7.200005926160668e-15
  5.78540884760829e-17
  4.573883974512219e-15
 -2.858537270733039e-14
  1.570066156412183e-14
 -2.4512827773735094e-14
  9.551606133757615e-15
 -1.1600666345571795e-14
  5.468429999665181e-15
  7.129369996859868e-15
  1.8724548290078468e-14
 -3.7759896354234164e-16
  1.8125874882377547e-14
 -3.574851120683754e-15
 -2.435780893289072e-15
  1.1847684231506737e-14
 -6.3674879231118985e-15
  4.764398278798727e-15
  3.6025714037562356e-15
 -3.010111919300141e-15

The density wave works in 1D.

julia> trixi_include("examples/1d/elixir_euler_density_wave.jl", 
           solver=DGSEM(5, flux_lax_friedrichs, VolumeIntegralLocalFluxComparison(flux_central, flux_ranocha)), 
           initial_refinement_level=4, 
           tspan=(0.0, 100.0))

However, the EOC is reduced significantly the more "discontinuous" the switch between the fluxes gets...

[gregorgassner]

Are you really sure about the „non-typical“ symbols? It makes the code hard to read for me.

[gregorgassner]

I observed, that the anti-dissipation flux comparison is better behaved for the value of the switch smoothing.
Also, 1e-12 basically is the hard switch...maybe relax as default?

[ranocha]

I changed the names of the fluxes and also increased the switch constant c = 1.0e-7. The EOC seems to be a bit nicer, but we still get linear instability for Burgers

julia> max_real_parts = Float64[]; for k in 1:20
           trixi_include("examples/1d/elixir_burgers_linear_stability.jl", solver=DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxComparison(flux_central, flux_ec)), k=k)
           J = jacobian_ad_forward(semi)
           λ = eigvals(J)
           push!(max_real_parts, maximum(real, λ))
       end; max_real_parts
...
20-element Vector{Float64}:
 -8.081296981009977e-15
  8.759236275368565e-5
  0.002672830885772526
  0.0047136084249901344
 -1.0506381259358614e-14
  4.320843044164221e-15
 -3.7196624897745393e-16
  1.0166290131881313e-14
 -4.715186516430352e-15
  9.422661826721416e-15
 -1.1932961017170007e-14
 -1.9287309015117426e-15
  9.550440399427289e-15
 -4.6938725197635265e-15
  2.2519700026712824e-15
 -9.7154478061007e-15
  1.2910526426416925e-14
  3.4252871508758034e-15
  1.2244933605080075e-14
 -1.1586891245214021e-14

[ranocha]

Now in 1D, 2D, and 3D