Engee documentation

Solving the heat equation with diffusion-implicit time-stepping

In this tutorial, we’ll be solving the heat equation:

²

with boundary conditions: . We’ll solve these equations numerically using Finite Difference Method on cell faces. The same exercise could easily be done on cell centers.

Code loading and parameters

First, we’ll use / import some packages:

import Plots
using LinearAlgebra
using DiffEqBase
using OrdinaryDiffEq: SplitODEProblem, solve, IMEXEuler
import SciMLBase

Next, we’ll define some global problem parameters:

a, b, n = 0, 1, 10               # zmin, zmax, number of cells
n̂_min, n̂_max = -1, 1            # Outward facing unit vectors
α = 100;                        # thermal diffusivity, larger means more stiff
β, γ = 10000, π;                # source term coefficients
Δt = 1000;                      # timestep size
N_t = 10;                       # number of timesteps to take
FT = Float64;                   # float type
Δz = FT(b - a) / FT(n)
Δz² = Δz^2;
∇²_op = [1 / Δz², -2 / Δz², 1 / Δz²]; # interior Laplacian operator
∇T_bottom = 10;                 # Temperature gradient at the top
T_top = 1;                      # Temperature at the bottom
S(z) = β * sin(γ * z)               # source term, (sin for easy integration)
zf = range(a, b, length = n + 1);   # coordinates on cell faces
0.0:0.1:1.0

Derivation of analytic solution

Here, we’ll derive the analytic solution:

²²

Apply bottom boundary condition:

Apply top boundary condition:

And now let’s define this in a Julia function:

function T_analytic(z) # Analytic steady state solution
    c1 = ∇T_bottom - β * cos(γ * a) / (γ * α)
    c2 = T_top - (β * sin(γ * b) / (γ^2 * α) + c1 * b)
    return β * sin(γ * z) / (γ^2 * α) + c1 * z + c2
end
T_analytic (generic function with 1 method)

Derive the temporal discretization

Here, we’ll derive the matrix form of the temporal discretization we wish to use (diffusion-implicit and explicit Euler):

²²²²

Note that, since the ² reaches to boundary points, we’ll need to modify the stencils to account for boundary conditions.

Derive the finite difference stencil

For the interior domain, a central and second-order finite difference stencil is simply:

²²²²²²

At the boundaries, we need to modify the stencil to account for Dirichlet and Neumann BCs. Using the following index denotation:

  • i first interior index

  • b boundary index

  • g ghost index

the Dirichlet boundary stencil & source:

²²²

and Neumann boundary stencil & source:

²²

Define the discrete diffusion operator

# Initialize interior and boundary stencils:
∇² = Tridiagonal(ones(FT, n) .* ∇²_op[1],
                 ones(FT, n + 1) .* ∇²_op[2],
                 ones(FT, n) .* ∇²_op[3]);

# Modify boundary stencil to account for BCs

∇².d[1] = -2 / Δz²
∇².du[1] = +2 / Δz²

# Modify boundary stencil to account for BCs
∇².du[n] = 0  # modified stencil
∇².d[n + 1] = 0 # to ensure `∂_t T = 0` at `z=zmax`
∇².dl[n] = 0  # to ensure `∂_t T = 0` at `z=zmax`
D = α .* ∇²
11×11 Tridiagonal{Float64, Vector{Float64}}:
 -20000.0   20000.0        ⋅         ⋅   …        ⋅         ⋅         ⋅    ⋅
  10000.0  -20000.0   10000.0        ⋅            ⋅         ⋅         ⋅    ⋅
       ⋅    10000.0  -20000.0   10000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅    10000.0  -20000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅    10000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅   …        ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅       10000.0        ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅      -20000.0   10000.0        ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅       10000.0  -20000.0   10000.0   ⋅
       ⋅         ⋅         ⋅         ⋅            ⋅    10000.0  -20000.0  0.0
       ⋅         ⋅         ⋅         ⋅   …        ⋅         ⋅        0.0  0.0

Define boundary source

Here, we’ll compute the boundary source Missing close brace }{Δz²}\right)]

AT_b = zeros(FT, n + 1);
AT_b[1] = α * 2 / Δz * ∇T_bottom * n̂_min;
AT_b[end - 1] = α * T_top / Δz²;
9999.999999999998

Set initial condition

Let’s just initialize the solution to 1, and also set the top boundary condition:

T = zeros(FT, n + 1);
T .= 1;
T[n + 1] = T_top; # set top BC
1

Define right-hand side sources

Here, we define the right-hand side (RHS) sources:

function rhs!(dT, T, params, t)
    n = params.n
    i = 1:n # interior domain
    dT[i] .= S.(zf[i]) .+ AT_b[i]
    return dT
end;
rhs! (generic function with 1 method)

Next, we’ll package up parameters needed in the RHS function, define the ODE problem, and solve.

params = (; n)

tspan = (FT(0), N_t * FT(Δt))

prob = SplitODEProblem(SciMLBase.DiffEqArrayOperator(D),
                       rhs!,
                       T,
                       tspan,
                       params)
alg = IMEXEuler()
println("Solving...")
sol = solve(prob,
            alg,
            dt = Δt,
            saveat = range(FT(0), N_t * FT(Δt), length = 5),
            progress = true,
            progress_message = (dt, u, p, t) -> t);
retcode: Success
Interpolation: 1st order linear
t: 5-element Vector{Float64}:
     0.0
  2500.0
  5000.0
  7500.0
 10000.0
u: 5-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 [22.578070798888803, 23.56410096026957, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34220342338085, 15.737763609271497, 11.317321881651878, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]

Visualizing results

Now, let’s visualize the results of the solution and error:

T_end = sol.u[end]

p1 = Plots.plot(zf, T_analytic.(zf), label = "analytic", markershape = :circle,
                markersize = 6)
p1 = Plots.plot!(p1, zf, T_end, label = "numerical", markershape = :diamond)
p1 = Plots.plot!(p1, title = "T ∈ cell faces")

p2 = Plots.plot(zf, abs.(T_end .- T_analytic.(zf)), label = "error", markershape = :circle,
                markersize = 6)
p2 = Plots.plot!(p2, title = "T ∈ cell faces")

Plots.plot(p1, p2)