寫大規模三維並行計算流體力學(CFD)求解器時，有什麼經驗和心得？

01-07

祖傳 CFD 代碼二十多年，各路神仙都更改過代碼，使得求解器很不友好。計劃重寫組內求解器，徵求各位關於結構、內容以及個人在寫求解器時的一些看法、經驗與心得。

關注這個問題一段時間了，覺得自己有資格回答一下。

本人花了博三到博四一年半時間寫了一套三維流固耦合程序，目前已經基本完工。主要是想計算葉輪機械的葉片流固耦合問題。

主要模塊：

1. 二維/三維兼容的可壓縮CFD：最基本的部分是重構、通量、隱式/顯式時間推進等，然後是SST湍流模式，接著是基於ALE的動網格，在接著是openmp和mpi並行技術，最後是目前剛剛完成了旋轉非慣性系和滑移網格。

2.二維/三維結構動力學CSD，包括多種二階混合單元，靜力學計算，動力學計算，模態分析，旋轉非慣性系等等。

3.流固耦合FSI，包括流固耦合界面雙向傳值，隱式時間推進，顫振分析等等。

由於實驗室沒有計算流CFD的程序積累，所以我是完全從零開始，獨自開發，比較好的方面是沒有歷史包袱，可以按照自己的想法寫。

下面是我的一些經驗：

1. 編程語言選擇。一句話：初學原理Matlab，正式編程c++。現在cfd領域內有很多用Fortran，迷信fortran的高效，我的看法是，除非實驗室有大量歷史遺留Fortran代碼可用，否則用c++。我開始編程的前半年用Fortran，後面見到c++，相見恨晚，於是花一個多月用c++全部重寫了，帶來的是編程效率大幅提高和計算性能小幅提高。

2. c++編程技巧。c++編程一定要注意熟練地用一些技巧來提高效率，c++不像Fortran，隨手一寫效率就很高，c++需要程序員要有一定的經驗和技巧來寫，如果有熟練的技巧，c++的效率會比fortran高。比如a = a+b改成a+=b、矩陣連續定址的i j順序、多使用局部變數等，不一而足。

3. 矩陣模板準備。由於c++沒有原生的矩陣類型，而vector模板在多維的時候定址效率較低，所以建議在開始時就先編寫一個高性能矩陣的模板文件，而且，編寫模版也是對不熟悉c++的同學一個很好的入門和練手過程。我寫了一維到六維的矩陣和稀疏矩陣模板，重載了+-*/等基本的矩陣操作符，添加了很多有用的小功能，比如支持負指標（對虛擬網格編號很有用）、debug模式下數組越界檢查、矩陣輸出到txt文件等等。

4. CFD的編寫順序。先寫無粘歐拉方程，再加粘性湍流等。固體部分可以同步寫。歐拉方程求解器的完成是一個重要的節點，這時需要對後續程序進行好好規劃。如果要做MPI並行，要在歐拉方程求解器做完之後就開始改寫成MPI，否則後面全部完成之後改起來太麻煩。

5. 兼容問題。 編寫的時候一定要儘可能兼容二維三維，儘可能做到子程序通用、復用，不然維護很麻煩。我是採用一個宏變數#define dimension 3or2來控制程序是計算二維還是三維，如果某一段代碼，二維三維處理完全無法統一，用

#if（dimension == 2）

...

#endif

#if（dimension ==3）

...

#endif

來分別實現。這樣從編譯上兼容，並且二維程序的計算效率並不會因為對三維的兼容有明顯降低。同理，兩種並行模式MPI和openMP的切換也可以用這種宏變數的方式，定義一個參數，來控制編譯出來的版本採用哪種並行模式。

6. 動網格以及ALE。一開始就要想好程序需不需要涉及動網格，一旦有可能涉及動網格就要做好兼容ALE求解的準備。可以在完成歐拉方程之後就著手ALE改寫。一定要注意，很多文獻提到，動網格下的CFD要保證幾何守恆（GCL），實際上只有顯式時間推進才需要，對於雙時間步隱式時間推進不需要。

