python - 更快的多重网格泊位方程求解器？

Question

我正在尝试制作自己的 CFD 求解器，而计算成本最高的部分之一是求解压力项。一种更快求解泊松微分方程的方法是使用多重网格方法。基本的递归算法是：

function phi = V_Cycle(phi,f,h)
    % Recursive V-Cycle Multigrid for solving the Poisson equation (\nabla^2 phi = f) on a uniform grid of spacing h

    % Pre-Smoothing
    phi = smoothing(phi,f,h);
    
    % Compute Residual Errors
    r = residual(phi,f,h);
    
    % Restriction
    rhs = restriction(r);

    eps = zeros(size(rhs));

    % stop recursion at smallest grid size, otherwise continue recursion
    if smallest_grid_size_is_achieved
            eps = smoothing(eps,rhs,2*h);
    else        
            eps = V_Cycle(eps,rhs,2*h);        
    end
    
    % Prolongation and Correction
    phi = phi + prolongation(eps);
    
    % Post-Smoothing
    phi = smoothing(phi,f,h);    
end

我试图自己实现这个算法，但是它非常慢并且没有给出好的结果，所以很明显它做错了。我一直试图找出原因太久了，我认为值得看看是否有人可以帮助我。

如果我使用 2^5 x 2^5 点的网格大小，那么它可以解决它并给出合理的结果。然而，一旦我超过这个，无论执行多少 V-Loop，它都需要成倍地更长的时间来解决并且基本上会陷入某种程度的不准确性。在 2^7 x 2^7 点处，代码花费的时间太长而无法使用。

我认为我的主要问题是我对雅可比迭代的实现是使用线性代数来计算每一步的更新。一般来说，这应该很快，但是更新矩阵 A 是一个 n*m 大小的矩阵，计算 2^7 * 2^7 大小的矩阵的点积是昂贵的。由于大多数单元格只是零，我应该使用不同的方法计算结果吗？

如果有人对多重网格方法有任何经验，我将不胜感激！

谢谢

标签： pythonnumpymathdifferential-equations