DiffEqOperators and DifferentialEquations: solving a coupled system of PDEs

[ash]

I am trying to solve a system of coupled PDEs shown below:

and
.

I have written some simple code using DifferentialEquations.jl and DiffEqOperators.jl hoping to solve this system:

using DifferentialEquations
using DiffEqOperators

# 300x300 grid, x \in 0,30 and t \in 0,0.3

# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
bc = Dirichlet0BC(Float64)

#Parameters
f = 1.2
m = 190
q = 0.001
ϵ₁ = 0.01
ϵ = 0.01
K = 1000
Θ = 0.0025
d₀ = 0.01
D₁ = 1
p = [f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ]

function myModel!(du,u,p,t)
   f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ = p
   c₁, c₂ = u
   du[1] = 1/ϵ*(f*c₂.*(q .- c₁)./(q .+ c₁) .+ c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₁.^2) .+
           D₁*Δ*bc*c₁      
   du[2] = c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₂ .+ D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2).*Δ*bc*c₂
end

# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (1, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
solve(prob, Tsit5(), dt = dt)

However, I am greeted with this error:

ERROR: LoadError: MethodError: no method matching iterate(::DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64})

which makes me think I’m doing something wrong with DiffEqOperators, but it seems to work when you just have one equation. Any thoughts?

[ChrisRackauckas]

That looks like a broadcasting issue on your side. It shouldn’t be itereating the DiffEqOperator but it is, so somewhere you have a .* where it should be *.

[ash]

Thanks! Managed to hunt this down and fix it. To further debug the code I’ve simplified it as much as possible to the following:

using DiffEqOperators, OrdinaryDiffEq

# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
const bc = Dirichlet0BC(Float64)

function myModel!(du,u,p,t)
    du[1] = Δ*bc*u[1]
    #du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
    du[2] = Δ*bc*u[2]
end

# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (1, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
solve(prob, Tsit5(), dt = dt)

However I am getting the error:

ERROR: LoadError: MethodError: no method matching *(::GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64}, RobinBC{Float64, Vector{Float64}}}, ::Float64)

It seems to work perfectly fine in the single equation case though. I’ve already tried everything I can think of to debug this but will keep working on it.

[CollinsPeriodLee]

You could try this, it works

using DiffEqOperators, OrdinaryDiffEq

# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
const Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
const bc = Dirichlet0BC(Float64)

function myModel!(du, u,p,t)
    du[1,:] .= Δ* bc * u[1,:]
    #du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
    du[2,:] .= Δ* bc * u[2,:]
end

# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (2, Nₓ))
prob = ODEProblem(myModel!, u0, tspan)
dt = 0.001
solve(prob, Tsit5(), dt = dt)

Since you have 2 PDEs, the initial condition should be u0 = rand(Float64, (2, Nₓ)) and you also need to match du & u in the function mymodel!

function myModel!(du, u,p,t)
    du[1,:] .= Δ* bc * u[1,:]
    #du[2] .= D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ₂*bc*u[2]
    du[2,:] .= Δ* bc * u[2,:]
end

[ash]

That works, thank you!

The full code now is:

using OrdinaryDiffEq
using DiffEqOperators, LinearAlgebra
using Plots

#Parameters
f = 1.2
m = 190
q = 0.001
ϵ₁ = 0.01
ϵ = 0.01
K = 1000
Θ = 0.0001
d₀ = 0.01
D₁ = 1
p = [f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ]

# Define FD operator
Nₓ = 300
Δx = 15.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
bc = Neumann0BC(Δx, 1)


function myModel!(du,u,p,t)
    f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ = p
    c₁ = u[1,:]
    c₂ = u[2,:]
    du[1,:] .= 1/ϵ*(f*c₂.*(q .- c₁)./(q .+ c₁) .+ c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₁.^2) .+ D₁*Δ*bc*c₁      
    du[2,:] .= c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₂ .+ D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2)*Δ*bc*c₂
end

# Grid
# 300x300
# Spatial step: 0.1 s.u.
# Temporal step: 0.001 t.u.
# i.e.,
# Time runs from 0 to 0.3 and space from 0 to 30

# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (2, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
sol = solve(prob, Tsit5(), dt = dt, reltol=1e-8, abstol=1e-8) # Should probably use KenCarp4 but it makes no difference


