CN117217062A - 基于刚度矩阵的流体仿真方法及装置 - Google Patents

Abstract

本发明公开了一种基于刚度矩阵的流体仿真方法及装置。其中，该方法包括：确定目标流体在问题域上的有限元方程，其中，问题域包括多个有限元，每个有限元对应问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；配置静态化的坐标COO格式的稀疏张量类模板，其中，稀疏张量类模板的模板参数包括张量维数和非零元的位置；根据有限元方程和稀疏张量类模板计算单位刚度矩阵；采用单位刚度矩阵在处理器中装填生成问题域的总刚度矩阵；采用总刚度矩阵模拟计算目标流体在几何形状中的流动特征参数。通过本发明，通过提高张量解析算力，加快单元刚度矩阵和总刚度矩阵的计算速度，从而提升了仿真速度和仿真效率。

Description

基于刚度矩阵的流体仿真方法及装置

技术领域

本发明涉及数字建模领域，具体而言，涉及一种基于刚度矩阵的流体仿真方法及装置。

背景技术

相关技术中，有限元方法是一种数值分析方法，用于求解复杂物理和工程问题中的偏微分方程，通过将问题域离散化为许多小的、简单的子域，即"有限元"，用于求解复杂结构的物理问题。它在工程学和科学领域中得到广泛应用，尤其是在结构力学、流体力学、热传导和电磁场等领域。其中，刚度矩阵计算是有限元方法计算程序中的关键部分，而单元刚度矩阵计算过程中，通常会选取多个积分点处，进行矩阵的元素计算；而这一矩阵元素的计算过程，涉及多个不同维度的张量缩并运算，由于张量的阶数最大可以是4阶，这一运算开销大。

相关技术中在单元刚度矩阵计算过程中，采用稠密张量表示和计算单位刚度矩阵，导致计算量大，进而导致模型仿真的速度慢，效率低。

针对相关技术中存在的上述问题，目前尚未发现有效的解决方案。

发明内容

本发明实施例提供了一种基于刚度矩阵的流体仿真方法及装置。

根据本申请实施例的一个方面，提供了一种基于刚度矩阵的流体仿真方法，包括：确定目标流体在问题域上的有限元方程，其中，所述问题域包括多个有限元，每个有限元对应所述问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；配置静态化的坐标COO格式的稀疏张量类模板，其中，所述稀疏张量类模板的模板参数包括张量维数和非零元的位置；根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵；采用所述单位刚度矩阵在处理器中装填生成所述问题域的总刚度矩阵；采用所述总刚度矩阵模拟计算所述目标流体在几何形状中的流动特征参数。

进一步，配置静态化的COO格式的稀疏张量类模板包括：配置正数类型的第一参数包；将所述第一参数包中的一个第一参数选择为张量维数d，并将所述第一参数的所在位置之后的每d个数配置为一个非零元的位置，得到静态化的COO格式的第一稀疏张量类模板。

进一步，配置静态化的COO格式的稀疏张量类模板包括：配置正数类型的第二参数包；将所述第二参数包中的一个第二参数选择为张量维数d，并在所述第二参数的所在位置之后配置常量表达式std::array<int, d>，采用所述常量表达式配置每个非零元的位置，得到静态化的COO格式的第二稀疏张量类模板。

进一步，根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵包括：采用有限元方程的基函数和网格参数推导得到单元刚度矩阵元素方程和缩并维度；对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合;基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果；采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵。

进一步，对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合，包括：判断所述稀疏张量类模板中的非零元是否有解析表达式；若所述稀疏张量类模板中的非零元有解析表达式，利用静态for循环遍历所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合；若所述稀疏张量类模板中的非零元没有解析表达式，显式写出所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合。

进一步，基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果，包括：判断所述缩并维度是否为标量；若所述缩并维度是否为标量为标量，采用C++语法中的折叠表达式对所述稀疏张量集合各个非零元的值乘以相应的数后求和，得到中间结果。

进一步，采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵，包括：确定以下中间结果：、/>、/>、/>，其中，/>是形状为Field*Field的稀疏张量，/>为形状Field * Dim * Field的稀疏张量，/>是形状为Dim* Field * Field的稀疏张量，/>为形状Dim * Field * Dim * Field的稀疏张量，Dim为所述问题域的维数，Field为所述有限元方程的数量；采用所述中间结果配置以下单元刚度矩阵元素方程/>：，得到单位刚度矩阵；其中，∫表示数值积分，K为网格单元，/>为基函数，i、j分别为网格单元K内的行和列。

根据本申请实施例的另一个方面，还提供了一种基于刚度矩阵的流体仿真装置，包括：确定模块，用于确定目标流体在问题域上的有限元方程，其中，所述问题域包括多个有限元，每个有限元对应所述问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；配置模块，用于配置静态化的坐标COO格式的稀疏张量类模板，其中，所述稀疏张量类模板的模板参数包括张量维数和非零元的位置；计算模块，用于根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵；生成模块，用于采用所述单位刚度矩阵在处理器中装填生成所述问题域的总刚度矩阵；模拟模块，用于采用所述总刚度矩阵模拟计算所述目标流体在几何形状中的流动特征参数。

进一步，所述配置模块包括：第一配置单元，用于配置正数类型的第一参数包；第一选择单元，用于将所述第一参数包中的一个第一参数选择为张量维数d，并将所述第一参数的所在位置之后的每d个数配置为一个非零元的位置，得到静态化的COO格式的第一稀疏张量类模板。

进一步，所述配置模块包括：第二配置单元，用于配置正数类型的第二参数包；第二选择单元，用于将所述第二参数包中的一个第二参数选择为张量维数d，并在所述第二参数的所在位置之后配置常量表达式std::array<int, d>，采用所述常量表达式配置每个非零元的位置，得到静态化的COO格式的第二稀疏张量类模板。

进一步，所述计算模块包括：求解单元，用于采用有限元方程的基函数和网格参数推导得到单元刚度矩阵元素方程和缩并维度；运算单元，用于对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合；生成单元，用于基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果；配置单元，用于采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵。

进一步，所述运算单元包括：判断子单元，用于判断所述稀疏张量类模板中的非零元是否有解析表达式；第一生成子单元，用于若所述稀疏张量类模板中的非零元有解析表达式，利用静态for循环遍历所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合；第二生成子单元，用于若所述稀疏张量类模板中的非零元没有解析表达式，显式写出所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合。

进一步，所述生成单元包括：判断子单元，用于判断所述缩并维度是否为标量；缩并子单元，用于若所述缩并维度是为标量，采用C++语法中的折叠表达式对所述稀疏张量集合各个非零元的值乘以相应的数后求和，得到中间结果。

进一步，所述配置单元包括：确定子单元，用于确定以下中间结果：、、/>、/>，其中，/>是形状为Field*Field的稀疏张量，/>为形状Field * Dim * Field的稀疏张量，/>是形状为Dim * Field * Field的稀疏张量，为形状Dim * Field * Dim * Field的稀疏张量，Dim为所述问题域的维数，Field为所述有限元方程的数量；配置子单元，用于采用所述中间结果配置以下单元刚度矩阵元素方程/>：/>，得到单位刚度矩阵；其中，∫表示数值积分，K为网格单元，/>为基函数，i、j分别为网格单元K内的行和列。

根据本申请实施例的另一方面，还提供了一种存储介质，该存储介质包括存储的程序，程序运行时执行上述的步骤。

根据本申请实施例的另一方面，还提供了一种电子设备，包括处理器、通信接口、存储器和通信总线，其中，处理器，通信接口，存储器通过通信总线完成相互间的通信；其中：存储器，用于存放计算机程序；处理器，用于通过运行存储器上所存放的程序来执行上述方法中的步骤。

本申请实施例还提供了一种包含指令的计算机程序产品，当其在计算机上运行时，使得计算机执行上述方法中的步骤。

通过本发明，确定目标流体在问题域上的有限元方程，其中，问题域包括多个有限元，每个有限元对应问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵，配置静态化的坐标COO格式的稀疏张量类模板，其中，稀疏张量类模板的模板参数包括张量维数和非零元的位置，根据有限元方程和稀疏张量类模板计算单位刚度矩阵，采用单位刚度矩阵在处理器中装填生成问题域的总刚度矩阵，采用总刚度矩阵模拟计算目标流体在几何形状中的流动特征参数，通过配置并采用稀疏张量类模板计算目标流体在问题域上的单位刚度矩阵，采用稀疏化的张量减少了计算量，解决了相关技术采用稠密张量表示单位刚度矩阵导致模拟计算流体流动特性速度慢的技术问题，通过提高张量解析算力，加快单元刚度矩阵和总刚度矩阵的计算速度，从而提升了仿真速度和仿真效率。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明实施例的一种计算机的硬件结构框图；

图2是根据本发明实施例的一种基于刚度矩阵的流体仿真方法的流程图；

图3是本发明实施例中计算单位刚度矩阵的流程图；

图4是本发明实施例中静态化COO稀疏张量的模板实例化的示意图；

图5是根据本发明实施例的一种基于刚度矩阵的流体仿真装置的结构框图。

具体实施方式

为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分的实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

需要说明的是，本申请的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本申请的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

实施例1

本申请实施例一所提供的方法实施例可以在服务器、计算机、或者类似的运算装置中执行。以运行在计算机上为例，图1是本发明实施例的一种计算机的硬件结构框图。如图1所示，计算机10可以包括一个或多个（图1中仅示出一个）处理器102（处理器102可以包括但不限于微处理器MCU或可编程逻辑器件FPGA等的处理装置）和用于存储数据的存储器104，可选地，上述计算机还可以包括用于通信功能的传输设备106以及输入输出设备108。本领域普通技术人员可以理解，图1所示的结构仅为示意，其并不对上述计算机的结构造成限定。例如，计算机10还可包括比图1中所示更多或者更少的组件，或者具有与图1所示不同的配置。

存储器104可用于存储计算机程序，例如，应用软件的软件程序以及模块，如本发明实施例中的一种基于刚度矩阵的流体仿真方法对应的计算机程序，处理器102通过运行存储在存储器104内的计算机程序，即本发明实施例中的一种基于刚度矩阵的流体仿真方法对应的计算机程序，从而执行各种功能应用以及数据处理，即实现上述的方法。存储器104可包括高速随机存储器，还可包括非易失性存储器，如一个或者多个磁性存储装置、闪存、或者其他非易失性固态存储器，用于存储数据库，数据库中包括搜索空间中的高维数据，以及存储每次迭代周期的查找结果和最终的查找结果。在一些实例中，存储器104可进一步包括相对于处理器102远程设置的存储器，这些远程存储器可以通过网络连接至计算机10。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。

传输设备106用于经由一个网络接收或者发送数据。上述的网络具体实例可包括计算机10的通信供应商提供的无线网络。在一个实例中，传输设备106包括一个网络适配器（Network Interface Controller，简称为NIC），其可通过基站与其他网络设备相连从而可与互联网进行通讯。在一个实例中，传输设备106可以为射频（Radio Frequency，简称为RF）模块，其用于通过无线方式与互联网进行通讯，用于从人机交互设备向处理器102传输基于刚度矩阵的流体仿真指令，以及从处理器102向显示设备传输查找结果。

在本实施例中提供了一种基于刚度矩阵的流体仿真方法，图2是根据本发明实施例的一种基于刚度矩阵的流体仿真方法的流程图，如图2所示，该流程包括如下步骤：

步骤S202，确定目标流体在问题域上的有限元方程，其中，问题域包括多个有限元，每个有限元对应问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；

在本实施例中，目标流体可以是空气、水流等流体，目标流体一个非线性有限元问题可以表示成如下格式：，求解(公式 1)，其中，F为由控制方程决定的多元线性或非线性函数。/>是计算区域；/>是计算区域中的点；待求解的/>是关于x的函数。/>是梯度算子；/>是拉普拉斯算子；B为由边界条件决定的多元线性或非线性函数。/>是计算区域的边界。/>是有限元函数空间。

选取检验函数。/>是在/>内，且在第一类边值的边界上值为0的函数。将/>乘以公式（1）式两侧，使用分部积分后可以获得如下形式：(公式2)，其中为关于x, u的表达式，由F和B决定。其中“·”号表示缩并维数为1的张量缩并运算，具体包括(a)向量内积(b)向量乘以矩阵(c)向量与3阶以上张量对应维度的缩并运算。其中“:”符号表示缩并阶数为2的张量缩并运算，具体包括(a)矩阵与矩阵对应元素相乘求和（b）矩阵与3阶以上张量对应维度的缩并运算。

设有限元解函数,其中/>为有限元基函数，/>为待求解的系数，则取/>即得关于/>的线性或非线性方程组/>}。其中/>的表达式由/>以及有限元基函数决定。上述线性或非线性方程组的求解需要计算}关于/>的偏导，即雅可比矩阵/>，记为/>,也即总刚度矩阵，需要在每个单元网格上计算单元刚度矩阵并关于每个单元求和。下面以内部单元网格为例，边界单元网格的积分方式与内部单元网格类似，对于网格单元K,其单元刚度矩阵非零元素：/>（公式3），其中i,j可取的值由单元K决定，K上的积分可以采用高斯积分点计算数值积分，对于积分点x，由变分形式和基函数可以相应的计算出/>。

设问题域的维数为Dim, 即,问题域的原方程数为Field, 则/>为形状为/>的四维张量，/>为/>的矩阵。如采用稠密张量，设单元内自由度数为L，高斯积分点数为M,则张量缩并/>的计算复杂度为/>。

步骤S204，配置静态化的坐标COO格式的稀疏张量类模板，其中，稀疏张量类模板的模板参数包括张量维数和非零元的位置；

步骤S206，根据有限元方程和稀疏张量类模板计算单位刚度矩阵；

由于计算复杂度大，而实际问题中，有很多位置是0,因此，可以采用稀疏化的张量,以减少计算。

步骤S208，采用单位刚度矩阵在处理器中装填生成问题域的总刚度矩阵；

步骤S210，采用总刚度矩阵模拟计算目标流体在几何形状中的流动特征参数；

本实施例的方案应用在流体力学的仿真中，流体力学采用的有限元方程为的Navier-Stokes（纳维-斯托克斯）方程为例，本实施例的方案同样适用于固体力学、电磁学等其他领域的有限元方法。在流体力学中，有限元方法可以用于模拟和分析流体在复杂几何形状中的流动行为。流动特征参数可以是速度、压力、温度和浓度等参数，采用总刚度矩阵可以实时仿真或者预测流体的速度、压力、温度和浓度等分布，了解流体流动的特性和行为，具体如风场模拟、水力学问题、空气动力学和气候模拟等。

通过上述步骤，确定目标流体在问题域上的有限元方程，其中，问题域包括多个有限元，每个有限元对应问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵，配置静态化的坐标COO格式的稀疏张量类模板，其中，稀疏张量类模板的模板参数包括张量维数和非零元的位置，根据有限元方程和稀疏张量类模板计算单位刚度矩阵，采用单位刚度矩阵在处理器中装填生成问题域的总刚度矩阵，采用总刚度矩阵模拟计算目标流体在几何形状中的流动特征参数，通过配置并采用稀疏张量类模板计算目标流体在问题域上的单位刚度矩阵，采用稀疏化的张量减少了计算量，解决了相关技术采用稠密张量表示单位刚度矩阵导致模拟计算流体流动特性速度慢的技术问题，通过提高张量解析算力，加快单元刚度矩阵和总刚度矩阵的计算速度，从而提升了仿真速度和仿真效率。

本实施例的方案采用C++模板元编程实现，因此无需引入领域专用语言，使得有限元编程中可以直接调用现有C++库，增加了编程的灵活性；同时，避免了代码生成的过程和对第三方JIT技术的依赖，使得该方案可以轻松地部署到任意提供C++17以上标准编译器的计算平台。

在本实施例的一个实施方式中，以C++17为例，配置静态化的COO格式的稀疏张量类模板包括：配置正数类型的第一参数包；将第一参数包中的一个第一参数选择为张量维数d，并将第一参数的所在位置之后的每d个数配置为一个非零元的位置，得到静态化的COO格式的第一稀疏张量类模板。

在本实施例的另一个实施方式中，以C++20为例，配置静态化的COO格式的稀疏张量类模板包括：配置正数类型的第二参数包；将第二参数包中的一个第二参数选择为张量维数d，并在第二参数的所在位置之后配置常量表达式std::array<int, d>，采用常量表达式配置每个非零元的位置，得到静态化的COO格式的第二稀疏张量类模板。

定义静态化的COO格式稀疏张量类模板。其中，张量的维数和非零元的位置作为模板参数。本实施例除了上述基于C++17和C++20标准的实现方式之外，也可以采用其他将非零元位置作为模板参数的实现方式。

如使用C++17标准，使用正数类型的参数包(parameter pack)作为模板参数，将参数包中的一个参数作为张量维数（记为d），之后每d个数视作一个非零元的位置，如有nnz个非零元，则参数包的长度为nnz*d + 1。如使用C++20标准，则可以使用一个参数作为张量维数（记为d），之后是一个类型为std::array<int, d>的参数包，参数包长度为nnz，分别表示每个非零元的位置。类的成员变量主要是一个非常量数组value，长度为nnz，类型根据需要可为浮点数或其他类型，记录稀疏张量的元素值。value的第i个元素即为第i个非零元的值。非零元的位置允许重复，重复的非零元的数学值定义为多个值的和。

图3是本发明实施例中计算单位刚度矩阵的流程图，在本实施例的一个实施方式中，根据有限元方程和稀疏张量类模板计算单位刚度矩阵包括：

S31，采用有限元方程的基函数和网格参数推导得到单元刚度矩阵元素方程和缩并维度；

S32，对稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合；

在一个示例中，根据缩并维度对稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合，包括：判断稀疏张量类模板中的非零元是否有解析表达式；若稀疏张量类模板中的非零元有解析表达式，利用静态for循环遍历稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合；若稀疏张量类模板中的非零元没有解析表达式，显式写出稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合。

针对如上静态化的COO格式稀疏张量，可以定义张量缩并运算。如果缩并成标量，例如运算，则可以直接使用C++17中的折叠表达式（fold expression）, 将各个非零元的值乘以相应的数后求和。对于其他缩并的情况，则可以利用静态for循环遍历每个非零元进行相应操作。

在一些示例中，为便于静态化COO稀疏张量的定义，通过重载加法运算符实现两个静态化COO稀疏矩阵的静态加法，从而可以通过加法的方式定义和计算矩阵。

图4是本发明实施例中静态化COO稀疏张量的模板实例化的示意图，以二维矩阵为例，4个静态化COO稀疏矩阵为（0,0），（0,2），（1,1），（2,2）。

S33，基于缩并维度对稀疏张量集合执行张量缩并运算，得到缩并后的中间结果；

可选的，基于缩并维度对稀疏张量集合执行张量缩并运算，得到缩并后的中间结果，包括：判断缩并维度是否为标量；若缩并维度是为标量，采用C++语法中的折叠表达式对稀疏张量集合各个非零元的值乘以相应的数后求和，得到中间结果。

通过使用静态循环和C++17折叠表达式，能够以表达式的方式定义稀疏矩阵,而避免手动展开每个非零元。

S34，采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵。

在一个实例中，采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵，包括：确定以下中间结果：、/>、/>、/>，其中，/>是形状为Field*Field的稀疏张量，/>为形状Field * Dim * Field的稀疏张量，/>是形状为Dim * Field * Field的稀疏张量，/>为形状Dim * Field * Dim * Field的稀疏张量，Dim为所述问题域的维数，Field为所述有限元方程的数量；采用所述中间结果配置以下单元刚度矩阵元素方程/>：，得到单位刚度矩阵；其中，∫表示数值积分，K为网格单元，/>为基函数，i、j分别为网格单元K内的行和列。

在一个实施场景中，有限元方程如下：；

(公式4)，/>，其中u为速度,分为水平方向速度和竖直方向速度两个分量，p为压力。/>为边界值。/>为雷诺数。

将(公式4)式变换为变分形式，得到形式（公式2）式，其中这里的共同组成(2)式中的u,/>形式如下：

；

；

；

，其中，(;)表示向量的拼接，[;]表示矩阵的拼接。

对于单元刚度矩阵计算中需要用到的g_0,g_1, g_2, g_3, 即为静态化稀疏张量，具体表达如下：

g_0为2阶张量形状3*3，(0,0),(0,1),(1,0),(1,1)处的非零元分别为Re乘以u的梯度在该处的值。

g_1为3阶张量形状3*3*2。(0,2,0),(1,2,1)处的非零元为Re, (2,0,0),(2,1,1)处的非零元为-1.0; 对于 i=0,1,(i,0,0)和(i,1,0)处的非零元为Re*u[0], (i,0,1)和(i,1,1)处的非零元值为Re*u[1].其他位置值为0。

g_2为3阶张量形状3*2*3,所有元素为0。

g_3为4阶张量形状3*2*3*2, 对于i=0,1, d=0,1, (i,d,i,d)和(i,d,d,i)处的值为0.5, 其他位置为0。

其中，g_0是形状为3*3=9的张量，g_1为形状3 * 3* 2=18的张量，g_2是形状为3*2 * 3=18的张量，g_3为形状3 * 2 * 3 * 2=36的张量，而针对本问题，g_0 含有4个非零元，g_1含有12个非零元，g_2含有0个非零元, g_3含有6个非零元。相比采用稠密张量的方案，大幅减少了计算量。

通用计算平台Xeon(R) Gold 6240，8核心, 1000单元进行实验, 采用P1-P2混合元求解二维Navier-Stokes方程方腔流问题（问题域）, 单元刚度矩阵计算时间如采用稠密张量需要0.19s，如采用稀疏静态化张量，需要0.036s，加速约5倍。本实施例的可以正常运行并加速计算出正确结果。

采用本实施例的方案，实现了一种静态化稀疏张量方法，用于加速有限元方法中单元刚度矩阵的计算中。采用C++模板元编程实现，因此无需引入领域专用语言，使得有限元编程中可以直接调用现有C++库，增加了编程的灵活性；同时，避免了代码生成的过程和对其他第三方技术的依赖，可以轻松地部署到任意提供C++17以上标准编译器的计算平台。本实施例的方案能够在缩短有限元方法中计算单元刚度矩阵的时间，同时适用于多种计算平台，从而加速流体力学、固体力学、电磁学等问题中有限元方法仿真分析的过程，通过提高算力，加快单元刚度矩阵和总刚度矩阵的计算速度，从而提高仿真速度和仿真效率。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到根据上述实施例的方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对相关技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质（如ROM/RAM、磁碟、光盘）中，包括若干指令用以使得一台终端设备（可以是手机，计算机，服务器，或者网络设备等）执行本发明各个实施例所述的方法。

实施例2

在本实施例中还提供了一种基于刚度矩阵的流体仿真装置，用于实现上述实施例及优选实施方式，已经进行过说明的不再赘述。如以下所使用的，术语“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图5是根据本发明实施例的一种基于刚度矩阵的流体仿真装置的结构框图，如图5所示，该装置包括：

确定模块50，用于确定目标流体在问题域上的有限元方程，其中，所述问题域包括多个有限元，每个有限元对应所述问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；

配置模块52，用于配置静态化的坐标COO格式的稀疏张量类模板，其中，所述稀疏张量类模板的模板参数包括张量维数和非零元的位置；

计算模块54，用于根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵；

生成模块56，用于采用所述单位刚度矩阵在处理器中装填生成所述问题域的总刚度矩阵；

模拟模块58，用于采用所述总刚度矩阵模拟计算所述目标流体在几何形状中的流动特征参数。

可选的，所述配置模块包括：第一配置单元，用于配置正数类型的第一参数包；第一选择单元，用于将所述第一参数包中的一个第一参数选择为张量维数d，并将所述第一参数的所在位置之后的每d个数配置为一个非零元的位置，得到静态化的COO格式的第一稀疏张量类模板。

可选的，所述配置模块包括：第二配置单元，用于配置正数类型的第二参数包；第二选择单元，用于将所述第二参数包中的一个第二参数选择为张量维数d，并在所述第二参数的所在位置之后配置常量表达式std::array<int, d>，采用所述常量表达式配置每个非零元的位置，得到静态化的COO格式的第二稀疏张量类模板。

可选的，所述计算模块包括：求解单元，用于采用有限元方程的基函数和网格参数推导得到单元刚度矩阵元素方程和缩并维度；运算单元，用于对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合；生成单元，用于基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果；配置单元，用于采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵。

