并行化代码(DifferentialEquations.jl) -朱莉娅

Question

我有以下耦合的PDE系统(来自于离散积分微分PDE)：

​

习是我控制的x网格上的点。我可以用以下简单的代码来解决这个问题：

using DifferentialEquations

function ode_syst(du,u,p, t)
    N = Int64(p[1])
    beta= p[2]
    deltax = 1/(N+1)
    xs = [deltax*i for i in 1:N]
    for j in 1:N
        du[j] = -xs[j]^(beta)*u[j]+deltax*sum([u[i]*xs[i]^(beta) for i in 1:N])
    end
end

N = 1000
u0 = ones(N)
beta = 2.0
p = [N, beta]
tspan = (0.0, 10^3);

prob = ODEProblem(ode_syst,u0,tspan,p);
sol = solve(prob);

然而，随着网格的细化，即增加N，计算时间迅速增长(我猜在N中是二次定标)。对于如何使用分布式并行或多线程来实现这一点，有什么建议吗？

附加信息：--我附加了可能有助于理解程序大部分时间花费在何处的剖析图

​

Answer 1

我查看了您的代码，发现了一些问题，例如由于重新计算和项而意外引入了O(N^2)行为。

我的改进版本使用Tullio包来进一步提高矢量化的速度。Tullio也有可调参数，如果系统变得足够大，就可以进行多线程处理。在这里查看可以调优在“选项”部分的哪些参数。您可能还会看到GPU支持在那里，我没有测试，但它可能会产生进一步的速度或崩溃。我也选择从穴位阵列得到长度，这应该使使用更经济，更容易出错。

using Tullio

function ode_syst_t(du,u,p, t)
    N = length(du)
    beta = p[1]
    deltax = 1/(N+1)
    @tullio s := deltax*(u[i]*(i*deltax)^(beta))
    @tullio du[j] = -(j*deltax)^(beta)*u[j] + s
    return nothing
end

你的代码：

 @btime sol = solve(prob);
  80.592 s (1349001 allocations: 10.22 GiB)

我的代码：

prob2 = ODEProblem(ode_syst_t,u0,tspan,[2.0]);
@btime sol2 = solve(prob2);
  1.171 s (696 allocations: 18.50 MiB)

结果大致一致：

julia> sum(abs2, sol2(1000.0) .- sol(1000.0))
1.079046922815598e-14

Lutz Lehmanns解决方案-我也设定了基准：

prob3 = ODEProblem(ode_syst_lehm,u0,tspan,p);

@btime sol3 = solve(prob3);
  1.338 s (3348 allocations: 39.38 MiB)

然而，当我们用tspan (0.0，10.0)来标度N到1000000时，

prob2 = ODEProblem(ode_syst_t,u0,tspan,[2.0]);

@time solve(prob2);
  2.429239 seconds (280 allocations: 869.768 MiB, 13.60% gc time)

prob3 = ODEProblem(ode_syst_lehm,u0,tspan,p);

@time solve(prob3);
  5.961889 seconds (580 allocations: 1.967 GiB, 11.08% gc time)

由于在我的旧机器和生锈机器中使用了两个核，我的代码变得比以前快了一倍多。

Answer 2

分析公式。显然，原子术语重复。所以它们只能计算一次。

function ode_syst(du,u,p, t)
    N = Int64(p[1])
    beta= p[2]
    deltax = 1/(N+1)
    xs = [deltax*i for i in 1:N]
    term = [ xs[i]^(beta)*u[i] for i in 1:N]
    term_sum = deltax*sum(term)
    for j in 1:N
        du[j] = -term[j]+term_sum
    end
end

这应该只会线性地增加N。