splitting_lax_friedrichs for LinearScalarAdvection1D (#1546)

* splitting_lax_friedrichs for LinearScalarAdvection1D

* Update src/equations/linear_scalar_advection_1d.jl

Co-authored-by: Andrew Winters <[email protected]>

---------

Co-authored-by: Andrew Winters <[email protected]>

examples/tree_1d_fdsbp/elixir_advection_upwind.jl

@@ -0,0 +1,57 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear scalar advection equation equation
+
+equations = LinearScalarAdvectionEquation1D(1.0)
+
+function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation1D)
+    return SVector(sinpi(x[1] - equations.advection_velocity[1] * t))
+end
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order = 1,
+                         accuracy_order = 4,
+                         xmin = -1.0, xmax = 1.0,
+                         N = 16)
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral = SurfaceIntegralUpwind(flux_splitting),
+               volume_integral = VolumeIntegralUpwind(flux_splitting))
+
+coordinates_min = -1.0
+coordinates_max =  1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sin, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
+            ode_default_options()..., callback=callbacks);
+summary_callback() # print the timer summary

src/equations/inviscid_burgers_1d.jl

@@ -132,7 +132,7 @@ end
                              equations::InviscidBurgersEquation1D)
 
 Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)`
-and `f⁻ = 0.5 (f - λ u)` where λ = abs(u).
+and `f⁻ = 0.5 (f - λ u)` where `λ = abs(u)`.
 
 Returns a tuple of the fluxes "minus" (associated with waves going into the
 negative axis direction) and "plus" (associated with waves going into the

src/equations/linear_scalar_advection_1d.jl

@@ -172,6 +172,44 @@ end
     return abs.(equation.advection_velocity)
 end
 
+"""
+    splitting_lax_friedrichs(u, orientation::Integer,
+                             equations::LinearScalarAdvectionEquation1D)
+    splitting_lax_friedrichs(u, which::Union{Val{:minus}, Val{:plus}}
+                             orientation::Integer,
+                             equations::LinearScalarAdvectionEquation1D)
+
+Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)`
+and `f⁻ = 0.5 (f - λ u)` where `λ` is the absolute value of the advection
+velocity.
+
+Returns a tuple of the fluxes "minus" (associated with waves going into the
+negative axis direction) and "plus" (associated with waves going into the
+positive axis direction). If only one of the fluxes is required, use the
+function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`.
+
+!!! warning "Experimental implementation (upwind SBP)"
+    This is an experimental feature and may change in future releases.
+"""
+@inline function splitting_lax_friedrichs(u, orientation::Integer,
+                                          equations::LinearScalarAdvectionEquation1D)
+    fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation, equations)
+    fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation, equations)
+    return fm, fp
+end
+
+@inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer,
+                                          equations::LinearScalarAdvectionEquation1D)
+    a = equations.advection_velocity[1]
+    return a > 0 ? flux(u, orientation, equations) : zero(u)
+end
+
+@inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer,
+                                          equations::LinearScalarAdvectionEquation1D)
+    a = equations.advection_velocity[1]
+    return a < 0 ? flux(u, orientation, equations) : zero(u)
+end
+
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u
 

test/test_tree_1d_fdsbp.jl

@@ -7,6 +7,24 @@ include("test_trixi.jl")
 
 EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_fdsbp")
 
+@testset "Linear scalar advection" begin
+  @trixi_testset "elixir_advection_upwind.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_upwind.jl"),
+      l2   = [1.7735637157305526e-6],
+      linf = [1.0418854521951328e-5],
+      tspan = (0.0, 0.5))
+
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+      t = sol.t[end]
+      u_ode = sol.u[end]
+      du_ode = similar(u_ode)
+      @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+  end
+end
+
 @testset "Inviscid Burgers" begin
   @trixi_testset "elixir_burgers_basic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_basic.jl"),