全文約8000字，仔細閱讀約20分鐘

在偏微分方程求解的幾種數值方法中，有限元方法(FEM）最為通用，在工程領域應用也最為廣泛，沒有之一。

FEM有時候作Finite Element Method，

有時候Finite Element Model，

注意區分。

本文將重點討論有限元方法中的隱式方法，也就是需要求解線性方程組，且線性方程組的剛度矩陣一般是稀疏對稱矩陣。還是老規矩，盡可能少用專業術語和方程，軟硬件偏軟件介紹。

將從以下幾方面討論：

1. 有限元基本原理和特點

2. 數據的生產，存儲和管理

3. 業務計算流程

4. 并行處理和分布式處理

5. 線性方程組求解方法

6. 后處理

7. 優化改進

從顯式，隱式和準靜態說起

在力學分析中常用的求解方式有兩種：顯式(Explicit)和隱式(Implicit)。在動力學分析中，尤其是對時間比較敏感的場合，普遍采用顯式方法，其特點是不需要組裝整體剛度矩陣求解線性方程組，在時間上分別對單個單元進行連續分析，時間步驟可以取得非常小，不存在收斂性問題，容易編寫程序進行并行分布式計算。

而隱式則相反，適合處理處理大跨度時間，需要組裝整體剛度矩陣，求解線性方程組，同時存在收斂性問題，但是在時間選擇上需要考慮兩步之間的特定狀況，太大誤差偏大，太小浪費計算資源。比如在熱瞬態分析中，針對兩步溫度推導方程，歐拉參數取值不同，分別使用Crank-Nicolson和Backward Euler方法。且每一步計算過后都需要動態評估時間步長是否合理。

以上圖片來自網絡

準靜態(Quasi-Static)是一種簡化的瞬態分析方法，在某些特定場合可以用較少的資源獲得不錯的計算結果。比如在瞬態電磁場分析中，如果電磁場隨時間變化可以忽略；在動力學分析中，不考慮加速度等隨時間變化明顯的因素，或者把加速度進行適當轉換。把一個連續隨時間變化的過程分解成相對獨立的多步驟計算。每一個單步都可以方便獨立設置邊界荷載，降低了瞬態分析的復雜度和變量因素。

1. 有限元基本原理和特點

這里簡要介紹一下有限元方法基本原理和特點，有限元本質上是求解偏微分方程的一種方法。針對偏微分方程整體求解解析解困難的情況，有限元方法是將分析目標對象本身進行細化，劃分成小的有限個網格單元，比如三角形，四面體，六面體等，有時也會根據求解特點將高緯度幾何降為低緯度幾何單元，將構造的基函數和形函數應用在每一個網格單元，利用加權余量方法求解未知量。

有限元方法最早應用于結構力學分析，后來逐步推廣到熱，聲，電磁，流體等領域，是解決工程數值分析的通用方法。最早基于三角面片，后來推廣到各種形狀網格單元，主流網格單元為三角形，四邊形，四面體和六面體。目前大部分介紹有限元的教程集中在結構領域，對多物理場以及多物理場耦合介紹的非常少，這是國內教材可以加強的一個方向。

有限元方法自誕生以來，在理論基礎，物理，數學，工程應用等都得到了驗證，是目前求解偏微分方程以及多物理場問題最完備，最有效，最廣泛的數值計算方法。

需要提及的是我國的胡海昌于1954年提出了廣義變分原理，錢偉長最先研究了拉格朗日乘子法和廣義變分原理之間關系。而馮康更是在20世紀60年代于國內獨立提出了有限元的方法和概念，并研究了有限元分析的精度，邊界條件以及收斂性問題，指導了國內多項重大工程實踐。（1965年馮康發表了論文“基于變分原理的差分格式”，這篇論文是國際學術界承認我國獨立發展有限元方法的主要依據）。知Zienkiewicz而不知馮康，是我們科普和教材需要改進的地方。

關于有限元更詳細的介紹，點擊鏈接查看

一篇文章入門多物理場有限元（全篇）

2. 數據的生產，存儲和管理

大模型處理首先要考慮的是數據結構和數據存儲方式，我們以1億自由度double實數存儲為例，1億為1e8，矩陣規模為1e8*1e8，即1e16，假設稀疏矩陣非0數據為總數的0.01%(1e-4)，即1e12。double類型數據占用8字節，內存占用數據為8e12，即8000G，考慮到編譯器設置，頁面對齊，硬件管理等因素，實際占用數據只會多不會少。

在以往實際用例中，一般自由度越大，其稀疏矩陣中非零數據比例越低。這個比例并沒有固定規律，和仿真類型，幾何形狀，邊界，荷載數量，網格大小密度質量等都有關系。假設實際情況其比例在1e-6到1e-4之間，即使取最小1e-6，其1億自由度的內存消耗也在80G以上，普通臺式機無法存儲。這還僅僅是一個總剛矩陣數據的存儲，沒有涉及到任何計算內存消耗。

稀疏矩陣的存儲方式

對于一般矩陣，我們常使用一個二維數組來定義，通過行列和索引來存儲數據，如：

const int NUMSIZE = 1e6;double value[NUMSIZE][NUMSIZE];

棧存儲不了太多的元素，一般在堆上分配數據，即

const int NUMSIZE = 1e6;double **value;value = new double*[NUMSIZE];int i;for(i=0;i<NUMSIZE;i++){value[i]=new double[NUMSIZE];}

但是稀疏矩陣使用該方法非常浪費空間，也無必要，實際分析中稀疏矩陣有其自定義的優化存儲方式。

稀疏矩陣常用的存儲方式有：

1. Coordinate Format(COO)

該方法采用了三個一維數組，分別存儲行索引，列索引和對應的數值

int x[NUMSIZE];int y[NUMSIZE];double value[NUMSIZE];

C++STL標準模板庫中提供了map容器，該容器可以將一對數據作為key存儲。也就是

map<pair<int,int>,double> values;

表面上看起來，通過pair來存儲行和列索引，通過map來取對應數值，很適合存儲矩陣數據，但是map使用的是紅黑樹結構，也就是插入數據時會對pair排序，這就導致數據越多，后期插入性能損耗和內存開銷越嚴重。

下圖展示了分別使用二維數組和map存儲double型數據隨元素個數N，也就是NUMSIZE值增加的內存開銷。下圖可以看出到后期，相比數組，std::map的內存開銷隨N的增加呈指數級增加。

實際存儲double數據個數為N*N

在對矩陣數據操作中，一般不會使用STL或第三方庫，而是直接使用數據指針。使用過矩陣求解庫的朋友應該也有體會，函數參數以及數據傳輸一般都是直接使用地址操作，矩陣操作首先要保證性能。這也是為什么很多庫都會使用C風格的編碼方式，即使使用C++，也主要是用在管理數據上，實際數據操作都是偏向C風格。

2. Compressed Sparse Row (CSR)和

Compressed Sparse Column(CSC)

COO方法雖然降低了內存需求，但從數據壓縮角度看，并不是最優。CSR和CSC數據結構可以進一步壓縮數據，是最常用的稀疏矩陣存儲方式。

以上圖片來源于網絡

以CSR為例， 在上圖4*4矩陣中：

row offsets表示行偏移值：

1. 第一行從0開始，偏移值為0；
   

2. 第二行，之前有兩個值1和7，因此偏移值為2；

3. 同理，第三行，前面有1，7，2，8四個數，第三個偏移為4
   

column indices和values則分別表示列索引和對應值：

比如按順序，第0列的第一個數為1， 第1列的第一個為7，第二個為2。以此類推。

可以看出CSR是一種全局編碼，相比COO，進一步減少了行索引值占用數據。CSC數據結構類似。

3. DIA(Diagonal Sparse Matrix)

即對角稀疏矩陣格式

對于矩陣元分布在主對角線以及主對角線附近的稀疏矩陣，采取這種存儲方式比較有效，但通用性不強。

還有其它很多稀疏矩陣存儲方式，一般是以上幾種方式的衍生，擴充或者組合。

稀疏矩陣的數據結構定義

矩陣數據結構是最基礎的數據結構。從計算角度看，稀疏矩陣一般單獨定義，單獨實現接口；但是從平臺和系統架構看，最好是將數據結構融入到系統設計中。

看幾個關于矩陣計算的需求：

1. 普通矩陣和稀疏矩陣進行加法計算

2. 將自定義的稀疏矩陣數據傳到不同的線性方程組求解庫求解

3. 稀疏矩陣求逆

4.按照業務封裝成模塊供其它模塊調用

這些需求從計算角度看是要盡量避免的，但是從業務和程序邏輯看，又是合理和很可能出現的。

一般矩陣結構定義

C++中提供了template模板功能，可以使用相同的數據結構形式，處理不同的數據類型(double, float, int, complex等)，在矩陣定義中可以靈活使用。

裸指針 or 智能指針？

C++中，我們把直接使用關鍵字new在堆上分配內存創建的指針對象稱為裸指針，顧名思義也就是沒有任何裝飾，裸指針在使用結束后需要顯式調用delete刪除對象。

使用裸指針最大的挑戰在于如何管理對象資源。如果數據的生命周期非常短，且作用域非常明確，那問題不大；相反，如果數據貫穿程序大部分生命周期，且容易出現拷貝，賦值，在第三方數據庫中操作使用，如果不清楚指針的生命周期和所有權，就非常容易出現空指針，野指針以及內存泄漏等各種問題。

為了從根本上解決指針操作出現的各種問題，C++先后引入了各種智能指針，包括auto_ptr，unique_ptr，shared_ptr等。而一般的開源庫軟件也都有自己的智能指針，

比如VTK的vtkSmartPointer，OCC的handle等。

在一般業務流程操作和小規模數據操作中，推薦使用智能指針；求解器中對大數據的操作仍然建議使用裸指針。

在以往非求解器和非UI開發代碼規范中，會強制使用unique_ptr

3. 業務計算流程

典型的一次業務計算流程：

1. 解析網格數據

2. 創建單剛矩陣

3. 組裝總剛矩陣

4. 組裝各種邊界荷載激勵

5. 求解線性方程組

6. 提取結果數據后處理

這樣的一次流程在軟件開發中需要進行封裝，便于網格加密迭代，優化算法，掃頻，或者多步瞬態方法等外部程序調用。盡可能做成獨立通用的模塊，即不依賴具體的業務，不依賴特定的算法，在使用過程中按照外部輸入參數定制。

不同物理場由于基函數取值不同，實際組裝內容是不同的，不同的仿真類型和網格單元，其單剛矩陣具體內容也非常不同，計算單剛矩陣和組裝總剛矩陣是有限元計算方法的核心業務。在常用的三種最基本類型邊界和荷載基礎上，不同物理場會衍生出多種具體形式的邊界荷載，進一步結合幾何模型，仿真類型，實際操作方法等，會固定出業務上的常用邊界荷載，這也是商業軟件的一大突出特點。

不同物理場計算其剛度矩陣組裝過程大同小異，在求解器架構設計和開發中是可以流程化和平臺化的。

4. 并行和分布式處理

拋開算法，從軟件角度看，并行和分布計算體現在以下幾個層次：

1.業務并行

這個是最容易理解的，業務并行的最大特點是每個任務之間沒有互斥的數據，前后順序依賴關系也比較弱。比如在DOE設計，掃頻，優化設計中，各個模型之間共享模型數據，只是某些求解參數設置不同。所以很容易將每個任務派發到不同機器節點上執行。這里沒有涉及到數據交互，網絡帶寬，負載均衡等，是最容易的并行計算。

2.Culster集群式并行

這種主要是針對大型機，服務器以及HPC集群。典型的是利用MPI在不同核心，節點上進行計算。節點之間的數據高度相關，順序依賴關系強。需要對模型進行分治，數據分塊；節點之間數據需要進行通信。根據模型數據和硬件特點選擇合適的分治策略就非常重要。常見例子，某些大線性方程組求解在節點達到一定數目后，其加速比呈下降趨勢，原因在節點多到一定程度后，分治數據之間的通信開銷開始成為性能瓶頸。

3.單機并行

在單機上，CPU多核已經非常普遍，利用多線程，多進程可以輕松做到在不同CPU核上并行計算。而單機的線程工具已經非常普及，比如OpneMP，C++11新增線程類，QT提供的線程進程操作工具。對于類似共享內存式并行，網絡以及通信開銷可以忽略不計。

4.指令優化

單指令流多數據流機器（SIMD）

SIMD是采用一個指令流處理多個數據流。這類機器在數字信號處理、圖像處理、以及多媒體信息處理等領域非常有效。

Intel處理器實現的MMXTM、SSE（Streaming SIMD Extensions）、SSE2及SSE3擴展指令集，都能在單個時鐘周期內處理多個數據單元。也就是說我們現在用的單核計算機基本上都屬于SIMD機器。

多指令流多數據流機器（MIMD）

MIMD機器可以同時執行多個指令流，這些指令流分別對不同數據流進行操作。最新的多核計算平臺就屬于MIMD的范疇，例如Intel和AMD的雙核處理器等都屬于MIMD。

指令優化主要針對具體硬件的指令支持情況，程序中合理使用指令集，計算可以起到很好的加速效果。

5.異構架構式并行
DPU（Data Processing Unit）是以數據為中心構造的專用處理器，采用軟件定義技術路線支撐基礎設施層資源虛擬化，支持存儲、安全、服務質量管理等基礎設施層服務。2020年NVIDIA公司發布的DPU產品戰略中將其定位為數據中 心繼CPU和GPU之后的“第三顆主力芯片”，掀起了一波行業熱潮。DPU的出 現是異構計算的一個階段性標志。與GPU的發展類似，DPU是應用驅動的體系 結構設計的又一典型案例；但與GPU不同的是，DPU面向的應用更加底層。

從工程學角度看，并行計算的本質是提升性價比，GPU的應用對于數值計算而言是異構架構一大成功。NVIDIA提供的CUDA計算語言成為GPU計算的主流工具。

有限元求解器主要操作對象是矩陣，對于大矩陣計算，涉及到的是矩陣的分塊和分解以及并行計算。矩陣的分塊和分解是兩個不同的概念，注意區分。

求解線性方程組，不管是迭代法還是直接法，最多的運算是矩陣和矩陣乘法，矩陣和向量乘法。乘法的并行計算是各個基礎庫的計算核心功能之一。

圖畫分

對于一個超大矩陣，通常需要分塊，分塊的目的是在已有硬件資源基礎上，讓每塊獲得合理的計算資源，達到并行計算的目的。而并行計算一個核心問題就是負載均衡，即整體上每個計算資源要與計算內容匹配，避免部分硬件資源長期閑置。這就要用到圖畫分工具。

圖論：一般將網格單元抽象成一張圖的頂點，頂點的權重代表網格單元的計算量，網格單元之間的數據依賴關系抽象為圖的邊，邊的權重代表網格單元之間需要交換的數據量，這張圖反映了問題的基本計算量和數據依賴關系，可以作為該問題的負載模型。對于稀疏矩陣，該圖可以描述向量乘法的基本計算量和數據依賴關系，稀疏矩陣向量乘法的負載均衡可以轉化為圖的多路劃分問題。在一些線性方程求解庫中，可以看到相應的圖劃分工具。

5. 線性方程組求解

線性方程組求解是隱式有限元方法求解的核心。考慮到大規模稀疏對稱矩陣線性方程組求解方法的復雜性以及求解庫繁多，作為系列文章，后面專門分多文介紹。先簡要介紹兩個解線性方程組的工具庫，MUMPS和OpenBLAS，具體使用方法不展開。

1.MUMPS網址：

http://mumps-solver.org/

最新版是5.5.0(2022年4月)

MUMPS全稱

MUltifrontal Massively Parallel sparse direct Solver

是一款經歷了長期實踐檢驗，專業可靠的求解大規模稀疏矩陣線性方程組的工具庫，在很多商業以及開源CAE/EDA/CFD等數值仿真軟件中都有采用。

MUMPS采用的是直接法，支持串行和并行計算，底層編碼語言Fortran和C。

MUMPS求解過程：

1.主進程執行排序和符號分解；

2.進行LU分解或LDLT分解，將矩陣A分布到多個處理器上，每個處理器對每個波前矩陣進行數值分解

3.主進程將B和x分布到各個處理器進行因子計算后

在以往超大規模的熱，結構動力和高頻電磁有限元仿真計算中，相比其它求解庫，MUMPS表現出了內存消耗少，速度快，負載均衡合理等優勢。

2.OpenBLAS網址：

http://www.openblas.net/

OpenBLAS 是一個優化的 BLAS 庫，基于 GotoBLAS2 1.13 BSD 版本(該版本不再維護更新)

BLAS（Basic Linear Algebra Subprograms 基礎線性代數程序集）是一個應用程序接口（API）標準（可以理解為C++里函數做了聲明，但沒有定義，大家都可以來完成定義實現），BLAS本身也有一個實現。BLAS用以規范發布基礎線性代數操作的數值庫（如矢量或矩陣乘法），其程序集最初發布于 1979 年，并用于建立更大的數值程序包（如 LAPACK）。在高性能計算領域，BLAS 被廣泛使用，例如，LINPACK 的運算成績則很大程度上取決于 BLAS 中子程序 DGEMM 的表現。為提高性能，各軟硬硬件廠商針對其產品對 BLAS 接口實現進行高度優化，大的的廠商譬如INETL的MKL，NVDIA的CUBLAS。

參考之前的文章，點擊鏈接查看

一篇文章入門大規模線性方程組求解

6. 后處理

數據的后處理涉及到的也是關于大數據的公式計算，查詢，提取以及存儲。有限元方法計算的結果通常放在節點，邊或者高斯積分點上，根據實際需要將其映射到對應單元節點上，即可方便實現繪圖。

1.后處理的計算一般依據已有公式，沒有太復雜的算法。對于大數據，單機可方便使用共享內存OpenMP即可實現并行加速。圖形渲染上可以參考Paraview的客戶機/服務器模式。

2.在后處理中，盡量避免全局計算：用戶通常感興趣的是物理量最大最小值以及一些特殊特定幾何位置的數值。合理設置探測點，同時通過業務邏輯過濾，減少數據的計算，避免全局任何n方以及以上的操作。

3.之前介紹的一篇文章入門仿真軟件性能優化(點擊鏈接查看)對后處理仍然有效。利用分治原則，大數據盡可能通過數量多的小文件存儲，數據存儲可以使用加密明文，提升壓縮效率。對于三維網格，可以優先只計算表面網格結果，在實際需要的時候再展開內部三維結果數據。

7. 優化改進

1.設計模式

點擊鏈接查看

為什么求解器開發需要強勢使用設計模式

設計模式本質上是一種代碼規范，它要求開發人員針對某類問題能使用相同或類似的設計思路和方法。這樣做的好處是：

1.能反映出問題的本質，幫助做出更好的設計；

2.方便開發人員之間協作和溝通；

3.有利于提升軟件質量和代碼維護。

求解器開發并不僅僅是底層算法實現，還包括了大量業務和應用層面的需求。比如線性方程組的庫的包裝，切換，更新；網格自適應迭代集成；優化設計模塊集成；第三方平臺調用；硬件參數自適應；程序穩定后的功能維護擴展或者重構等等。

這些非常需要從軟件工程和系統工程角度考慮問題和設計，而非僅僅有限元算法本身。這也是長期以來強調軟件開發要重設計架構，輕語言代碼的體現。

2.調試

大模型的求解器調試一直是非常有挑戰性的難題，主要原因還是數據量太大，導致平常一些不起眼的小功能也能導致性能瓶頸，比如文件讀取，網格解析，有效性檢查，中間數據導出等。很多的這種性能瓶頸集合在一起，最終的結果就是一天可能只能調試2-3次，嚴重影響開發效率。可以從以下幾個方面改進調試：

1. 完善模型數據導出。在任意一步計算中，可以方便地將數據導出為常用的數據格式，比如CSV，MATLAB文件格式。這個操作和模型大小無關，是基礎性的功能開發，需要保證高效和正確性。

2. 詳細的日志系統。日志系統需要詳細記錄程序運行狀態，包括時間，硬件資源內存，CPU，硬盤，網絡使用狀況，每一步程序是否運行成功，按照狀態分級給出信息，警告，錯誤等具體信息。提供實時硬件運行狀態是大模型仿真基本功能之一。一般服務器，大型機都會提供接口。

3. 模塊化分級。將整個計算仿真流程模塊化，按照功能有重要程度分級，每個模塊能夠通過文件提供接口。這樣做的好處是流程中每一步可以做有效性檢查，一旦出現錯誤可以快速定位出問題的步驟。也方便更新golden模型。

4. 獨立的數據分析模塊。該模塊可以獨立提供整體模型的數據特點，包括網格質量分析(求解角度)，敏感度分析，各種矩陣特征分析等。有些功能第三方庫會提供。

3.硬件使用

前面講過超大模型的計算和硬件緊密關聯，有些工具庫甚至需要在運行機器上編譯，根據硬件實際情況優化后使用。根據業務合理的選擇硬件和軟件資源是加速求解的關鍵。

本文介紹了超大模型有限元求解器計算方法的一些研發知識，可以作為有限元方法工業級應用開發的入門參考。需要說明的是，超大規模的有限元模型求解方法非常依賴模型數據的特點，并沒有一個黃金標準方法，需要在實踐中選擇合適的方法。

后話：發現一個有趣的現象，現在很多研究高性能矩陣計算和線性方程組求解的，已經不是在工業自動化，工業軟件數值計算這塊，而是在最近幾年火熱的AI計算，深度學習方面。正如工業設計仿真軟件與數學和多物理場與數學(1)--求解器開發的數學基礎(點擊鏈接查看)中所說的：學科之間都是相通的，而這個“通”就是靠數學聯系起來，很多看起來不相關的應用學科，底層技術都是類似或相同。