In summary, the steady state heat equation example from the documentation works great at dx == dy == 0.1 but produces incorrect results for every smaller discretization step that I have tried.
Example with dx == dy == 0.09:
Code to produce this:
using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve
using CairoMakie
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
bcs = [u(0, y) ~ x * y,
u(1, y) ~ x * y,
u(x, 0) ~ x * y,
u(x, 1) ~ x * y]
# Space and time domains
domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])
dx, dy = 0.1, 0.1 # works great
#dx, dy = 0.09, 0.09 # instability near (1, 1)
#dx, dy = 0.08, 0.08 # unstable region grows
#dx, dy = 0.05, 0.05 # all NaN on the interior of the domain
iters_power = 5
maxiters = Int(10^iters_power)
# Note that we pass in `nothing` for the time variable `t` here since we
# are creating a stationary problem without a dependence on time, only space.
discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)
prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson(), maxiters=maxiters)
begin
fig = Figure(resolution=(1000, 800))
ax = Axis(fig[1, 1], title="dx = $dx, dy = $dy, maxiters = 1e$iters_power")
h = heatmap!(ax, sol[x], sol[y], sol[u(x, y)], colorrange=(0, 1))
cb = Colorbar(fig[1, 2], h, label="u(x, y)")
fig
end
save("heat_ss_$(dx)_1e$(iters_power).png", fig)
Started from a discourse thread here: https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651
In summary, the steady state heat equation example from the documentation works great at
dx == dy == 0.1
but produces incorrect results for every smaller discretization step that I have tried.Example with
dx == dy == 0.09
:Code to produce this:
Similar issues seem to occur with
chebyspace
:Code (modifications) for chebyspace example: