使用分布式离散化的多重网格预处理器求解微分方程

此示例说明如何使用预条件迭代求解器和分布式数组来求解泊松方程。通过使用分布式数组，您可以使用计算机集群的内存（而不仅仅是单台计算机的内存）来扩展计算。您无需更改代码即可进行扩展。

此示例延续了使用分布式数组通过迭代方法求解线性方程组中涵盖的主题。基于 [1]，该示例使用泊松方程（其形式称为均质稳态热方程）仿真房间内的热量分布。稳定状态意味着热量不随时间变化，均匀意味着没有外部热源。

$- Δ u = - (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}}) = 0$

在等式中， $u$ 代表房间每个点 $(x, y, z)$ 的温度。为了求解该方程，首先需要使用有限差分离散化方法，通过线性方程组对其进行近似。然后，使用预条件共轭梯度法（pcg）来解决该系统。预条件转换问题以提高数值求解器的性能。通过使用分布式数组，您可以利用计算机集群的组合内存并实现更精细的离散化。

了解如何：

建立离散的三维网格和边界条件。
定义离散化和多重网格预处理器。
应用预条件数值求解器来求解三维体积中的热方程。

离散化空间维度

在这个示例中，边长为 1 的立方体仿真了房间。第一步是使用 3-D 网格将其离散化。

本示例中的预条件方法使用了几个具有不同粒度级别的网格。每个级别都会使网格在每个维度上粗化 2 倍。定义多重网格层的数量。

multigridLevels = 2;

定义最精细网格每个维度中的点数，X、Y 和 Z。预条件方法要求该网格中的点数能够被 2^multigridLevels 整除。在这种情况下，点的数量必须可以被 4 整除，因为多重网格层数是 2。

numPoints.X = 32;
numPoints.Y = 32;
numPoints.Z = 32;

使用 meshgrid 函数将空间维度离散化为三维网格。使用 linspace 按照点数均匀划分各个维度。请注意，要包含立方体的边界，必须添加两个额外的点。

[X,Y,Z] = meshgrid(linspace(0,1,numPoints.X+2), ...
    linspace(0,1,numPoints.Y+2), ...
    linspace(0,1,numPoints.Z+2));

定义边界条件

假设房间有一扇窗户和一扇门。墙壁和天花板的温度恒定为 0 度，窗户的温度恒定为 16 度，门的温度恒定为 15 度。地板为 0.5 度。目标是确定房间内部的温度分布。

使用关系运算符定义地板、窗户和门的坐标，并定义这些边界元素上的温度。边界是立方体的各个面，因此 X、Y 或 Z 中的一个必须是 0 或 1。将立方体的其余边界和内部设置为 0。

floor  = (0.0 <= X & X <= 1.0) & (0.0 <= Y & Y <= 1)   & (Z == 0.0);
window = (X == 1)              & (0.2 <= Y & Y <= 0.8) & (0.4 <= Z & Z <= 0.6);
door   = (0.4 <= X & X <= 0.6) & (Y == 1.0)            & (0.0 <= Z & Z <= 0.6);

u = zeros(size(X));
u(floor)  = 0.5;
u(window) = 16;
u(door)   = 15;

这些边界条件指定了解决方案沿域边界必须采用的常数值。这种类型的边界条件称为狄利克雷边界条件。

使用 slice 函数将边界条件可视化。使用位于立方体边界的切片来显示非零边界条件。

xSlices = 1;
ySlices = 1;
zSlices = 0;
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Constant nonzero boundary conditions'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;
set(f,'EdgeColor',[0 0 0]);

离散化并求解微分方程

此示例使用有限差分近似方法将微分方程离散化为线性系统，并使用多重网格预处理器来提高迭代求解器的性能。对于此示例，使用 discretizePoissonEquation 和 multigridPreconditioner 支持函数中的离散化和预处理器。discretizePoissonEquation 函数作为支持文件附加到此示例。要访问此文件，请将示例作为实时脚本打开。

对于其他问题，请选择适合您应用的离散化和预处理器。

在这个示例中，discretizePoissonEquation 使用七点模板有限差分方法将泊松方程离散化为具有不同粒度的多个网格。该函数创建了离散化的多重网格结构，具有预先计算的三角分解和在粗略级别和精细级别之间映射的运算符。预处理器利用多重网格信息以有效的方式近似求解。

除其他技术外，该预处理器还应用平滑技术，通过一系列近似值来最小化误差。定义平滑步骤的数量。使用更多的步骤可以使近似值更加准确，但计算量也会更大。然后，离散化微分方程并建立预处理器。

numberOfSmootherSteps = 1;
[A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,u);

Level 0: The problem is of dimension 32768 with 223232 nonzeros.
Level 1: The problem is of dimension 4096 with 27136 nonzeros.
Level 2: The problem is of dimension 512 with 3200 nonzeros.

preconditioner = setupPreconditioner(multigridData);

使用预条件共轭梯度法求解线性系统。

tol = 1e-12;
maxit = numel(A);
pcg(A,b,tol,maxit,preconditioner);

pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

使用分布式数组进行扩展

如果您需要更多的计算资源，例如内存，则可以使用分布式数组进行扩展，而无需更改代码。分布式数组将您的数据分布在多个工作单元上，它们可以利用计算机集群的计算性能和内存。

启动一个并行工作单元池。默认情况下，parpool 使用您的默认集群。在 MATLAB 主页选项卡上的环境区域中的并行 > 选择默认集群中检查您的默认集群配置文件。

parpool;

