Stokes-Cahn-Hilliard方程的二阶并行解耦算法

Second-Order Parallel Decoupled Algorithms for the Stokes-Cahn-Hilliard Equation

摘要: 本文为Stokes-Cahn-Hilliard相场模型开发了一种具有二阶时间精度的时间并行解耦算法，在时间上基于谱延迟校正方法和时间并行Parareal算法。此外，我们使用稳定化方法处理Cahn-Hilliard方程，运用了压力投影方法处理Stokes方程中压力p和速度u的耦合。相较于解耦的方案，进一步提高了算法的计算效率。最后，通过数值实验证明了它的稳定性和高效性。

Abstract: In this paper, we develop a time-parallel decoupled algorithm with second-order temporal accuracy for the Stokes-Cahn-Hilliard phase-field model, based in time on the spectral-deferred correction method and the time-parallel Parareal algorithm. In addition, we use a stabilisation method for the Cahn-Hilliard equation and employ a pressure projection method for the coupling of pressure p and velocity u in the Stokes equation. Compared with the decoupled scheme, the computational efficiency of the algorithm is further improved. Finally, numerical experiments demonstrate its stability and efficiency.

1. 引言

流体、固体和气体混合物在多孔介质中的流动发生在各种物理、生物和工业过程中。因此，了解和准确预测多孔介质中的多相多组分流动对许多科学和工程应用都具有重大意义。相场方法是建模和模拟多相流问题的常用工具，相关概述见[1]-[4]。Stokes-Cahn-Hilliard (SCH)模型[5]由Cahn-Hilliard方程和不可压缩的Stokes方程构成，是相场模型中的一种，该模型描述了两种不可压缩不相容流体的混合物的相场演化过程。但是想要建立有效的数值格式来求解该耦合相场模型仍然面临刚性大和耦合项的问题。凸分裂方法和稳定线性隐式方法[6]是构造无条件能量稳定格式的两种常用方法。Eyre [7]推广的凸分裂格式是一种能保持能量耗散定律的数值格式。该格式隐式地处理化学势的凸部分和凹部分，并建立了唯一可解的和无条件能量稳定的格式。Kim等[8]提出了具有对数自由能的Cahn-Hilliard方程的非线性凸分裂傅里叶谱格式，这是对多项式自由能能量泛函凸分裂思想的适当推广。不幸的是，凸分裂方法通常会导致非线性格式。同时，如何将算法推广到更高阶仍是一个问题。

谱延迟校正(SDC)方法是一种具有良好的稳定性和高精度的高阶数值方法，它在2000年由Dutt等人[9]首次提出。然而，SDC方法由于本身的数值积分节点和迭代校正次数导致了较高的计算代价。为了降低计算成本，Minion [10]提出了一种半隐式含刚性和非刚性的常微分方程的谱延迟校正方法(SISDC)，该方法对非刚性项显式处理，刚性项隐式处理。在这个半隐式框架内，利用一阶欧拉法求解修正方程。假设数值积分是精确的，那么修正方程每迭代一次，数值解的精度就会提高一阶。之后，为了进一步提高精度和迭代速度，郭瑞晗[11]开发了一类基于二阶时间积分法的半隐式SDC方法，每增加一次迭代，精度就会提高两阶。这为更快地构建高阶的数值方案提供了一种方法。

Parareal方法作为时间并行计算领域[12]的重要突破，自2001年由Lions等人提出以来，受到了学术界和工业界的广泛关注。该方法巧妙地将时间区间剖分为若干个子区间，然后利用多重处理器，在每个时间子区间上并行求解，是一种在时间方向上并行计算常微分方程或时间离散的偏微分方程数值解的迭代方法。这种方法也被称作多重射击法，或者时间多重网格法，Gander等人[13]曾针对该方法展开过探讨。自从Parareal方法问世，其再度激发了学界对于构建时间并行方法的研究热情[14]-[16]。在此背景下，一系列更适配高性能并行计算机的数值算法相继涌现并不断发展，进而在实际应用中诸如航空航天领域，加快了飞行器的设计和优化进程，降低研发成本等。此外，文献[17] [18]中提出了不同的改进算法来提高算法收敛阶。特别地，由于Parareal方法和SDC方法都是预测校正方法，Minion等人[19]将两者结合起来。在文献[20]中给出该混合算法的详细分析，该融合算法的优势主要体现在两个方面：首先，将时间并行的Parareal算法与高阶SDC算法相结合，可以有效地平摊SDC算法在每个时间步上较高的计算成本；另一方面，相较于其它高阶算法，例如高阶龙格库塔方法等，选取高阶SDC算法作为Parareal算法中的高阶格式可以利用低阶数值方法迭代求解得到高阶精度，减少了Parareal算法中每次迭代的计算成本，提高了并行效率。因此，在本文中，我们基于谱延迟校正方法和时间并行Parareal算法，提出了SCH模型的二阶时间并行算法，该算法在提高计算效率和精度的同时，保证离散算法的稳定性。

