I am trying to run the reaction diffusion example in this blog post. Running the code as-is, I notice that I need to include using LinearAlgebra to use the Tridiagonal method. I still get an error UndefVarError: full not defined. It appars that I’m missing some package to run full. I’ve tried searching how to fix this, but most stack-exchange posts arent explicit about package usage. For example SparseArrays doesnt supply full.
If I just take out the full, then I get an error when solving the ODE: sol = solve(prob,BS3()). I have tried other solvers and the issue remains. I am able to run the examples at the end of the readme in OrdinaryDiffEq.jl, so I think the problem is with the problem definition. I am new to julia and trying to learn from these examples.
All packages are up to date
Julia Version 1.2.0
Commit c6da87ff4b (2019-08-20 00:03 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
#######################################################
### Solve the PDE
#######################################################
using OrdinaryDiffEq, RecursiveArrayTools, LinearAlgebra
###
# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const D = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 100
const X = reshape([i for i in 1:100 for j in 1:100],N,N)
const Y = reshape([j for i in 1:100 for j in 1:100],N,N)
const α₁ = 1.0.*(X.>=80)
#const Mx = full(Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1]))
## I don't know which package supplies "full" so I tried running without it
const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0
# Define the initial condition as normal arrays
A = zeros(N,N); B = zeros(N,N); C = zeros(N,N)
u0 = ArrayPartition((A,B,C))
const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N)
# Define the discretized PDE as an ODE function
function f(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(MyA,My,A)
A_mul_B!(AMx,A,Mx)
@. DA = D*(MyA + AMx)
@. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end
# Solve the ODE
prob = ODEProblem(f,u0,(0.0,100.0))
@time sol = solve(prob,BS3(),progress=true,save_everystep=false,save_start=false)
using Plots; pyplot()
p1 = surface(X,Y,sol[end].x[1],title = "[A]")
p2 = surface(X,Y,sol[end].x[2],title = "[B]")
p3 = surface(X,Y,sol[end].x[3],title = "[C]")
plot(p1,p2,p3,layout=grid(3,1))
#######################################################
### Solve the PDE using CLArrays
#######################################################
using CLArrays
gA = CLArray(A); gB = CLArray(B); gC = CLArray(C)
const gMx = CLArray(Mx)
const gMy = CLArray(My)
const gα₁ = CLArray(α₁)
gu0 = ArrayPartition((gA,gB,gC))
const gMyA = CLArray(MyA)
const gAMx = CLArray(AMx)
const gDA = CLArray(DA)
function gf(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
A_mul_B!(gMyA,gMy,A)
A_mul_B!(gAMx,A,gMx)
@. gDA = D*(gMyA + gAMx)
@. dA = gDA + gα₁ - β₁*A - r₁*A*B + r₂*C
@. dB = α₂ - β₂*B - r₁*A*B + r₂*C
@. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end
prob2 = ODEProblem(gf,gu0,(0.0,100.0))
GPUArrays.allowslow(false)
@time sol = solve(prob2,BS3(),progress=true,dt=0.003,adaptive=false,save_everystep=false,save_start=false)
prob2 = ODEProblem(gf,gu0,(0.0,100.0))
sol = solve(prob2,BS3(),progress=true,save_everystep=false,save_start=false)
# Adaptivity currently fails due to https://github.com/JuliaGPU/CLArrays.jl/issues/10
#######################################################
### Solve the SPDE
#######################################################
using StochasticDiffEq
function g(t,u,du)
A,B,C = u.x
dA,dB,dC = du.x
@. dA = γ₁*A
@. dB = γ₂*A
@. dC = γ₃*A
end
prob3 = SDEProblem(f,g,u0,(0.0,100.0))
sol = solve(prob3,SRIW1(),progress=true,save_everystep=false,save_start=false)
p1 = surface(X,Y,sol[end].x[1],title = "[A]")
p2 = surface(X,Y,sol[end].x[2],title = "[B]")
p3 = surface(X,Y,sol[end].x[3],title = "[C]")
plot(p1,p2,p3,layout=grid(3,1))
# Exercise: Do SPDE + GPU
type NullParameters has no field x
Stacktrace:
[1] getproperty(::Any, ::Symbol) at ./Base.jl:20
[2] f(::ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}}, ::DiffEqBase.NullParameters, ::Float64) at ./In[16]:39
[3] (::ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}}, ::DiffEqBase.NullParameters, ::Vararg{Any,N} where N) at /home/user/.julia/packages/DiffEqBase/jMOwn/src/diffeqfunction.jl:219
[4] initialize!(::OrdinaryDiffEq.ODEIntegrator{BS3,false,ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1},ODESolution{Float64,2,Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1},Nothing,Nothing,Array{Float64,1},Array{Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1},1},ODEProblem{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},BS3,OrdinaryDiffEq.InterpolationData{ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1},Array{Float64,1},Array{Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1},1},OrdinaryDiffEq.BS3ConstantCache{Float64,Float64}},DiffEqBase.DEStats},ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.BS3ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Float64,Nothing}, ::OrdinaryDiffEq.BS3ConstantCache{Float64,Float64}) at /home/user/.julia/packages/OrdinaryDiffEq/ROQW3/src/perform_step/low_order_rk_perform_step.jl:4
[5] #__init#329(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Float64, ::Float64, ::Bool, ::Bool, ::Rational{Int64}, ::Nothing, ::Nothing, ::Rational{Int64}, ::Int64, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::BS3, ::Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/user/.julia/packages/OrdinaryDiffEq/ROQW3/src/solve.jl:356
[6] __init(::ODEProblem{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::BS3, ::Array{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/user/.julia/packages/OrdinaryDiffEq/ROQW3/src/solve.jl:66 (repeats 4 times)
[7] #__solve#328 at /home/user/.julia/packages/OrdinaryDiffEq/ROQW3/src/solve.jl:4 [inlined]
[8] __solve at /home/user/.julia/packages/OrdinaryDiffEq/ROQW3/src/solve.jl:4 [inlined]
[9] #solve_call#442(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.solve_call), ::ODEProblem{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::BS3) at /home/user/.julia/packages/DiffEqBase/jMOwn/src/solve.jl:40
[10] solve_call at /home/user/.julia/packages/DiffEqBase/jMOwn/src/solve.jl:37 [inlined]
[11] #solve#443 at /home/user/.julia/packages/DiffEqBase/jMOwn/src/solve.jl:57 [inlined]
[12] solve(::ODEProblem{ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::BS3) at /home/user/.julia/packages/DiffEqBase/jMOwn/src/solve.jl:45
[13] top-level scope at In[18]:3