Skip to content
This repository has been archived by the owner on Jun 24, 2022. It is now read-only.

Why this doesn't work #40

Open
kanav99 opened this issue Jul 9, 2020 · 0 comments
Open

Why this doesn't work #40

kanav99 opened this issue Jul 9, 2020 · 0 comments

Comments

@kanav99
Copy link

kanav99 commented Jul 9, 2020

using OrdinaryDiffEq, LinearAlgebra, SparseArrays, SparsityDetection

# 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 _DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 32
const X = reshape([i for i in 1:N for j in 1:N],N,N)
const Y = reshape([j for i in 1:N for j in 1:N],N,N)
const α₁ = 1.0.*(X.>=4*N/5)

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 discretized PDE as an ODE function
function f(u,p,t)
    A = @view  u[:,:,1]
    B = @view  u[:,:,2]
    C = @view  u[:,:,3]
    MyA = My*A
    AMx = A*Mx
    DA = @. _DD*(MyA + AMx)
    dA = @. DA + α₁ - β₁*A - r₁*A*B + r₂*C
    dB = @. α₂ - β₂*B - r₁*A*B + r₂*C
    dC = @. α₃ - β₃*C + r₁*A*B - r₂*C
    cat(dA,dB,dC,dims=3)
end

u0 = zeros(N,N,3)
MyA = zeros(N,N);
AMx = zeros(N,N);
DA = zeros(N,N);
prob = ODEProblem(f,u0,(0.0,10.0))
inp = similar(u0)
out = similar(u0)
function f2(du,u,p,t)
    du .= f(u,p,t)
end
jacobian_sparsity(f2,out,inp,SparsityDetection.Fixed(nothing),0.0)

For some reason, this code doesn't work. P.S - The model is badly designed delibrately
@shashi @ChrisRackauckas

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant