工业软件研发中处理超大模型(6)--有限元求解器

全文约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)--求解器开发的数学基础(点击链接查看)中所说的：学科之间都是相通的，而这个“通”就是靠数学联系起来，很多看起来不相关的应用学科，底层技术都是类似或相同。