可选的，所述运算单元包括：判断子单元，用于判断所述稀疏张量类模板中的非零元是否有解析表达式；第一生成子单元，用于若所述稀疏张量类模板中的非零元有解析表达式，利用静态for循环遍历所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合；第二生成子单元，用于若所述稀疏张量类模板中的非零元没有解析表达式，显式写出所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合。

可选的，所述生成单元包括：判断子单元，用于判断所述缩并维度是否为标量；缩并子单元，用于若所述缩并维度是为标量，采用C++语法中的折叠表达式对所述稀疏张量集合各个非零元的值乘以相应的数后求和，得到中间结果。

可选的，所述配置单元包括：确定子单元，用于确定以下中间结果：、、/>、/>，其中，/>是形状为Field*Field的稀疏张量，/>为形状Field * Dim * Field的稀疏张量，/>是形状为Dim * Field * Field的稀疏张量，为形状Dim * Field * Dim * Field的稀疏张量，Dim为所述问题域的维数，Field为所述有限元方程的数量；配置子单元，用于采用所述中间结果配置以下单元刚度矩阵元素方程/>：/>，得到单位刚度矩阵；其中，∫表示数值积分，K为网格单元，/>为基函数，i、j分别为网格单元K内的行和列。

需要说明的是，上述各个模块是可以通过软件或硬件来实现的，对于后者，可以通过以下方式实现，但不限于此：上述模块均位于同一处理器中；或者，上述各个模块以任意组合的形式分别位于不同的处理器中。

实施例3

本发明的实施例还提供了一种存储介质，该存储介质中存储有计算机程序，其中，该计算机程序被设置为运行时执行上述任一项方法实施例中的步骤。

可选地，在本实施例中，上述存储介质可以被设置为存储用于执行以下步骤的计算机程序：

S1，确定搜索空间中的初始值，其中，所述初始值包括初始位置值和初始元素值，所述初始位置值由N个向量构成，每个向量长度为高维数组对应阶的维度大小，N为大于2的正整数；

S2，基于所述初始值查找所述搜索空间在单一维度的单维最大值，以及所述单维最大值的数组位置；

S3，将所述单维最大值和所述数组位置输出为所述搜索空间在当前迭代周期的当前最大元素。

可选地，在本实施例中，上述存储介质可以包括但不限于：U盘、只读存储器（Read-Only Memory，简称为ROM）、随机存取存储器（Random Access Memory，简称为RAM）、移动硬盘、磁碟或者光盘等各种可以存储计算机程序的介质。

本发明的实施例还提供了一种电子设备，包括存储器和处理器，该存储器中存储有计算机程序，该处理器被设置为运行计算机程序以执行上述任一项方法实施例中的步骤。

可选地，上述电子设备还可以包括传输设备以及输入输出设备，其中，该传输设备和上述处理器连接，该输入输出设备和上述处理器连接。

可选地，在本实施例中，上述处理器可以被设置为通过计算机程序执行以下步骤：

S1，确定搜索空间中的初始值，其中，所述初始值包括初始位置值和初始元素值，所述初始位置值由N个向量构成，每个向量长度为高维数组对应阶的维度大小，N为大于2的正整数；

S2，基于所述初始值查找所述搜索空间在单一维度的单维最大值，以及所述单维最大值的数组位置；

S3，将所述单维最大值和所述数组位置输出为所述搜索空间在当前迭代周期的当前最大元素。

可选地，本实施例中的具体示例可以参考上述实施例及可选实施方式中所描述的示例，本实施例在此不再赘述。

上述本申请实施例序号仅仅为了描述，不代表实施例的优劣。

在本申请的上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。

在本申请所提供的几个实施例中，应该理解到，所揭露的技术内容，可通过其它的方式实现。其中，以上所描述的装置实施例仅仅是示意性的，例如所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，单元或模块的间接耦合或通信连接，可以是电性或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本申请各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本申请的技术方案本质上或者说对相关技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备（可为个人计算机、服务器或者网络设备等）执行本申请各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、只读存储器（ROM，Read-Only Memory）、随机存取存储器（RAM，Random Access Memory）、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述仅是本申请的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本申请原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本申请的保护范围。

Claims (10)

1.一种基于刚度矩阵的流体仿真方法，其特征在于，包括：