本文的主要内容安排如下。在第2节中，首先介绍了SCH模型和Parareal时间并行方法。然后给出求解SCH系统的二阶SDC解耦算法，然后基于SDC算法和Parareal算法构造了SCH模型的二阶时间并行解耦算法。在第3节中，通过数值示例来验证了二阶时间并行解耦算法的准确性、稳定性和高效性。第4节做出总结。

2. 模型与算法

2.1. Stokes-Cahn-Hilliard模型

两种不可压缩不相容流体的混合物的相场演化模型是由在时间区域$\left(0,T\right)$ 和有界的Lipschitz空间域$\Omega \subset {ℝ}^d\left(d=2,3\right)$ 的Stokes-Cahn-Hilliard方程所描述的：

$\{\begin{array}{ll}{c}_t+\nabla \cdot \left(cu\right)=M\Delta \omega ,\hfill & \left(0,T\right)\times \Omega \hfill \\ \omega =\lambda \left(-\Delta c+f\left(c\right)\right),\hfill & \left(0,T\right)\times \Omega \hfill \\ u_t+\nabla p-\nu \Delta u=-c\nabla \omega ,\hfill & \left(0,T\right)\times \Omega \hfill \\ \nabla \cdot u=0,\hfill & \left(0,T\right)\times \Omega \hfill \end{array}$ (1)

其中，$M$ 是与松弛时间尺度相关的迁移率常数，$\nu$ 是粘度系数，$\lambda$ 与表面张力相对应，$u$ 是速度场，$p$ 是修正压力，$\omega$ 是化学势。相场函数$c$ 起着观察微观浓度的作用。比如在纯流体的区域内，相场函数

$c$ 可以达到数值−1和1，而在扩散界面内，$\left|c\right|<1$ 。对于非线性项$f\left(c\right)=F^{\prime }\left(c\right)$ ，其中$F\left(c\right)=\frac{1}{4{\eta }^2}{\left(c-1\right)}^2$

表示双阱势能，用于惩罚偏离物理约束$\left|c\right|\le 1$ 的行为。方程(1)的前两项是Cahn-Hilliard方程，$\nabla \cdot \left(cu\right)$ 项模拟平流效应。方程(1)的后两项是不可压缩的Stokes方程，其中的非线性项$c\nabla \omega$ 模拟了表面张力效应。

下面将讨论模型(1)的边界条件和初始条件，我们假设

$\left(u\left(0\right),c\left(0\right)\right)=\left(u^0,c^0\right),$ (2)

和

${u|}_{\partial \Omega }=0,{{\partial }_n{c}^{n+1}|}_{\partial \Omega }=0,{{\partial }_n{\omega }^{n+1}|}_{\partial \Omega }=0.$ (3)

此外，若是假设该系统除重力外没有其他外力，则Stokes-Cahn-Hilliard方程满足以下能量定律：

$\frac{\text{d}}{\text{d}t}\left\{W\left(\varphi ,\nabla \varphi \right)+{\displaystyle {\int }_{\Omega }\frac{1}{2}{\left|u\right|}^2\text{d}x}\right\}=-{\displaystyle {\int }_{\Omega }M{\left|\nabla \mu \right|}^2+\nu {\left|\nabla u\right|}^2\text{d}x}.$ (4)

2.2. 时间并行Parareal算法

基于Parareal算法的理念，当时间区间$\left[0,T\right]$ 被分割为$N$ 个互不重叠的子区间时，每个子区间将被分配给不同的处理器。为简化表述，我们将这$N$ 个处理器记为$P_1,\cdots ,P_N$ 。在Parareal时间并行算法[20]里，会用到两个数值逼近算法，分别记为传播函数$\mathcal{G}$ 与$ℱ$ 。为确保算法在并行计算时具备高效性，$\mathcal{G}$ 的计算成本必须低于$ℱ$ 。另外，在时间区间$\left[t^n,t^{n+1}\right]$ 内，若以$t^n$ 时刻的$\widehat{u}$ 作为初始值，通过$\mathcal{G}$ 和$ℱ$ 计算得到$t^{n+1}$ 时刻的数值解，则分别用$\mathcal{G}\left(t^{n+1},t^n,\widehat{u}\right)$ 与$ℱ\left(t^{n+1},t^n,\widehat{u}\right)$ 来表示。

在Parareal方法中，$\mathcal{G}$ 常被用来串行计算所有节点的临时逼近解，即

$u_1^{n+1}=\mathcal{G}\left(t^{n+1},t^n,u_1^n\right),n=0,1,\cdots ,N-1.$

只要每个处理器$P_{n+1}$ 获得初始值$u_1^n$ 和临时逼近解$u_1^{n+1}$ ，它们就能并行计算出更精确的近似解 $ℱ\left(t^{n+1},t^n,u_1^n\right)$ 。然后使用如下的校正步串行得到最终的逼近解

$u_2^{n+1}=\mathcal{G}\left(t^{n+1},t^n,u_2^n\right)+ℱ\left(t^{n+1},t^n,u_1^n\right)-\mathcal{G}\left(t^{n+1},t^n,u_1^n\right).$ (5)

2.3. SCH方程的二阶算法

在本节中，我们首先给出了SCH方程的二阶解耦SDC算法。然后，基于Parareal并行算法，建立了一种求解SCH系统的二阶解耦时间并行算法。首先我们通过添加稳定项(7)得到SCH方程的一阶解耦算法。然后在一阶解耦算法的基础上通过谱延迟校正方法构造二阶解耦SDC算法。

算法1. (SCH方程的二阶SDC解耦算法)

1) (一阶解耦算法)给定初值$\left(c^0,{\omega }^0,u^0,p^0\right)$ ，求解$\left(c_1^{n+1},{\omega }_1^{n+1},u_1^{n+1},p_1^{n+1}\right),n\ge 0$ 满足

$\{\begin{array}{l}\frac{c_1^{n+1}-c_1^n}{\Delta t}+\nabla \cdot \left(c_1^n{u}_{*1}^n\right)-M\Delta {\omega }_1^{n+1}=0,\\ {\omega }_1^{n+1}-\frac{\lambda }{{\eta }^2}\left(c^{n+1}-c^n\right)+\lambda \left(\Delta c^{n+1}-f\left(c^n\right)\right)=0;\end{array}$ (6)

$u_{*1}^n=u_1^n-\Delta t{c}_1^n\nabla {\omega }_1^{n+1};$ (7)

$\frac{{\tilde{u}}_1^{n+1}-u_1^n}{\Delta t}=\nu \Delta {\tilde{u}}_1^{n+1}-\nabla p_1^n-c_1^n\nabla {\omega }_1^{n+1};$ (8)

$\{\begin{array}{l}-\Delta \left(p_1^{n+1}-p_1^n\right)=-\frac{1}{\Delta t}\nabla \cdot {\tilde{u}}_1^{n+1},\\ {{\partial }_n\left(p^{n+1}-p^n\right)|}_{\partial \Omega }=0,\end{array}{u}_1^{n+1}={\tilde{u}}_1^{n+1}-\Delta t\nabla \left(p_1^{n+1}-p_1^n\right).$ (9)

2) (二阶解耦SDC算法)由一阶解耦算法求得$\left(c_1^{n+1},{\omega }_1^{n+1},u_1^{n+1},p_1^{n+1}\right)$ ，求解$\left(c_2^{n+1},{\omega }_2^{n+1},u_2^{n+1},p_2^{n+1}\right)$ 满足

$\{\begin{array}{l}\frac{c_2^{n+1}-c_2^n}{\Delta t}+\nabla \cdot \left(c_2^n{u}_2^n\right)-M\Delta {\omega }_2^{n+1}=\frac{1}{2}\left(\nabla \cdot \left(c_1^n{u}_1^n\right)-\nabla \cdot \left(c_1^{n+1}{u}_1^{n+1}\right)\right)+\frac{1}{2}M\left(\Delta {\omega }_1^n-\Delta {\omega }_1^{n+1}\right),\\ {\omega }_2^{n+1}-\frac{\lambda }{{\eta }^2}\left(c_2^{n+1}-c_2^n\right)+\lambda \left(\Delta c_2^{n+1}-f\left(c_2^n\right)\right)\\ =\frac{1}{2}\left({\omega }_1^n-{\omega }_1^{n+1}\right)+\frac{\lambda }{{\eta }^2}\left(c_1^{n+1}-c_1^n\right)+\frac{1}{2}\lambda \left(\Delta c_1^n-\Delta c_1^{n+1}\right)-\frac{1}{2}\lambda \left(f\left(c_1^n\right)-f\left(c_1^{n+1}\right)\right),\end{array}$ (10)

$\{\begin{array}{l}\frac{{\tilde{u}}_2^{n+1}-u_2^n}{\Delta t}-\nu \Delta {\tilde{u}}_2^{n+1}=\frac{1}{2}\nu \left(\Delta {\tilde{u}}_1^n-\Delta {\tilde{u}}_1^{n+1}\right)-\nabla p_2^n\\ +\frac{1}{2}\nabla \left(p_1^n-\nabla p_1^{n+1}\right)-\frac{1}{2}\left(c_2^n\nabla {\omega }_2^n+c_2^{n+1}\nabla {\omega }_2^{n+1}\right),\end{array}$ (11)

$\begin{array}{l}\{\begin{array}{l}-\Delta \left(p_2^{n+1}-p_2^n-p_1^{n+1}+p_1^n\right)=-\frac{1}{\Delta t}\nabla \cdot {\tilde{u}}_2^{n+1},\\ {{\partial }_n\left(p_2^{n+1}-p_2^n-p_1^{n+1}+p_1^n\right)|}_{\partial \Omega }=0,\end{array}\\ u_2^{n+1}={\tilde{u}}_2^{n+1}-\Delta t\nabla \left(p_2^{n+1}-p_2^n-p_1^{n+1}+p_1^n\right).\end{array}$ (12)

为了进一步提高数值求解SCH系统的计算效率，我们将基于SDC方法和并行算法，提出一种具有时间二阶精度且能量稳定的时间并行算法。该算法是对上述解耦算法的进一步优化，相较于一般的SDC算法，降低了计算规模，大大提高了计算效率。在本文中，我们采用一阶解耦算法作为传播函数$\mathcal{G}$ ，并采用类似(10)~(12)式的SDC方法作为传播函数$ℱ$ ，在确保精度和能量稳定性的前提下，利用时间步并行技巧提高了计算效率。

算法2. (SCH方程的二阶解耦并行算法)

1) 串行求解$\mathcal{G}\left(t^{n+1},t^n,c_1^n,{\omega }_1^n,u_1^n,p_1^n\right)$ ，$n=0,\cdots ,N-1$ ，即在时间方向依次计算(10)式得到 $\left(c_1^{n+1},{\omega }_1^{n+1},u_1^{n+1},p_1^{n+1}\right)$ 。

2) 利用二阶解耦SDC算法的校正步并行求解$ℱ\left(t^{n+1},t^n,c_1^n,{\omega }_1^n,u_1^n,p_1^n\right)$ ，$n=0,\cdots ,N-1$ ，以第$n+1$ 个处理器$P_{n+1}$ 为例：求解$\left(c_F^{n+1},{\omega }_F^{n+1},u_F^{n+1},p_F^{n+1}\right)$ 满足

$\{\begin{array}{l}\frac{c_F^{n+1}-c_1^n}{\Delta t}+\nabla \cdot \left(c_1^n{u}_{*1}^n\right)-M\Delta {\omega }_F^{n+1}=\frac{1}{2}\left(\nabla \cdot \left(c_1^n{u}_1^n\right)-\nabla \cdot \left(c_1^{n+1}{u}_1^{n+1}\right)\right)+\frac{1}{2}M\left(\Delta {\omega }_1^n-\Delta {\omega }_1^{n+1}\right),\\ {\omega }_F^{n+1}-\frac{\lambda }{{\eta }^2}\left(c_F^{n+1}-c_1^n\right)+\lambda \left(\Delta c_F^{n+1}-f\left(c_1^n\right)\right)\\ =\frac{1}{2}\left({\omega }_1^n-{\omega }_1^{n+1}\right)+\frac{\lambda }{{\eta }^2}\left(c_1^{n+1}-c_1^n\right)+\frac{1}{2}\lambda \left(\Delta c_1^n-\Delta c_1^{n+1}\right)-\frac{1}{2}\lambda \left(f\left(c_1^n\right)-f\left(c_1^{n+1}\right)\right),\end{array}$ (13)

$u_{*1}^n=u_1^n-\Delta t{c}_1^n\nabla {\omega }_1^{n+1};$ (14)

$\frac{{\tilde{u}}_F^{n+1}-u_1^n}{\Delta t}-\nu \Delta {\tilde{u}}_F^{n+1}=\frac{1}{2}\nu \left(\Delta {\tilde{u}}_1^n-\Delta {\tilde{u}}_1^{n+1}\right)-\frac{1}{2}\left(\nabla p_1^n+\nabla p_1^{n+1}\right)-c_1^n\nabla {\omega }_1^{n+1},$ (15)

$\begin{array}{l}\{\begin{array}{l}-\Delta \left(p_F^{n+1}-p_1^n-p_1^{n+1}+p_1^n\right)=-\frac{1}{\Delta t}\nabla \cdot {\tilde{u}}_F^{n+1},\\ {{\partial }_n\left(p_F^{n+1}-p_1^n-p_1^{n+1}+p_1^n\right)|}_{\partial \Omega }=0,\end{array}\\ u_F^{n+1}={\tilde{u}}_F^{n+1}-\Delta t\nabla \left(p_F^{n+1}-p_1^n-p_1^{n+1}+p_1^n\right)\end{array}$ (16)

3) 串行求解$\mathcal{G}\left(t^{n+1},t^n,c_2^n,{\omega }_2^n,u_2^n,p_2^n\right)$ ，$n=0,\cdots ,N-1$ ，即在时间方向依次求解$\left(c_G^{n+1},{\omega }_G^{n+1},u_G^{n+1},p_G^{n+1}\right)$ ，满足

$\{\begin{array}{l}\frac{c_G^{n+1}-c_2^n}{\Delta t}+\nabla \cdot \left(c_2^n{u}_2^n\right)-M\Delta {\omega }_G^{n+1}=0,\\ {\omega }_G^{n+1}-\frac{\lambda }{{\eta }^2}\left(c_G^{n+1}-c_2^n\right)+\lambda \left(\Delta c_G^{n+1}-f\left(c_2^n\right)\right)=0;\end{array}$ (17)

$\{\begin{array}{l}{c}_{2,h}^{n+1}=c_G^{n+1}+c_F^{n+1}-c_{1,h}^{n+1},\\ {\omega }_{2,h}^{n+1}={\omega }_G^{n+1}+{\omega }_F^{n+1}-{\omega }_{1,h}^{n+1};\end{array}$ (18)

$\frac{{\tilde{u}}_G^{n+1}-u_2^n}{\Delta t}=\nu \Delta {\tilde{u}}_G^{n+1}-\nabla p_2^n-\frac{1}{2}\left(c_2^n\nabla {\omega }_2^n+c_2^{n+1}\nabla {\omega }_2^{n+1}\right);$ (19)

$\{\begin{array}{l}-\Delta \left(p_G^{n+1}-p_2^n\right)=-\frac{1}{\Delta t}\nabla \cdot {\tilde{u}}_G^{n+1},\\ {{\partial }_n\left(p^{n+1}-p^n\right)|}_{\partial \Omega }=0;\end{array}$ (20)

$u_G^{n+1}={\tilde{u}}_G^{n+1}-\Delta t\nabla \left(p_G^{n+1}-p_2^n\right);$ (21)

$\{\begin{array}{l}{\tilde{u}}_{2,h}^{n+1}={\tilde{u}}_G^{n+1}+{\tilde{u}}_F^{n+1}-{\tilde{u}}_{1,h}^{n+1},\\ u_{2,h}^{n+1}=u_G^{n+1}+u_F^{n+1}-u_{1,h}^{n+1},\\ p_{2,h}^{n+1}=p_G^{n+1}+p_F^{n+1}-p_{1,h}^{n+1}.\end{array}$ (22)

随后将介绍数值实验，以证明二阶解耦并行算法算法在求解SCH方程时确实具有高阶精确性、能量稳定性和高效性。

3. 数值实验

在这一部分，我们将借助数值实验来检验针对SCH方程的二阶时间并行解耦算法的有效性。所有数值运算均基于MATLAB软件开展。就空间离散化而言，对于Stokes方程，我们运用Taylor-Hood (P2-P1)有限元进行近似处理；对于Cahn-Hillard方程，则采用二次多项式有限元(P2)来加以近似。

3.1. 收敛性实验

我们通过数值实验验证算法2的时间收敛阶。我们考虑在计算区域为$\Omega =\left[0,1\right]\times \left[0,1\right]$ ，最终时间为$T = 0.1$ ，物理参数为$\eta =0.04,\lambda ={0.04}^2,M=0.1,\nu =\frac{1}{100}$ 。初始条件为

${\varphi }_0=0.24\cos \left(2\pi x\right)\cos \left(2\pi y\right)+0.4\cos \left(\pi x\right)\cos \left(3\pi y\right),$

$u_0=\left(-{\sin }^2\left(\pi x\right)\sin \left(2\pi y\right),{\sin }^2\left(\pi y\right)\sin \left(2\pi x\right)\right).$

在文献[21]中，给出了以下计算收敛阶的方法：假定

$v_h^{\Delta t}\left(x,t_m\right)\approx v\left(x,t_m\right)+C_1\left(x,t_m\right)\Delta t^{\gamma }+C_2\left(x,t_m\right)h^{\mu },$

定义

${\rho }_{v,\Delta t,i}=\frac{{\left|v_h^{\Delta t}\left(x,t_m\right)-v_h^{\frac{\Delta t}{2}}\left(x,t_m\right)\right|}_i}{{\left|v_h^{\frac{\Delta t}{2}}\left(x,t_m\right)-v_h^{\frac{\Delta t}{4}}\left(x,t_m\right)\right|}_i}\approx \frac{4^{\gamma }-2^{\gamma }}{2^{\gamma }-1},$

Table 1. The temporal convergence order of Algorithm 2 with the fixed grid size $h=\frac{1}{32}$ , the final moment $T=0.1$ and different time steps $\Delta t=\frac{T}{80},\frac{T}{160},\frac{T}{320},\frac{T}{640}$

表1. 算法2的时间收敛阶，其中固定网格大小为$h=\frac{1}{32}$ ，最终时刻$T=0.1$ ，并选取不同的时间步长 $\Delta t=\frac{T}{80},\frac{T}{160},\frac{T}{320},\frac{T}{640}$

$\Delta t$	$\Vert {\delta }_{u_1}^{\Delta t}\Vert$	order	${\Vert {\delta }_{u_1}^{\Delta t}\Vert }_{H^1\left(\Omega \right)}$	order	$\Vert {\delta }_{u_2}^{\Delta t}\Vert$	order	${\Vert {\delta }_{u_1}^{\Delta t}\Vert }_{H^1\left(\Omega \right)}$	order
T/80	2.9932e−05	1.9036	4.8890e−04	1.8604	3.0472e−05	1.9030	4.4057e−04	1.8399
T/160	7.9998e−06	1.9518	1.3464e−04	1.9275	8.1478e−06	1.9518	1.2307e−04	1.9175
T/320	2.0678e−06	1.9760	3.5393e−05	1.9631	2.1062e−06	1.9760	3.2578e−05	1.9581
T/640	5.2565e−07		9.0774e−06		5.3539e−07		8.3846e−06
$\Delta t$	$\Vert {\delta }_{\varphi }^{\text{Δ}t}\Vert$	order	${\Vert {\delta }_{\varphi }^{\Delta t}\Vert }_{H^1\left(\Omega \right)}$	order	$\Vert {\delta }_p^{\text{Δ}t}\Vert$	order
T/80	7.0617e−04	1.8723	1.4863e−02	1.8732	1.8003e−04	1.8902
T/160	1.9288e−04	1.9349	4.0571e−03	1.9357	4.8567e−05	1.9436
T/320	5.0447e−05	1.9671	1.0605e−03	1.9676	1.2626e−05	1.9715
T/640	1.2903e−05		2.7114e−04		3.2194e−06

${\rho }_{v,h,i}=\frac{{\left|v_h^{\Delta t}\left(x,t_m\right)-v_{\frac{h}{2}}^{\Delta t}\left(x,t_m\right)\right|}_i}{{\left|v_{\frac{h}{2}}^{\Delta t}\left(x,t_m\right)-v_{\frac{h}{4}}^{\Delta t}\left(x,t_m\right)\right|}_i}\approx \frac{4^{\mu }-2^{\mu }}{2^{\mu }-1}.$

这里，$v$ 可以为$u$ 、$p$ 或$\phi$ ，且$i=0,1$ 。如果${\log }_2{\rho }_{v,\Delta t,0}$ ，${\log }_2{\rho }_{v,h,0}\approx 2$ ，则表明$L^2$ 误差的空间收敛阶和时间收敛阶分别为$O\left(\Delta t^2\right)$ 和$O\left(h^2\right)$ 。同样，通过计算${\log }_2{\rho }_{v,\Delta t,1}$ 和${\log }_2{\rho }_{v,h,1}$ 可以得到$H^1$ 误差的空间和

时间收敛阶。为方便记述，我们将定义${\delta }_v^{\Delta t}=v_h^{\Delta t}-v_h^{\frac{\Delta t}{2}}$ 和${\delta }_v^h=v_h^{\Delta t}-v_{\frac{h}{2}}^{\Delta t}$ 。

表1列出了算法2的时间收敛阶。我们固定空间步长为$h=\frac{1}{32}$ ，并选择不同的时间步长 $\Delta t=\frac{T}{80},\frac{T}{160},\frac{T}{320},\frac{T}{640}$ 。表1中的结果显示，对于$c,p$ 和$u$ ，算法2的时间收敛阶在$L^2$ 范数下和$H^1$ 范数下均达到了二阶精度。

3.2. 稳定性实验

在这一部分中，我们将验证算法2的能量稳定性，首先，先定义离散时间$t = n\Delta t$ 的离散能量函数

$E_t=\lambda \left(F\left({\varphi }^n\right),1\right)+\frac{\lambda }{2}{\Vert \nabla {\varphi }^n\Vert }^2+\frac{1}{2}{\Vert u^n\Vert }^2.$

在本实验中，取空间步长$h=\frac{1}{64}$ ，其他参数与上述收敛性实验中的参数相同。图1给出了算法2在

$\Delta t =0.1,0.05,0.01$ 三个时间步长下的能量变化情况。图1显示，算法2在不同的时间步长下，离散能量函数$E_t$ 随时间是递减的。因此，我们可以得到，算法2具有良好的稳定性。

(A) $\Delta t=0.1$

(B) $\Delta t=0.05$

(C) $\Delta t=0.01$

Figure 1. Evolution of discrete energy $E_t$ for Algorithm 2 are shown at $\Delta t=0.1,0.05,0.01$

图1. 算法2在$\text{Δ}t=0.1,0.05,0.01$ 时，离散能量$E_t$ 变化情况

3.3. 运行时间比较

为了验证算法2和算法3的并行有效性，我们比较了时间并行算法与串行二阶SDC算法，即算法1的计算时间。在本实验中，计算区域和参数选择与实验3.1中相同。表2和表3中的数值结果表明，算法2和3的运算时间相近，且与串行算法1相比，都提升了运行效率，但通过实验2表明了算法2的稳定性更好。因此，在同等条件下，算法2在稳定性和计算效率均具有优势。

Table 2. The CPU time of the time-parallel and serial decoupling schemes with the fixed mesh $h=\frac{1}{32}$

表2. 时间并行方案和串行解耦方案在固定网格$h=\frac{1}{32}$ 时的CPU运行时间

数值算法	CPU (s)
	$\text{Δ}t=1/40$	$\text{Δ}t=1/80$	$\text{Δ}t=1/160$	$\text{Δ}t=1/320$
二阶解耦算法	11.598	23.5438	45.2505	89.1838
二阶解耦并行算法	6.2058	10.9925	22.3059	44.272
加速比	1.9	2.1	2.0	2.0

Table 3. The CPU time of the time-parallel and serial decoupling schemes with the fixed mesh $h=\frac{1}{64}$

表3. 时间并行方案和串行解耦方案在固定网格$h=\frac{1}{64}$ 时的CPU运行时间

数值算法	CPU (s)
	$\text{Δ}t=1/40$	$\text{Δ}t=1/80$	$\text{Δ}t=1/160$	$\text{Δ}t=1/320$
二阶解耦算法	62.912	121.499	229.6978	461.2383
二阶解耦并行算法	28.3711	56.9073	114.9398	234.079
加速比	2.2	2.1	2.0	2.0

4. 结论

为提高数值求解具有匹配密度的Stokes-Cahn-Hilliard相场模型的计算效率，本文基于半隐的SDC算法和Parareal算法构造了一个具有时间二阶精度的时间并行算法。我们在数值上验证了算法的收敛性和稳定性，同时通过不同算法的比较说明了我们的时间并行算法的高效性。