在有限差分格式上使用 PARFOR

Question

我正在尝试使用有限差分法求解非线性扩散方程。它本质上是两个嵌套的for 循环，内层是网格上的循环，外层在时间上推进方程，有点像这样：

for i=1:(Nt-1)
    for j=2:(Nx-1)

        %%% simple diffusion equation
        c(2,i+1) = c(1,i) + D*dt/dx/dx*(cL(1,j-1)2*cL(1,j)+cL(1,j+1));

    end

    %re-impose BCs
    c(2,1)=1; 
    c(2,Nx)=1; 

    %save into solution matrix
    if mod(i,floor(Nt/sol_no))==0
        sol_c(sol_i+1,:) = c(2,:);
        t_axis(sol_i+1,1) = i*dt;
        sol_i=sol_i+1;
    end

    %flip rows 1 and 2
    c = flipud(c);

end

为了节省内存，c 只是一个 2×Nx 矩阵，我在每次增量时都会不断翻转它，例如，将每 100 个实例保存到一个解矩阵中。

请注意，为了清楚起见，我在这里使用了简单的扩散方程。问题是，实际上我必须使用一个非常精细的网格（理想情况下约为 200,000 个点），以考虑一些快速变化的系数，这些系数是空间的函数。这一点，再加上稳定性准则，对时间步长和总时间的大小施加了严格的限制。

现在，我想优化我的代码，我想出一种方法是在内部循环中使用parfor 而不是for。 c(2,:) 行中的每个元素仅依赖于 c(1,:) 行中的元素，所以这应该可以工作，对吧？但是当我尝试运行它时，它会永远运行并且不会产生任何东西。我没有在 sn-p 中显示它，但是当循环开始时所有变量都已经初始化，所以我不确定是什么导致了问题。我通读了parfor 文档，但没有发现任何有用的东西。我怀疑它一定与切片有关，但我不知道如何解决它。另外，也许我错过了一些其他更有效的方法来解决这个问题，所以任何建议都值得赞赏。

编辑：我正在添加我的代码的要点，精简为裸露的骨头。 parfor 的执行速度比 for 慢得多。

tic
gridpoints = 400000;
time_iter = 200;

c1 = ones(1,gridpoints);
c2 = ones(1,gridpoints);
x = linspace(0,1,gridpoints);
k = 1.0000001;


for j=1:time_iter

    parfor i=2:gridpoints-1
       c2(i) = c1(i-1) + k*c1(i) - c1(i+1);
    end

c1=c2;
end

toc

Answer 1

您可以完全消除内部循环，并使用 Matlab 的矢量化（使用您的第二个示例），这要快得多：

c2(2:end-1) = c1(1:end-2) + k * c1(2:end-1) - c1(3:end);

话虽如此，你甚至可以完全消除 c2：

c(2:end-1) = c(1:end-2) + k * c(2:end-1) - c(3:end);

并且可能想要添加边界条件：

c = [1; c(1:end-2) + k * c(2:end-1) - c(3:end); 1];

您还可以在外循环之外创建一个向量 (idx = [2:gridpoints-1];) 并适当地替换索引以澄清代码。

此外，如果您遇到稳定性问题，隐式方案可能更可取。