Finite differences with Julia

General Usage Performance
1

Dear Julia users,
Is there a way to improve performance of this kind of code without devectorizing (use loops). All variables are 2d arrays. I solve a PDE problem with finite differences method. Thanks for your help

996988748                divV .=  diff(Vx,dims=1)/dx
992096000                divV .+= diff(Vy,dims=2)/dy
        -
992096000                Exxc .= diff(Vx,dims=1)/dx
495600000                Exxc .-= 1/2*divV
992096000                Eyyc .= diff(Vy,dims=2)/dy
495600000                Eyyc .-= 1/2*divV
        -
507248000                Exyv  .= diff(Vx_exp,dims=2)
        0                Exyv ./= dy
507248000                Exyv .+= diff(Vy_exp,dims=1)
        0                Exyv ./= dx
        0                Exyv .*= 0.5
        -
1984192000         @views Exyc  .= 0.25 * ( Exyv[1:end-1,1:end-1]
        -                                + Exyv[2:end  ,1:end-1]
        -                                + Exyv[1:end-1,2:end  ]
        -                                + Exyv[2:end  ,2:end  ])
2

Don’t use diff, that creates temporary arrays. Don’t use multiple lines of .+=, ./=, etcetera — you want a single fused loop for cache efficiency. Don’t use * and + and / when you want .* and .+ and ./ in order to fuse the loops.

(But honestly, this is the sort of thing that is probably more readable if you just write a loop, and there is no performance reason not to do so in Julia.)

7 Likes
3

Thanks for your answer, i rewrite the code following your advices and it is better. With loops the code is faster. If i understood well, in julia the fortran-like syntax is better than the
matlab-like syntax. Unfortunately, the code i wanted to translate is in matlab :slight_smile:
I really like Julia as a language but fortran is still faster for what i am doing (for now).

  1024000         @views divV .= ((f.Vx[2:end,:] .- f.Vx[1:end-1,:])./dx
        -                       .+(f.Vy[:,2:end] .- f.Vy[:,1:end-1])./dy)
        -
  1024000         @views Exyv .= 0.5 .*( (Vx_exp[:,2:end] .- Vx_exp[:,1:end-1])./dy
        -                          .+  (Vy_exp[2:end,:] .- Vy_exp[1:end-1,:])./dx)
        -
  1024000         @views Exyc  .= 0.25 .* (  Exyv[1:end-1,1:end-1]
        -                                 .+ Exyv[2:end  ,1:end-1]
        -                                 .+ Exyv[1:end-1,2:end  ]
        -                                 .+ Exyv[2:end  ,2:end  ])
        -
        0                Eii2  .= 0.5 .* (Exxc.^2 .+ Eyyc.^2) .+ Exyc.^2
4

Check out this for more information. One example is on a reaction-diffusion equation.

https://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb

2 Likes
5

Thank you very much for help, the last link will help
My code is here GitHub - pnavaro/M2Dpt.jl
I will try to beat fortran and matlab during next few days.

6

I agree this may be clearer & faster with loops, or functions containing the loops. But sometimes it’s good to have something which matches the old code almost line-for-line.

You may want to try @uviews from UnsafeArrays.jl for creating lots of views you’ll only use once.