甄玉潔，胥 康，蔣 濤，任金蓮

（揚(yáng)州大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，江蘇 揚(yáng)州 225002）

引言

含高階導(dǎo)數(shù)的多元Cahn-Hilliard (C-H)方程常被用來(lái)描述多相流體力學(xué)或熱力學(xué)等領(lǐng)域中的復(fù)雜非線性相分離現(xiàn)象[1-4].近些年來(lái)，多元C-H 的高效性數(shù)值模擬方法研究受到許多國(guó)際學(xué)者的關(guān)注，目前，已有多種基于網(wǎng)格類(lèi)方法被成功地用于單變量C-H 方程下非線性擴(kuò)散或高維多元C-H 方程下相分離過(guò)程的數(shù)值預(yù)測(cè)，如有限元法[5]、有限差分法[6-8]和多重網(wǎng)格法等[9-10].上述方法對(duì)多元C-H 描述的多維復(fù)雜非規(guī)則相界面的模擬實(shí)現(xiàn)比較復(fù)雜，且不易實(shí)施多CPU 或異構(gòu)圖形處理器(GPU)的并行計(jì)算.于是，近年來(lái)不完全依賴(lài)網(wǎng)格的純無(wú)網(wǎng)格法在國(guó)際科學(xué)計(jì)算領(lǐng)域中受到普遍關(guān)注.相較于網(wǎng)格類(lèi)方法，無(wú)網(wǎng)格方法[11-15]具有易推廣到復(fù)雜區(qū)域問(wèn)題求解和高維問(wèn)題并行計(jì)算實(shí)施的優(yōu)點(diǎn)[15-20]，如有限點(diǎn)集法(finite pointset method,FPM)[21-24]和光滑粒子流體動(dòng)力學(xué)(簡(jiǎn)稱(chēng)SPH)法等[15,20-23].

然而，上述純無(wú)網(wǎng)格FPM 和SPH 法推廣應(yīng)用到C-H 方程下多相分離過(guò)程的模擬在國(guó)際上還鮮有報(bào)道.筆者在2020年首次將FPM 拓展到二維單分量C-H 方程的求解，并發(fā)展了基于局部加密的LR-FPM 法(見(jiàn)文獻(xiàn)[24])，將其直接推廣應(yīng)用到多元C-H 方程下三維復(fù)雜多相分離過(guò)程的模擬涉及巨大計(jì)算耗時(shí)問(wèn)題，且計(jì)算精度有待進(jìn)一步提高.FPM 數(shù)值模擬多元C-H 下三維多相分離過(guò)程較網(wǎng)格類(lèi)方法的主要優(yōu)點(diǎn)在于：① 易處理高維復(fù)雜相界面和易離散求解四階空間導(dǎo)數(shù)；② 易并行計(jì)算實(shí)施.

基于上述分析，拓展FPM，采用Wendland 權(quán)函數(shù)，引入基于CUDA 的CPU-GPU 異構(gòu)并行計(jì)算，提出了一種能夠高效、準(zhǔn)確地模擬多元C-H 方程下高維復(fù)雜多相界面分離過(guò)程的改進(jìn)FPM GPU 并行算法(CFPMGPU).首先采用Wendland 權(quán)函數(shù)的FPM 對(duì)四階導(dǎo)數(shù)分裂成的兩個(gè)二階空間導(dǎo)數(shù)依次離散；其次準(zhǔn)確施加Neumann 邊界條件；最后，采用基于CUDA 語(yǔ)言的單個(gè)CPU-GPU 異構(gòu)協(xié)同的并行計(jì)算方式給出上述離散過(guò)程的GPU并行算法以加速計(jì)算模擬.本文并行計(jì)算模擬中采用一個(gè)型號(hào)為T(mén)ITAN RTX 的GPU，其包含24 GB顯存和4608 個(gè)CUDA 核心數(shù).數(shù)值算例中首先展現(xiàn)了方法的準(zhǔn)確性和并行計(jì)算效果；其次對(duì)二維復(fù)雜二或三相分離過(guò)程，以及三維復(fù)雜二相分離過(guò)程進(jìn)行了模擬預(yù)測(cè).

1 三元C-H 方程

本文主要研究多元C-H 描述的相分離現(xiàn)象[1,10]，考慮如下三元C-H 方程在區(qū)域? ?Rd(維數(shù)d=2 或3)內(nèi)的形式：

在邊界 ??上滿足齊次Neumann 邊界條件

其中，?i(x,t)(i=1,2,3)為三元混合物中每個(gè)分量的摩爾數(shù)，n為邊界上的單位法向量.根據(jù)質(zhì)量守恒律，摩爾數(shù)總和須滿足

在本文數(shù)值模擬中，利用?3表示一個(gè)任意復(fù)雜區(qū)域(相分離現(xiàn)象發(fā)生在區(qū)域內(nèi))，區(qū)域上?3≈0，其余區(qū)域，在區(qū)域邊界上考慮?3≈0.5(可詳見(jiàn)文獻(xiàn)[10]).此外，根據(jù)式(3)可得.因此，方程(1)描述的相分離過(guò)程可以簡(jiǎn)化為求解關(guān)于?1的C-H 方程：

2 CFPM-GPU 并行算法

本文主要考慮多元C-H 方程下高維復(fù)雜多相界面分離過(guò)程的純無(wú)網(wǎng)格法模擬，基于CUDA 的單個(gè)GPU 并行計(jì)算，結(jié)合Wendland 權(quán)函數(shù)的優(yōu)點(diǎn)[17]，拓展文獻(xiàn)[24]中二維FPM 離散過(guò)程得到一種能夠準(zhǔn)確、高效模擬多元高維C-H 方程的CFPM-GPU 并行算法.

2.1 三維CFPM 離散格式

對(duì)帶Neumann 邊界條件(2)含四階導(dǎo)數(shù)三維C-H 方程(5)的求解，拓展應(yīng)用文獻(xiàn)[24]里對(duì)二維情況下單個(gè)C-H 方程連續(xù)應(yīng)用兩次FPM 的離散過(guò)程，時(shí)間更新采用二階預(yù)估校正格式.基于Taylor 展開(kāi)和加權(quán)最小二乘思想的FPM 主要用來(lái)離散函數(shù)空間一階/二階導(dǎo)數(shù)，其詳細(xì)離散過(guò)程參見(jiàn)文獻(xiàn)[25-27].相對(duì)于網(wǎng)格類(lèi)方法，純無(wú)網(wǎng)格FPM 更易實(shí)施非均勻或局部加密節(jié)點(diǎn)分布計(jì)算的優(yōu)點(diǎn)也已在文獻(xiàn)[24]中對(duì)二維問(wèn)題進(jìn)行了詳細(xì)分析.文獻(xiàn)[24]中FPM 離散過(guò)程采用樣條權(quán)函數(shù)，而根據(jù)文獻(xiàn)[17]知，Wendland 權(quán)函數(shù)較樣條函數(shù)具有更好的光滑性和較高的數(shù)值精度，本文的離散過(guò)程采用Wendland 權(quán)函數(shù).

以三維空間上均勻布點(diǎn)為例，采用的Wendland 權(quán)函數(shù)[17]形式如下：

其中，r=|xj?x|為支持域半徑，j代表位置點(diǎn)x處支持域內(nèi)的相鄰點(diǎn)，wd為權(quán)函數(shù)系數(shù)，其在二維、三維空間中分別為7/(64πh2)和21/(211πh3)，h為光滑長(zhǎng)度，在離散節(jié)點(diǎn)不動(dòng)的情況下常取h≈0.95d0(d0為離散點(diǎn)初始間距).

由文獻(xiàn)[24-27]知FPM 離散過(guò)程中，對(duì)一個(gè)函數(shù)f(x)在計(jì)算區(qū)域的任意分布離散點(diǎn)x處的一階/二階導(dǎo)數(shù)的近似值，通過(guò)x點(diǎn)在其支持域內(nèi)所有相鄰節(jié)點(diǎn)xj(j=1,2,···,n)(n為相鄰節(jié)點(diǎn)數(shù))上進(jìn)行Taylor 級(jí)數(shù)展開(kāi)和加權(quán)最小二乘思想過(guò)程得到的一個(gè)線性方程組求解而獲得，即需求的局部線性方程組為

其中，上標(biāo)符號(hào)“T”為轉(zhuǎn)置，MTWM為9×9的局部系數(shù)矩陣；a9×1為x處函數(shù)一階/二階導(dǎo)數(shù)形成的未知矢量，為所有相鄰節(jié)點(diǎn)處與x處函數(shù)差和施加邊界條件構(gòu)成的已知矢量，對(duì)角線上元素為“w1,w2,···,wn,1” ；矩陣M(n+1)×9為

單位外法向量n=(nx,ny,nz).

本文數(shù)值模擬中函數(shù)f(x)為方程(5)中函數(shù)? 或μ.值得注意的是，上述線性方程組的求解過(guò)程需要求解兩次來(lái)得到四階空間導(dǎo)數(shù)的近似值，再結(jié)合二階精度的時(shí)間離散格式可得方程(5)中 ?1新時(shí)間層里的數(shù)值解，然后可得 ?2，該過(guò)程稱(chēng)為高維多元C-H 方程的CFPM 離散格式.

2.2 單個(gè)GPU 并行加速算法

