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
dx == dy == 0.1 (as in the original tutorial):
dx == dy == 0.09:
dx == dy == 0.05 (all
NaNs on the interior):
The return code in every
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
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.