I am pretty new to Julia and ultimately I want to transfer my python code for solving a nonlinear PDE into Julia (and see whether I can get the speedup I hope).
The PDE I want to solve is similar to a nonlinear heat equation. In python I have a function which computes dh/dt. I use the method of lines to discretize the vector h(x,t) in spatial domain, and then use a stiff ODE solver scipy.integrate.solve_ivp(… method=‘BDF’) to solve the PDE. I am trying to adopt a similar approach in Julia with the DifferentialEquations.jl package.
The first step I need to do is to figure out a way to implement numerical differentiation using finite difference (I need first and third derivative, with either 2nd or 4th order accuracy) with periodic boundary condition applied to both ends of the vector. I have the following four versions of codes:
function ndiff1(u::AbstractVector, n::Integer, order::Integer, dx::Float64) Nx = length(u) if n==1 if order==2 du = 0.5*(circshift(u,-1) - circshift(u,1)) elseif order==4 du = (circshift(u,2) - circshift(u,-2) + 8*(circshift(u,-1) - circshift(u,1)))/12.0 end elseif n==3 # similar to above with just different stencils end dudx = du/dx^n return dudx end function ndiff2(u::AbstractVector, n::Integer, order::Integer, dx::Float64) Nx = length(u) dudx = similar(u) if n==1 if order==2 for i = 2:Nx-1 dudx[i] = 0.5*(u[i+1] - u[i-1])/dx end dudx = 0.5*(u-u[end])/dx dudx[end] = 0.5*(u-u[end-1])/dx elseif order==4 # similar to above with different stencils end elseif n==3 # similar to above with different stencils end return dudx end using DiffEqOperators const periodic_bc = PeriodicBC(Float64) const D1_O2 = CenteredDifference(1, 2, dx, Nx)*periodic_bc const D1_O4 = CenteredDifference(1, 4, dx, Nx)*periodic_bc const D3_O2 = CenteredDifference(3, 2, dx, Nx)*periodic_bc const D3_O4 = CenteredDifference(3, 2, dx, Nx)*periodic_bc function ndiff3(u::AbstractVector, n::Integer, order::Integer) if n==1 if order==2 dudx = D1_O2*u elseif order==4 dudx = D1_O4*u end elseif n==3 if order==2 dudx = D3_O2*u elseif order==4 dudx = D3_O4*u end end return dudx end function ndiff4(u::AbstractVector, n::Integer, order::Integer, dx::Float64) Nx = length(u) dudx = CenteredDifference(n, order, dx, Nx)*periodic_bc*u return dudx end
When I used @btime to benchmark the above four functions on an array of size 256 (typical for my problem),
ndiff3 is a bit faster than
ndiff2, which is a bit faster than
ndiff4 is close to an order of magnitude slower than the rest. My questions are:
Are there anything noticeable areas for improvement for
ndiff2? That’s similar to how I implemented the finite difference in python (except I try to ‘vectorize’ it in python using numpy without the for loop for
ndiff2) and that’s already close to 2 orders of magnitude of speedup which I’m pretty happy about, but I’m not sure if I there is anything significant I’m missing in Julia.
I try to use the
ndiff4so I don’t have to write everything explicitly from scratch myself. However I’m not sure what’s the best way to implement
CenteredDifference(n, order, dx, Nx)*periodic_bc. If I place them inside the function as in
ndiff4, does Julia need to initialize it every time the function is called (and hence leads to slower speed)? If I place them outside the function as I did for
ndiff3, then obviously it’s very inflexible as I will need to define the problem size
any comment/help is appreciated. thanks