Merge bd3b6a5 into 99534c8

src/HyperbolicDiffEq.jl

@@ -126,6 +126,7 @@ include("balance_laws/sextic.jl")
 include("balance_laws/septic.jl")
 include("balance_laws/octic.jl")
 include("balance_laws/buckley_leverette.jl")
+include("balance_laws/quartic_nonconvex.jl")
 include("balance_laws/shallow_water.jl")
 include("balance_laws/euler.jl")
 include("balance_laws/keyfitz_kranzer.jl")
@@ -136,7 +137,8 @@ include("utils.jl")
 # models
 export ConstantLinearAdvection, Burgers, IntegralQuantitiesBurgers,
        Cubic, Quartic, Quintic, Sextic, Septic, Octic,
-       BuckleyLeverette
+       BuckleyLeverette,
+       QuarticNonconvex
 export ShallowWater, ShallowWaterVar1D
 export Euler, EulerVar1D, EulerVar2D, EulerVar3D, IntegralQuantitiesEuler
 export KeyfitzKranzer, KeyfitzKranzerVar, IntegralQuantitiesKeyfitzKranzer

src/balance_laws/quartic_nonconvex.jl

@@ -0,0 +1,76 @@
+
+"""
+    QuarticNonconvex{T}
+
+The scalar conservation law
+``
+  \\partial_t u + \\partial_x ( u(t,x)^4 - 10 u(t,x)^2 + 3 u(t,x) ) = 0
+``
+in one space dimensions using `T` as scalar type.
+"""
+struct QuarticNonconvex{T} <: ScalarBalanceLaw{T,1} end
+
+function QuarticNonconvex(T=Float64)
+  QuarticNonconvex{T}()
+end
+
+function Base.show(io::IO, model::QuarticNonconvex{T}) where {T}
+  print(io, "Scalar conservation law {T=", T, "} with flux f(u) = u^4 - 10 u^2 + 3 u")
+end
+
+
+"""
+    flux(u, model::QuarticNonconvex)
+
+Compute the flux of `u` for `model`.
+"""
+@inline flux(u, model::QuarticNonconvex) = u^2 * (u^2 - 10) + 3u
+
+"""
+    speed(u::Real, model::QuarticNonconvex)
+
+Compute the speed f'(`u`) for `model`.
+"""
+@inline speed(u, model::QuarticNonconvex) = 4*u^3 - 20*u + 3
+
+
+function min_max_speed(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
+  min_speed = zero(T)
+  max_speed = zero(T)
+
+  if ur > 1.2909944487358056284 && -2.5819888974716112568 < ul < 1.2909944487358056284
+    min_speed = T(-14.213259316477408379)
+  elseif (ur > 1.2909944487358056284 && (ul >= 1.2909944487358056284 || ul <= -2.5819888974716112568)) || (ur <= 1.2909944487358056284 && 2 * ul + ur + sqrt(20 - 3 * ur^2) <= 0) || 1.2909944487358056284 + ur <= 0
+    min_speed = 3 + 4 * ul * (-5 + ul^2)
+  else
+    min_speed = 3 + 4 * ur * (-5 + ur^2)
+  end
+
+  if (1.2909944487358056284 + ul <= 0 && 1.2909944487358056284 + ur > 0 && ur < 1.2909944487358056284) || (1.2909944487358056284 + ul < 0 && 1.2909944487358056284 <= ur < 2.5819888974716112568)
+    max_speed = T(20.213259316477408379)
+  elseif (1.2909944487358056284 + ul == 0 && ur == 2.5819888974716112568) || (1.2909944487358056284 + ul >= 0 && (ur ≈ 1.2909944487358056284 || (sqrt(20 - 3 * ur^2) >= 2 * ul + ur && 1.2909944487358056284 < ur < 2.5819888974716112568))) || (ur < 1.2909944487358056284 && 1.2909944487358056284 + ul > 0)
+    max_speed = 3 + 4 * ul * (-5 + ul^2)
+  else
+    max_speed = 3 + 4 * ur * (-5 + ur^2)
+  end
+
+  min_speed, max_speed
+end
+
+function max_abs_speed(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
+  min_speed, max_speed = min_max_speed(ul, ur, model)
+  max(abs(min_speed), abs(max_speed))
+end
+
+
+################################################################################
+
+
+"""
+    (::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::QuarticNonconvex{T}) where {T}
+
+Compute the energy (L₂ entropy) conservative flux between `uₗ` and `uᵣ` for `model`.
+"""
+function (flux::EnergyConservativeFlux)(ul::T, ur::T, model::QuarticNonconvex{T}) where T<:Real
+  (ul^4 + ur*ul^3 + ur^2*ul^2 + ur^3*ul + ur^4) / 5 - 10 * (ul^2 + ul*ur + ur^2) / 3 + 3 * (ul + ur) / 2
+end

test/balance_laws/quartic_nonconvex_test.jl

@@ -0,0 +1,34 @@
+using Test, OrdinaryDiffEq, HyperbolicDiffEq, DiffEqCallbacks
+
+function runtest()
+    xmin = -1.0
+    xmax = +1.0
+    p = 1
+    N = 2^6
+    tspan = (0., 0.02)
+    ode_solver = OrdinaryDiffEq.SSPRK33()
+
+    mesh = UniformPeriodicMesh1D(xmin, xmax, N)
+    basis = LobattoLegendre(p)
+    balance_law = QuarticNonconvex()
+    fvol = EnergyConservativeFlux()
+    fnum = FluxPlusDissipation(EnergyConservativeFlux(), LocalLaxFriedrichsDissipation(max_abs_speed))
+    semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, mesh, basis, fvol, fnum)
+    ode = semidiscretise(semidisc, sinpi, tspan)
+
+    save_func = (u,t,integrator) -> integrate(u->IntegralQuantitiesBurgers(u,balance_law), u, integrator.f.f)
+    saved_values = SavedValues(Float64, IntegralQuantitiesBurgers{Float64})
+    saving = SavingCallback(save_func, saved_values, saveat=range(tspan..., length=10^2))
+    stepsize = StepsizeLimiter((u,p,t)->max_dt(t,u,semidisc), safety_factor=0.5)
+
+    sol = solve(ode, ode_solver, dt=0.5*max_dt(tspan[1], ode.u0, semidisc), save_everystep=false, callback=CallbackSet(stepsize,saving))
+
+    mass = first.(saving.affect!.saved_values.saveval)
+    entropy = last.(saving.affect!.saved_values.saveval)
+    @test mass[1] ≈ mass[end] atol=10*eps()
+    @test entropy[1] >= entropy[end]
+
+    nothing
+end
+
+runtest()

test/runtests.jl

@@ -9,6 +9,7 @@ using HyperbolicDiffEq
   @time @testset "Flux Difference" begin include("flux_difference_test.jl") end
   @time @testset "ENO Reconstruction" begin include("modified_eno_test.jl") end
   @time @testset "TeCNO" begin include("tecno_test.jl") end
+  @time @testset "Quartic Nonconvex Flux" begin include("balance_laws/quartic_nonconvex_test.jl") end
   @time @testset "Keyfitz-Kranzer System" begin include("balance_laws/keyfitz_kranzer_test.jl") end
 
 end