CN116258042B - 一种基于ddm的大规模传热异构并行仿真方法 - Google Patents
一种基于ddm的大规模传热异构并行仿真方法 Download PDFInfo
- Publication number
- CN116258042B CN116258042B CN202310047349.XA CN202310047349A CN116258042B CN 116258042 B CN116258042 B CN 116258042B CN 202310047349 A CN202310047349 A CN 202310047349A CN 116258042 B CN116258042 B CN 116258042B
- Authority
- CN
- China
- Prior art keywords
- calculation
- matrix
- scale
- temperature
- heat transfer
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 145
- 238000012546 transfer Methods 0.000 title claims abstract description 53
- 238000004088 simulation Methods 0.000 title claims abstract description 32
- 238000004364 calculation method Methods 0.000 claims abstract description 91
- 230000008569 process Effects 0.000 claims abstract description 67
- 238000005516 engineering process Methods 0.000 claims abstract description 16
- 238000004458 analytical method Methods 0.000 claims abstract description 8
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 81
- 238000005192 partition Methods 0.000 claims description 61
- 238000000638 solvent extraction Methods 0.000 claims description 17
- 239000000463 material Substances 0.000 claims description 16
- 230000006870 function Effects 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 10
- 238000006243 chemical reaction Methods 0.000 claims description 8
- 230000000903 blocking effect Effects 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000004907 flux Effects 0.000 claims description 3
- 238000012805 post-processing Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims 1
- 238000000354 decomposition reaction Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 238000004043 dyeing Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000005272 metallurgy Methods 0.000 description 1
- 238000004377 microelectronic Methods 0.000 description 1
- 238000007639 printing Methods 0.000 description 1
- 239000004753 textile Substances 0.000 description 1
- 238000002076 thermal analysis method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Computer Graphics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于DDM的大规模传热异构并行仿真方法，属于高性能计算仿真技术领域，所述的方法首先建立传热分析的几何模型，其次对要分析的大规模几何模型或者几何装配体模型进行有限元网格划分，然后利用区域分解算法(DomainDecompositionMethod，DDM)将网格计算模型划分成几块计算区域(包括模型、边界条件、物理参数)，利用MPI技术对每块计算区域实现单机或者集群并行计算，在多进程下(OpenMP)通过利用多线程进行循环拆段来加速计算[K]e，并在计算后处理上利用GPU众核优势计算节点温度热流等导出量，从而实现大规模传热异构并行仿真。
Description
技术领域
本发明属于高性能计算仿真技术领域，具体是一种基于DDM的大规模传热异构并行仿真方法。
背景技术
传热不仅是常见的自然现象，而且广泛存在于工程技术领域。在能源动力、化工制药、材料冶金、机械制造、电气电信、建筑工程、文通运输、航空抗天、芯片、纺织印染、农业林业、生物工程、环境保护和气象预报等部门中存在大量的热量传递问题，常常还起着关键作用。
目前，虽然传热学的理论已逐渐完善，但是由于传热问题处理的问题多为非线性的实际问题，一般通过线性分析法进行近似再进行计算，但是大规模的求解使得求解的速度和正确性会受到影响。随着电子计算机的软硬件技术的飞速发展以及数值仿真理论和算法的不断优化，数值仿真技术越来越受到工业界的重视。近年来随着航空航天、微电子技术，计算机技术等，数值仿真方法等的迅猛发展，在热分析方面计算求解的规模不断扩大，单台计算机的计算资源已经不能够满足大规模的计算要求。对于大规模的传热仿真计算有了更急迫的需求。
因此，建立一种基于DDM的大规模传热异构并行仿真方法，应用发展成熟的有限元数值计算理论和传热学理论，并以大规模计算为切入点，能够有效的计算求解大规模的热传导问题，不仅可以应用于航空航天等高端装备大区域传热分析，还可以在为芯片设计、微纳制造等高精端领域提供更加精确的求解，能够极好的促进材料和结构设计的热力学性能优化。
发明内容
本发明公开了一种基于DDM的大规模传热异构并行仿真方法，该方法结合将传热学理论和有限元法结合起来、同时基于DDM、MPI、OpenMP、GPU技术利用计算机卓越的数值计算能力来解决传热学的相关问题，可以快速实现大规模传热的仿真计算。
本发明解决上述技术问题所提供的技术方案是：一种基于DDM的大规模传热异构并行仿真方法，包括以下步骤：
S1、建立任意大规模传热分析对象的几何模型或者几何装配体模型；
S2、对要分析的大规模模型进行空间离散化，即网格划分，生成计算需要的网格数据；
S3、对划分的网格模型进行边界条件、约束、材料、导热系数等参数设定；
S4、使用分区算法自动根据硬件情况确定分块数和分块方法，将求解网格模型划分成多个计算的子区域同时使得每个区域携带自己的基本信息参数；
S5、使用CPU多进程并行计算，对每个区域进行单独的迭代计算，实现对计算负载均衡划分和计算资源的充分利用；
S6、在每个进程下采用CPU多线程技术（OpenMP）计算单元传导矩阵，从而充分利用每个处理机的计算资源，并使用CRS方式存储稀疏矩阵和求解迭代温度场；
S7、各分区完成计算，进程0收集温度场的计算结果；
S8、通过GPU计算单元的温度梯度、热流梯度，快速实现计算的后处理；
S9、通过GPU计算节点的平均热流等从而实现大规模传热异构并行仿真；
进一步的技术方案是，所述步骤S1具体实现方法为：建立宏观尺度大规模的几何模型或者大规模的几何装配体模型，从而建立一个连续的求解域；
进一步的技术方案是，所述步骤S2具体实现方法为：
S21、离散参数的设定；
结合几何连续模型的几何特征，设定几何模型的全局网格种子、部分边的局部种子、部分区域的几何物理组，在部分变化较为剧烈的区域/>设置局部细化参数；
S22、对传热现象进行数学建模，会产生偏微分方程（PDE），数值分析即有限元方法计算使得PDE能够被近似求解原来模型域的数值解，其离散方程
其中n是离散化的自由度数，F是线性的，可以使用线性代数的方法求解u；
S23、建立大规模几何模型的离散模型；
通过读取网格划分参数，使用LiToMesh网格划分工具开展几何模型的网格划分，实现对不同形状的几何模型的空间离散，将模型离散为以四面体为单元的空间连续离散体，同时将模型的离散数据写出到本地文件中；
进一步的技术方案是，所述步骤S3具体实现方法为：
S31、设置材料参数，设置材料的导热系数λ、材料泊松比ν、材料弹性模量E、热膨胀系数CTE等；
S32、施加边界条件，指定约束和温度区域载荷来指定求解模型的边界；
进一步的技术方案是，所述步骤S4具体实现方法为：
S41、根据大规模模型的大小，同时结合要开启的硬件的进程的数量N，决定划分区域的数量，一般将大规模区域划分成N个区域，同时当N小于8时使用多级递归平分法，当N大于等于8时使用多级K-way 划分方法，将有限元的单元非结构化图离散区域换分成N块（N视模型和计算规模而定）区域；
其中线性传热系统被划分为N个集，使用如下计算矩阵形式：
式中为传导矩阵，/>为每个块的解集，/>为每个块的温度载荷与约束集；
S42、分区完成后将各个分区的节点及单元重新编号，并保存各分区重新编号前后的节点信息及单元信息，完成计算模型区域划分工作；
进一步的技术方案是，所述步骤S5具体实现方法为：
S51、在程序运行开始运行前，通过指定运行该SIMD程序被执行的进程数来决定进程的数量，执行如下的指令
其中mpiexec为MPI进程启动指令，-host为指定运行节点主机情况，-n为指定执行的进程数命令后面跟着运行的进程数，task为被执行并行程序的名字；
S52、在每个MPI进程中获取分区传导矩阵、分区边界条件/>；
S53、在每个并行MPI进程中迭代计算各划分区域的温度场并判断时候满足残差要求，具体计算过程如下；
计算温度场，迭代温度场的计算公式：
其中Ri为N到Ni的限制算子，转置算子Ri T是从Ni到N的扩展算子，；
以两个分区的具体计算形式：
其中，代表第n选代步的第1分区/>非重叠部分，/>代表第n选代步的第1分区重叠部分，/>代表第n选代步的第2分区/>重叠部分，/>代表第n选代步的第2分区/>非重叠部分，/>代表第n选代步的第1分区残差/>非重叠部分，/>代表第n迭代步的第2分区残差/>非重叠部分，/>代表第n选代步的第1分区残差与第2分区残差重叠部分，所以计算公式为：
两个MPI进程求解各自接收到的线性方程组，得到结果后按照重合部分乘以二分之一并相加，这种形式适用于所有分区数情况，当分区数大于 2 的区域划分情况的可以参照式(1.6)进行改写，只需要把重叠的节点所在分区情况体现在常数项矩阵即可，即可完成多分区的多MPI进程的分布式并行策略；
进一步的技术方案是，所述步骤S6具体实现方法为：
S61、首先计算单元的形函数，具体计算步骤如下；
使用的线性四面体单元网格，单元温度分布计算公式为：
Ni~l为单元形函数，i，j，k，l是每个单元的节点编号；
然后计算梯度矩阵；
其中[B] 是结构力学中的应变矩阵， {T} 是温度矢量；可以得到导热系数矩阵变为：
S62、在每个MPI进程并行计算区域内使用OpenMP多线程计算单元传导矩阵[K]e；计算方法为将 [D] 和 [B] 矩阵代入如下公式，即可得到单元传导矩阵[K]e；
S63、单元的方程由以下方程导出，从而建立单元的传热方程组；
S64、由局部温度和全局温度分量关系可以得到分区区域的整体的传导方矩阵，局部温度和全局温度分量的转换计算公式为：
其中为局部坐标系下的温度节点温度值，/>为分区下全局坐标系下的温度节点温度值；
S65、通过上面的转换矩阵得到单元在整体坐标下的传导矩阵，计算公式为：
其中，T为转换矩阵，为分区下一个单元的全局传导矩阵，/>为分区下一个单元的局部传导矩阵；
S66、组装区域的整体传导矩阵，上面得到区域全局下的传导矩阵，然后保存其非0元素的位置索引及其值，使用稀疏矩阵的COO存储方式进行稀疏矩阵的存储；
进一步的技术方案是，所述步骤S7具体实现方法为：
S71、实现进程的同步；各MPI进程计算各分区温度场，由于每块区域的计算耗时存在差异部分，所以需要进程同步，使用MPI_Barrier()使得先抵达执行这个接口的进程进入阻塞状态，以等待其他区域温度场计算的计算；
S72、温度场结果的收集；程序阻塞直到所有进程都开始执行MPI_Barrier()接口，当完成进程同步后，使用MPI_Gather()接口实现对其他进程数据的收集工作，将结果数据放在0号根进程中；
进一步的技术方案是，所述步骤S8具体实现方法为：
S81、计算求解温度梯度，整个求解域的温度场计算完成后温度场值，将温度场数据拷贝到CUDA显存中，计算好GPU的线程数划分线程块及线程网格，再使用kernal<<<block,grid>>>(double* t,double N)配置并调用CUDA的核函数，并在核函数内将数据分块使用众核线程的分块运算，得到温度梯度；计算公式如下，
其中Ti~j为单元节点温度，中间矩阵为温度梯度算子；
S82、计算热流梯度，整个求解域的温度梯度计算完成后，接着使用GPU计算单元热流，热传导满足傅里叶定律，公式如下，
其中qx，qy，qz是 x ，y，z方向的热通量； k 是热导率，材料的一种固有特性，dT /dx，dT /dy，dT /dz是温度梯度；
根据热通量计算平均热流，平均热流矢量计算公式为：
式中：为平均热流矢量；
本发明的有益效果是：本发明基于LiToMesh完成几何模型的离散建立有限元模型，然后使用DDM技术将大规模离散模型（含1.6亿四面体网格）进行区域分解，将分解后的各区域使用MPI技术对每块计算区域实现单机或者集群并行计算，在多进程下通过利用OpenMP进行循环拆段的多线程技术加速计算[K]e，并在后处理上利用GPU技术通过节点位移计算应力应变等导出量。采用现代大规模计算技术提出一种基于区域分解并结合多进程、CPU多线程、GPU计算的大规模异构并行仿真计算方法，使用“分而治之”的思想能够实现单台计算系统由于计算资源限制无法完成的计算传热任务，或者对求解速度有极高要求的传热计算场景，从而实现大规模传热异构并行仿真。
附图说明
图1是本发明的一种基于DDM的大规模传热异构并行仿真方法的流程图；
图2是本发明的几何计算模型示意图；
图3是离散单元的示意图；
图4是几何模型离散后的有限元网格图；
图5是基于DDM的离散后的整体和区域图；
图6是整体网格和区域分解后的信息图；
图7是进程数据收集示意图；
图8是计算出的温度结果；
图9是三个方向的平均热流矢量。
实施方式
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明提供一种基于DDM的大规模传热异构并行仿真方法，包括以下步骤：
S1、建立宏观尺度的几何模型，建立一个连续的求解域；
具体的实现步骤为：导入通过其他三维建模软件建立求解的大规模几何模型，将其设置为传热问题的几何连续求解域，即建得到连续的宏观模型，见图2；
S2、对要分析的几何传热模型进行空间离散化，即网格划分，生成计算需要的大规模网格数据；
具体包含以下步骤：
S21、离散参数的设定；
结合几何连续模型的几何特征，再由计算求解的规模和精度要求设置网格划分的全局网格种子数量或者距离、设置部分边的局部种子、并设置部分区域的几何物理组，这样可以根据几何物理组设置单元的材料参数，在部分变化较为剧烈的区域/>设置局部细化参数，防止应力集中导致的区域计算出现不正确的情况出现，该部分数据可以通过计算程序的数据库直接输入或者通过界面获取用户的输入；
S22、对传热现象进行数学建模，会产生偏微分方程（PDE），数值分析即有限元方法计算使得PDE能够被近似求解原来模型域的数值解，其离散方程
其中n是离散化的自由度数，F是线性的，可以使用线性代数的方法求解u；
S23、建立大规模几何模型的离散模型；
通过读取网格划分参数，使用LiToMesh网格划分工具开展几何模型的网格划分，实现对不同形状的几何模型的空间离散，将模型离散为以四面体为单元的空间连续离散体，单元模型见图3，离散的模型见图4，同时将模型的离散数据写出到本地文件中，文件包含节点的全局坐标情况和单元的节点信息；
S3、对划分的网格模型进行边界条件、约束、材料、导热系数等参数设定；
具体包括如下步骤：
S31、设置材料参数，设置材料的导热系数λ、材料泊松比ν、材料弹性模量E、热膨胀系数CTE等；该部分数据可以通过计算程序的数据库直接输入或者通过界面获取用户的输入，为后面的计算提供数据；
S32、施加边界条件，指定约束和温度区域载荷来指定求解模型的边界；从操作界面或者指定上选取节点作为边界作用位置然后设置其对应的温度边界条件和位移边界条件；
S4、使用DDM方法将求解网格模型划分成多个计算的子区域同时使得每个区域携带自己的基本信息参数；
具体包括以下步骤：
S41、通过和自研算法并结合由Karypis Lab开发的Mites工具包，根据大规模模型的大小，包括物理大小和网格数量大小，同时结合要开启的硬件的进程的数量N，决定划分区域的数量，一般将大规模区域划分成N个区域，同时当N小于8时使用多级递归平分法，当N大于等于8时使用多级K-way 划分方法，将有限元的单元非结构化图离散区域换分成N块（N视模型和计算规模而定）区域，当将求解离散模型划分为8个子区域时，其效果见图5，划分前的整个模型和划分后各分区的文件大小及节点、单元信息见图6；
其中线性传热系统被划分为N个集，使用如下计算矩阵形式：
式中为传导矩阵，/>为每个块的解集，/>为每个块的温度载荷与约束集；
S42、分区完成后将各个分区的节点及单元重新编号，并保存各分区重新编号前后的节点信息及单元信息，完成计算模型区域划分工作；
存储各个单元的编号、分区后的重新编号、单元类型、单元包含的节点及各节点坐标、单元体积、单元旋转矩阵、单元弹性矩阵、单元几何矩阵及单元传导矩阵；保存包括模型尺寸、单元信息及节点信息等所有模型信息；
S5、引入MPI多进程并行计算，对每个区域进行单独的迭代计算；
具体包括以下步骤：
S51、在程序运行开始运行前，通过指定运行该SIMD程序被执行的进程数来决定进程的数量，执行如下的指令
其中mpiexec为MPI进程启动指令，-host为指定运行节点主机情况，-n为指定执行的进程数命令后面跟着运行的进程数，task为被执行并行程序的名字；
S52、在每个MPI进程中获取分区传导矩阵、分区边界条件/>；
S53、在每个并行MPI进程中迭代计算各划分区域的温度场并判断时候满足残差要求，具体计算过程如下；
计算温度场，迭代温度场的计算公式：
其中Ri为N到Ni的限制算子，转置算子Ri T是从Ni到N的扩展算子，残差的计算公式如下：
；
以两个分区的具体计算形式：
其中，代表第n选代步的第1分区/>非重叠部分，/>代表第n选代步的第1分区重叠部分，/>代表第n选代步的第2分区/>重叠部分，/>代表第n选代步的第2分区/>非重叠部分，/>代表第n选代步的第1分区残差/>非重叠部分，/>代表第n迭代步的第2分区残差/>非重叠部分，/>代表第n选代步的第1分区残差与第2分区残差重叠部分，所以计算公式为：
两个MPI进程求解各自接收到的线性方程组，得到结果后按照重合部分乘以二分之一并相加，这种形式适用于所有分区数情况，当分区数大于 2 的区域划分情况的可以参照式上式进行改写，只需要把重叠的节点所在分区情况体现在常数项矩阵即可，即可完成多分区的多MPI进程的分布式并行策略；
S6、引入OpenMP在每个进程下采用多线程技术计算单元传导矩阵，和求解迭代温度场，运用OpenMP 的相关指令设置线程数及获取当前最大线程数n ，将计算传导矩阵需要用到的相关变量都设置为动态数组，已经是动态数组的提升一个维度，并将数组第一维大小设置为n ，以此防止调用多线程访问变量时发生冲突，在计算程序的配置属性中启用OpenMP 支持，运用OpenMP 的parallel for 指令标识原代码中用来计算全部单元的传导矩阵的for 循环，以此并行化原来的计算；
具体包括以下步骤：
S61、首先计算单元的形函数，具体计算步骤如下；
使用的线性四面体单元网格，单元温度分布计算公式为：
Ni~l为单元形函数，i，j，k，l是每个单元的节点编号；
然后计算梯度矩阵；
其中[B] 是结构力学中的应变矩阵， {T} 是温度矢量；可以得到导热系数矩阵变为：
S62、在每个MPI进程并行计算区域内使用OpenMP多线程计算单元传导矩阵[K]e；计算方法为将 [D] 和 [B] 矩阵代入如下公式，即可得到单元传导矩阵[K]e；
S63、单元的方程由以下方程导出，从而建立单元的传热方程组；
S64、由局部温度和全局温度分量关系可以得到分区区域的整体的传导方矩阵，局部温度和全局温度分量的转换计算公式为：
其中为局部坐标系下的温度节点温度值，/>为分区下全局坐标系下的温度节点温度值；
S65、通过上面的转换矩阵得到单元在整体坐标下的传导矩阵，计算公式为：
其中，T为转换矩阵，为分区下一个单元的全局传导矩阵，/>为分区下一个单元的局部传导矩阵；
S66、组装区域的整体传导矩阵，上面得到区域全局下的传导矩阵，然后保存其非0元素的位置索引及其值，使用稀疏矩阵的CSR方式中的COO存储方式进行稀疏矩阵的存储；分区传导矩阵的集成还需考虑各分区边界面上的节点受到其它分区包含该类节点的单元的影响，故首先需找出各分区与其他分区的边界面，并提取其节点，搜寻一个区域的边界节点，在除该区域外的其它分区中搜寻包含这些节点的单元，按照包含该类节点的个数及顺序提取该类单元的传导矩阵中的相关元素及对应的节点自由度编码，将节点自由度编码转换为该区域的节点自由度编码后可与对应的单元传导矩阵元素组成三元数组列表，将列表压缩存储进第一步的初始传导矩阵后最终完成了分区传导矩阵的组装；
S7、各分区完成计算，进程0收集温度场的计算结果；
具体包括以下步骤：
S71、实现进程的同步；各MPI进程计算各分区温度场，由于每块区域的计算耗时存在差异部分，所以需要进程同步，使用MPI_Barrier()使得先抵达执行这个接口的进程进入阻塞状态，以等待其他区域温度场计算的计算；
S72、温度场结果的收集；程序阻塞直到所有进程都开始执行MPI_Barrier()接口，当完成进程同步后，使用MPI_Gather()接口实现对其他进程数据的收集工作，将结果数据放在0号根进程中，如图7所示p0收集其他进程的数据，计算的温度结果汇总后效果如图8所示；
S8、通过GPU计算单元的温度梯度、热流梯度；
具体包括以下步骤：
S81、计算求解温度梯度，整个求解域的温度场计算完成后温度场值，将温度场数据拷贝到CUDA显存中，设计好GPU的线程数划分线程块及线程网格即block与grid的值，不同的配置会使得CUDA计算的速度有区别，再使用kernal<<<block,grid>>>(double* t,double N)配置并调用CUDA的核函数，并在核函数内将数据分块使用众核线程的分块运算，得到温度梯度；计算公式如下，
其中Ti~j为单元节点温度，中间矩阵为温度梯度算子；
S82、计算热流梯度，整个求解域的温度梯度计算完成后，接着使用GPU计算单元热流，热传导满足傅里叶定律，公式如下，
其中qx，qy，qz是 x ，y，z方向的热通量； k 是热导率，材料的一种固有特性，dT /dx，dT /dy，dT /dz是温度梯度；
9、通过GPU计算节点的平均热流；根据热通量计算平均热流，平均热流矢量计算公式为：
式中：为平均热流矢量，计算机结果见图9；
本发明将传热学理论和有限元法结合起来、同时基于DDM、MPI、OpenMP、GPU技术利用计算机卓越的数值计算能力来解决传热学的相关问题，可以快速实现大规模传热的仿真计算。
以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，可利用上述揭示的技术内容做出变动或修饰为等同变化的等效实施例。但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
Claims (10)
1.一种基于DDM的大规模传热异构并行仿真方法，其特征在于，包括以下步骤：
S1、建立任意大规模传热分析对象的几何模型或者几何装配体模型；
S2、对要分析的大规模模型进行空间离散化，即网格划分，生成计算需要的网格数据；
S3、对划分的网格模型进行边界条件、约束、材料、导热系数参数设定；
S4、使用分区算法自动根据硬件情况确定分块数和分块方法，将求解网格模型划分成多个计算的子区域同时使得每个区域携带自己的基本信息参数；
S5、使用多进程并行计算技术，对每个区域进行单独的迭代计算，实现对计算负载均衡划分和计算资源的充分利用；
S6、在每个进程下采用CPU多线程技术OpenMP计算单元传导矩阵，从而充分利用每个处理机的计算资源，并使用CRS方式存储稀疏矩阵和求解迭代温度场；
S7、各分区完成计算，进程0收集温度场的计算结果；
S8、通过GPU计算单元的温度梯度、热流梯度，快速实现计算的后处理；
S9、通过GPU计算节点的平均热流从而实现大规模传热异构并行仿真。
2.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S1具体实现方法为：建立宏观尺度大规模的几何模型或者大规模的几何装配体模型，从而建立一个连续的求解域。
3.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S2具体实现方法为：
S21、离散参数的设定；
结合大规模的几何连续模型的几何特征，设定几何模型的全局网格种子、部分边的局部种子、部分区域的几何物理组，在部分变化较为剧烈的区域/>设置局部细化参数；
S22、对传热现象进行数学建模，会产生偏微分方程PDE，数值分析即有限元方法计算使得PDE能够被近似求解原来模型域的数值解，其离散方程
其中n是离散化的自由度数，F是线性的，使用线性代数的方法求解u；
S23、建立大规模几何模型的离散模型；
通过读取网格划分参数，使用网格划分工具开展大规模几何模型的网格划分，实现对不同形状的大规模几何模型的空间离散，将模型离散为以四面体为单元的空间连续离散体，同时将模型的离散数据写出到本地文件中。
4.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S3具体实现方法为：
S31、设置材料参数，设置材料的导热系数λ、材料泊松比ν、材料弹性模量E、热膨胀系数CTE；
S32、施加边界条件，指定约束和温度区域载荷来指定求解模型的边界。
5.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S4具体实现方法为：
S41、根据大规模模型的大小，同时结合要开启的硬件的进程的数量N，决定划分区域的数量，将大规模区域划分成N个区域，同时当N小于8时使用多级递归平分法，当N大于等于8时使用多级K-way 划分方法，将有限元的单元非结构化图离散区域换分成N块区域；
其中线性传热系统被划分为n个集，使用如下计算矩阵形式：
式中为传导矩阵，/>为每个块的解集，/>为每个块的温度载荷与约束集；
S42、分区完成后将各个分区的节点及单元重新编号，并保存各分区重新编号前后的节点信息及单元信息，完成计算大规模模型区域划分工作。
6.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S5具体实现方法为：
S51、在程序开始运行前，通过指定运行SIMD程序被执行的进程数来决定进程的数量，执行如下的指令
其中mpiexec为MPI进程启动指令，-host为指定运行节点主机情况，-n为指定执行的进程数命令后面跟着运行的进程数，task为被执行并行程序的名字；
S52、在每个MPI进程中获取分区传导矩阵、分区边界条件/>；
S53、在每个并行MPI进程中迭代计算各划分区域的温度场并判断时候满足残差要求，具体计算过程如下；
计算温度场，迭代温度场的计算公式：
其中Ri为N到Ni的限制算子，转置算子Ri T是从Ni到N的扩展算子，；
以两个分区的具体计算形式：
其中，代表第n选代步的第1分区/>非重叠部分，/>代表第n选代步的第1分区/>重叠部分，/>代表第n选代步的第2分区/>重叠部分，/>代表第n选代步的第2分区/>非重叠部分，/>代表第n选代步的第1分区残差/>非重叠部分，/>代表第n迭代步的第2分区残差/>非重叠部分，/>代表第n选代步的第1分区残差与第2分区残差重叠部分，所以计算公式为：
两个MPI进程求解各自接收到的线性方程组，得到结果后按照重合部分乘以二分之一并相加，这种形式适用于所有分区数情况，当分区数大于 2 的区域划分情况的可以参照上式进行改写，只需要把重叠的节点所在分区情况体现在常数项矩阵即可，即可完成多分区的多MPI进程的分布式并行策略。
7.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S6具体实现方法为：
S61、首先计算单元的形函数，具体计算步骤如下；
使用的线性四面体单元网格，单元温度分布计算公式为：
Ni~l为单元形函数，i，j，k，l是每个单元的节点编号；
然后计算梯度矩阵；
其中[B] 是结构力学中对应的应变矩阵，{T} 是温度矢量；可以得到导热系数矩阵变为：
S62、在每个MPI进程并行计算区域内使用OpenMP多线程计算单元传导矩阵[K]e；计算方法为将 [D] 和 [B] 矩阵代入如下公式，即可得到单元传导矩阵[K]e：
S63、单元/>的方程由以下方程导出，从而建立单元的传热方程组；
S64、由局部温度和全局温度分量关系可以得到分区区域的整体的传导方矩阵，局部温度和全局温度分量的转换计算公式为：
其中/>为局部坐标系下的温度节点温度值，/>为分区下全局坐标系下的温度节点温度值；
S65、通过上面的转换矩阵得到单元在整体坐标下的传导矩阵，计算公式为：其中，T为转换矩阵，/>为分区下一个单元的全局传导矩阵，/>为分区下一个单元的局部传导矩阵；
S66、组装区域的整体传导矩阵，上面得到区域全局下的传导矩阵，然后保存其非0元素的位置索引及其值，使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储。
8.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S7具体实现方法为：
S71、实现进程的同步；各MPI进程计算各分区温度场，由于每块区域的计算耗时存在差异部分，所以需要进程同步，使用MPI_Barrier()使得先抵达执行这个接口的进程进入阻塞状态，以等待其他区域温度场计算的计算；
S72、温度场结果的收集；程序阻塞直到所有进程都开始执行MPI_Barrier()接口，当完成进程同步后，使用MPI_Gather()接口实现对其他进程数据的收集工作，将结果数据放在0号根进程中。
9.根据权利要求1所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于，所述步骤S8具体实现方法为：
S81、计算求解温度梯度，整个求解域的温度场计算完成后温度场值，将温度场数据拷贝到CUDA显存中，计算好GPU的线程数划分线程块及线程网格，再使用kernal<<<block,grid>>>(double* t,double N)配置并调用CUDA的核函数，
并在核函数内将数据分块使用众核线程的分块运算，得到温度梯度；计算公式如下，
其中Ti~j为单元节点温度，中间矩阵为温度梯度算子；
S82、计算热流梯度，整个求解域的温度梯度计算完成后，接着使用GPU计算单元热流，热传导满足傅里叶定律，公式如下，
其中qx，qy，qz是 x ，y，z方向的热通量； k 是热导率，材料的一种固有特性，dT /dx，dT/dy，dT /dz是温度梯度。
10.根据权利要求9所述的一种基于DDM的大规模传热异构并行仿真方法，其特征在于根据热通量计算平均热流，平均热流矢量计算公式为：
式中：/>为平均热流矢量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310047349.XA CN116258042B (zh) | 2023-01-31 | 2023-01-31 | 一种基于ddm的大规模传热异构并行仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310047349.XA CN116258042B (zh) | 2023-01-31 | 2023-01-31 | 一种基于ddm的大规模传热异构并行仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116258042A CN116258042A (zh) | 2023-06-13 |
| CN116258042B true CN116258042B (zh) | 2023-11-17 |
Family
ID=86678683
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310047349.XA Active CN116258042B (zh) | 2023-01-31 | 2023-01-31 | 一种基于ddm的大规模传热异构并行仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116258042B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116737384B (zh) * | 2023-06-21 | 2024-02-27 | 上海玫克生储能科技有限公司 | 电化学模型仿真计算的加速方法、存储介质及电子设备 |
| CN116911146B (zh) * | 2023-09-14 | 2024-01-19 | 中南大学 | 三维重力场全息数值模拟及cpu-gpu加速方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2007201062A1 (en) * | 2006-03-15 | 2007-10-04 | Alver Pty Ltd | A heater |
| CN102033985A (zh) * | 2010-11-24 | 2011-04-27 | 南京理工大学 | 基于*-矩阵算法的高效时域电磁仿真方法 |
| WO2015095785A1 (en) * | 2013-12-19 | 2015-06-25 | University Of Louisville Research Foundation, Inc. | Multi-scale mesh modeling software products and controllers |
| CN109492317A (zh) * | 2018-11-20 | 2019-03-19 | 中冶赛迪工程技术股份有限公司 | 基于连铸机二维温度场仿真方法及监控模型的运行方法 |
| CN111859766A (zh) * | 2020-07-28 | 2020-10-30 | 深圳拳石科技发展有限公司 | 可变计算域的拉格朗日积分点有限元数值仿真系统及方法 |
| CN114117864A (zh) * | 2021-12-03 | 2022-03-01 | 厦门大学 | 自适应时间步长有限元法在电子器件热仿真中的应用方法 |
-
2023
- 2023-01-31 CN CN202310047349.XA patent/CN116258042B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2007201062A1 (en) * | 2006-03-15 | 2007-10-04 | Alver Pty Ltd | A heater |
| CN102033985A (zh) * | 2010-11-24 | 2011-04-27 | 南京理工大学 | 基于*-矩阵算法的高效时域电磁仿真方法 |
| WO2015095785A1 (en) * | 2013-12-19 | 2015-06-25 | University Of Louisville Research Foundation, Inc. | Multi-scale mesh modeling software products and controllers |
| CN109492317A (zh) * | 2018-11-20 | 2019-03-19 | 中冶赛迪工程技术股份有限公司 | 基于连铸机二维温度场仿真方法及监控模型的运行方法 |
| CN111859766A (zh) * | 2020-07-28 | 2020-10-30 | 深圳拳石科技发展有限公司 | 可变计算域的拉格朗日积分点有限元数值仿真系统及方法 |
| CN114117864A (zh) * | 2021-12-03 | 2022-03-01 | 厦门大学 | 自适应时间步长有限元法在电子器件热仿真中的应用方法 |
Non-Patent Citations (2)
| Title |
|---|
| 多物理场耦合软件GTEA开发及应用;明平剑;张文平;;计算机辅助工程(第06期);15-21 * |
| 高超声速飞行器并行仿真方法研究;孙学功;龚春叶;;系统仿真学报(05);32-42 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116258042A (zh) | 2023-06-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116258042B (zh) | 一种基于ddm的大规模传热异构并行仿真方法 | |
| Hoisie et al. | Performance and scalability analysis of teraflop-scale parallel architectures using multidimensional wavefront applications | |
| Slattery et al. | The Data Transfer Kit: A geometric rendezvous-based tool for multiphysics data transfer | |
| Karatarakis et al. | GPU-acceleration of stiffness matrix calculation and efficient initialization of EFG meshless methods | |
| CN104182571B (zh) | 基于Delaunay和GPU的Kriging插值方法 | |
| Liu et al. | JAUMIN: a programming framework for large-scale numerical simulation on unstructured meshes | |
| Du et al. | Model parallelism optimization for distributed inference via decoupled CNN structure | |
| Fu et al. | Auto-NBA: Efficient and effective search over the joint space of networks, bitwidths, and accelerators | |
| Kanov et al. | The Johns Hopkins turbulence databases: an open simulation laboratory for turbulence research | |
| CN104111871A (zh) | 一种用于执行电路仿真中动态负载平衡的方法及装置 | |
| CN112347638B (zh) | 一种基于对偶单元法的三维集成微系统电热耦合分析方法 | |
| CN105677995A (zh) | 一种基于全网格配点理论的模糊稳态热传导问题数值求解方法 | |
| Khimich et al. | Numerical study of the stability of composite materials on computers of hybrid architecture | |
| Davis et al. | Paradigmatic shifts for exascale supercomputing | |
| Deng et al. | CPU/GPU computing for an implicit multi-block compressible Navier-Stokes solver on heterogeneous platform | |
| González García et al. | Load balancing for parallel computations with the finite element method | |
| Chavarría-Miranda et al. | High-performance computing (HPC): Application & use in the power grid | |
| Oancea et al. | Developing a high performance software library with MPI and CUDA for matrix computations | |
| Gasparini et al. | A linear solver framework for flow and geomechanics reservoir simulation | |
| Ding et al. | An automatic performance model-based scheduling tool for coupled climate system models | |
| Khan et al. | Analyzing the Implementation of the Newton Raphson Based Power Flow Formulation in CPU+ GPU Computing Environment | |
| Caron et al. | Performance prediction and analysis of parallel out-of-core matrix factorization | |
| Chang et al. | Spark distributed system and GPU parallel computing | |
| Hossain et al. | A flexible-blocking based approach for performance tuning of matrix multiplication routines for large matrices with edge cases | |
| Saltz et al. | Programming tools and environments |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |