几何多重网格的系数推导

参考文献:A Fast Unsmoothed Aggregation Algebraic Multigrid Framework for the Large-Scale Simulation of Incompressible Flow

考虑求解方程-\nabla^2 p=b,或写成矩阵形式Ap=b,加负号是为了让矩阵正定。向量p的每一项对应某个n_x\times n_y\times n_z里面的一个格子。我们采用转换积分的方式:

\int_V\nabla^2 p\mathrm{d}V=\oint_S\nabla p\cdot\bm{n}\mathrm{d}S.

这里面V是一个格子,S是格子的六个面,实际就是在六个面上把p的相应偏导乘以面积加起来。假设格子的长度是h,那么相应的系数矩阵就是(忽略各种边界条件):

A_{\alpha,\beta}=\begin{cases}
6h,\text{ }\alpha=\beta,\\
-h,\text{ }\alpha\neq \beta\text{ but adjacent,}\\
0,\text{ otherwise.}
\end{cases}

这个h是这么来的:想象p_\alphap_\beta\alpha对应方程里面的形式,它会产生一项

\frac{-p_\beta+p_\alpha}{h}\cdot h^2=h(-p_\beta+p_\alpha).

系数就是h。然后,我们用RP矩阵把它变换到粗网格上,格子的长度为2h。这里面取

R_{\alpha,\gamma}=\begin{cases}
1/8,\text{ }\text{if coarse voxel }\alpha\text{ covers fine voxel }\gamma,\\
0,\text{ }\text{otherwise.}
\end{cases}

我们用A^h代表细网格,A^H代表粗网格,那么A^H=RA^hP,定义P=8R^T。计算这个矩阵乘法(以下使用爱因斯坦求和约定):

A^H_{\alpha,\beta}=R_{\alpha,\gamma}A^h_{\gamma,\eta}P_{\eta,\beta}=8R_{\alpha,\gamma}A^h_{\gamma,\eta}R_{\beta,\eta}.

只有当\alpha覆盖\gamma\beta覆盖\eta时两个R项非零。也就是说,在\alpha\beta的八个儿子里面各选一格,计算相应的乘积。假设\alpha,\beta相邻,那么只有边界上四对儿子的A^h非零,合起来就是4/8\cdot(-h)\cdot 8/8=-(2h)/4。如果\alpha=\beta,那么有8\gamma=\eta的情况和24\gamma\eta相邻的情况,合起来是8/8\cdot(6h)\cdot 8/8+24/8\cdot(-h)\cdot 8/8=3h=6\cdot(2h)/4,是\alpha\neq \beta的六倍,符合泊松方程的性质。

随后考虑restriction和prolongation操作。Restriction应当为:

r^H=Rr^h.

r^H_{\alpha}=R_{\alpha,\gamma}r^h_{\gamma}.

这实际上就是对\alpha的八个儿子取平均。求解方程A^Hp^H=r^H得到p^H,然后做prolongation:

p^h_{\text{new}}=p^h_{\text{old}}+\lambda P p^H.

这里面\lambda是松弛系数,取1或者2。计算

(P p^H)_{\gamma}=P_{\gamma,\alpha}p^H_{\alpha}=8R_{\alpha,\gamma}p^H_\alpha.

这就是直接取\gamma父亲的值。但是有一个问题是,现在我们得到的A^H=-\nabla^2/4,即粗网格上求解的是

-\nabla^2p/4=r^H.

而我们希望在粗网格上求解-\nabla^2。这时候可以把右端r^H乘以4,求解

-\nabla^2p = 4r.

这个额外的scaling是采用积分形式导致的。假如求解的是逐点的微分,那么A^h里面的相邻项是-1/h,而A^H里面的相邻项是-1/2h,正好符合-\nabla^2的形式。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注