In [1]:
# import Pkg; Pkg.add("Flux")
using OrdinaryDiffEq, RecursiveArrayTools, LinearAlgebra, Test, SparseArrays, SparseDiffTools, Sundials
using Plots; pyplot()
using BenchmarkTools
using DifferentialEquations, DiffEqSensitivity
using Flux, DiffEqFlux, Optim
using CuArrays
using DataDrivenDiffEq
using ModelingToolkit
using DiffEqSensitivity
Gray Scott Reaction Diffusion¶
$$A'=A+(D_A \Delta^2A-AB^2+F(1-A)\Delta t$$$$B'=B+(D_B \Delta^2B+AB^2-(K+F)B)\Delta t$$
$F(1-A)$ feed at rate f scaled by (1-A) so it doesn't exceend 1.0
$D_A$ & $D_B$ diffusion rates
$\Delta^2A$ & $\Delta^2B$ 2d lablacian functions
$-AB^2$ & $+AB^2$ Reaction the chance one A and two Bs come together A is converted to B
$+F(1-A)$ & $-(k+f)$ Kill
Basic Version¶
In [2]:
# Generate the constants
# p = (0.14, 0.06, 0.035, 0.065) # a,α,ubar,β,D1,D2
p = (0.16, 0.08, 0.035, 0.060) # a,α,ubar,β,D1,D2
N = 100
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
# Ax[2,1] = 2.0
# Ax[end-1,end] = 2.0
# Ay[1,2] = 2.0
# Ay[end,end-1] = 2.0
function basic_version!(dr,r,p,t)
D1, D2, F, K = p
u = r[:,:,1]
v = r[:,:,2]
Du = D1*(Ay*u + u*Ax)
Dv = D2*(Ay*v + v*Ax)
dr[:,:,1] .= Du .- u.*v.*v .+ F .* (1.0 .-u)
dr[:,:,2] .= Dv .+ u.*v.*v .- (F+K) .*v
end
Du, Dv, F, K = p
r0 = zeros(N,N,2)
r0[:,:,1]=ones(N,N)
r0[:,:,1] .+= 0.02 .*rand(N,N)
r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=16
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.5
r0[s:e, s:e,2] .= 0.25
tend=100;
prob = ODEProblem(basic_version!,r0,(0.0,tend),p,saveat=1:tend/20:tend)
sol= solve(prob,Tsit5());
sol.t;
In [3]:
anim = @animate for i=1:length(sol.t)
heatmap(sol[:,:,2,i],c = :redsblues,legend = false)
end
gif(anim, fps = 5)
Out[3]:
Vectorized Version¶
In [5]:
N = 256
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
# Ax[2,1] = 2.0
# Ax[end-1,end] = 2.0
# Ay[1,2] = 2.0
# Ay[end,end-1] = 2.0
Ayu = zeros(N,N)
uAx = zeros(N,N)
Du = zeros(N,N)
Ayv = zeros(N,N)
vAx = zeros(N,N)
Dv = zeros(N,N)
function gm4!(dr,r,p,t)
D1, D2, F, K,Ayu,uAx,Du,Ayv,vAx,Dv = p
u = @view r[:,:,1]
v = @view r[:,:,2]
du = @view dr[:,:,1]
dv = @view dr[:,:,2]
mul!(Ayu,Ay,u)
mul!(uAx,u,Ax)
mul!(Ayv,Ay,v)
mul!(vAx,v,Ax)
@. Du = D1*(Ayu + uAx)
@. Dv = D2*(Ayv + vAx)
@. du = Du - u*v*v + F * (1.0 -u)
@. dv = Dv + u*v*v - (F+K) *v
end
D1, D2, F, K = (0.16, 0.08, 0.06, 0.062)#corals
# D1, D2, F, K = (0.14, 0.06, 0.035, 0.065)#bacteria
# D1, D2, F, K = (0.12, 0.08, 0.020, 0.050)#spirals
# D1, D2, F, K = ( 0.16, 0.08, 0.035, 0.060)#zebrafish
p = (D1, D2, F, K,Ayu,uAx,Du,Ayv,vAx,Dv)
r0 = zeros(N,N,2)
r0[:,:,1]=ones(N,N)
r0[:,:,1] .+= 0.02 .*rand(N,N)
r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=16
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.5
r0[s:e, s:e,2] .= 0.25
tend=32000;
numSaves=20;
prob = ODEProblem(basic_version!,r0,(0.0,tend),p,saveat=0.0:tend/numSaves:tend,progress=true)
# prob = ODEProblem(gm4!,r0,(0.0,tend),p,progress=true,save_everystep=false,save_start=true)
# sol= solve(prob,BS3());
sol= solve(prob,ROCK2());
# sol= solve(prob);
sol.t;
In [6]:
heatmap(sol[:,:,2,end],c = :redsblues)
Out[6]:
In [7]:
dim=2
mi=minimum(sol[:,:,dim,:])
ma=maximum(sol[:,:,dim,:])
anim = @animate for i=1:length(sol.t)
heatmap(sol[:,:,2,i],c = :redsblues,legend = false)
end
gif(anim,"grayscott.gif" ,fps = 5)
Out[7]:
GPU¶
In [10]:
N = 256
r0 = zeros(N,N,2)
r0[:,:,1]=ones(N,N)
r0[:,:,1] .+= 0.02 .*rand(N,N)
r0[:,:,2] .+= 0.02 .*rand(N,N)
N2=convert(Int,N/2)
r=16
s=convert(Int,N2-r)
e=convert(Int,N2+r)
r0[s:e, s:e,1] .= 0.5
r0[s:e, s:e,2] .= 0.25
Ax = Array(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
Ay = copy(Ax)
gAx = CuArray(Float32.(Ax))
gAy = CuArray(Float32.(Ay))
gAyu = CuArray(zeros(Float32,N,N))
guAx = CuArray(zeros(Float32,N,N))
gDu = CuArray(zeros(Float32,N,N))
gAyv = CuArray(zeros(Float32,N,N))
gvAx = CuArray(zeros(Float32,N,N))
gDv = CuArray(zeros(Float32,N,N))
gr0 = CuArray(Float32.(r0))
D1, D2, F, K = (0.16, 0.08, 0.06, 0.062)#corals
D1, D2, F, K = (0.14, 0.06, 0.035, 0.065)#bacteria
# D1, D2, F, K = (0.12, 0.08, 0.020, 0.050)#spirals
# D1, D2, F, K = ( 0.16, 0.08, 0.035, 0.060)#zebrafish
# D1, D2, F, K = ( 0.16, 0.08, 0.025, 0.06)#PULSATING SOLITONS
# D1, D2, F, K = ( 0.16, 0.08, 0.078, 0.061)#WORMS
# D1, D2, F, K = ( 0.16, 0.08, 0.014, 0.054)#MOVING SPOTS
# D1, D2, F, K = ( 0.16, 0.08, 0.018, 0.051)#SPOTS AND LOOPS
# D1, D2, F, K = ( 0.16, 0.08, 0.014, 0.045)#waves
# D1, D2, F, K = ( 0.16, 0.08, 0.062, 0.061)#ukate
# p = (D1, D2, F, K,Ayu,uAx,Du,Ayv,vAx,Dv)
function gff(dr,r,p,t)
# D1, D2, F, K = p
u = @view r[:,:,1]
v = @view r[:,:,2]
du = @view dr[:,:,1]
dv = @view dr[:,:,2]
mul!(gAyu,gAy,u)
mul!(guAx,u,gAx)
mul!(gAyv,gAy,v)
mul!(gvAx,v,gAx)
@. gDu = D1*(gAyu + guAx)
@. gDv = D2*(gAyv + gvAx)
@. du = gDu - u*v*v + F * (1.0 -u)
@. dv = gDv + u*v*v - (F+K) *v
end
CuArrays.allowscalar(false)
tend=32000;
numSaves=20;
prob = ODEProblem(gff,gr0,(0.0,tend),saveat=0.0:tend/numSaves:tend,progress=true)
# sol = solve(prob,ROCK2())
sol = solve(prob,BS3())
sol=Array(sol)
heatmap(sol[:,:,2,end],c = :redsblues)
Out[10]:
In [11]:
dim=2
mi=minimum(sol[:,:,dim,:])
ma=maximum(sol[:,:,dim,:])
anim = @animate for i=1:numSaves
heatmap(sol[:,:,2,i],c = :redsblues,legend = false)
end
gif(anim,"grayscott3.gif" ,fps = 5)
Out[11]:
In [ ]: