CN109359397A - 一种有限元区域分解改进ssorpcg求解渗流场的并行方法 - Google Patents

Abstract

本发明公开了一种有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，包括如下步骤：建立求解模型，对该模型进行有限元离散，得到模型的结构信息；利用区域分解法对模型进行划分，得到P个子区域和每个子区域相对应的数据文件；消息传递接口MPI初始化，通过变分原理形成有限元支配方程；根据达西定律进而得到每个分区的渗流场；步骤5，合并P个子区域的结果，得到整体模型的计算结果。优点：所有的处理器在独立进行各自分区的计算，每个进程是平等的，从而可以有效避免了主进程的限制，可以最大限度地发挥计算节点的能力，进行更大规模的并行计算；在有限元计算中有效提高了求解大型稀疏系数矩阵方程组的计算效率。

Description

一种有限元区域分解改进SSORPCG求解渗流场的并行方法

技术领域

本发明涉及一种有限元区域分解改进SSORPCG求解渗流场的并行方法，属于水利工程有限元工程数值模拟计算和大规模并行计算技术领域。

背景技术

近年来，计算机正在以惊人的速度不断发展进步，在计算性能方面取得了长足的进步和卓越的成果。但是随着计算规模的提升和计算精度的提高，尤其是在数据密集和计算密集的数值分析应用领域，利用有限元分析的方式越频繁，普通计算机在性能上已经不能满足这些领域的需求，寻找效率高的并行算法成为当前计算不可缺少的方式之一。

在实际水利工程中，经常要求解地基及土坝的渗流问题。由于工程结构外形不规则，边界条件有比较复杂，渗流场问题使用古典方法是很难求解的，但是有限单元法可以得到良好的计算结果。目前对于有限元并行技术一般采用物理区域分割并行方法，计算网格和计算区域分解是实现粗粒度并行的最有效直接的任务分配方式。区域分解方法最大特点是适用于构造并行算法，逐步推广到各个计算领域，已经成为并行计算中常见的方法之一，在有限元分析领域发挥着重要作用。

对于有限元需要求解的大规模稀疏系数矩阵方程组问题，目前大多数采用对称逐步超松弛预处理共轭梯度(SSORPCG)求解有限元方程组，可以有效提高计算效率。现有技术中的求解方法并不适合并行求解。

发明内容

本发明所要解决的技术问题是克服现有技术的缺陷，提供一种有限元区域分解改进SSORPCG求解渗流场的并行方法，提高求解大规模稀疏系数矩阵方程组问题的能力。

为解决上述技术问题，本发明提供一种有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，包括如下步骤：

步骤1，建立求解模型，对该模型进行有限元离散，得到有限元模型，对单元和结点全部从1开始编号，得到整体编号，统计模型的结构信息，包括单元总数、单元信息和结点信息；

步骤2，模型区域分解，根据进程数目或者需要用于计算的进程数P及处理器的拓扑结构，利用区域分解法对模型进行划分，得到P个分区和每个分区相对应的数据文件；

步骤3，消息传递接口MPI初始化，确定进程的总数目，使得子区域与进程号一一对应，分别分配到P个进程中，在各个进程中形成分区内单元渗流传导矩阵和单元结点荷载列阵，将分区内所有单元集合组装两个矩阵，通过变分原理形成有限元支配方程；

步骤4，调用改进的SSORPCG求解器并行求解各个未知结点的水头值，根据达西定律进而得到每个分区的渗流场；

步骤5，合并P个分区的结果，得到整体模型的计算结果。

进一步的，所述单元总数、单元信息和结点信息具体如下：

式中，φ^e(x,y,z)和φ(ξ,η,ζ)分别为整体笛卡尔坐标系(x,y,z)和局部笛卡尔坐标系(ξ,η,ζ)下的单元内任一点的水头，m为单元的结点个数，N_i为形函数，i＝1,2,…,m,N为形函数矩阵，N＝[N₁N₂…N_m]，{φ}^e为单元结点水头列阵，{φ}^e＝[φ₁φ₂…φ_m]^T,T表示矩阵的转置。

进一步的，所述步骤2，区域分解法在串行下完成对有限元网格进行分区，依据所需进程数P的需求完成分区，具体包括步骤：

21)对步骤1中有限元离散后的模型进行单元分区和结点分区：整体有限元模型调用Metis函数库将单元分成P个区域，根据单元的分区确定其结点所属分区，将结点属于多个分区的点指定只属于一个分区；

22)确定交界面单元和in/out点：如果其中一个单元的全部结点属于某个分区，则该单元才属于该分区，属于多个分区的单元就是交界面单元，交界面单元的内边界点为out点，外边界点为in点，对于两个相邻分区，交界面公共点互为in点和out点，对于in点来说，是没有自由度的，后期并行计算时需要进行in点和out点传值；

23)对单元和结点重新编号：将上述每个分区内的单元和结点在该分区内进行重新编号，将每个分区的结点和单元编号映射到从1开始的正整数集上，得到在该分区内的局部编号，使得局部编号与整体编号具有相应的映射关系；

24)分配边界和荷载信息。

进一步的，所述单元内渗流水头传导矩阵是根据各分区内单元形状以及插值函数形成；

所述结点荷载列阵是根据分区内单元结点的体积力和边界表面力作用形成；

所述有限元支配方程是通过组装上述分区内的单元渗流水头传导矩阵和单元结点荷载列阵通过变分原理形成。

进一步的，所述改进的SSORPCG求解器是将SSORPCG求解过程并行化，修改其串行下的方程组求解方式。

进一步的，所述交界面单元的形成主要用于数据传递交换的求解，在初始迭代计算时各个分区在交界面单元处通过相互传递计算值，多次迭代后计算出各个分区的结点温度，针对大型稀疏系数方程组Ax＝b，其中A为方程组的系数矩阵、x为求解的量和b为常数项，具体实现的并行方式为：

D为A的对角阵，L为A的严格下三角，即A＝L+D+L^T，L＝L′+L_c，L_c为耦合部分，L'为各分区矩阵的严格下三角，上标T表示为矩阵的转置，当不分区时，L'＝L，Lc＝0，W＝D/ω+L′，V＝(2-ω)D/ω，则ω为松弛因子，0＜ω＜2，W、V为分解系数矩阵A过程中引进的两个算子；

初值：x⁰，g⁰＝Ax⁰-b，y⁰＝W-¹g⁰，其中x⁰为迭代求解初值，大多数情况下为0；g⁰表示初始剩余向量；y⁰初始处理x⁰用于后期迭代计算的向量算子；

第k次迭代,从k＝0开始,k表示迭代次数：

δ_k＝(y^k,Vy^k)，如果δ_k＜ε停止，否则继续如下循环：

对δ_k规约求和，并广播给每一个处理器

β_k＝δ_k/δ_k-1，如果k＝0则β_k＝0

z^k＝-Vy^k+β_kz^k-1

d^k＝W-^Tz^k

u^k＝z^k-Vd^k

如果myid＞id，传递out点d值给其他分区，从其他分区接收d值给in点

如果myid＜id，从其他分区接收d值给in点，传递out点d值给其他分区，如果myid＝id，不通讯，其中myid表示为当前进程号，id表示为除当前进程之外的其他进程号，

s^k＝W-¹u^k

t^k＝(d^k,z^k+u^k)

对t^k规约求和，并广播给每一个处理器

τ_k＝δ_k/t^k

x^k+1＝x^k+τ_kd^k

y^k+1＝y^k+τ_k(d^k+s^k)

k＝k+1表示k+1赋于k；

上述迭代中，δ_k为第k次迭代计算的y与V·y的内积，用于判断是否达到给定精度ε要求；在共轭梯度法(PCG)中，β_k为方向共轭系数，初始迭代时为0；d^k为选取的接近精确值的搜寻方向；y^k为处理每步x值的向量算子，τk为y的修正系数，z^k、u^k、s^k和t^k为迭代过程中第k次迭代引进的各项简化计算的向量算子，

计算过程中，除了传递，每个分区都有计算出的δ_k和τ_k，这是整体的误差算子，无法传递，将每次迭代的δ_k和τ_k进行规约求和，并广播到各个处理器中，保证整体迭代精度，计算出的x即为每个结点水头值，进而根据达西定律求解流速得到区域内的渗流场。

本发明所达到的有益效果：

所有的处理器在独立进行各自分区的计算，计算过程中每次迭代到给定精度，通讯一次，采用对等式方式并行，每个进程是平等的，从而可以有效避免了主进程的限制，可以最大限度地发挥计算节点的能力，进行更大规模的并行计算；在有限元计算中有效提高了求解大型稀疏系数矩阵方程组的计算效率。

附图说明

图1本发明实施例中总体计算流程图；

图2本发明实施例中区域分解4个分区图；

图3本发明实施例中交界面单元信息传递简化示意图；

图4本发明实施例中SSORPCG并行求解流程图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明提出的一种有限元区域分解改进SSORPCG求解渗流场的并行方法。本方法通过将有限元模型串行进行区域分解和合并结果，中间过程进行并行计算，具体的整个计算流程图如附图1所示。主要包括5个步骤：⑴建立求解模型，对该模型进行有限元离散，得到模型的结构信息，包括单元总数、单元信息和结点信息；⑵模型区域分解，根据P个进程数及处理器的拓扑结构，利用区域分解法对模型进行划分，得到P个子区域和每个子区域相对应的数据文件；⑶消息传递接口MPI初始化，确定进程的总数目，使得子区域与进程号一一对应，分别分配到P个进程中，在各个进程中形成分区内渗流传导矩阵和结点荷载列阵，从而形成有限元支配方程；⑷调用改进的SSORPCG求解器并行求解各个未知结点的水头值，再根据达西定律得到每个分区的渗流场；⑸合并P个子区域的结果，得到整体模型的计算结果。本实施例中根据4个进程分解为4个分区。

下面按照流程图进行详细介绍说明：

第一步，建立求解模型，对该模型进行有限元离散，对单元和结点全部从1开始编号，得到整体编号，统计模型的结构信息，包括单元总数、单元信息和结点信息；

将求解域划分为有限元网格，各单元的渗流函数用该单元结点渗流水头插值得到，设单元的结点数为m，则单元的渗流水头可以表示为：

式中，φ^e(x,y,z)和φ(ξ,η,ζ)分别为整体笛卡尔坐标系(x,y,z)和局部笛卡尔坐标系(ξ,η,ζ)下的单元内任一点的水头，m为单元的结点个数，N_i为形函数，i＝1,2,…,m,N为形函数矩阵，N＝[N₁N₂…N_m]，{φ}^e为单元结点水头列阵，{φ}^e＝[φ₁φ₂…φ_m]^T,T表示矩阵的转置。

第二步，串行下利用区域分解法对有限元模型进行分区，根据所需要计算的P个进程数，将有限元模型分解成P个子区域。如图2所示调用metis函数库(Metis函数库是由Karypis Lab开发的一个具有强大功能的串行图切分软件包，用户使用Metis软件包时，可以根据需要选择相应的切分方式)分4个区，获得相应的分区文件。具体过程和分区包括以下4个方面内容和步骤：

⑴采用metis函数库进行分区或者手动分区成P个区域，先根据单元分区，再进行结点分区。定义每个区域所属的进程号，同时根据单元和结点所属的分区，将结点属于多个分区且所属单元向内扩展形成交界面单元，最后将属于多个分区的结点指定只属于一个分区。

⑵根据第⑴步，处理各个分区的交界面单元，如图3所示，交界面单元的外边界点设置为in点，内边界点为out点，交界面公共点互为in点、out点。对于in点来说，是没有自由度的。

⑶根据第⑴步，将每个分区中的单元，结点全部进行重新编号，保证在各个分区的局部编号与整体编号具有一定的映射关系，将分区的单元信息、结点信息写入分区文件中。

⑷根据第⑶步，调用改进的SSORPCG求解器并行求解各个未知结点的水头值。

第三步，消息传递接口MPI初始化，确定进程的总数目，使得子区域与进程号一一对应，分别分配到P个进程中，在各个进程中形成分区内单元渗流水头传导矩阵[H]^e和结点荷载列阵{F}^e，进而将分区内所有单元集合组装两个矩阵，通过变分原理形成有限元支配方程。具体内容如下：

根据广义达西定律，流速用矩阵表示{v}＝-[k]{φ'}，{v}＝[v_x v_y v_z]^T，其中[k]为渗透系数，三维问题有9个系数，由于对称性k_ij＝k_ji，独立系数有6个，对于平面问题，其有3个独立系数。根据步骤1所述，各分区将此值代入到上式中得到：

{v}＝[v_x v_y v_z]^T＝-[k]{φ'}＝-[k][B]{φ}

其中

渗流条件还满足不可压缩流体的连续方程：其中Q为内源。

有限元离散后的单元渗流场泛函I^e(φ)为

其中，V表示求解区域，S为边界，只有靠近边界S的单元才会出现这项。

将上述泛函对φ_i(i＝1,2,…,m)求微分，得到

化简得到

两者的具体元素为

式中[H]^e和{F}^e即为单元传导矩阵和单元结点荷载，对于单元结点荷载{F}^e，右端第一项为体积力引起的结点荷载，后一项为边界上表面力引起的结点荷载。

对于上述方程，每个进程将其内分区各个单元的加以集合，那么上述元素中的i、j表示分区中重新编号后每个结点号，对于求解区域内的全部结点，得到方程组

即[H]{φ}＝{F}，

第四步，将串行程序改写成并行程序时，只需要修改其串行下的方程组求解方式。交界面单元的形成主要用于下步所述的数据传递交换的求解。通过改进的SSORPCG求解器，并行求解各个未知结点的水头值，再根据达西定律求解流速进而得到分区内的渗流场。针对上步中的有限元方程，统一写成Ax＝b形式。林绍忠在河海大学学报，1998(03):114-117上发表的文章“用预处理共轭梯度法求解有限元方程组及程序设计”中修改的SSORPCG的迭代方式如下：

初值：x⁰，g⁰＝Ax⁰-b，y⁰＝W^-1g⁰

第k次迭代(从k＝0开始)：

δ_k＝(y^k,Vy^k)，如果δ_k＜ε停止，否则继续如下循环：

β_k＝δ_k/δ_k-1(如果k＝0则β_k＝0)

z^k＝-Vy^k+β_kz^k-1

d^k＝W-^Tz^k

u^k＝z^k-Vd^k

s^k＝W-¹u^k

t^k＝(d^k,z^k+u^k)

τ_k＝δ_k/t^k

x^k+1＝x^k+τ_kd^k

y^k+1＝y^k+τ_k(d^k+s^k)

k＝k+1

其中，D为A的对角阵，L为A的严格下三角，即A＝L+D+L^T，W＝D/ω+L，V＝(2-ω)D/ω，ω为松弛因子(0＜ω＜2)，ε为给定的迭代精度。

以下针对上述迭代进行并行算法的修改，主要包括以下内容：

根据图3所示，采用矩阵分块思想描述有限元的矩阵形式：EAi表示A分区的第i个单元，Ai表示A分区的第i个结点，EBi表示B分区的第i个单元，Bi表示B分区第i个结点。A₃和B₃作为in点，是没有自由度的，但A₂和B₃以及A₃和B₂实际上是同一个点。所以写出的总体方程为：

其中各分区方程为：A_ij和B_ij分别表示A和B分区内，非交界面位置处形成的单元传导矩阵，即为j结点对i结点的水头传递的贡献；X_Ai和X_Bi、F_Ai和F_Bi分别表示A和B分区的i结点水头和i结点的荷载；

如果单独考虑A、B两个分区，也不考虑边界点，即忽略in点和out点，针对上述SSORPCG计算式d^k＝W^-Tz^k，转化为W^Td^k＝z^k，写出的两个分区的计算矩阵，并转化为方程组型式：

分区A：

分区B：

因为W^T＝D^T/ω+L^T，为了书写方便和理解，将两个分区中对角元素的/ω舍去，但计算时是包含此项的。d_Ai和d_Bi分别表示分区A、B中第i个点选取的接近精确值的搜寻方向。

同时，根据图3和(1)式，结合考虑in点(无自由度)和out点，A₂₃变成了A₂B₂，B₂₃变成了B₂A₂，上述等式变化为：

分区A：

分区B：

根据上述迭代过程可以看出，当进行计算u^k＝z^k-Vd^k时，计算时需要用到计算时需要用到 和为在交界面点A₂、B₂处理搜索方向d_A2和d_B2后的计算算子。每次迭代过程中需要将传递给分区A、传递给分区B，且过程中是通过in点和out点进行数据传递交换完成并行计算。

由此得到矩阵的的耦合项为：

在上述原SSORPCG中u^k＝z^k-Vd^k，此时并行计算需要加此耦合项，故得

最终改进的SSORPCG求解有限元方程组Ax＝b的并行算法流程如图4所示，具体说明如下：

针对Ax＝b，D为A的对角阵，L为A的严格下三角，即A＝L+D+L^T。L＝L′+L_c，L_c为耦合部分；L'为各分区矩阵的严格下三角，当不分区时，L'＝L，L_c＝0。W＝D/ω+L′，V＝(2-ω)D/ω，则其余各项参数与上述SSORPCG参数含义一致。用myid表示为当前进程号，id表示为除当前进程之外的其他进程号，对所有处理器myid(myid＝0,…,p-1)同时执行如下过程：

初值：x⁰，g⁰＝Ax⁰-b，y⁰＝W^-1g⁰

第k次迭代(从k＝0开始)：

对δ_k规约求和，并广播给每一个处理器

β_k＝δ_k/δ_k-1(如果k＝0则β_k＝0)

z^k＝-Vy^k+β_kz^k-1

d^k＝W^-Tz^k

u^k＝z^k-Vd^k

如果myid＞id，传递out点d值给其他分区，从其他分区接收d值给in点

如果myid＜id，从其他分区接收d值给in点，传递out点d值给其他分区

如果myid＝id，不通讯

s^k＝W^-1u^k

t^k＝(d^k,z^k+u^k)

对t^k规约求和，并广播给每一个处理器

τ_k＝δ_k/t^k

x^k+1＝x^k+τ_kd^k

y^k+1＝y^k+τ_k(d^k+s^k)

k＝k+1

在上述迭代过程中，除了传递，因为每个分区都有计算出的δ_k和τ_k，这是整体的误差算子，无法传递，所以还需要将每次迭代的δ_k和τ_k进行规约求和，并广播到各个处理器中，保证整体迭代的精度，计算出的x即为每个结点水头值，再将此值代入达西定律中求解结点流速即可得到4个分区的渗流场。

第五步，合并4个子区域的结果，得到整体模型的结算结果。由于每个分区都有交界面单元，只要合并交界面单元的结点和单元即可，最终合并成整体模型的结果。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (6)

1.一种有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，包括如下步骤：

步骤1，建立求解模型，对该模型进行有限元离散，得到有限元模型，对单元和结点全部从1开始编号，得到整体编号，统计模型的结构信息，包括单元总数、单元信息和结点信息；

步骤2，模型区域分解，根据进程数目或者需要用于计算的进程数P及处理器的拓扑结构，利用区域分解法对模型进行划分，得到P个分区和每个分区相对应的数据文件；

步骤3，消息传递接口MPI初始化，确定进程的总数目，使得子区域与进程号一一对应，分别分配到P个进程中，在各个进程中形成分区内单元渗流传导矩阵和单元结点荷载列阵，将分区内所有单元集合组装两个矩阵，通过变分原理形成有限元支配方程；

步骤4，调用改进的SSORPCG求解器并行求解各个未知结点的水头值，根据达西定律进而得到每个分区的渗流场；

步骤5，合并P个分区的结果，得到整体模型的计算结果。

2.根据权利要求1所述的有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，所述单元总数、单元信息和结点信息具体如下：

式中，φ^e(x,y,z)和φ(ξ,η,ζ)分别为整体笛卡尔坐标系(x,y,z)和局部笛卡尔坐标系(ξ,η,ζ)下的单元内任一点的水头，m为单元的结点个数，N_i为形函数，i＝1,2,…,m,N为形函数矩阵，N＝[N₁ N₂ … N_m]，{φ}^e为单元结点水头列阵，{φ}^e＝[φ₁ φ₂ … φ_m]^T,T表示矩阵的转置。

3.根据权利要求1所述的有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，所述步骤2，区域分解法在串行下完成对有限元网格进行分区，依据所需进程数P的需求完成分区，具体包括步骤：

21)对步骤1中有限元离散后的模型进行单元分区和结点分区：整体有限元模型调用Metis函数库将单元分成P个区域，根据单元的分区确定其结点所属分区，将结点属于多个分区的点指定只属于一个分区；

22)确定交界面单元和in/out点：如果其中一个单元的全部结点属于某个分区，则该单元才属于该分区，属于多个分区的单元就是交界面单元，交界面单元的内边界点为out点，外边界点为in点，对于两个相邻分区，交界面公共点互为in点和out点，对于in点来说，是没有自由度的，后期并行计算时需要进行in点和out点传值；

23)对单元和结点重新编号：将上述每个分区内的单元和结点在该分区内进行重新编号，将每个分区的结点和单元编号映射到从1开始的正整数集上，得到在该分区内的局部编号，使得局部编号与整体编号具有相应的映射关系；

24)分配边界和荷载信息。

4.根据权利要求1所述的有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，所述单元内渗流水头传导矩阵是根据各分区内单元形状以及插值函数形成；

所述结点荷载列阵是根据分区内单元结点的体积力和边界表面力作用形成；

所述有限元支配方程是通过组装上述分区内的单元渗流水头传导矩阵和单元结点荷载列阵通过变分原理形成。

5.根据权利要求1所述的有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，所述改进的SSORPCG求解器是将SSORPCG求解过程并行化，修改其串行下的方程组求解方式。

6.根据权利要求3所述的有限元区域分解改进SSORPCG求解渗流场的并行方法，其特征在于，所述交界面单元的形成主要用于数据传递交换的求解，在初始迭代计算时各个分区在交界面单元处通过相互传递计算值，多次迭代后计算出各个分区的结点温度，针对大型稀疏系数方程组Ax＝b，其中A为方程组的系数矩阵、x为求解的量和b为常数项，具体实现的并行方式为：

D为A的对角阵，L为A的严格下三角，即A＝L+D+L^T，L＝L′+L_c，L_c为耦合部分，L'为各分区矩阵的严格下三角，上标T表示为矩阵的转置，当不分区时，L'＝L，L_c＝0，W＝D/ω+L′，V＝(2-ω)D/ω，则ω为松弛因子，0＜ω＜2，W、V为分解系数矩阵A过程中引进的两个算子；

