Unstable steady-state solutions for smaller grid sizes in heat equation example

[thomasteisberg]

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:

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)

Similar issues seem to occur with chebyspace:

Code (modifications) for chebyspace example:

cheb_n = 14
discx = chebyspace(cheb_n, domains[1])
discy = chebyspace(cheb_n, domains[2])
discretization = MOLFiniteDifference([discx, discy], 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="chebyspace points per axis = $cheb_n, 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_cheb$(cheb_n)_1e$(iters_power).png", fig)

[PeterARBork]

+1 For this problem to be solved. Would really like to use MethodOfLines.jl for my problem, which is very similar to the one described here.