flux plus dissipation and global Lax-Friedrichs dissipation

[ranocha]

I've added a global Lax-Friedrichs dissipation. Sadly, it looks like the Jacobians of the dissipation and the central scheme do not commute. Thus, we can in general not use this tool to prove anything about linear stability with the global Lax-Friedrichs dissipation.

julia> using Revise, LinearAlgebra, Plots; using Trixi

julia> equations = InviscidBurgersEquation1D();

julia> mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level=4, n_cells_max=10^5);

julia> function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
         k = 1
         2 + sinpi(k * (x[1] - 0.7)) |> SVector
       end
initial_condition_linear_stability (generic function with 1 method)

julia> semi_central = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(3, flux_central));

julia> semi_llf =  SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(3, flux_lax_friedrichs));

julia> semi_glf = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(3, FluxPlusDissipation(flux_central, DissipationGlobalLaxFriedrichs(3.0))));

julia> u0 = compute_coefficients(0.0, semi_central); du = similar(u0); u = u0;

julia> J_central = jacobian_ad_forward(semi_central, u0_ode=u); λ_central = eigvals(J_central); maximum(real, λ_central)
1.4210854715202004e-14

julia> J_llf = jacobian_ad_forward(semi_llf, u0_ode=u); λ_llf = eigvals(J_llf); maximum(real, λ_llf)
-1.2189806126760645e-14

julia> J_glf = jacobian_ad_forward(semi_glf, u0_ode=u); λ_glf = eigvals(J_glf); maximum(real, λ_glf)
1.6996289275834217e-14

julia> commutator(A, B) = A*B - B*A
commutator (generic function with 1 method)

julia> J_glf_diss = J_central - J_glf; commutator(J_glf_diss, J_central) |> norm
240235.5068510898

I could find random initial conditions where the global Lax-Friedrichs surface dissipation was not linearly stable.

julia> u = 2 .+ rand(size(u0)...); extrema(u)
(2.0006933714471224, 2.9590000282032705)

julia> Trixi.rhs!(du, u, semi_central, 0.0); Trixi.integrate(du .* u, semi_central) |> first
0.007953079349215408

julia> Trixi.rhs!(du, u, semi_llf, 0.0); Trixi.integrate(du .* u, semi_llf) |> first
-1.9638326213941046

julia> Trixi.rhs!(du, u, semi_glf, 0.0); Trixi.integrate(du .* u, semi_glf) |> first
-2.1171609819396595

julia> J_central = jacobian_ad_forward(semi_central, u0_ode=u); λ_central = eigvals(J_central); maximum(real, λ_central)
1.7763568394002505e-14

julia> J_llf = jacobian_ad_forward(semi_llf, u0_ode=u); λ_llf = eigvals(J_llf); maximum(real, λ_llf)
1.531368783355667e-14

julia> J_glf = jacobian_ad_forward(semi_glf, u0_ode=u); λ_glf = eigvals(J_glf); maximum(real, λ_glf)
0.006843788081297042

I couldn't find such an example for the local Lax-Friedrichs dissipation for polydeg = 3 but e.g. for polydeg = 7:

julia> mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level=2, n_cells_max=10^5);

julia> semi_central = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(7, flux_central));

julia> semi_llf =  SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(7, flux_lax_friedrichs));

julia> semi_glf = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, DGSEM(7, FluxPlusDissipation(flux_central, DissipationGlobalLaxFriedrichs(3.0))));

julia> u0 = compute_coefficients(0.0, semi_central); du = similar(u0); u = u0;

julia> for i in 1:1000
           u = 2 .+ rand(size(u0)...); extrema(u)
           # Trixi.rhs!(du, u, semi_central, 0.0); Trixi.integrate(du .* u, semi_central) |> first
           Trixi.rhs!(du, u, semi_llf, 0.0); Trixi.integrate(du .* u, semi_llf) |> first > 0 && (println("llf not globally entropy dissipative: ", Trixi.integrate(du .* u, semi_llf) |> first))
           # Trixi.rhs!(du, u, semi_glf, 0.0); Trixi.integrate(du .* u, semi_glf) |> first
           # J_central = jacobian_ad_forward(semi_central, u0_ode=u); λ_central = eigvals(J_central); maximum(real, λ_central)
           J_llf = jacobian_ad_forward(semi_llf, u0_ode=u); λ_llf = eigvals(J_llf); maximum(real, λ_llf) > 1.0e-3 && (println("found it for llf: ", maximum(real, λ_llf)); break)
           # J_glf = jacobian_ad_forward(semi_glf, u0_ode=u); λ_glf = eigvals(J_glf); maximum(real, λ_glf) > 1.0e-13 && (println("found it for glf"); break)
       end
llf not globally entropy dissipative: 0.015894818768715097
llf not globally entropy dissipative: 0.01729277608839308
found it for llf: 0.0017462663659202349

[sloede]

Is it possible that you are pulling in changes from main to es-vol-int via your hr/flux_plus_dissipation branch? I am very confused about changes that look familiar but are marked as new...

[ranocha]

Is it possible that you are pulling in changes from main to es-vol-int via your hr/flux_plus_dissipation branch? I am very confused about changes that look familiar but are marked as new...

Yeah, I'm sorry for that. I based my changes on main instead of es-vol-int and recognized that too late 😒 It's just FluxPlusDissipation and DissipationGlobalLaxFriedrichs

[ranocha]

You only need to look at src/equations/numerical_fluxes.jl and src/equations/equations.jl.

[sloede]

LGTM. I'd add a small comment to one of the includes, but in general this can be merged.