(2017/11/13更新如下）

7.數據類型double/float。double和float類型在x64浮點數計算上其實速度差別並不大，但是對於內存中的連續大矩陣的訪問提取，float類型由於小一半的內存，因此訪問耗時更小。串列計算的時候可能體會不到，並行計算（openmp/MPI）的時候，由於內存帶寬的限制，內存訪問的速度往往是非常耗時，經過測試，在我的工作站和台式機上，從堆里訪問數據（也就是new出來的數據）的並行加速倍數一般最多3倍，而在棧裡面的局部變數訪問和運算幾乎是線性加速，因此要儘可能減小內存訪問相對於純計算的比例。在某些精度本身要求不高的大矩陣上用float數據類型替換double能大大加速內存訪問數據，提高並行效率，比如湍流模式要用的壁面距離矩陣、隱式迭代的對角矩陣、梯度運算元矩陣等等，用float之後，計算精度並沒有下降。

（想到再寫）

我正在做這件事，甚至工作量更大: 我已經新寫了非線性固體求解器，tbb做的並行，還沒寫mpi；剛寫了不可壓縮流體求解器，還處在串列階段，剛算了幾個2D的benchmark cases，也準備做並行；未來還要寫流固耦合的部分。

許多人都說不可能或者值得，我倒不這麼認為。但首先得想清楚幾個問題：1. 老代碼有多爛； 2. 自己有多熟悉數值方法；3. 你有多少時間。

每個人都覺得老代碼爛，但也許你見了別人的代碼就不這麼認為了。我覺得爛有幾個級別： (a) 丑但能用 -- 結構不清晰，代碼重複，變數重複，幾十個參數傳來傳去，但代碼大體正確，數學上沒問題； (b) 效率低下 -- 在以上基礎之上，在不影響正確性的地方有公式錯誤，收斂緩慢，但結果在可接受的範圍；(c) 充滿苟且 -- 出現神秘參數不能改動，數學公式明顯和論文對不上，但運行熟悉的case的結果竟然是對的，換新的case就不一定了；(d) 忍無可忍 -- 在以上基礎上，有重大公式錯誤，造成求解不收斂，根本無法通過Benchmark測試，出於某種原因前人一直在拿這稀爛的結果忽悠人。

如果你的代碼的癥狀在a和b的範疇之內，只要你還看得懂就不要考慮重寫了；到c這一步就慎重考慮；到d這一步，要麼直接走人要麼就壯士斷腕。

其次要考慮你的基礎知識如何，包括但不限於要會推導空間和時間離散格式，理解線性系統左邊和右邊分別是什麼，稀疏矩陣的性質，懂一些線性方程迭代求解的方法，還要有一些編程基礎。如果你連個線彈性固體求解器都寫不出來就還是算了。

然後還要考慮時間因素，如果要求不是很高又有紮實基礎的話，幾個月時間是可以做一個能用的並行CFD程序的（不涉及湍流）。如果你是碩士或者本科生趁早放棄，如果是博士第一年或者第二年我說go for it! 第三第四年就看因人而異了。

我這個回答主要想記錄和分享一下我寫CFD代碼的經驗，我目前做的只是最基本的CFD：不可壓縮/沒有湍流模型/牛頓流體，方法是傳統有限元，求不鄙視。。

使用庫

絕對不要純手工打造一個日常使用的CFD代碼，否則你寫出來的東西很可能不堪重用而且比現有的代碼還要難看。如果你用有限元的話 dealii 是一個非常好的庫，不光提供了有限元需要的所有的標準的組件（比如shape function, Gauss integration, dof renumbering etc.）還集成了p4est/PETSc/Trilinons等，給你「一站式」體驗！這裡我不推薦FEnics/ngsolve等因為我覺得他們太「高層」，底層的東西控制不了，只適合用來快速測試你的formulation。（如果說錯了請指正。）假如你是做FDM的，如果沒有類似dealii這樣的庫那起碼你能用PETSc，這是一個專註於大規模矩陣運算的庫，用起來有點難就是了。

使用合適的數值格式

數值格式決定了求解的難度。CFD計算中最難搞的是對流項，因為它把矩陣變成非對稱的從而增加矩陣求解的難度，另外當它佔主導地位時還經常需要額外的穩定項。但如果你把這一項做顯式處理，把它挪到方程右邊，那你的方程左邊就跟stokes方程一樣了。這樣能降低求解的難度，但代價的CFL條件更苛刻。如果你要用純隱式方法你就要正確地線性化NS方程再使用牛頓法求解，總之要選擇並推導出合適的數值格式。

求解Ax = b

假設你不管使用什麼方法成功地組裝起來了一個稀疏、分塊矩陣。那麼恭喜你完成了任務的10%。剩下90%的功夫都得花在求解Ax = b上。如果你不僅僅要寫一個玩具求解器那麼迭代式求解是必須的，我用的是GMRES。dealii/PETSc什麼的有矩陣求解器，但遺憾的是對於Navier Stokes方程，沒有一個高效的preconditioner是不可能在有限時間裡求出解來的！關於怎樣precondition NS方程的文獻浩如煙海而且現在仍然是一個活躍的研究領域，如果你的目標是要解幾十億個自由度的系統，geometric multigrid是唯一現實的選擇。這個如果你不是數學系的PhD我估計你不會有興趣去做。但是做不了超大規模，中等規模的問題還是可以解決的。使用正確的preconditioner，在需要求逆的地方用ILU方法或許可以解決幾百萬自由度的系統。（當然是在超算上）

並行

PETSc本身提供了分散式存儲的矩陣，並行求解器等一系列好東西。千萬不要傻乎乎地去用MPI親自擼一整套東西，不要在2017年去干1957年的活。不過即使有了這些，你的計算效率還是不一定有你想像中的那麼好，multigrid preconditioner可以讓求解的複雜度控制在O(n)級別（這已經是完美了），但你很可能實現不了，即便給你很多計算資源，也還是算不了超大的3D case。但是，中等規模的網格還是能算的。

最後我想說做這個事情需要你老闆的支持，要做好期望控制，別想著一下子搞出一個比商業軟體牛逼的東西。。。

謝邀 @崔忱

其實沒什麼經驗，還在用老舊的fortran在老舊的code上做一些微小的工作。講點小體會吧…

1，要不要造輪子？先要想清楚是自己寫還是在開源code上面二次開發，開源的openfoam蠻強大的，做很多科學問題完全可以在上面二次開發，毫不遜色。(身邊有些組在上面做出很優秀的工作)

2，要懂cfd，而不是一個碼農。CFL3D已經開源了，也可以從很多渠道弄到DLR以及老美的很多大規模並行code，所以寫cfd程序已經不是什麼秘密。我覺得最重要的是懂cfd，cfd是藝術勝於技術。重構，限制器，通量格式，時間推進格式，最核心的東西一定要完全弄明白，而不是只會抄公式寫一堆子函數。

3，問題導向。看你研究的課題，我覺得沒必要所有cfder都在幾年內寫一個大code。如果攻演算法格式，那就從一維寫起，二維就足夠，一步都不能含糊，一個算例一個算例的過。我們組側重格式，師兄姐都做二維居多，能有很多優秀的工作。如果你必須算三維大規模並行計算，那也得把一二維的算好，一二維的benchmark能暴露很多問題的其實…

4，各種細節。首先是網格的問題，結構化網格plot3d蠻流行，還有更通用的cgns。網格並行分塊工具，pestc還比較強大吧，只聞大名著實不懂，但是這一塊處理好可以給後期並行帶來很大的便利(想在code上加上，但時間精力都不允許)。數據結構的問題，結構網格還是非結構網格要先明確。並行採用MPI還是OpenMP還是時尚的GPU之類的，依據自己研究方向選吧，一般來說MPI是足夠的。重構方法幾個經典的一階，二階，WENO之類的是要有的吧，通量格式幾類經典的FVSFDS混合通量是要有的(更重要是知道他們的原理優劣勢應用範圍)，也有一些更出色的通量格式…時間推進，顯格式，LUSGS是比較流行的，還有一些更高級的(有一些著實不懂)

5，湍流模式。為何把湍流單列？不是因為自己在這個坑裡，也不是為了解這個困難的問題，而是實際中如果不算帶湍流的非定常問題，那談大規模並行的cfd程序就沒有太大意義了。定常的時代是還交給做格式的cfder吧(表示很羨慕在desktop玩cfd的，很高級很高級)，非定常的繞不開湍流，一方程sa兩方程kw類還有一些研究特定問題的湍流模式是大部分code的標配。les和dns隨著超算機器的強大，也不是不可接受，現在工程問題里des類方法也是很流行的。各個模型水也是很深的，可能是除網格之外cfd裡面最玄的一塊了吧。(想弄明白也不是很容易，其實弄的越明白會越失落越傷感越無奈)

5，團隊合作。我覺得做大code的確需要良好的團隊，大科學的時代即使數學家也需要更多的交流和刺激，何況這種妥妥的體力活。非數十年累計，數十人開發而不得，如NASA之流。我這方面沒有太多經驗，我的體會是如果蒙頭寫code改你們的code，我覺得畢業都是問題…當然如果完全沒有後顧之憂，開著大奔玩cfd，那這類工作於課題組是很有益的。

總之，不能一時衝動就決定翻新code，我們是做research的，code是為我們研究服務的，不是我們研究為code服務的。

要用最小的體力活實現最多的腦力活！(自己深深被打臉了其實)

每一個cfder都是很可愛的人，code的頑皮只有自己最懂，人艱不拆，像四爺學習努力打敗八阿哥！

心得：

這種爛攤子，重寫的工作量之巨大超乎想像。很多地方都不知道為啥那麼寫，一改就算錯了。

如果不是導師叫你重寫，那就不要重寫，不然就會對「吃力不討好」有深刻的體會。

每個博士生都覺得：師兄師姐的代碼是屎啊讓我來重新造輪子，投入無盡的光陰和精力。而這是大部分導師所極其不喜的。

代碼上的建議：

把模塊分清楚了一點一點改，準備好一堆標模，改一個函數跑一遍，不然到後面跑出來不一樣根本不知道是哪裡導致的。

PhD是以跑模擬、理解流體物理為主，還是以研究新的計算格式為主？如果是做模擬，但不喜歡前輩的代碼，那麼最好去找一套合適的library來用。如果是研究計算格式，並且導師給你開發新解法器的時間是2-3年以上，那麼真的有太多的東西要考慮，因為開發一套massively parallel DNS solver這項任務本身就足夠完成一個博士的培養。

過去三年里我的主要工作就是寫了一個「大規模三維並行計算流體求解器」，挺多教訓的。我覺得最重要的經驗是only learn from the best。關於開發這套求解器具體實踐的細節，只聽取了解這套具體演算法、這份代碼的人的意見，只聽取能夠花時間幫你一個問題一個問題解決的人的意見，這個人會是你的導師或者是實驗室里資歷最老的前輩。People think they know CFD they know numerical methods, but they don"t know what they are talking about. "If your method is not stable, why not implement an implicit solver?", "why don"t you try that boundary conditions", "why don"t you use FVM/FEM/AMR etc?", "check out this paper, it seems useful to you", "this stackexchange post seems solve your problem".... In CFD, things are much easier said than done.

如果身邊沒有這樣一個特別靠譜的人，我建議還是用現成的代碼，有問題可以在mailing list裡面問developers。

如果做「大規模三維並行計算流體求解器」不是為了拿PhD或者你只有1-2年的時間來做這份工作，我的建議是，do something else，life is short...

先不要急於動手。認真思考一下：

1、現有程序哪些地方不符合自己的要求，自己期望修改後的程序是什麼樣的？

2、現有的程序怎麼進行模塊化？可以使修改工作量最小。

3、哪些模塊可以採用第三方庫來代替，哪些模塊需要進行重寫？

4、是否需要多人合作？採用什麼版本控制工具？

不要重複造輪子不要重複造輪子不要重複造輪子

除非日後以後你打算靠此賺錢

解釋一下：

關鍵是你重寫的目的是為啥？要知道你在某組工作讀博士本質上你所寫出的東西你只有署名權你畢業之後的代碼產權你是沒有的。除非你自己寫自己的玩自己的然後把它乾脆開源掉否則in house代碼那做出來就是便宜組裡面了啊。後人雖然是萬分感激你但是你畢業之後這代碼該賣就賣也落不到你頭上。當然你說我為了提高cfd水平增長cfd見識。那你還不如直接去學讀成熟的代碼直接再上面大改就好。整個開源世界比你自己造輪子玩要學起來快的多。除非你所做的領域根本沒有啥開源玩。但是單純cfd 你一個剛入門的小博士生要挑戰這麼牛逼的沒有人目前涉足的cfd領域。。勇氣可嘉。。

所以除非你老闆能保證你付出若干心血寫完代碼之後的收益權（比如留組繼續作為生產力存在接班，或者有機會分錢）否則還不如自己去玩開源了事。嗯

看樣子可能是結構化網格的吧

如果一定要搞，別從最底層搞起

用現成的框架，比如SUMRAI，人家把並行、IO對接面處理都整好了

基本上你只需要關心分塊內部的操作和上層邏輯

題主用的什麼編譯器？建議直接用 Intel 的 MKL 函數庫進行求解，不需要自己寫求解器，只需要按照要求提供係數矩陣（可能是一位壓縮形式的）、右端向量等數組/向量，然後調用 MKL 的方程求解函數，比如說大名鼎鼎的 Pardiso 稀疏矩陣求解器進行求解。又快又方便。

我曾經也開發過基於非結構化網格的CFD軟體，但是目前開源的優秀軟體很多，後來工作上就長期使用OpenFOAM。個人認為開發工作固然有趣，但是一定要有所創新，要差異化競爭，否則代碼開發出來也不會有用戶。

我後來選擇了開發基於拉格朗日體系的無網格粒子法CFD軟體，支持混合併行。粒子法擺脫了網格的束縛，是下一代CFD發展的方向之一。歡迎大家關注我的知乎文章和微信公眾號。

第一種情況，你有實驗室的代碼。

1、能參考前人代碼就參考前人代碼

2、 CFD求解器的思路一般比較簡單，關鍵是並行化的時候如何share ghost cell等問題，這些會涉及到對MPI語言的運用。

3、讀懂之後自己寫一遍。這樣自己記住了一輩子忘不了。而且可以舉一反三，推而廣之做更複雜的問題。

第二種情況：沒有代碼，自己從頭寫。

0、自己在黑板上描述一下演算法流程圖，CFD演算法應該是可以自己講一遍能講出來每個時間布做什麼的。

1、學一下高級MPI IO的庫，HDF5。

2、先做周期邊界條件的問題。

3、寫好一維的代碼，實現所有功能，所有格式。TIP：所有輸出都寫出成數據文件，由單獨寫PYTHON來畫圖。

4、進行典型CASE驗證，提升PROCESS 數目後結果不變。

5、開始將邊界條件做文章，考慮其他邊界條件。

6、擴充到2D，進行算例驗證。

7、擴充到3D，基本上不會有什麼問題了。

先寫一維問題，然後寫二維，然後再寫三維，然後分別選擇比較經典的算例進行驗證，千萬別直接上二維或者三維。version control目前比較容易做到了。如果有能力，一定要做好regression test。Optimization是個大學問，不好講。

過來人經驗哈如果是改模塊那你隨意從頭編的話真心難受。你這是典型的強迫症跟以前的我一樣後來我發現這真心吃力不討好，對於老闆而言他又不關你過程咋樣只看結果。。結果就是自己顛倒黑白的干最後老闆來一句，哦你重新抄了一下程序啊，真行。。

關鍵是你能不能寫的比原來的好,如果只是好一丁點,好一兩倍,就算了,考慮是付出的時間和精力,還不如將就用原來的算了,要是你預期最少能比以前好上10倍,我完全贊同你重寫.

你要想想，20多年了有這個想法甚至付出了行動的不是你一人，為什麼大家還在湊合用祖傳代碼？如果不是導師要求，就別干這種體力活了，畢竟代碼寫得再漂亮也不是你自己的

1. 搞清楚、弄明白你們祖傳的代碼，其空間格式、時間格式細節，或者搞不懂細節，但可以把它們單獨寫為一個個單獨的子程序；

2.搞清楚祖傳代碼的數據結構，明白其對於並行的優缺點；在此基礎上設計適用於大規模並行的數據結構，這很重要；

3 程序IO的並行化很重要，初次改的時候可以先用串列IO，後面可以慢慢改為並行IO；

4.並行分區工具，是大規模並行的必要條件，初次改的時候可以手動分區，後面再加，不會對solver產生影響；

5.編程語言方面，要非常熟悉，對編程語言的熟悉程度會嚴重影響工作進程。

個人建議在動手之前，先搞清楚幾個問題。

1）你所在學校對於學生代碼的版權規定

研究生畢業論文相關的代碼，有些學校規定版權歸學生，有些規定版權歸學校。

2）代碼性質和計算規模

驗證數值方法？DNS? 工程問題？

課題組有多少計算資源？計算資源是什麼形式/架構？

3）痛點

相比於前代代碼，你想要解決哪些關鍵問題？

計算效率？並行效率？二次開發和維護性？數值格式？易用性和前後處理？擴展應用領域？

4）開發周期和代碼量

幾個人？幾條槍？搞多久？開發者的個人能力和知識背景如何？（數值方法、HPC、物理模型、代碼工程能力、寫大規模軟體的經驗）

6）開源 / 不開源

開源對學生作者有利（增加影響力和引用）對課題組不利（不好拿研究經費）

以上這些問題搞清楚之後，再決定要不要動手和怎麼動手。如果決定要寫，第一個問題就是技術選型。寫代碼和產品設計一樣，技術選型（或者說概念設計）決定了績效的80%。技術選型要考慮的問題包括：編程語言，數值格式，結構/非結構網格，並行模式。這些選擇決定了在什麼層面上使用何種第三方庫。最後的最後才是具體實現上的技術問題如何解決。

先看看程序怎麼用的，然後寫測試啊，保證各種情況都覆蓋掉，然後一點一點改，每次改後跑測試。。。。