Solving quasi-linear heat diffusion using DiffEqTools.jl and DifferentialEquations.jl

Hi everyone.

I am new to Julia and decided to give a try to the SciML framework.
My goal is first to solve the quasi-linear heat equation using the method of Lines.

\partial_t T = \partial_x (\kappa(x) \,\, \partial_x T(x,t))

In the future, I would like to use ModelingToolkit.jl, but I wanted to start by using only DiffEqTools.jl and DifferentialEquations.jl to better understand what’s going on inside the black box.
I found a similar topic here with what I have in mind but for a constant diffusivity parameter: Solving heat diffusion PDE using DiffEqTools.jl and DifferentialEquations.jl

What I have so far is:

using DiffEqOperators, DifferentialEquations

#=
Solving heat equation with variable diffusity κT:
    ∂tT = ∂x(κT ∂x(T))
=#

const Lx = 1.0  # length of the rod in the x direction
const nx = 40  # number of points in the x direction
const Δx = Lx / (nx - 1)  # grid spacing in the x direction
const x = range(0, length=nx, stop=Lx)  # Define the locations along a gridline.
const u0 = sin.(pi*x/Lx)  # initial condition

Q = Dirichlet0BC(eltype(u0))  # homogeneous dirichlet on the sides
u0_padded = Q * u0

κTx_padded = ones((length(u0_padded))) .* 1.22e-3  # has to be same length as T_padded

∂ut = similar(u0)  # pre-allocate array
nonlinear_diffusion!(∂ut, 1, 1, 4,
κTx_padded, u0_padded, Δx, nx)  # finite differences in space

tspan = (0.0, 10.0)

prob = ODEProblem(∂ut, u0, tspan)
solve(prob)

I got an error when using ODEProblem. What did I miss?

Thank you very much in advance :smiley:

1 Like

∂ut isn’t a function. The f should be a function.

Thank you Chris, it was indeed the problem.

It seems that I didn’t quite understand this f parameter :sweat_smile:

Everything is working now, and here is the solution for people coming across the same problem in the future:

using DiffEqOperators, DifferentialEquations

#=
Solving heat equation with variable diffusity κT:
    ∂tT = ∂x(κT ∂x(T))
=#

const Lx = 1.0  # length of the rod in the x direction
const nx = 40  # number of points in the x direction
const Δx = Lx / (nx - 1)  # grid spacing in the x direction
const x = range(0, length=nx, stop=Lx)  # Define the locations along a gridline.
const u0 = sin.(pi*x/Lx)  # initial condition

Q = Dirichlet0BC(eltype(u0))  # homogeneous dirichlet on the sides
κTx_padded = ones((length(u0) + 2)) .* 1.22e-3  # has to be same length as T_padded

tspan = (0.0, 100.0)

function f(u,p,t)
    nonlinear_diffusion(1, 1, 4,
    κTx_padded, Q * u, Δx, nx)
end

prob = ODEProblem(f, u0, tspan)
sol = solve(prob)

Thanks again!

1 Like