眾所周知，運(yùn)用純無(wú)網(wǎng)格方法或粒子法對(duì)高維問(wèn)題的計(jì)算模擬中，為得到較高精度數(shù)值結(jié)果常需要上百萬(wàn)離散節(jié)點(diǎn)或粒子數(shù)，這將會(huì)導(dǎo)致巨大存儲(chǔ)或較長(zhǎng)時(shí)間計(jì)算的問(wèn)題[15-20].根據(jù)文獻(xiàn)[15,17-18]知，對(duì)不含全局大型系數(shù)矩陣的純無(wú)網(wǎng)格法或粒子法易于并行實(shí)施以提高計(jì)算效率.為此，已有許多基于MPI 多CPU 并行技術(shù)或算法被提出以加速純無(wú)網(wǎng)格方法的計(jì)算[17-18]，為提高加速比就需要更多的CPU，會(huì)使得計(jì)算硬件和服務(wù)器的管理成本增加.近些年來(lái)，隨著計(jì)算機(jī)GPU 的發(fā)展和NVIDA 公司不斷更新的CUDA 編程模型[16,18-20]，使得GPU 能夠進(jìn)行大規(guī)模的并行計(jì)算，比如不涉及系數(shù)矩陣的純無(wú)網(wǎng)格粒子法(SPH 或耗散粒子動(dòng)力學(xué)法(DPD))的計(jì)算模擬已被成功地與GPU 并行計(jì)算技術(shù)結(jié)合發(fā)展了高效的并行算法 (詳見(jiàn)文獻(xiàn)[16,18-19,23])，且單個(gè)GPU 的并行計(jì)算加速比可以達(dá)到100 及以上(與單個(gè)CPU 串行計(jì)算相比)，但一個(gè)GPU 小工作站的管理成本遠(yuǎn)遠(yuǎn)低于有上百個(gè)CPU 工作站的管理.

純無(wú)網(wǎng)格FPM 在高維多相分離過(guò)程的模擬還處在發(fā)展探索階段，其并行算法實(shí)施的研究更是鮮有文獻(xiàn)報(bào)道.本文首次考慮了基于CUDA 的單個(gè)GPU 來(lái)實(shí)現(xiàn)上述多元C-H 方程CFPM 離散格式的并行計(jì)算以提高效率，發(fā)展了CFPM-GPU 并行算法.在CFPM 離散格式的實(shí)施過(guò)程中涉及較大計(jì)算量或長(zhǎng)時(shí)間計(jì)算的主要有兩個(gè)方面：① 相鄰節(jié)點(diǎn)的標(biāo)定[23]，在時(shí)間層更新前每個(gè)節(jié)點(diǎn)i(對(duì)應(yīng)位置xi，i=1,2,···,N，N為離散節(jié)點(diǎn)數(shù))都需要通過(guò)搜索所有離散節(jié)點(diǎn)確定其支持域內(nèi)相鄰節(jié)點(diǎn)(每個(gè)點(diǎn)的相鄰節(jié)點(diǎn)數(shù)不盡相同，其與離散點(diǎn)分布情況和位置有關(guān)系，本文數(shù)值模擬過(guò)程中離散點(diǎn)位置不變，故模擬實(shí)施中只需要循環(huán)一次相鄰節(jié)點(diǎn)的標(biāo)定)，隨著離散節(jié)點(diǎn)數(shù)的增加該標(biāo)定的計(jì)算時(shí)間也會(huì)變大(值得注意的是，本文數(shù)值模擬中節(jié)點(diǎn)或粒子位置在計(jì)算中固定)；② 在每個(gè)時(shí)間層的物理量更新中，N個(gè)離散點(diǎn)處空間導(dǎo)數(shù)的計(jì)算以及函數(shù)值的更新在串行編程實(shí)施過(guò)程中需要循環(huán)N次，包括局部矩陣的形成和求解，隨著離散點(diǎn)數(shù)N的增大將會(huì)導(dǎo)致物理量每更新一次將需要較長(zhǎng)時(shí)間的計(jì)算，但該循環(huán)可以同時(shí)并發(fā)進(jìn)行(完全不依賴(lài)網(wǎng)格方法的特點(diǎn)).鑒于此，并行算法思想是重點(diǎn)對(duì)相鄰節(jié)點(diǎn)編號(hào)標(biāo)定和每次物理量更新中每個(gè)離散點(diǎn)處導(dǎo)數(shù)模型的循環(huán)實(shí)施并發(fā)性的計(jì)算.為實(shí)施CPFM 離散格式的GPU并行計(jì)算(并行算法較單個(gè)CPU 的加速比可見(jiàn)4.3 小節(jié)的討論)，本文采用基于CPU 和GPU 的異構(gòu)計(jì)算框架，其中以一個(gè)CPU 為主處理器(程序中常用host 表示)和一個(gè)GPU 為并發(fā)性計(jì)算實(shí)施的協(xié)處理器(程序中常用device 表示)，具體實(shí)施以CUDA-C + + 編程語(yǔ)言為主.另外需要注意的是，CFPM 格式中含9×9局部矩陣的形成和求解，其局部矩陣形成與相鄰節(jié)點(diǎn)有關(guān)，在單個(gè)CPU 上串行計(jì)算實(shí)施時(shí)可以采用向量vector 的動(dòng)態(tài)方式定義二維數(shù)組，準(zhǔn)確地標(biāo)記每個(gè)離散點(diǎn)處的相鄰節(jié)點(diǎn)編號(hào)及相鄰節(jié)點(diǎn)數(shù)，則在物理量更新計(jì)算的循環(huán)語(yǔ)句中形成局部矩陣時(shí)也需要用動(dòng)態(tài)指針定義方式.但在GPU 并行計(jì)算實(shí)施中發(fā)現(xiàn)，將上述串行實(shí)施采用的動(dòng)態(tài)指針?lè)绞絹?lái)定義計(jì)算循環(huán)里的局部矩陣會(huì)使得GPU 內(nèi)存不足，導(dǎo)致程序計(jì)算運(yùn)行的阻塞.因此，在GPU 并行算法編程中需要首先用常量的方式來(lái)給出最大相鄰節(jié)點(diǎn)數(shù)(估計(jì)值)，然后動(dòng)態(tài)定義一個(gè)一維或二維數(shù)組來(lái)標(biāo)記相鄰節(jié)點(diǎn)編號(hào)，同時(shí)在物理量更新計(jì)算的循環(huán)語(yǔ)句中可用靜態(tài)方式來(lái)定義局部矩陣的形成.

