Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Quartic nonconvex #55

Merged
merged 4 commits into from Mar 13, 2020
Merged

Quartic nonconvex #55

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6c96757
WIP> quartic nonconvex flux
ranocha Mar 13, 2020
7661426
test for quartic nonconvex periodic semidiscretisation
ranocha Mar 13, 2020
107dc95
test printing of some semidiscretisations
ranocha Mar 13, 2020
223cd4d
add tests for quartic nonconvex flux
ranocha Mar 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 3 additions & 1 deletion src/SummationByPartsOperators.jl
View file
Open in desktop
Expand Up @@ -65,6 +65,7 @@ include("conservation_laws/burgers.jl")
include("conservation_laws/cubic.jl")
include("conservation_laws/variable_linear_advection.jl")
include("second_order_eqs/wave_eq.jl")
include("conservation_laws/quartic_nonconvex.jl")


# exports
Expand Down Expand Up @@ -97,6 +98,7 @@ export Tadmor1989, MadayTadmor1989, Tadmor1993,
export BurgersPeriodicSemidiscretisation, BurgersNonperiodicSemidiscretisation,
CubicPeriodicSemidiscretisation, CubicNonperiodicSemidiscretisation,
VariableLinearAdvectionNonperiodicSemidiscretisation,
WaveEquationNonperiodicSemidiscretisation
WaveEquationNonperiodicSemidiscretisation,
QuarticNonconvexPeriodicSemidiscretisation

end # module
99 changes: 99 additions & 0 deletions src/conservation_laws/quartic_nonconvex.jl
View file
Open in desktop
@@ -0,0 +1,99 @@

"""
QuarticNonconvexPeriodicSemidiscretisation{T,Derivative,Dissipation}

A semidiscretisation of the quartic nonconvex conservation law
\$\\partial_t u(t,x) + \\partial_x ( u(t,x)^4 - 10 u(t,x)^2 + 3 u(t,x) ) = 0\$
with periodic boundary conditions.
"""
struct QuarticNonconvexPeriodicSemidiscretisation{T,Derivative<:AbstractDerivativeOperator{T},
Dissipation,
SplitForm<:Union{Val{false}, Val{true}}} <: AbstractSemidiscretisation
derivative::Derivative
dissipation::Dissipation
tmp1::Vector{T}
tmp2::Vector{T}
split_form::SplitForm

function QuarticNonconvexPeriodicSemidiscretisation(derivative::Derivative, dissipation::Dissipation, split_form::SplitForm=Val{false}()) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}}
if dissipation != nothing
@argcheck size(derivative) == size(dissipation) DimensionMismatch
@argcheck grid(derivative) == grid(dissipation) ArgumentError
end
N = size(derivative, 2)
tmp1 = Array{T}(undef, N)
tmp2 = Array{T}(undef, N)
new{T,Derivative,Dissipation,SplitForm}(derivative, dissipation, tmp1, tmp2, split_form)
end
end


function Base.show(io::IO, semidisc::QuarticNonconvexPeriodicSemidiscretisation)
print(io, "Semidiscretisation of the quartic nonconvex conservation law\n")
print(io, " \$ \\partial_t u(t,x) + \\partial_x ( u(t,x)^4 - 10 u(t,x)^2 + 3 u(t,x) ) = 0 \$ \n")
print(io, "with periodic boundaries using")
if semidisc.split_form == Val{true}()
print(io, " a split form and: \n")
else
print(io, " no split form and: \n")
end
print(io, semidisc.derivative)
print(io, semidisc.dissipation)
end


function (disc::QuarticNonconvexPeriodicSemidiscretisation)(du, u, p, t)
@unpack tmp1, tmp2, derivative, dissipation, split_form = disc
@boundscheck begin
@argcheck length(u) == length(tmp1)
@argcheck length(u) == length(tmp2)
@argcheck length(u) == length(du)
end

# volume terms
if typeof(split_form) <: Val{true}
m_2_5 = -2*one(eltype(u)) / 5
m_20_3 = -20*one(eltype(u)) / 3
m_3 = -3*one(eltype(u))

# # D u^4
@. tmp2 = u^4
mul!(tmp1, derivative, tmp2)
@. du = m_2_5 * tmp1

@. tmp2 = u^3
mul!(tmp1, derivative, tmp2)
@. du += m_2_5 * u * tmp1

@. tmp2 = u^2
mul!(tmp1, derivative, tmp2)
@. du += m_2_5 * u^2 * tmp1

mul!(tmp1, derivative, u)
@. du += m_2_5 * u^3 * tmp1

