Fail at discretization

[guanming-zhang]

Hello, I am a fresher on Julia working on active matter physics. I am trying to simulate a famous equation, conserved model B, in active matter physics(see https://arxiv.org/pdf/1904.01330.pdf, page 4 for more details) using method of lines.However, I got an error in discretization. Here is the error message.

ArgumentError: Differential w.r.t. multiple variables Any[2.608695652173913, 3.4782608695652173, 4.782608695652174, 1.7391304347826086, 4.3478260869565215, 6.956521739130435, 3.9130434782608696, t, 8.26086956521739, 3.0434782608695654, 7.391304347826087, 1.3043478260869565, 0.8695652173913043, 8.695652173913043, 0.43478260869565216, 2.1739130434782608, 6.086956521739131, 6.521739130434782, 5.6521739130434785, 7.826086956521739, 9.130434782608695, 9.565217391304348, 10.0, 5.217391304347826] are not allowed.

I also paste my codes below. I find that if I set K=0 and A,B,T>0, then there is no error. However, there is always an error when K>0. I would appreciate it if you could help me solve the problem.

using Symbolics,ModelingToolkit
using SymbolicUtils
using DomainSets,MethodOfLines, OrdinaryDiffEq

@parameters t,x,y,z
@parameters T,A,B,K
@syms rho(x,y,t)

rhs = (2B*rho(x, y, t)*Differential(x)(rho(x, y, t)) - A*Differential(x)(rho(x, y, t)) - K*(Differential(x)(Differential(x)(Differential(x)(rho(x, y, t)))) + Differential(x)(Differential(y)(Differential(y)(rho(x, y, t))))))*Differential(x)(rho(x, y, t)) + (2B*rho(x, y, t)*Differential(y)(rho(x, y, t)) - A*Differential(y)(rho(x, y, t)) - K*(Differential(y)(Differential(x)(Differential(x)(rho(x, y, t)))) + Differential(y)(Differential(y)(Differential(y)(rho(x, y, t))))))*Differential(y)(rho(x, y, t)) + (2B*(Differential(x)(rho(x, y, t))^2) + 2B*rho(x, y, t)*Differential(x)(Differential(x)(rho(x, y, t))) - A*Differential(x)(Differential(x)(rho(x, y, t))) - K*(Differential(x)(Differential(x)(Differential(x)(Differential(x)(rho(x, y, t))))) + Differential(x)(Differential(x)(Differential(y)(Differential(y)(rho(x, y, t)))))))*rho(x, y, t) + (2B*(Differential(y)(rho(x, y, t))^2) + 2B*rho(x, y, t)*Differential(y)(Differential(y)(rho(x, y, t))) - A*Differential(y)(Differential(y)(rho(x, y, t))) - K*(Differential(y)(Differential(y)(Differential(x)(Differential(x)(rho(x, y, t))))) + Differential(y)(Differential(y)(Differential(y)(Differential(y)(rho(x, y, t)))))))*rho(x, y, t) + T*Differential(x)(Differential(x)(rho(x, y, t))) + T*Differential(y)(Differential(y)(rho(x, y, t)))

rhs = substitute(rhs, Dict([T=>0.5,A=>0.1,B=>0.1,K=>0.01]))
x_min = y_min = t_min = 0.0
x_max = y_max = 10.0
t_max = 0.5
rho0(x,y,t) = sin(x) * cos(y)
domains = [x ∈ Interval(x_min, x_max),y ∈ Interval(y_min, y_max),t ∈ Interval(t_min, t_max)]
bcs = [rho(x,y,0) ~ rho0(x,y,0),
       rho(0,y,t) ~ rho(x_max,y,t),
       rho(x,0,t) ~ rho(x,y_max,t)] 

@named pdesys = PDESystem(Dt(rho(x,y,t))~rhs,bcs,domains,[x,y,t],[rho(x,y,t)])

N = 24

order = 2 

discretization = MOLFiniteDifference([x=>N, y=>N], t, approx_order=order,advection_scheme=WENOScheme())
print(discretization)
println("Discretization:")
@time prob = discretize(pdesys,discretization)

[xtalax]

Hi, I will take a look at this soon, I think it's related to this term

 Differential(y)(Differential(y)(Differential(y)(rho(x, y, t))))))*Differential(y)(rho(x, y, t))

We have only a rudimentary scheme for 3rd order derivatives with high dispersion, and it is possible that the upwind scheme is interacting poorly here, but you could try wrapping this 3rd order deriv with an auxiliary variable.