确定目标流体在问题域上的有限元方程，其中，所述问题域包括多个有限元，每个有限元对应所述问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；

配置静态化的坐标COO格式的稀疏张量类模板，其中，所述稀疏张量类模板的模板参数包括张量维数和非零元的位置；

根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵；

采用所述单位刚度矩阵在处理器中装填生成所述问题域的总刚度矩阵；

采用所述总刚度矩阵模拟计算所述目标流体在几何形状中的流动特征参数。

2.根据权利要求1所述的方法，其特征在于，配置静态化的COO格式的稀疏张量类模板包括：

配置正数类型的第一参数包；

将所述第一参数包中的一个第一参数选择为张量维数d，并将所述第一参数的所在位置之后的每d个数配置为一个非零元的位置，得到静态化的COO格式的第一稀疏张量类模板。

3.根据权利要求1所述的方法，其特征在于，配置静态化的COO格式的稀疏张量类模板包括：

配置正数类型的第二参数包；

将所述第二参数包中的一个第二参数选择为张量维数d，并在所述第二参数的所在位置之后配置常量表达式std::array<int, d>，采用所述常量表达式配置每个非零元的位置，得到静态化的COO格式的第二稀疏张量类模板。

4.根据权利要求1所述的方法，其特征在于，根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵包括：

采用有限元方程的基函数和网格参数推导得到单元刚度矩阵元素方程和缩并维度；

对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合;

基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果；

采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵。

5.根据权利要求4所述的方法，其特征在于，对所述稀疏张量类模板中的非零元进行运算，生成静态化的稀疏张量集合，包括：

判断所述稀疏张量类模板中的非零元是否有解析表达式；

若所述稀疏张量类模板中的非零元有解析表达式，利用静态for循环遍历所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合；若所述稀疏张量类模板中的非零元没有解析表达式，显式写出所述稀疏张量类模板中的每个非零元，并依次执行预设运算操作，生成静态化的稀疏张量集合。

6.根据权利要求4所述的方法，其特征在于，基于所述缩并维度对所述稀疏张量集合执行张量缩并运算，得到缩并后的中间结果，包括：

判断所述缩并维度是否为标量；

若所述缩并维度是为标量，采用C++语法中的折叠表达式对所述稀疏张量集合各个非零元的值乘以相应的数后求和，得到中间结果。

7.根据权利要求4所述的方法，其特征在于，采用所述中间结果计算单元刚度矩阵元素，得到单位刚度矩阵，包括：

确定以下中间结果：、/>、/>、/>，其中，/>是形状为Field*Field的稀疏张量，/>为形状Field * Dim * Field的稀疏张量，/>是形状为Dim * Field * Field的稀疏张量，/>为形状Dim * Field * Dim * Field的稀疏张量，Dim为所述问题域的维数，Field为所述有限元方程的数量；

采用所述中间结果配置以下单元刚度矩阵元素方程：，得到单位刚度矩阵；

其中，∫表示数值积分，K为网格单元，为基函数，i、j分别为网格单元K内的行和列。

8.一种基于刚度矩阵的流体仿真装置，其特征在于，包括：

确定模块，用于确定目标流体在问题域上的有限元方程，其中，所述问题域包括多个有限元，每个有限元对应所述问题域的一个离散化子域，每个离散化子域对应一个单位刚度矩阵；

配置模块，用于配置静态化的坐标COO格式的稀疏张量类模板，其中，所述稀疏张量类模板的模板参数包括张量维数和非零元的位置；

计算模块，用于根据所述有限元方程和所述稀疏张量类模板计算单位刚度矩阵；

生成模块，用于采用所述单位刚度矩阵在处理器中装填生成所述问题域的总刚度矩阵；

模拟模块，用于采用所述总刚度矩阵模拟计算所述目标流体在几何形状中的流动特征参数。

9.一种存储介质，其特征在于，所述存储介质包括存储的程序，其中，所述程序运行时执行上述权利要求1至7中任一项所述的方法步骤。

10.一种电子设备，包括处理器、通信接口、存储器和通信总线，其中，处理器，通信接口，存储器通过通信总线完成相互间的通信；其中：

存储器，用于存放计算机程序；

处理器，用于通过运行存储器上所存放的程序来执行权利要求1至7中任一项所述的方法步骤。