2D time-independent heat equation#
From the tutorial
Using MethodOfLines.jl (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}
\]
using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using NonlinearSolve
using Plots
Setup variables and differential operators
@independent_variables x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Differential(y) ∘ Differential(y)
PDE equation
eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
\[ \begin{equation}
\frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y \right) = 0
\end{equation}
\]
Boundary conditions
bcs = [u(0, y) ~ x * y,
u(1, y) ~ x * y,
u(x, 0) ~ x * y,
u(x, 1) ~ x * y
]
\[\begin{split} \begin{align}
u\left( 0, y \right) &= x y \\
u\left( 1, y \right) &= x y \\
u\left( x, 0 \right) &= x y \\
u\left( x, 1 \right) &= x y
\end{align}
\end{split}\]
Space and time domains
domains = [ x ∈ Interval(0.0, 1.0), y ∈ Interval(0.0, 1.0)]
2-element Vector{Symbolics.VarDomainPairing}:
Symbolics.VarDomainPairing(x, 0.0 .. 1.0)
Symbolics.VarDomainPairing(y, 0.0 .. 1.0)
PDE system
@named pdesys = PDESystem(eq, bcs, domains, [x, y], [u(x, y)])
\[ \begin{align}
\frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y \right) &= 0
\end{align}
\]
Discretize the PDE system into an Nonlinear system
Pass nothing to the time parameter
@time prob = let dx=0.1
discretization = MOLFiniteDifference([x=>dx, y=>dx], nothing, approx_order=2, grid_align = MethodOfLines.EdgeAlignedGrid())
prob = discretize(pdesys, discretization)
end
22.309808 seconds (44.76 M allocations: 2.150 GiB, 1.84% gc time, 95.24% compilation time: 18% of which was recompilation)
NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 24-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
⋮
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
Solve the PDE
@time sol = NonlinearSolve.solve(prob, NewtonRaphson())
19.080613 seconds (51.81 M allocations: 2.072 GiB, 2.10% gc time, 80.73% 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.BasicSymbolic{Real}}:
y
xdomain:(-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.0025265 … -0.0718899 0.0; -0.00252433 0.00252612 … 0.04403…
Extract data
discrete_x = sol[x]
discrete_y = sol[y]
u_sol = sol[u(x,y)]
12×12 Matrix{Float64}:
0.0 -0.0025265 -0.0075784 … -0.0411707 -0.0718899 0.0
-0.00252433 0.00252612 0.00757853 0.0421997 0.0440374 0.0625449
-0.00757664 0.00757679 0.0227332 0.128218 0.143295 0.162896
-0.0126244 0.0126246 0.0378801 0.21395 0.238027 0.262007
-0.0176659 0.017666 0.0530093 0.300248 0.332857 0.362694
-0.0226961 0.022696 0.0681054 … 0.388621 0.430461 0.468921
-0.0277091 0.0277087 0.0831491 0.481381 0.531444 0.498799
-0.0326996 0.032699 0.0981221 0.586601 0.715133 1.02745
-0.0376651 0.0376646 0.113016 0.64849 0.71504 0.694732
-0.0426088 0.0426085 0.12784 0.721795 0.801804 0.871943
-0.0475379 0.0475378 0.142619 … 0.797319 0.898437 1.00343
0.0 0.0524621 0.157381 0.857583 0.991195 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.