一個(gè)離散模型需消耗CPU 計(jì)算的實(shí)際時(shí)間，除了與編程語(yǔ)言和算法流程實(shí)施有一定關(guān)系外，受計(jì)算硬件性能指標(biāo)的影響也很明顯，本文數(shù)值模擬中采用的并行服務(wù)器工作站的主要硬件參數(shù)為：Intel Xeon Silver 4214 型號(hào)的CPU 處理器兩顆，每顆CPU 處理器12 核心24 線程，主頻為2.2 GHz，動(dòng)態(tài)加速頻率為3.2 GHz；三星DDR4-2666 的內(nèi)存共128 GB；一個(gè)NVIDIA Geforce Titan RTX 型號(hào)的GPU，其顯存容量為24 GB，核心頻率1770 MHz，CUDA 核心數(shù)4608 個(gè).本文給出的CFPM-GPU 并行算法流程可簡(jiǎn)化為如下幾個(gè)步驟：

1)選定光滑長(zhǎng)度h和時(shí)間步長(zhǎng)dt，并用常量定義方式給計(jì)算中需要的一些參數(shù)變量賦值.值得注意的是，考慮計(jì)算區(qū)域 ?為長(zhǎng)方形或長(zhǎng)方體 (相分離發(fā)生在一個(gè)復(fù)雜非規(guī)則區(qū)域內(nèi)且對(duì) ?進(jìn)行離散得節(jié)點(diǎn)數(shù)為N，其中上離散點(diǎn)采用均勻分布情況，上離散節(jié)點(diǎn)可根據(jù)其形狀特點(diǎn)進(jìn)行非均勻布點(diǎn)方式，但相鄰兩個(gè)節(jié)點(diǎn)的距離近似等于均勻布點(diǎn)的情況(比如圓形或球體上可以采用圓形或球形方式布點(diǎn)，也體現(xiàn)了本文采用的純無(wú)網(wǎng)格法可以任意布點(diǎn)的優(yōu)點(diǎn)).

2)動(dòng)態(tài)方式定義所需的變量分配相關(guān)變量的host 內(nèi)存及數(shù)據(jù)初始化，其中對(duì)函數(shù) ?3賦初值(計(jì)算中不變)，函數(shù) ?1和 ?2賦初始時(shí)刻的值；分配device 內(nèi)存，并把host 上數(shù)據(jù)拷貝到device 上(CPU 到GPU)；調(diào)用CUDA 核函數(shù)在device 上并行方式運(yùn)算標(biāo)記出相鄰節(jié)點(diǎn)的編號(hào)，便于后面物理量隨時(shí)間更新循環(huán)運(yùn)算.

3)當(dāng)Nnum=1 時(shí)(設(shè)Nnum為計(jì)算時(shí)間層數(shù)，實(shí)際時(shí)間t=dt×Nnum)，調(diào)用兩次CUDA 核函數(shù)在device 上并行實(shí)現(xiàn)方程(5)右端項(xiàng)CFPM 格式的運(yùn)算，并結(jié)合時(shí)間項(xiàng)格式的第一步，得到每個(gè)離散點(diǎn)處函數(shù) ?1,i的第一個(gè)中間值其中包括局部矩陣形成和調(diào)用在device 與host 上都編譯的求解線性方程組的子函數(shù).

4)結(jié)合時(shí)間項(xiàng)格式的第二步重復(fù)步驟3)的運(yùn)算，將得到每個(gè)離散點(diǎn)處函數(shù) ?1,i的第二個(gè)中間值.5)通過(guò)上述得到的兩個(gè)中間值結(jié)合時(shí)間項(xiàng)上采用的預(yù)估校正格式的第三步，調(diào)用CUDA 核函數(shù)在device(GPU)上并行得到Nnum=1 時(shí)間層上每個(gè)離散點(diǎn)處第一個(gè)函數(shù)值；然后把device 上的計(jì)算結(jié)果傳遞到host 上(GPU 上信息傳到CPU 上).6)運(yùn)用式(3)在host(CPU)上計(jì)算得到Nnum=1 時(shí)間層上每個(gè)離散點(diǎn)處第二個(gè)函數(shù)值，然后再把host 上結(jié)果傳遞到device 上.

7)重復(fù)步驟3)～6)可以得到下一個(gè)時(shí)間層上每個(gè)離散點(diǎn)上的函數(shù)值 ?1,i和 ?2,i.

8)根據(jù)需要輸出指定時(shí)間層上的計(jì)算結(jié)果.

3 CFPM-GPU 并行算法有效性驗(yàn)證

