Skip to content

Instantly share code, notes, and snippets.

@SteffenPL

SteffenPL/Project.toml

Last active September 20, 2023 12:00
Show Gist options
  • Save SteffenPL/478a0545cb7067a2854a887a2adf7f36 to your computer and use it in GitHub Desktop.
Save SteffenPL/478a0545cb7067a2854a887a2adf7f36 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
Project.toml
[deps]
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Raw
RDSolver optimised.jl
using SparseArrays, LinearAlgebra, OrdinaryDiffEq
@assert dimensions == 2
m = 100
L = 30.0
a = 0.01
b = 1.0
Du = 1.0
Dv = 100.0
dx = L / (m-1)
tspan = (0, 100)
function laplace!(dU, U, D, dx, m)
@inbounds for i in 1:m, j in 1:m
i₋ = clamp(i-1,1,m)
i₊ = clamp(i+1,1,m)
j₋ = clamp(j-1,1,m)
j₊ = clamp(j+1,1,m)
dU[i,j] += D * ( (U[i₊,j] + U[i₋,j] - 2*U[i,j])/dx^2 + (U[i,j₊] + U[i,j₋] - 2*U[i,j])/dx^2 )
end
return dU
end
function rhs_opt!(dz, z, p, t)
(;a, b, Du, Dv, m, dx) = p
u = view(z, :, :, 1)
v = view(z, :, :, 2)
du = view(dz, :, :, 1)
dv = view(dz, :, :, 2)
@. du = a - u + u^2 * v
@. dv = b - u^2 * v
laplace!(du, u, Du, dx, m)
laplace!(dv, v, Dv, dx, m)
return nothing
end
u0 = [fill(a+b, m, m) ;;; fill(b/(b+a)^2, m, m)] + 1e-2*rand(m, m, 2)
p = (;a, b, Du, Dv, N, dx, Lap, m)
prob = ODEProblem(rhs_opt!, u0, tspan, p)
@time sol = solve(prob, ROCK2())
f = ODEFunction(rhs_opt!, jac_prototype = pattern)
prob = ODEProblem(f, u0, tspan, p)
@time sol = solve(prob, QNDF(autodiff=false))
# create video
using GLMakie
using GLMakie.Makie.LaTeXStrings: latexstring
t = Observable(0.0)
u = @lift sol($t)[:,:,1]
title = @lift latexstring("t = $(round(Int64,100*$t/tspan[end]))%")
fig = Figure(resolution = (800, 600))
ax = Axis(fig[1, 1], title = lift(t -> latexstring("t = $(round(Int64,100*t/tspan[end]))%"), t))
hm = heatmap!(ax, u, colormap = :jet, colorrange = (0.0, 10.0))
Colorbar(fig[1,2], hm)
record(fig, "RD.gif", LinRange(tspan..., 100)) do i
t[] = i
end
Raw
RDSolver.jl
using SparseArrays, LinearAlgebra, OrdinaryDiffEq
dimensions = 2
m = 100
N = dimensions == 1 ? m : m^2
L = 30.0
a = 0.01
b = 1.0
Du = 1.0
Dv = 100.0
dx = L / (m-1)
tspan = (0, 100)
e = ones(m)
Lap = spdiagm(-1 => e[2:end], 0 => -2*e, 1 => e[2:end])
Lap[1,1] = -1.0
Lap[end,end] = -1.0
if dimensions == 1
Lap = Lap / dx^2
elseif dimensions == 2
Lap = kron(Lap / dx^2, I(m)) + kron(I(m), Lap / dx^2)
end
function rhs!(dz, z, p, t)
(a, b, Du, Dv, N, Lap) = p
u = view(z, 1:N)
v = view(z, N+1:2*N)
du = view(dz, 1:N)
dv = view(dz, N+1:2*N)
@. du = a - u + u^2 * v
@. dv = b - u^2 * v
du .+= Du * (Lap * u)
dv .+= Dv * (Lap * v)
return nothing
end
u0 = [fill(a+b, N); fill(b/(b+a)^2, N)] + 1e-2*rand(2*N)
p = (;a, b, Du, Dv, N, Lap)
pattern = [Lap I(N); I(N) Lap]
f = ODEFunction(rhs!, jac_prototype = pattern)
prob = ODEProblem(f, u0, tspan, p)
@time sol = solve(prob, QNDF(autodiff=false)) # ode15s: QNDF or FBDF
# create video
using GLMakie
using GLMakie.Makie.LaTeXStrings: latexstring
t = Observable(0.0)
u = @lift reshape(sol($t)[1:N], m, m)
title = @lift latexstring("t = $(round(Int64,100*$t/tspan[end]))%")
fig = Figure(resolution = (800, 600))
ax = Axis(fig[1, 1], title = lift(t -> latexstring("t = $(round(Int64,100*t/tspan[end]))%"), t))
hm = heatmap!(ax, u, colormap = :jet, colorrange = (0.0, 10.0))
Colorbar(fig[1,2], hm)
record(fig, "RD.gif", LinRange(tspan..., 100)) do i
t[] = i
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment