Issue with (2+1)D Diffusion Reaction PDE with DiffEqOperators.jl

Specific Domains Numerics
1

I am studying a (2+1)D diffusion-reaction PDE. I have asked about this before and got plenty of help and that resolved my issues. However, I just realized I was getting the wrong results because the model should be in two spatial dimensions, not one, so I started working on modifying it.

Here’s a minimal non-working example, with a much simpler PDE:

using DifferentialEquations, DiffEqOperators, LinearAlgebra, Symbolics, SparseArrays

# Define FD operator
N = 60
δ = 30.0/N
x = range(0, step=δ, length=N)
y = x
ord_deriv = 2
ord_approx = 4
Δx = CenteredDifference{1}(ord_deriv, ord_approx, δ, N)
Δy = CenteredDifference{2}(ord_deriv, ord_approx, δ, N)
bc = Neumann0BC(δ, 2)
Qx, Qy = MultiDimBC(bc, size(rand(N,N)))

# Define ODE
function myModel!(du,u,p,t)
    c₁ = u[1,:,:]
    c₂ = u[2,:,:]
    du[1,:,:] = - c₂.^2 .+ 2.0*(Δx*Qx*c₁ + Δy*Qy*c₁)
    du[2,:,:] = - c₁ .+ 3.0*(Δx*Qx*c₂ + Δy*Qy*c₂)
end

# Prepare the ODE integrator
tspan = (0.0,0.03)
u0 = rand(Float64, (2, N, N))

# Prepare the Jacobian stuff
du0 = copy(u0)
jac_sparsity = Symbolics.jacobian_sparsity((du,u)->myModel!(du,u,p,0.0),du0,u0)
f = ODEFunction(myModel!;jac_prototype=float.(jac_sparsity))
prob = ODEProblem(f,u0,tspan,p)

# Solve!
sol = solve(prob,TRBDF2(linsolve = KLUFactorization()), reltol=1e-8, abstol=1e-10)

which leads to the following error:

MethodError: no method matching *(::MultiDimDirectionalBC{Float64, RobinBC{Float64, Vector{Float64}}, 1, 2, 1}, ::Matrix{Num})

which comes from 3rd line in myModel!. It seems like a simple issue, it probably has something to do with how I’m indexing u but I can’t quite figure it out. Any advice would be appreciated.

If I just do:

prob = ODEProblem(myModel!,u0,tspan,p)

and solve it (i.e., without the fancy Jacobian stuff), it just takes forever or crashes (with the actual PDE under study, not this minimal one).