文獻(xiàn)[24]主要對(duì)FPM 離散格式求解二維單分量C-H 方程的數(shù)值收斂性進(jìn)行了分析，但未對(duì)多分量CH 描述的明顯相分離過(guò)程以及三維C-H 方程進(jìn)行數(shù)值研究.本節(jié)采用文獻(xiàn)[8]中具有二維圓環(huán)和三維球殼收縮相分離過(guò)程的兩個(gè)基準(zhǔn)算例，對(duì)給出的CFPM-GPU 算法求解帶Neumann 邊界二維/三維C-H 方程的準(zhǔn)確性和高效性進(jìn)行分析，從而驗(yàn)證本文算法捕捉相分離過(guò)程的有效性.

文獻(xiàn)[8]中兩個(gè)相分離的基準(zhǔn)算例沒(méi)有理論解析解，其通過(guò)坐標(biāo)變換給出了高精度的數(shù)值解作為基準(zhǔn)解，并列出了某一時(shí)刻64 個(gè)網(wǎng)格點(diǎn)上的基準(zhǔn)解(詳見(jiàn)文獻(xiàn)[8]中表1 和表3)用來(lái)進(jìn)行誤差分析.為了體現(xiàn)本文算法的準(zhǔn)確性，參照文獻(xiàn)[8]中誤差分析方式，則有：

二范數(shù)誤差

和最大范數(shù)誤差

其中，d為空間維數(shù)，表示本文方法數(shù)值解和文獻(xiàn)[8]中的基準(zhǔn)解.

本節(jié)數(shù)值模擬二維圓環(huán)和三維球殼收縮過(guò)程對(duì)應(yīng)的C-H方程可詳見(jiàn)文獻(xiàn)[8]，其中勢(shì)函數(shù)F(?)=0.25(?2?1)2，梯度能量系數(shù)值得注意的是，為與文獻(xiàn)[8]中的基準(zhǔn)解做對(duì)比，本節(jié)數(shù)值模擬中的離散點(diǎn)均采用的是均勻分布方式，且離散點(diǎn)間距(空間步長(zhǎng))d0=1/64.

3.1 二維圓環(huán)收縮過(guò)程模擬

考慮二維空間區(qū)域?=(?1,1)2下圓環(huán)收縮現(xiàn)象，并將區(qū)域均勻離散成N=129×129個(gè)點(diǎn)，對(duì)應(yīng)節(jié)點(diǎn)間距d0=1/64，時(shí)間步長(zhǎng)取運(yùn)算結(jié)束時(shí)間取T=8 000 000dt≈2.384.選取初始時(shí)刻函數(shù)：

經(jīng)物理量更新循環(huán)8 000 000時(shí)間步運(yùn)算，本文CFPM-GPU 法得到的數(shù)值解與基準(zhǔn)解(見(jiàn)文獻(xiàn)[8]中表1)的二范數(shù)誤差E2=0.003 39，最大范數(shù)誤差E∞=0.004 96.上述表明，本文方法求解二維圓環(huán)收縮過(guò)程對(duì)應(yīng)的C-H 方程是準(zhǔn)確穩(wěn)定的，得到的圓環(huán)收縮現(xiàn)象也是可靠的.圖1 給出了三個(gè)時(shí)刻二維圓環(huán)收縮現(xiàn)象的CFPMGPU 數(shù)值結(jié)果，其與文獻(xiàn)[8]中圖6 數(shù)值結(jié)果吻合.從圖1 中可以看出圓環(huán)外圓與內(nèi)圓隨時(shí)間延長(zhǎng)逐漸縮小形成圓環(huán)的收縮過(guò)程，結(jié)合上述算法誤差精度可知，圓環(huán)的內(nèi)外圓半徑隨時(shí)間延長(zhǎng)的變小過(guò)程與文獻(xiàn)[8]的基準(zhǔn)解一致，從而體現(xiàn)本文給出的CFPM-GPU 法模擬二維情況下相分離過(guò)程是準(zhǔn)確有效的.

圖1 二維圓環(huán)收縮現(xiàn)象CFPM-GPU 數(shù)值模擬結(jié)果：(a)t=0；(b)t=1.31；(c)t=2.38Fig.1 The 2D shrinking annuli obtained with the CFPM-GPU:(a)t=0;(b)t=1.31;(c)t=2.38

3.2 三維球殼收縮過(guò)程模擬

三維球殼收縮過(guò)程對(duì)應(yīng)的是球?qū)ΨQ(chēng)的C-H 方程，數(shù)值模擬中考慮三維空間域?=(?1,1)3，并將區(qū)域均勻離散成N=129×129×129個(gè)點(diǎn)，對(duì)應(yīng)節(jié)點(diǎn)間距d0=1/64，時(shí)間步長(zhǎng)取運(yùn)算結(jié)束時(shí)間取T=400 000dt≈1.192.選取初始時(shí)刻函數(shù)：

