2D time-independent heat equation

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)
)
_images/bfa2eb1baf05ae6ed59e077792159a46ebe3a6fe2c56745b1062fce2903c0fa1.png

This notebook was generated using Literate.jl.