The steady-state heat equation example from the MethodOfLines.jl documentation seems to break for smaller step sizes. The example I’m referring to is here: Steady State Heat Equation - No Time Dependence - NonlinearProblem · MethodOfLines.jl
Output for dx == dy == 0.1 (as in the original tutorial):
Output for dx == dy == 0.09:
Output for dx == dy == 0.05 (all NaNs on the interior):
The return code in every PDESolution is MaxIters (including the 0.1 case), so I also tried significantly increasing maxiters from the default. This doesn’t seem to do much:
Is there some reason why this is expected behavior?
Minimal example (essentially identical to the example in the docs, except that I’m using Makie for plotting and varying dx/dy and maxiters):
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)
Thanks in advance for taking a look!
Edit: Updated with more examples.