經(jīng)400000 時(shí)間步運(yùn)算，CFPM-GPU 法得到t=1.192 時(shí)刻的三維球殼收縮現(xiàn)象，并與初始時(shí)刻結(jié)果做對(duì)比(見(jiàn)圖2).圖2 給出了兩個(gè)時(shí)刻三維球殼收縮形狀和二維截面圓環(huán)，其展現(xiàn)出明顯的三維球殼收縮現(xiàn)象與文獻(xiàn)[8]中圖10 結(jié)果吻合.為體現(xiàn)本文方法數(shù)值求解的準(zhǔn)確性，將沿半圓環(huán)線上的數(shù)值結(jié)果與基準(zhǔn)解(見(jiàn)文獻(xiàn)[8]中表3)進(jìn)行對(duì)比，得到二范數(shù)誤差E2=0.005 4和最大范數(shù)誤差E∞=0.007 81.通過(guò)上述數(shù)值誤差和三維球殼收縮現(xiàn)象表明本文方法模擬三維相分離過(guò)程準(zhǔn)確.

圖2 三維球殼收縮現(xiàn)象CFPM-GPU 數(shù)值模擬結(jié)果：(a)t=0；(b)t=1.192Fig.2 The 3D spherical shells obtained with the CFPM-GPU:(a)t=0;(b)t=1.192

3.3 GPU 并行算法加速比分析

為體現(xiàn)本文基于CUDA 的單個(gè)GPU 并行算法的高效性，將其與單個(gè)CPU 的串行計(jì)算時(shí)間進(jìn)行對(duì)比，關(guān)于GPU 并行服務(wù)器的性能參數(shù)見(jiàn)2.2 小節(jié).為體現(xiàn)并行算法的效果，采用3.2 小節(jié)三維C-H 方程的算例，并定義加速比：

其中，tC表示單個(gè)CPU 計(jì)算消耗的時(shí)間，tC-G表示單個(gè)CPU-GPU 異構(gòu)并行計(jì)算消耗的時(shí)間，δSur為加速比.

根據(jù)2.2 小節(jié)里并行算法思想的討論，需要并行加速提高計(jì)算效率的有相鄰節(jié)點(diǎn)標(biāo)定和物理量更新循環(huán)兩個(gè)部分，表1 和表2 分別給出了相鄰節(jié)點(diǎn)標(biāo)記消耗計(jì)算時(shí)間和每個(gè)時(shí)間層里物理量更新循環(huán)所需平均計(jì)算時(shí)間及加速比.由表1 和表2 知，隨節(jié)點(diǎn)數(shù)增加，單個(gè)CPU 計(jì)算時(shí)間和GPU 并行計(jì)算時(shí)間也都隨之非線性增加，但GPU 并行消耗計(jì)算時(shí)間遠(yuǎn)少于CPU 計(jì)算時(shí)間，表明基于CPU-GPU 異構(gòu)并行計(jì)算較單個(gè)CPU 串行提高運(yùn)算效率是顯著的，且GPU 并行算法的加速比約為160.

表1 不同節(jié)點(diǎn)數(shù)下，相鄰節(jié)點(diǎn)標(biāo)定消耗計(jì)算時(shí)間對(duì)比Table 1 The consumed computing time for the calibration of neighbor nodes under different node numbers

表2 不同節(jié)點(diǎn)數(shù)下，每個(gè)時(shí)間層里物理量更新循環(huán)所需平均計(jì)算時(shí)間對(duì)比Table 2 The average computing time costs at each time step under different node numbers

4 二維/三維多相分離過(guò)程數(shù)值預(yù)測(cè)

為進(jìn)一步體現(xiàn)本文給出的CFPM-GPU 并行算法數(shù)值模擬高維復(fù)雜區(qū)域上多相分離現(xiàn)象(m相分離情況對(duì)應(yīng)m+ 1 元C-H 方程)，本節(jié)選取5 個(gè)代表性的復(fù)雜區(qū)域上的相分離算例[10]，其中包含3 個(gè)二維算例(圓盤(pán)域內(nèi)和腦切面域內(nèi)二相分離，星形域內(nèi)三相分離)和2 個(gè)三維算例.

4.1 二維多相分離現(xiàn)象數(shù)值預(yù)測(cè)

算例1二維圓盤(pán)域上二相分離現(xiàn)象

其中rand(x,y)表示介于?1 和1 之間的隨機(jī)數(shù).選取節(jié)點(diǎn)間距為d0=1/512，時(shí)間步長(zhǎng)為dt=10?6.運(yùn)算到1 05時(shí)間步(對(duì)應(yīng)時(shí)間t=105dt)結(jié)束，圖3 給出了幾個(gè)時(shí)刻下圓盤(pán)域上二相分離現(xiàn)象的CFPM-GPU 模擬結(jié)果，其中藍(lán)色為 ?1，紅色為 ?2.通過(guò)圖3 可觀察到在較短時(shí)間里(約t=0.05)就出現(xiàn)了明顯的相分離現(xiàn)象，且隨時(shí)間延長(zhǎng)相分離越來(lái)越明顯.本文結(jié)果與文獻(xiàn)[10]中圖3 的多重網(wǎng)格法得到的二相分離過(guò)程趨勢(shì)一致，表明本文CFPMGPU 法模擬預(yù)測(cè)二相分離現(xiàn)象是可靠的(注意，兩者間存在的差異主要由隨機(jī)數(shù)rand(x,y)生成的初始值不完全相同導(dǎo)致).