Starting parallel pool (parpool) using the 'MyCluster' profile ...
Connected to the parallel pool (number of workers: 12).

使用 u 函数将温度变量 distributed 分布在集群中各个工作单元的内存中。

distU = distributed(u);

您可以使用与之前相同的代码；不需要进行任何更改，因为如果输入是分布式数组，则离散化和预条件子函数会创建分布式数组。许多 MATLAB 函数针对分布式数组进行了增强，因此您可以按照使用内存数组的方式使用它们。

请注意，discretizePoissonEquation 返回包含分布式数据的结构体。要以分布式方式使用结构体内的分布式数据，必须在 spmd 代码块内创建结构体。您还必须调用在 spmd 代码块内使用它的任何函数。

spmd
    [A,b,multigridData] = discretizePoissonEquation(numPoints,multigridLevels,numberOfSmootherSteps,distU);
    preconditioner = setupPreconditioner(multigridData);
end

Analyzing and transferring files to the workers ...done.
Lab  1: 
  Level 0: The problem is of dimension 32768 with 223232 nonzeros.
  Level 1: The problem is of dimension 4096 with 27136 nonzeros.
  Level 2: The problem is of dimension 512 with 3200 nonzeros.

使用 pcg 代码块内的 spmd 以分布式方式解决线性系统。

spmd
    x = pcg(A,b,tol,maxit,preconditioner);
end

Lab  1: 
  pcg converged at iteration 45 to a solution with relative residual 5.4e-13.

绘制结果

求解器的解决方案是一个适合内存的向量。使用 gather 将数据从工作单元发送到客户端。将数据重新整形为三维数组并重新排序维度以产生最终解决方案。将 u 的内部部分设置为该解决方案。外部，即边界，已经包含了边界条件的值。

x3D = reshape(gather(x),numPoints.X,numPoints.Y,numPoints.Z);
u(2:end-1,2:end-1,2:end-1) = permute(x3D, [2, 1, 3]);

使用 slice 函数将解决方案可视化。添加额外的切片来绘制立方体内部的温度。您可以使用 Rotate 工具或改变切片的位置来探索解决方案。

xSlices = [.5,1];
ySlices = [.5,1];
zSlices = [0,.5];
f = slice(X,Y,Z,u,xSlices,ySlices,zSlices,'nearest');
title('Heat distribution'), xlabel('x'), ylabel('y'), zlabel('z');
colorbar, colormap cool;
shading interp;

您可以尝试此示例中的 numPoints 的不同值来测试不同级别的离散化。使用较大的值可以提高分辨率，但需要更多的内存。此外，multigridLevels 越大，预处理器的内存效率越高。然而，较大的 multigridLevels 意味着预处理器的精度较低，因为粗化会降低每个级别的精度。因此，求解器可能需要更多次迭代才能达到相同的精度水平。

定义预条件子

定义一个多重网格预处理器以便与预条件共轭梯度法一起使用。这种类型的预处理器使用具有不同粒度的多个离散化网格来更有效地近似线性方程组的解。本示例中的预处理方法基于 [2]，并遵循以下主要阶段：

使用高斯-赛德尔逼近方法进行预平滑。
在更粗略的层面上计算残差解。
如果是在较粗的层次上，则递归地进行预处理；如果是在最粗的层次上，则直接求解。
使用更粗的网格解决方案更新解决方案。
使用高斯-赛德尔逼近方法进行后平滑。

function x = multigridPreconditioner(mgData,r,level)

if(level < mgData(level).MaxLevel)
    x = zeros(size(r),like=r);
    
    % Presmooth using Gauss-Seidel
    for i=1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
    
    % Compute residual on a coarser level
    Axf = mgData(level).Matrices.A*x;
    rc = r(mgData(level).Fine2CoarseOperator)- Axf(mgData(level).Fine2CoarseOperator);
    
    % Recursive call until coarsest level is reached
    xc = multigridPreconditioner(mgData, rc, level+1);
    
    % Update solution with prolonged coarse grid solution
    x(mgData(level).Fine2CoarseOperator) = x(mgData(level).Fine2CoarseOperator)+xc;
    
    % Postsmooth using Gauss-Seidel
    for i = 1:mgData(level).NumberOfSmootherSteps
        x = mgData(level).Matrices.LD \ (-mgData(level).Matrices.U*x + r);
        x = mgData(level).Matrices.DU \ (-mgData(level).Matrices.L*x + r);
    end
else
    % Obtain exact solution on the coarsest level
    x = mgData(level).Matrices.A \ r;
end

end

创建一个函数，该函数采用多重网格数据并返回将预条件子应用于输入数据的函数句柄。在此示例中，此函数句柄是 pcg 的预条件器输入。您必须创建此函数，因为无法在 spmd 代码块内定义匿名函数。

function preconditioner = setupPreconditioner(multigridData)

if ~isempty(multigridData)
    preconditioner = @(x,varargin) multigridPreconditioner(multigridData,x,1);
else
    preconditioner = [];
end

end

参考资料

[1] Dongarra, J., M. A. Heroux, and P. Luszczek.HPCG Benchmark:A New Metric for Ranking High Performance Computing Systems.Knoxville, TN:University of Tennessee, 2015.

[2] Elman, H. C., D. J. Silvester, and A. J. Wathen.Finite Elements and Fast Iterative Solvers:With Applications in Incompressible Fluid Dynamics.Oxford, UK:Oxford University Press, 2005, Section 2.5.

另请参阅

pcg | distributed | spmd | gather