Discretization error due to boundary conditions

[dave7895]

While trying to replicate a textbook example of a driven cavity in fluid dynamics with this package I encountered a Discretization Error. When further simplifying the governing equations, the error remained.

Erroring code

println("Imports")
@time using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, DifferentialEquations

@parameters x y t ν vel0
@variables ω(..) ψ(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

∆(u) = Dxx(u) + Dyy(u)

xmin = ymin = tmin = 0.0
xmax = ymax = 1.0

tmax = 10.0

eq = [Dt(ω(x, y, t)) ~ ∆(ω(x, y, t)),
   ∆(ψ(x, y, t)) ~ ω(x, y, t)
]

domains = [x ∈ Interval(xmin, xmax),
   y ∈ Interval(ymin, ymax),
   t ∈ Interval(tmin, tmax)
]

boundary = [
   # stationary walls
#   Dy(ψ(xmin, y, t)) ~ 0.0,
   Dx(ψ(xmin, y, t)) ~ 0.0,
#   Dy(ψ(xmax, y, t)) ~ 0.0,
   Dx(ψ(xmax, y, t)) ~ 0.0,
   Dy(ψ(x, ymax, t)) ~ 0.0,
#   Dx(ψ(x, ymax, t)) ~ 0.0,
   # moving wall
   Dy(ψ(x, ymin, t)) ~ 0.0,
#   Dx(ψ(x, ymin, t)) ~ vel0,
   # time boundary
   ψ(x, y, tmin) ~ 0.0,
   ω(x, y, tmin) ~ 0.0
]

@named pde = PDESystem(eq, boundary, domains, [x, y, t], [ω(x, y, t), ψ(x, y, t)], [vel0 => 1.0, ν => 0.1])

Nx = Ny = 6
dt = 0.1
Nt = round(Int, 1/dt)

order = 2

discretization = MOLFiniteDifference([x => Nx, y => Ny], t, approx_order=order)
@show discretization
println("Discretization:")
@time problem = discretize(pde, discretization)

I narrowed it down to the boundary conditions. They apparently may not contain a derivative w.r.t. to a variable that is not fixed. e.g. Dx(u(xmin, y, t)) ~ 0 is allowed while Dx(u(xmin, y, t)) ~ 0 is not. Is there any way around this or is there no way to specify sth like a No-slip boundary using a stream function where the relevant velocities/BCs are defined tangentially to a boundary?