MethodOfLines issue with coupled non-linear PDE system

[RayanJobs]

Hi. I am trying to solve a highly-coupled system of nonlinear PDEs. I have attached the code, and I am happy to answer any questions. I have a 2D domain (x,z) and time dependence. My two variables are c (x,z,t), a 2D concentration, and h(x,t), a height. Periodicity is assumed in x and the z domain is transformed onto a fixed domain (such that it is fixed between 0 and 1 and the interface does not move). I keep having issues (current issue is ERROR: LoadError: AssertionError: islinear
Stacktrace). I am not sure why this is happening, in particular because my system of equations will include an additional pde later on.

hb = 0.5
s = 3.0


Ma = 0.01
Ma_s = 0.01
eps = 0.1
cppp = 0.01
vv = 0.1
Pe = 10000.0
α = 0.3

using Plots, ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets


@parameters x z t 
@variables h(..) c(..)
Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2
Dxxx = Differential(x)^3
Dxxxx = Differential(x)^4
Dxz = Differential(x) * Differential(z)  # Mixed derivative, d^2(dx dz)

brusselator_f(x, z) =  (1 + cppp * cos(α * x) * exp(-(z.-hb)^2/(2*vv^2))) * ((z .<= 0.5 - 1/s) * 1.0 +  (0.5 .- 1/s .< z .< 0.5 + 1/s) * exp(-1/(1-((z .- (hb - 1.0/s))/(.0/s)) ))/(exp(-1/(1-((z .- (hb - 1.0/s))/(.0/s)))) + exp(-1/((z .- (hb - 1.0/s))/(.0/s)))) + (z .>= 0.5 + 1/s) * 0.0)




#plot(z, brusselator_f(0,z,0), label="T(z)", xlabel="z", ylabel="f(z)", title="Plot of initial condition")






#∇²(u) = Dxx(u) + Dyy(u)

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 10.0


eqs = [Dt(c(x,z,t)) ~ 1/Pe*((2*z*(1/h(x,t)*Dx(h(x,t)))^2-z/h(x,t)*Dxx(h(x,t)))*Dz(c(x,z,t))+(z/h(x,t)*Dx(h(x,t)))^2*Dzz(c(x,z,t))-2*z/h(x,t)*Dx(h(x,t))*Dxz(c(x,z,t))+ Dxx(c(x,z,t))) + 1/(eps^2 * Pe * (h(x,t))^2) * Dzz(c(x,z,t)), 
Dt(h(x,t)) ~ -(h(x,t))^2 * Dx(h(x,t))* Dxxx(h(x,t)) - 1/3 * h(x,t)^3 * Dxxxx(h(x,t)) + h(x,t) * Dx(h(x,t)) * Dx(c(x,1.0, t)) + 1/2 * (h(x,t))^2 *  Dxx(c(x,1.0, t))]

domains = [x ∈ Interval(x_min, x_max),
              z ∈ Interval(y_min, y_max),
              t ∈ Interval(t_min, t_max)]

# Periodic BCs
bcs = [h(x,0.0) ~ 1.0,
       h(0.0,t) ~ h(1.0,t),
       
       c(x,z,0.0) ~ brusselator_f(x, z),
       c(0.0,z,t) ~ c(1.0,z,t),
       Dz(c(x,1.0,t)) ~ 0.0,
       Dz(c(x,0.0,t)) ~ 0.0] 

# Initial condition above models a disturbance

@named pdesys = PDESystem(eqs,bcs,domains,[x,z,t],[c(x,z,t),h(x,t)])

N = 100
order = 2 # This may be increased to improve accuracy of some schemes


discretization = MOLFiniteDifference([x=>N, z=>N], t, approx_order=order)





println("Discretization:")
@time prob = discretize(pdesys,discretization)

println("Solving RN:")
@time sol = solve(prob, TRBDF2(), saveat=0.1)

using JLD2
@save "new_solutionn.jld2" sol