ArgumentError: Differential w.r.t. multiple variables Any[t, ...] are not allowed. about methodoflines.jl HOT 2 OPEN

using ModelingToolkit, MethodOfLines, DomainSets;

inch = 2.54e-2
lbs = 0.4536
gravity = 9.80643738977
L = 34inch
T0 = 34lbs * gravity
d = 0.095inch
ρ = 7.9e3 # 7.85-8.05
λ = 0.4
η0 = ρ * π * (d/2)^2 * λ
m = η0 * L
k = 3.5e4
ρ_air = 1.29
c = 0.41
h = 0.01
b = 2L / 3

@parameters t x
@variables u(..) v(..) w(..) η(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs = [ Dt(v(x, t)) ~ ((T0 + k * L * (η0 / η(x, t) - 1)) * Dxx(u(x, t)) / (1 + Dx(u(x, t))^2)^2 - c * ρ_air * d * (v(x, t)^2 + w(x, t)^2)) / η0,
        Dt(u(x, t)) ~ v(x, t),
        Dx(Dt(η(x, t))) ~ -w(x, t) * Dx(η(x, t)),
        Dx(u(x, t))^2 ~ (η0 / η(x, t))^2 - 1];

f(x; b = L/2) = h * min(x/b, (L-x)/(L-b))
bcs = [ u(0, t) ~ 0.0,
        u(L, t) ~ 0.0,
        u(x, 0) ~ f(x; b=b),
        v(x, 0) ~ 0.0,
        η(x, 0) ~ m / (√(b^2+h^2) + √((L-b)^2+h^2) - L),
        Dt(η(x, 0)) ~ 0.0,
        w(0, t) ~ 0.0,
        w(L, t) ~ 0.0,
        w(x, 0) ~ 0.0]

domains = [t in Interval(0.0, 1.0),
           x in Interval(0.0, L)]

@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [u(t,x), v(t,x), w(t,x), η(t,x)])

dx = 0.1;
order = 2;
discretization = MOLFiniteDifference([x=>6], t);
prob = discretize(pdesys, discretization)

1|debug> 
In symbolic_discretize(pdesys, discretization) at /home/cmhyett/.julia/packages/PDEBase/Aqj4G/src/symbolic_discretize.jl:7
 63  
 64      ####
 65      # Loop over equations, Discretizing them and their dependent variables' boundary conditions
 66      ####
 67      for pde in pdeeqs
 68          # Read the dependent variables on both sides of the equation
 69          depvars_lhs = get_depvars(pde.lhs, v.depvar_ops)
 70          depvars_rhs = get_depvars(pde.rhs, v.depvar_ops)
 71          depvars = collect(depvars_lhs ∪ depvars_rhs)
 72          depvars = filter(u -> !any(map(x -> x isa Number, arguments(u))), depvars)
 73  
 74          eqvar = get_eqvar(vareqmap, pde)
 75  
 76          # * Assumes that all variables in the equation have same dimensionality except edgevals
 77          args = ivs(eqvar, v)
 78          indexmap = Dict([args[i] => i for i in 1:length(args)])
 79              # Generate the equations for the interior points
 80          discretize_equation!(disc_state, pde, vareqmap, eqvar, bcmap,
 81                               depvars, s, derivweights, indexmap, discretization)
>82      end

1|julia> disc_state.eqs[11:end]
4-element Vector{Equation}:
 ifelse((w(t))[2] > 0, (w(t))[2]*(-5.789717461787865(η(t))[1] + 5.789717461787865(η(t))[2]), (w(t))[2]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3])) + Differential(0.17271999999999998)(Differential(t)((η(t))[2])) ~ 0
 Differential(0.34543999999999997)(Differential(t)((η(t))[3])) + ifelse((w(t))[3] > 0, (w(t))[3]*(-5.789717461787865(η(t))[2] + 5.789717461787865(η(t))[3]), (w(t))[3]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4])) ~ 0
 ifelse((w(t))[4] > 0, (w(t))[4]*(-5.789717461787865(η(t))[3] + 5.789717461787865(η(t))[4]), (w(t))[4]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5])) + Differential(0.51816)(Differential(t)((η(t))[4])) ~ 0
 Differential(0.6908799999999999)(Differential(t)((η(t))[5])) + ifelse((w(t))[5] > 0, (w(t))[5]*(-5.789717461787865(η(t))[4] + 5.789717461787865(η(t))[5]), (w(t))[5]*(-5.789717461787865(η(t))[5] + 5.789717461787865(η(t))[6])) ~ 0