SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
160 stars 28 forks source link

Discretization error due to boundary conditions #305

Open dave7895 opened 1 year ago

dave7895 commented 1 year ago

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 ```julia 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?

xtalax commented 1 year ago

You need to add bcs for \omega

e.g. Dx(u(xmin, y, t)) ~ 0 is allowed while Dx(u(xmin, y, t)) ~ 0 is not.

these read identical to me and should work, what error are you getting?

dave7895 commented 1 year ago

so for omega not only IC via time constraint but also for arbitrary time?

ah yes they read identical, that's my mistake. I meant to write that Dx(psi(xmin, y, t)) works while Dy(psi(xmin, y, t)) produces an error

xtalax commented 1 year ago

I'm afraid that only derivs in the constant direction are supported in BCs, sorry about this.