# Plot
c1 = []
c2 = []
for i in 1:length(sol.t)
    append!(c1, sol.u[i][1,:])
    append!(c2, sol.u[i][2,:])
end
c1r = reshape(c1, (Nₓ, length(sol.t)))
c2r = reshape(c2, (Nₓ,length(sol.t)))
p1 = heatmap(x, sol.t, c1r')
p2 = heatmap(x, sol.t, c2r')
plot(p1, p2, layout = (2,1), fmt=:png)

which runs perfectly fine. Now I just have to figure out why it produces no patterns

Expected solution (for c_1, presumably, they don’t say):

(depending on the value of \Theta, see here).

I get:

I tried playing around with dt and dx, approximation orders, etc. I have double-checked the model from two papers. I guess this will be an adventure.

[ChrisRackauckas]

@ash

sol = solve(prob, Tsit5(), dt = dt, reltol=1e-8, abstol=1e-8) # Should probably use KenCarp4 but it makes no difference

You are still using adaptive time stepping here since you did not set adaptive=false and instead only set the initial dt.

[ash]

Thank you. I had already tried that as well, it requires a much finer dt (~0.00001) to avoid instability and gives the same solution. I am not quite sure what is going on.

[ash]

I also tried enforcing the no flux (Neumann 0) boundary condition on the initial condition but the result didn’t change much.

u0[:,2] .= u0[:,1]
u0[:, Nₓ-1] .= u0[:, Nₓ]

[ChrisRackauckas]

Are you sure they are discretizing the same way? Did you look through their code?

[ash]

They were discretizing quite differently but the issue here specifically is that I was using (1+1)D and they were using (2+1)D, so a very silly mistake.

Even after using (2+1)D however I haven’t been able to fully reproduce the results of this model. I am not sure what the issue is (could be discretization, they use a fully explicit scheme?). With a little effort, I moved my code over to the Gray-Scott model and it works flawlessly; the results perfectly agree with the literature, so I think I’m done with this specific model for now.

Thanks again for the help!

[ChrisRackauckas]

@ash

Even after using (2+1)D however I haven’t been able to fully reproduce the results of this model. I am not sure what the issue is (could be discretization, they use a fully explicit scheme?)

How accurate is their method? Euler with a big dt? If someone uses a very inaccurate method you can still get qualitatively okay results with quantitatively wrong boundaries for the different behaviors. Then no matter how well you try to match you might not, because being more accurate may mean matching their results less . I’ve seen that in bio PDE results before with some silly "run once with a big Euler and no error estimation.

[ash]

They use a fully explicit three-level Du Fort-Frankel scheme with dt = 0.001, which is too large for the explicit integrators I tried in DifferentialEquations.jl (with adaptivity off of course), so I think something might be going on there. The other unknown is the exact nature of the initial condition, they just say “random noise” without giving more details.

Thanks for the advice, this is important to know!

[ziolai]

@ash: thank you so much for sharing a wonderful example! Could you please share a reference for the D_1 / D_2 expression? The Guiu-Souto-2010 paper does not seem to contain it. Thx!

[edit]: I found the expression for D_1 / D_2 is a more recent (2012) paper by the same authors. Apologies for the noice.

[ziolai]

@ash: I am developing a notebook for my student on Turing patterns using your implementation as point of departure. I would like to acknowledge your contribution. Your contact details would be good to have. Thank you.

[ash]

@ziolai My (slightly outdated) website is in my discourse bio, it has all my contact information, thanks!

Several things to note about this model

1. It’s (2+1) D, i was treating it as (1+1)D in this thread which is wrong. I was just confused by the paper.
2. I haven’t been able to reproduce the patterns in the paper exactly. They use an explicit method with a large dt and I was using a far more accurate method but still couldn’t reproduce the thresholds for the phase transitions between dots/labyrinth etc.
3. The Gray-Scott model is much more fun for learning about Turing patterns in my opinion. See here. I have implemented it in this unlisted and very preliminary package, but it is very heavily optimized (see my post history). There is an example in the root of the repository. Use the :NoisePatches init_cond to reproduce mrob’s results.

The code is a little messy and is not documented at the moment, so let me know if you have any questions!