I am currently bug-fixing a more complicated set of equations I am trying to solve using MethodOfLines.jl and I have narrowed down the current issue I am facing to a particular term which involves a perpendicular Neumann boundary condition.
Here is a simple example using a diffusion equation which reproduces my issue:
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
@parameters t x y
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dy = Differential(y)
Dyy = Differential(y)^2
# 2D PDE and boundary conditions
eq = Dt(u(t, x, y)) ~ 0.1 * (Dxx(u(t, x, y)) + Dyy(u(t, x, y)))
bcs = [u(0, x, y) ~ 0.0,
u(t, 0, y) ~ 0.0,
u(t, 1, y) ~ 0.0,
u(t, x, 0) ~ 0.0,
Dx(u(t, x, 1)) ~ sin(2 * pi * x)] #<-- This term produces the issue
# Space and time domains
domains = [t ∈ Interval(0.0, 10.0),
x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
# Method of lines discretization
discretization = MOLFiniteDifference([x => 0.1, y => 0.1], t, approx_order=2)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)
This results in a large printout of the system ending by saying the Discretization failed and then gives the following error:
ERROR: LoadError: ArgumentError: Differential w.r.t. multiple variables Any[0.4, 0.7, 0.3, 0.5, 0.2, t, 0.9, 0.1, 0.8, 0.6] are not allowed.
Have I made a mistake in my implementation here, or is there otherwise a better way to use such a boundary condition? Thanks!