圖3 二維圓盤(pán)域上二相分離現(xiàn)象CFPM-GPU 模擬結(jié)果Fig.3 The 2-phase separation obtained with the CFPM-GPU in the 2D disk domain

算例2二維腦剖面域上二相分離現(xiàn)象

為進(jìn)一步體現(xiàn)本文方法模擬較復(fù)雜區(qū)域上二相分離現(xiàn)象的能力，將文獻(xiàn)[10]中二維腦剖面輪廓形成區(qū)域視為發(fā)生二相分離的復(fù)雜區(qū)域，并嵌入到單位計(jì)算區(qū)域?=(?0.5,0.5)×(?0.5,0.5)內(nèi).圖4 給出了二維腦剖面原圖及復(fù)雜腦剖面輪廓作為區(qū)域，并通過(guò)圖像處理技術(shù)和迭代計(jì)算得到模擬的初始值具體步驟如下：第一步，將腦剖面輪廓圖4(b)拷貝到畫(huà)圖工具里，設(shè)置像素大小為512 × 512；第二步，用MATLAB 處理圖像命令讀取第一步里的RGB 像素矩陣，并輸出到一個(gè)txt 文件里；第三步，通過(guò)編程對(duì)第二步得到的數(shù)據(jù)進(jìn)行0-1 化處理，得到另一個(gè)txt 文件，將其作為模擬初始數(shù)據(jù)的估計(jì)值；第四步，將第三步里得到的數(shù)據(jù)代入到C-H 方程中，經(jīng)過(guò)10 次迭代運(yùn)算得到的值作為模擬初始的數(shù)據(jù).模擬中初始時(shí)刻函數(shù)值除了 ?3，其他兩個(gè)與上述算例1 的相同.對(duì)單位計(jì)算區(qū)域采用離散節(jié)點(diǎn)513×513，時(shí)間步長(zhǎng)dt=10?6，經(jīng)1 05時(shí)間步運(yùn)算結(jié)束.圖5 給出了CFPM-GPU 并行算法模擬的二維復(fù)雜腦剖面域上隨時(shí)間演化的二相分離過(guò)程，其中藍(lán)色為 ?1，紅色為 ?2.通過(guò)圖5 可知，本文方法能夠模擬復(fù)雜區(qū)域上的二相分離現(xiàn)象，且與文獻(xiàn)[10]中圖5 結(jié)果趨勢(shì)一致(兩者存在差異的主要原因在于隨機(jī)數(shù)產(chǎn)生的初始值不完全相同)，進(jìn)一步表明CFPM-GPU 模擬預(yù)測(cè)復(fù)雜區(qū)域上二相分離現(xiàn)象是有效的.

圖4 二維腦剖面原圖及其復(fù)雜腦剖面輪廓[10]：（a）腦剖面原圖；（b）腦剖面計(jì)算區(qū)域圖Fig.4 The 2D brain section picture and the complex brain shape domian[10]:(a)the brain section picture;(b)the brain section calculation domain

圖5 二維腦剖面域上二相分離現(xiàn)象CFPM-GPU 模擬結(jié)果Fig.5 The 2-phase separation obtained with the CFPM-GPU in the 2D brain shape domain

算例3二維星形域上三相分離現(xiàn)象

許多物理問(wèn)題系統(tǒng)中存在多組分的三相分離現(xiàn)象，比如多相合金和三元納米線制備中[28].因此，本節(jié)考慮了在單位計(jì)算區(qū)域?=(?0.5,0.5)×(?0.5,0.5)內(nèi)星形域上三相分離過(guò)程的CFPM-GPU 法模擬，對(duì)應(yīng)四元CH 系統(tǒng)的求解.四元混合物中各個(gè)分量?i(x,y,t)(i=1,2,3,4)同樣滿足質(zhì)量守恒律?1+?2+?3+?4=1，且星形域各個(gè)分量初始條件為

其中θ=tan?1(y/x)且x≠0，隨機(jī)數(shù)rand(x,y)同上.值得注意的是，?1和 ?2的初始值表達(dá)式相同，但編程計(jì)算時(shí)使其取的隨機(jī)數(shù)值不完全相同，且計(jì)算中需要求解關(guān)于 ?1和 ?2的C-H 方程組.采取的計(jì)算節(jié)點(diǎn)數(shù)與時(shí)間步長(zhǎng)同上，經(jīng)1 06時(shí)間步運(yùn)算得到幾個(gè)不同時(shí)刻的三相分離現(xiàn)象(見(jiàn)圖6).觀察圖6 可知，在t=0.1 時(shí)刻出現(xiàn)了明顯的三相分離現(xiàn)象，隨時(shí)間演化三相分離現(xiàn)象越來(lái)越清晰(紅色區(qū)域?yàn)??1，藍(lán)色區(qū)域?yàn)??2，綠色區(qū)域?yàn)??3)，表明本文方法模擬預(yù)測(cè)非規(guī)則區(qū)域上三相分離過(guò)程是有效的.

圖6 二維星形域上三相分離現(xiàn)象CFPM-GPU 模擬結(jié)果Fig.6 The 3-phase separation obtained with the CFPM-GPU in the 2D star domain

