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 prob, x, y, u = heat_equation()

 65.751945 seconds (65.63 M allocations: 4.330 GiB, 1.50% gc time, 98.66% compilation time: 10% of which was recompilation)

Solve the PDE

@time sol = NonlinearSolve.solve(prob, NewtonRaphson())

 12.156046 seconds (17.23 M allocations: 1.068 GiB, 1.65% gc time, 96.09% compilation time: 2% 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)]

12×12 Matrix{Float64}:
  0.0     -0.0025  -0.0075  -0.0125  …  -0.0375  -0.0425  -0.0475  0.0
 -0.0025   0.0025   0.0075   0.0125      0.0375   0.0425   0.0475  0.0525
 -0.0075   0.0075   0.0225   0.0375      0.1125   0.1275   0.1425  0.1575
 -0.0125   0.0125   0.0375   0.0625      0.1875   0.2125   0.2375  0.2625
 -0.0175   0.0175   0.0525   0.0875      0.2625   0.2975   0.3325  0.3675
 -0.0225   0.0225   0.0675   0.1125  …   0.3375   0.3825   0.4275  0.4725
 -0.0275   0.0275   0.0825   0.1375      0.4125   0.4675   0.5225  0.5775
 -0.0325   0.0325   0.0975   0.1625      0.4875   0.5525   0.6175  0.6825
 -0.0375   0.0375   0.1125   0.1875      0.5625   0.6375   0.7125  0.7875
 -0.0425   0.0425   0.1275   0.2125      0.6375   0.7225   0.8075  0.8925
 -0.0475   0.0475   0.1425   0.2375  …   0.7125   0.8075   0.9025  0.9975
  0.0      0.0525   0.1575   0.2625      0.7875   0.8925   0.9975  0.0

Visualize the solution¶

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.