Skip to content

Instantly share code, notes, and snippets.

@jw3126 jw3126/advection.jl
Created Aug 5, 2019

Embed
What would you like to do?
advection problems
# ]add AbstractPlotting
# ]add Makie
using Makie
function step!(u_new, u, o)
# u_new[i] = u[i] + dt*v*(u[i] - u[i-1])/dx
#
for i in reverse(eachindex(u))
s = o.v*o.dt/o.dx
u_new[i] = u[i] + s * (u[i] - get(u, i-1, zero(eltype(u))))
end
u_new
end
function step(u, o)
u_new = similar(u)
step!(u_new, u, o)
end
function advectionplot(u0, o)
t = Node(0.0)
u = Node(copy(u0))
display(Makie.plot(u))
map(t) do _
step!(u.val,u.val,o)
push!(u, u.val)
end
for _ in 1:length(u0)
sleep(o.dt/length(u0))
push!(t, 0.)
end
end
function delta(i, n=1000)
ret = zeros(n)
ret[i] = 1
ret
end
function gauss(i, n=1000)
ret = map(1:n) do idx
x = (idx - i) / 50
exp(-x^2/2)
end
end
# optimal parameters
o = (dt=1, v=-1, dx=1)
u0 = delta(50)
u0 = gauss(500)
advectionplot(u0, o)
# upwind/downwind wrong
o = (dt=5, v=1, dx=5)
u0 = gauss(500)
advectionplot(u0, o)
# numerical diffusion:
# dx too large
o = (dt=1, v=-1, dx=1.1)
u0 = gauss(500)
u0 = delta(500)
advectionplot(u0, o)
# instability:
# dx too small
o = (dt=1, v=-1, dx=0.9)
u0 = gauss(500)
advectionplot(u0, o)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.