2D time-independent heat equation

From the tutorial

Using MethodOfLines.jl (https://github.com/SciML/MethodOfLines.jl/) to symbolically define the PDE system and use the finite difference method (FDM) to solve the following PDE:

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}

(1)

using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using NonlinearSolve
using Plots

Setup variables and differential operators

function heat_equation(;dx=0.1)
    @independent_variables 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
    ]
    domains = [ x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]
    @named pdesys = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
    # Discretize the PDE system into an Nonlinear system
    # Pass `nothing` to the time parameter
    discretization = MOLFiniteDifference([x=>dx, y=>dx], nothing, approx_order=2, grid_align = MethodOfLines.EdgeAlignedGrid())
    prob = discretize(pdesys, discretization)
    return (; prob, x, y, u)
end

@time "Build problem" prob, x, y, u = heat_equation()
@time "Solve problem" sol = NonlinearSolve.solve(prob, NewtonRaphson())

Build problem: 115.609616 seconds (230.90 M allocations: 11.105 GiB, 1.63% gc time, 99.25% compilation time: 5% of which was recompilation)
Solve problem: 10.611527 seconds (20.83 M allocations: 951.809 MiB, 0.92% gc time, 96.90% compilation time: <1% of which was recompilation)

retcode: Stalled
Interpolation: Dict{Symbolics.Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}
ivs: 2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}:
 x
 ydomain:(-0.05:0.1:1.05, -0.05:0.1:1.05)
u: Dict{Symbolics.Num, Matrix{Float64}} with 1 entry:
  u(x, y) => [0.0 -0.0025 … -0.0475 0.0; -0.0025 0.0025 … 0.0475 0.0525; … ; -0…

Extract data

discrete_x = sol[x]
discrete_y = sol[y]
u_sol = sol[u(x,y)]

heatmap(
    discrete_x, discrete_y, u_sol,
    xlabel="x values", ylabel="y values", aspect_ratio=:equal,
    title="Steady State Heat Equation", xlims=(0, 1), ylims=(0, 1)
)

This notebook was generated using Literate.jl.