# D u^2
@. tmp2 = u^2
mul!(tmp1, derivative, tmp2)
@. du -= m_20_3 * tmp1

mul!(tmp1, derivative, u)
@. du -= m_20_3 * u * tmp1

# D u
mul!(tmp1, derivative, u)
@. du += m_3 * tmp1
else
@. tmp2 = u^4 - 10 * u^2 + 3 * u
mul!(tmp1, derivative, tmp2)
@. du = -tmp1
end

# dissipation
if dissipation != nothing
mul!(tmp1, dissipation, u)
@. du += tmp1
end

nothing
end
6 changes: 5 additions & 1 deletion test/conservation_laws/burgers_test.jl
View file
Open in desktop
Expand Up @@ -6,11 +6,12 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
N = 2^6
u0func = sinpi
tspan = (zero(T), one(T))

# Fourier
let D = fourier_derivative_operator(xmin, xmax, N)
Di = dissipation_operator(TadmorWaagan2012Standard(), D)
semidisc = BurgersPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -21,6 +22,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
D = periodic_derivative_operator(1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = BurgersPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -30,6 +32,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
let D = legendre_derivative_operator(xmin, xmax, N)
Di = nothing
semidisc = BurgersNonperiodicSemidiscretisation(D, Di, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -40,6 +43,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
D = derivative_operator(MattssonSvärdNordström2004(), 1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = BurgersNonperiodicSemidiscretisation(D, Di, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand Down
4 changes: 4 additions & 0 deletions test/conservation_laws/cubic_test.jl
View file
Open in desktop
Expand Up @@ -12,6 +12,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
let D = fourier_derivative_operator(xmin, xmax, N)
Di = dissipation_operator(TadmorWaagan2012Standard(), D)
semidisc = CubicPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -25,6 +26,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
D = periodic_derivative_operator(1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = CubicPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -37,6 +39,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
let D = legendre_derivative_operator(xmin, xmax, N)
Di = nothing
semidisc = CubicNonperiodicSemidiscretisation(D, Di, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -50,6 +53,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
D = derivative_operator(MattssonSvärdNordström2004(), 1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = CubicNonperiodicSemidiscretisation(D, Di, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand Down
37 changes: 37 additions & 0 deletions test/conservation_laws/quartic_nonconvex_test.jl
View file
Open in desktop
@@ -0,0 +1,37 @@
using Test, SummationByPartsOperators
using OrdinaryDiffEq, DiffEqCallbacks

for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
xmin = T(-1)
xmax = T(1)
N = 2^6
u0func = sinpi
tspan = (zero(T), one(T)/50)

# Fourier
let D = fourier_derivative_operator(xmin, xmax, N)
Di = dissipation_operator(TadmorWaagan2012Standard(), D)
semidisc = QuarticNonconvexPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))

saving = SavingCallback(semidisc, saveat=range(tspan..., length=10))
sol = solve(ode, SSPRK104(), dt=1/20N, save_everystep=false, callback=saving)
end

# Periodic FD
for acc_order in 2:2:8
D = periodic_derivative_operator(1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = QuarticNonconvexPeriodicSemidiscretisation(D, Di, split_form)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))

saving = SavingCallback(semidisc, saveat=range(tspan..., length=10))
sol = solve(ode, SSPRK104(), dt=1/20N, save_everystep=false, callback=saving)
end
end
4 changes: 3 additions & 1 deletion test/conservation_laws/variable_linear_advection_test.jl
View file
Open in desktop
Expand Up @@ -7,11 +7,12 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
u0func = sinpi
tspan = (zero(T), one(T))
afunc = one

# Legendre
let D = legendre_derivative_operator(xmin, xmax, N)
Di = nothing
semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand All @@ -22,6 +23,7 @@ for T in (Float32, Float64), split_form in (Val{true}(), Val{false}())
D = derivative_operator(MattssonSvärdNordström2004(), 1, acc_order, xmin, xmax, N)
Di = dissipation_operator(D)
semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, zero, zero)
println(devnull, semidisc)
ode = semidiscretise(u0func, semidisc, tspan)
du = similar(ode.u0)
semidisc(du, ode.u0, nothing, first(tspan))
Expand Down
1 change: 1 addition & 0 deletions test/runtests.jl
View file
Open in desktop
Expand Up @@ -12,6 +12,7 @@ using Test
include("conservation_laws/burgers_test.jl")
include("conservation_laws/cubic_test.jl")
include("conservation_laws/variable_linear_advection_test.jl")
include("conservation_laws/quartic_nonconvex_test.jl")
end
@time @testset "Second Order Equations" begin
include("second_order_eqs/wave_eq_test.jl")
Expand Down