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

（揚州大學 數(shù)學科學學院，江蘇 揚州 225002）

引言

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

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

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

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)為三元混合物中每個分量的摩爾數(shù)，n為邊界上的單位法向量.根據(jù)質(zhì)量守恒律，摩爾數(shù)總和須滿足

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

2 CFPM-GPU 并行算法

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

2.1 三維CFPM 離散格式

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

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

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

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

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

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

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

2.2 單個GPU 并行加速算法

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

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

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

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

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

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

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

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

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

3 CFPM-GPU 并行算法有效性驗證

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

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

二范數(shù)誤差

和最大范數(shù)誤差

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

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

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

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

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

圖1 二維圓環(huán)收縮現(xiàn)象CFPM-GPU 數(shù)值模擬結果：(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 三維球殼收縮過程模擬

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

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

圖2 三維球殼收縮現(xiàn)象CFPM-GPU 數(shù)值模擬結果：(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 的單個GPU 并行算法的高效性，將其與單個CPU 的串行計算時間進行對比，關于GPU 并行服務器的性能參數(shù)見2.2 小節(jié).為體現(xiàn)并行算法的效果，采用3.2 小節(jié)三維C-H 方程的算例，并定義加速比：

其中，tC表示單個CPU 計算消耗的時間，tC-G表示單個CPU-GPU 異構并行計算消耗的時間，δSur為加速比.

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

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

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

4 二維/三維多相分離過程數(shù)值預測

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

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

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

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

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

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

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

圖4 二維腦剖面原圖及其復雜腦剖面輪廓[10]：（a）腦剖面原圖；（b）腦剖面計算區(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 模擬結果Fig.5 The 2-phase separation obtained with the CFPM-GPU in the 2D brain shape domain

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

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

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

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

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

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

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

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

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

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

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

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

其中

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

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

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

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

5 結論

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

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

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

3)給出的CFPM-GPU 能準確預測二維/三維復雜區(qū)域上二相或三相分離現(xiàn)象隨時間的演化過程.

致謝本文作者衷心感謝揚州大學本科專業(yè)品牌建設與提升工程項目(ZYPP2018B007)對本文的資助.