Блог

Я пытаюсь численно решить уравнение теплопроводности в 1d:

Я использую конечные разности, и у меня возникли некоторые проблемы с использованием инструкции @threads в Julia. В частности, ниже приведены две версии одного и того же кода: первая — single thread, в то время как другая использует @threads (они идентичны, кроме инструкции @thread)

 function heatSecLoop(;T::Float64)

    println("start")
    L = 1
    ν = 0.5
    Δt = 1e-6
    Δx = 1e-3

    Nt = ceil(Int, T/Δt )
    Nx = ceil(Int,L/Δx   2)
    u = zeros(Nx)    
    u[round(Int,Nx/2)] = 1
    
    println("starting loop")
    for t=1:Nt-1
        u_old = copy(u)
        for i=2:Nx-1
            u[i] = u_old[i]   ν * Δt/(Δx^2)*(u_old[i.-1]-2u_old[i]   u_old[i. 1])
        end

        if t % round(Int,Nt/10) == 0
            println("time = " * string(round(t*Δt,digits=4)) )
        end
    end
    println("done")
    return u
end

function heatParLoop(;T::Float64)

    println("start")
    L = 1
    ν = 0.5
    Δt = 1e-6
    Δx = 1e-3

    Nt = ceil(Int, T/Δt )
    Nx = ceil(Int,L/Δx   2)
    u = zeros(Nx)    
    u[round(Int,Nx/2)] = 1
    
    println("starting loop")
    for t=1:Nt-1
        u_old = copy(u)
        Threads.@threads for i=2:Nx-1
            u[i] = u_old[i]   ν * Δt/(Δx^2)*(u_old[i.-1]-2u_old[i]   u_old[i. 1])
        end

        if t % round(Int,Nt/10) == 0
            println("time = " * string(round(t*Δt,digits=4)) )
        end
    end
    println("done")
    return u
end

  

Проблема в том, что последовательный быстрее, чем многопоточный. Вот временные рамки (после одного запуска для компиляции)

 julia> Threads.nthreads()
2

julia> @time heatParLoop(T=1.0)
start
starting loop
time = 0.1
time = 0.2
time = 0.3
time = 0.4
time = 0.5
time = 0.6
time = 0.7
time = 0.8
time = 0.9
done
  5.417182 seconds (12.14 M allocations: 9.125 GiB, 6.59% gc time)

julia> @time heatSecLoop(T=1.0)
start
starting loop
time = 0.1
time = 0.2
time = 0.3
time = 0.4
time = 0.5
time = 0.6
time = 0.7
time = 0.8
time = 0.9
done
  3.892801 seconds (1.00 M allocations: 7.629 GiB, 8.06% gc time)

  

Конечно, уравнение теплопроводности — это всего лишь пример для более сложной задачи. Я также пробовал использовать другие библиотеки, такие как SharedArrays, вместе с Distributed, но с худшими результатами.

Любая помощь приветствуется.

Похоже, это все еще верно, вероятно, из-за

1. накладные расходы на Threads.@threads
2. возможно, в меньшей степени тот факт, что сборка мусора в Julia является однопоточной, и исходная версия здесь генерирует изрядное количество мусора.

Более того, основываясь на совете из связанного потока discourse, возможно, стоит отметить, что теперь существует многопоточная версия макроса @avx (now @turbo ) от LoopVectorization.jl, которая использует очень легкую потоковую обработку от Polyester.jl, и ей удается добиться немного лучшей производительности, несмотря на все еще нетривиальные накладные расходы на потоковую обработку:

 function heatSecLoop(;T::Float64)

    println("start")
    L = 1
    ν = 0.5
    Δt = 1e-6
    Δx = 1e-3

    Nt = ceil(Int, T/Δt )
    Nx = ceil(Int,L/Δx   2)
    u = zeros(Nx)    
    u[round(Int,Nx/2)] = 1
    u_old = similar(u)

    println("starting loop")
    for t=1:Nt-1
        u_old, u = u, u_old
        for i=2:Nx-1
            u[i] = u_old[i]   ν * Δt/(Δx^2)*(u_old[i.-1]-2u_old[i]   u_old[i. 1])
        end

        if t % round(Int,Nt/10) == 0
            println("time = " * string(round(t*Δt,digits=4)) )
        end
    end
    println("done")
    return u
end
function heatVecLoop(;T::Float64)
    println("start")
    L = 1
    ν = 0.5
    Δt = 1e-6
    Δx = 1e-3

    Nt = ceil(Int, T/Δt )
    Nx = ceil(Int,L/Δx   2)
    u = zeros(Nx)
    u[round(Int,Nx/2)] = 1
    u_old = similar(u)

    println("starting loop")
    for t=1:Nt-1
       u_old, u = u, u_old
       @tturbo for i=2:Nx-1
           u[i] = u_old[i]   ν * Δt/(Δx^2)*(u_old[i-1]-2u_old[i]   u_old[i 1])
       end

       if t % round(Int,Nt/10) == 0
           println("time = " * string(round(t*Δt,digits=4)) )
       end
    end
    println("done")
    return u
end

function heatTVecLoop(;T::Float64)
    println("start")
    L = 1
    ν = 0.5
    Δt = 1e-6
    Δx = 1e-3

    Nt = ceil(Int, T/Δt )
    Nx = ceil(Int,L/Δx   2)
    u = zeros(Nx)
    u[round(Int,Nx/2)] = 1
    u_old = similar(u)

    println("starting loop")
    for t=1:Nt-1
       u_old, u = u, u_old
       @tturbo for i=2:Nx-1
           u[i] = u_old[i]   ν * Δt/(Δx^2)*(u_old[i-1]-2u_old[i]   u_old[i 1])
       end

       if t % round(Int,Nt/10) == 0
           println("time = " * string(round(t*Δt,digits=4)) )
       end
    end
    println("done")
    return u
end
  

 julia> @time heatSecLoop(T=1.0)
start
starting loop
time = 0.1
time = 0.2
time = 0.3
time = 0.4
time = 0.5
time = 0.6
time = 0.7
time = 0.8
time = 0.9
done
  1.786011 seconds (114 allocations: 22.094 KiB)

julia> @time heatVecLoop(T=1.0)
start
starting loop
time = 0.1
time = 0.2
time = 0.3
time = 0.4
time = 0.5
time = 0.6
time = 0.7
time = 0.8
time = 0.9
done
  0.314305 seconds (114 allocations: 22.094 KiB)

julia> @time heatTVecLoop(T=1.0)
start
starting loop
time = 0.1
time = 0.2
time = 0.3
time = 0.4
time = 0.5
time = 0.6
time = 0.7
time = 0.8
time = 0.9
done
  0.300656 seconds (114 allocations: 22.094 KiB)
  

Производительность однопоточной @turbo векторизованной версии также, по-видимому, значительно улучшилась с тех пор, как этот вопрос был задан впервые, и производительность многопоточной @tturbo версии, вероятно, продолжит улучшаться для задач большего размера.