4.2 三維多相分離現(xiàn)象數(shù)值預(yù)測(cè)

算例4三維環(huán)形體域上相分離現(xiàn)象

為體現(xiàn)CFPM-GPU 法模擬三維復(fù)雜域上相分離現(xiàn)象的有效性，考慮三維環(huán)形體域上二相分離，取環(huán)體區(qū)域嵌入在計(jì)算區(qū)域?=(?0.5,0.5)×(?0.5,0.5)×(?0.25,0.25)內(nèi).三元C-H 方程中每個(gè)分量的初始條件為

其中隨機(jī)數(shù)rand(x,y,z)同上.

模擬中采取離散節(jié)點(diǎn)129×129×65，初始間距d0=1/128和時(shí)間步長(zhǎng)為dt=10?6.圖7 展示了本文數(shù)值方法得到的關(guān)于 ?1三個(gè)時(shí)刻的相分離現(xiàn)象，在較短時(shí)間t=0.01 時(shí)就出現(xiàn)相分離，隨時(shí)間演化相分離現(xiàn)象越發(fā)明顯，且與文獻(xiàn)[10]中圖7 結(jié)果趨勢(shì)一致.

圖7 三維環(huán)體域上相分離現(xiàn)象CFPM-GPU 模擬結(jié)果Fig.7 The phase separation obtained with the CFPM-GPU in the 3D torus domain

算例5三維Schwarz-D 域上相分離現(xiàn)象

為進(jìn)一步體現(xiàn)CFPM-GPU 法模擬三維復(fù)雜域上相分離現(xiàn)象的可靠性，考慮Schwarz-D 形域[10]嵌入在計(jì)算域?=(0,1)×(0,1)×(0,1)的二相分離，其中發(fā)生相分離的復(fù)雜Schwarz-D 形域三元C-H 方程中每個(gè)分量的初始條件為

其中

C(x,y,z)=cos(2πx)cos(2πy)cos(2πz)?sin(2πx)sin(2πy)sin(2πz).

模擬中采取離散節(jié)點(diǎn)129×129×129，初始間距和時(shí)間步長(zhǎng)與算例4 相同.圖8 展示了三維復(fù)雜Schwarz-D 形域上CFPM-GPU 法模擬的相分離現(xiàn)象，可觀察到在較短時(shí)間t=0.01 時(shí)就出現(xiàn)明顯的相分離現(xiàn)象，本文數(shù)值結(jié)果與文獻(xiàn)[10]中的相分離趨勢(shì)一致.

圖8 三維Schwarz-D 區(qū)域上相分離現(xiàn)象CFPM-GPU 模擬結(jié)果Fig.8 The phase separation obtained with the CFPM-GPU in the 3D Schwarz-D domain

通過(guò)算例1～算例5 二維/三維復(fù)雜域上相分離現(xiàn)象的CFPM-GPU 模擬，并與文獻(xiàn)[10]中圖8(b)結(jié)果比較可知：本文提出的CFPM-GPU 求解三元或四元C-H 方程組是有效的；CFPM-GPU 數(shù)值預(yù)測(cè)二維/三維復(fù)雜域上的二相或三相分離現(xiàn)象是高效可靠的.

5 結(jié)論

本文針對(duì)二維/三維復(fù)雜域上多相分離現(xiàn)象的高效性與準(zhǔn)確性模擬，首次結(jié)合CPU-GPU 異構(gòu)并行計(jì)算、拓展求解含四階導(dǎo)數(shù)二維C-H 方程的FPM 和采用具有較好光滑性的Wendland 權(quán)函數(shù)，給出了一種基于CUDA 語(yǔ)言單個(gè)GPU 的CFPM-GPU 并行算法.為體現(xiàn)CFPM-GPU 模擬預(yù)測(cè)復(fù)雜域上多相分離的有效性，首先運(yùn)用CFPM-GPU 分別對(duì)兩個(gè)相分離的基準(zhǔn)算例(二維/三維)進(jìn)行了模擬，并對(duì)誤差和并行加速效率進(jìn)行了分析；其次，用CFPM-GPU 對(duì)無(wú)解析解的多個(gè)二維/三維復(fù)雜域上的二相或三相分離現(xiàn)象進(jìn)行了數(shù)值預(yù)測(cè)，并與多重網(wǎng)格法結(jié)果比較.所有數(shù)值結(jié)果表明：

1)本文的CFPM-GPU 并行算法求解二維/三維C-H 方程描述的相分離過(guò)程是準(zhǔn)確、高效的；

2)基于CUDA 的單個(gè)GPU 并行較CPU 串行運(yùn)算的加速比約為160 倍；

3)給出的CFPM-GPU 能準(zhǔn)確預(yù)測(cè)二維/三維復(fù)雜區(qū)域上二相或三相分離現(xiàn)象隨時(shí)間的演化過(guò)程.

致謝本文作者衷心感謝揚(yáng)州大學(xué)本科專(zhuān)業(yè)品牌建設(shè)與提升工程項(xiàng)目(ZYPP2018B007)對(duì)本文的資助.