初值：x⁰，g⁰＝Ax⁰-b，y⁰＝W-¹g⁰，其中x⁰为迭代求解初值，大多数情况下为0；g⁰表示初始剩余向量；y⁰表示初始处理x⁰用于后期迭代计算的向量算子；

第k次迭代,从k＝0开始,k表示迭代次数：

δ_k＝(y^k,Vy^k)，如果δ_k＜ε停止，否则继续如下循环：

对δ_k规约求和，并广播给每一个处理器

β_k＝δ_k/δ_k-1，如果k＝0则β_k＝0

z^k＝-Vy^k+β_kz^k-1

d^k＝W-^Tz^k

u^k＝z^k-Vd^k

如果myid＞id，传递out点d值给其他分区，从其他分区接收d值给in点

如果myid＜id，从其他分区接收d值给in点，传递out点d值给其他分区，如果myid＝id，不通讯，其中myid表示为当前进程号，id表示为除当前进程之外的其他进程号，

s^k＝W-¹u^k

t^k＝(d^k,z^k+u^k)

对t^k规约求和，并广播给每一个处理器

τ_k＝δ_k/t^k

x^k+1＝x^k+τ_kd^k

y^k+1＝y^k+τ_k(d^k+s^k)

k＝k+1表示k+1赋于k；

上述迭代中，δ_k为第k次迭代计算的y与V·y的内积，用于判断是否达到给定精度ε要求；在共轭梯度法(PCG)中，β_k为方向共轭系数，初始迭代时为0；d^k为选取的接近精确值的搜寻方向；y^k为处理每步x值的向量算子，τ_k为y的修正系数，z^k、u^k、s^k和t^k为迭代过程中第k次迭代引进的各项简化计算的向量算子，

将每次迭代的δ_k和τ_k进行规约求和，并广播到各个处理器中，保证整体迭代精度，计算出的x即为每个结点水头值，进而根据达西定律求解流速得到区域内的渗流场。