CN110275733A - 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 - Google Patents
基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 Download PDFInfo
- Publication number
- CN110275733A CN110275733A CN201910564493.4A CN201910564493A CN110275733A CN 110275733 A CN110275733 A CN 110275733A CN 201910564493 A CN201910564493 A CN 201910564493A CN 110275733 A CN110275733 A CN 110275733A
- Authority
- CN
- China
- Prior art keywords
- phonon
- gpu
- module
- temperature
- distribution
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 62
- 230000001133 acceleration Effects 0.000 title claims abstract description 20
- 238000005315 distribution function Methods 0.000 claims abstract description 6
- 238000002939 conjugate gradient method Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 46
- 230000005540 biological transmission Effects 0.000 claims description 7
- 230000006870 function Effects 0.000 claims description 6
- 239000006185 dispersion Substances 0.000 claims description 4
- 230000001617 migratory effect Effects 0.000 claims description 4
- 238000012805 post-processing Methods 0.000 claims description 4
- 238000004088 simulation Methods 0.000 claims description 4
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000001052 transient effect Effects 0.000 claims description 2
- 238000004134 energy conservation Methods 0.000 claims 1
- 239000000126 substance Substances 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 10
- 238000004886 process control Methods 0.000 abstract description 2
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 12
- 238000012360 testing method Methods 0.000 description 7
- 239000000047 product Substances 0.000 description 6
- 238000012804 iterative process Methods 0.000 description 5
- 238000013508 migration Methods 0.000 description 5
- 230000005012 migration Effects 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000004513 sizing Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 239000012466 permeate Substances 0.000 description 2
- 241000638935 Senecio crassissimus Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Classifications
-
- 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
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
- G06F9/38—Concurrent instruction execution, e.g. pipeline or look ahead
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种基于有限体积法求解声子玻尔兹曼方程的GPU并行加速方法，通过划分非结构网格，确定边界条件和计算参数并初始化能量密度分布后，从CPU内存向GPU显存传输每个网格单元之间的影响系数；然后计算声子散射项，并使用稳定双共轭梯度法(BiCGSTAB)求解线性方程组，对每个网格单元的能量密度分布进行更新并通过GPU对声子模式温度分布和平衡态分布函数进行更新，最后通过比较每个网格单元的能量密度分布更新前后的变化，并当满足收敛条件时停止计算并输出结果。本发明在GPU上并行计算求解过程中的主要迭代部分，CPU负责整个计算过程的数据读取、数据输出以及计算流程控制，从而显著提高了计算效率。
Description
技术领域
本发明涉及的是一种传热的数值模拟与高性能计算交叉领域的技术，具体是一种基于有限体积法求解声子玻尔兹曼方程的GPU并行加速方法。
背景技术
在宏观尺度下，使用傅立叶定律的扩散方程准确描述物体的导热现象。然而，随着当今技术的快速发展，使得研究能够到达的时间与空间尺度越来越小。宏观尺度的基于傅立叶定律的分析方法就不再适用了，声子玻尔兹曼方程(BTE)有效地模拟介观尺度下的导热问题。
发明内容
本发明针对现有技术存在的上述不足，提出一种基于有限体积法求解声子玻尔兹曼方程的GPU并行加速方法，在GPU上并行计算求解过程中的主要迭代部分，CPU负责整个计算过程的数据读取、数据输出以及计算流程控制，从而显著提高了计算效率。
本发明是通过以下技术方案实现的：
本发明涉及一种基于GPU的有限体积法求解玻尔兹曼偏微分方程的并行加速方法，包括以下步骤：
步骤1)划分非结构网格，确定边界条件和计算参数。
步骤2)初始化能量密度分布，求得所有声子模式及其所有方向上的系数矩阵。
步骤3)从CPU内存向GPU显存传输每个网格单元之间的影响系数。
步骤4)计算声子散射项，并使用稳定双共轭梯度法(BiCGSTAB)求解线性方程组，对每个网格单元的能量密度分布进行更新。
步骤5)通过GPU对声子模式温度分布和平衡态分布函数进行更新。
步骤6)通过比较每个网格单元的能量密度分布更新前后的变化，并当满足收敛条件时停止计算并输出结果，否则返回步骤4。
步骤2中所述的初始化是指：将声子玻尔兹曼方程进行角度的离散，声子色散曲线的离散以及计算域的离散：其中：eω，p为能量密度分布，为平衡态能量分布，p为声子模式，vg为声子群速度，为角度的方向，为网格单元法向量，τ为弛豫时间，ΔΩ为每个方向的控制角，A为计算域；因此对于每一个声子模式的每一个角方向及其每一个空间上的网格单元来说，以方向为例，声子玻尔兹曼方程的离散方程为： 其中e为能量密度分布，l为长度，V为控制体体积。方程离散之后对所有角度、声子模式下的能量密度分布赋予初值，所有能量密度初值被赋为零。
步骤2中所述的声子模式及其所有方向上的系数矩阵，利用声子玻尔兹曼方程模拟导热问题得到，具体为：在假定当前的平衡态分布函数后，该方程会离散为一个线性方程，而对每个角方向以及整体的区域来说，有线性方程组：该系数矩阵中Kij表示第j个网格单元对第i个网格单元的影响系数。
步骤4中所述的声子散射项是指：通过结合声子散射项和基于能量的声子玻尔兹曼方程，得到在弛豫时间近似下的瞬态项为零的基于能量的声子玻尔兹曼方程： eω，p为能量密度分布，为平衡态能量分布，vg为声子群速度，τω，p为弛豫时间，该等式右边项即为声子散射项。
步骤5中所述的更新是指:通过该能量密度获得给定频率和分支的声子的总的能量
Eω，p＝∫4πeω，pdΩ和热流量Ω为控制体；由此定义整体能量的平均度量值，与声子模式相关的温度Tω，p满足：Eω，p＝Cω，p(Tω，p-Tref)，其中：Cω，p是声子比热容，Tref为决定能量基准线的参考温度，Tω，p为声子模式温度。依此可得整体的平衡温度晶格温度TL满足：由于能量守恒，则散射项在所有声子模式的积包括0，则可知晶格温度TL与Tω，p满足：
其中：TL和T ω，p为相关晶格温度和声子模式的温度。通过上述公式完成对温度以及平衡态分布的更新。
本发明涉及一种实现上述方法的系统，包括：数据读取模块、初始化模块、主机端设备端数据传输模块、声子散射项计算模块、线性方程组求解模块、声子模式温度求解模块、晶格温度及平衡态分布求解模块以及后处理模块，其中：数据读取模块根据参数文件解析出网格、速度以及弛豫时间信息并输出至初始化模块，初始化模块与主机端设备端数据传输模块相连并传输稀疏矩阵信息，声子散射项计算模块与线性方程求解模块相连并传输声子散射项信息，线性方程组求解模块分别与声子模式温度求解模块、晶格温度及平衡态分布求解模块相连并传输能量密度分布信息，晶格温度及平衡态分布求解模块与后处理模块相连并传递温度及平衡密度分布信息。
技术效果
与现有技术相比，本发明采用CPU-GPU协作的异构计算模式，采用CUDA编程框架，实现了对采用有限体积法求解声子玻尔兹曼方程数值解法的并行加速。整个迭代过程的主要计算部分都在GPU上执行，一方面对热点计算部分进行加速，另一方面是使得整个迭代过程的数值更新都在GPU内进行，减少CPU与GPU之间的数据传输，进一步提高GPU计算效率。通过对网格大小为2.4万，离散的角度个数为256，band个数为1的算例进行测试，GPU迭代速度相较CPU的迭代计算有75倍的提高。
附图说明
图1为有限体积法求解声子玻尔兹曼方程基本流程图；
图2为实施例中GPU上BiCGSTAB算法求解流程；
图3为实施例中ELL稀疏矩阵存储格式图；
图4为实施例中GPU规约算法示意图；
图5为实施例中K80平台上针对不同大小算例的加速比图。
具体实施方式
如图1所示，本实施例包括以下步骤：
步骤1)使用COMSOL软件划分非结构网格，具体为：通过COMSOL对计算域计算模型进行非结构化网格划分，生成相应的网格文件inputmesh.dat。手动生成边界条件文件inputbc.dat以及计算参数文件inputdata.dat，通过读入所需数据确定求解网格单元的坐标、边界信息以及个网格单元之间的邻接关系，然后对环境变量、角度离散以及边界条件数据进行初始化。
所述的边界条件包括固定温度边界、漫反射性边界以及镜面反射型边界；求解计算域为二维矩形，因此边界条件文件inputbc.dat中，采用0～3来表示边界，每个边界数值可对应不同的边界条件；计算参数inputdata.dat文件中包含求解所需的迁移速度，比热容，弛豫时间以及误差系数。
步骤2)在对文件中读入参数进行计算前首先需要对计算过程中的声子模式和角度进行离散。
所述的离散包括：①角度的离散，对声子迁移方向的极角和方位角进行离散，确定每个网格单元对应的角度方向和该方向上的控制角的个数；②声子模式的离散，即非灰模型下对色散曲线谱的离散，该离散个数取决于输入文件计算参数中的声子模式个数。
本实施例中将声子迁移方向的极角和方位角共包括256个方向，求解过程中方向离散函数get_direction可获得离散角度的方向和每个方向的控制角。在求解过程中，每个声子模式和每个角方向都会对应一个线性方程组，并且该线性方程组的系数在整个迭代过程中是固定不变的，因此可根据不同声子模式和不同角方向，采用多维矩阵存储所需系数矩阵。
从网格文件inputmesh.dat中读取网格信息，从边界条件文件inputbc.dat中读取边界条件信息以及从计算参数文件inputdata.dat中读取声子迁移速度、比热容和弛豫时间参数。初始化能量密度分布为0，求得所有声子模式及其所有方向上的系数矩阵，具体为：通过读入的网格参数、迁移速度、弛豫时间数据、离散后的声子迁移方向的极角和方位角计算出如所示的每个声子模式及其每个角方向上的离散方程组的系数；对于非灰模型来说，每个声子模式及其角方向都会有一个线性方程组，该线性方程组的系数在整个迭代过程中固定不变且系数矩阵为稀疏矩阵，因此采用多维压缩稀疏矩阵存储所需系数矩阵。
步骤1和步骤2在整个计算过程中仅执行一次，且涉及到数据的读入，因此优选在CPU中执行，不进行GPU加速。
步骤3)在CPU内存与GPU显存之间传输每个网格单元之间的影响系数，具体为：将整个迭代过程需要的所有数据由CPU内存复制到GPU显存，主要为与线性方程组Ax＝B计算有关的参数，其中：A为影响系数矩阵，x为能量密度分布，B为声子散射项参数。
所述的将整个迭代过程需要的所有数据包括：迁移速度、弛豫时间、声子比热容、边界条件、网格单元离散角度、初始能量密度、网格单元信息、声子散射项中间值、参数系数中间值以及影响系数矩阵。
优选地，为保证数据结构读取和计算的一致性，在GPU上实现了同CPU类似的数据访存函数，即使用一维向量存储多维数据，并且使用多维数据坐标进行访问，使得在GPU计算过程中对数据的访存与CPU一致。
优选地，为提高GPU数据读取效率，将更新过程中固定不变的系数矩阵绑定到GPU纹理内存，以提高访存效率，进一步提升计算速度。
步骤4)中，声子散射项的计算需在GPU中执行。每一次迭代，每一个声子模式及角度都会生成一个声子散射项向量。每一个网格单元的散射项仅与上一次迭代结果有关，与其他网格单元之间无数据依赖，因此在GPU上，以网格单元为单位进行并行计算。每一个CUDA线程计算一个网格单元的散射值。为了兼顾计算效率，GPU调度以及GPU寄存器分配，每个CUDA块(block)分配256个线程。则总的CUDA网格(Grid)大小为(numberofcell+255)/256。这样设置是防止有网格单元未被分配。整个声子散射项的计算完全在GPU中进行，且所需数据均可由GPU显存获得，其更新的数据也均保存在GPU显存中，该方式较少了主机端(CPU)与设备端(GPU)的数据传输，提高了整体计算效率。整个声子散射项计算过程并不完全统一，不能充分利用GPU单指令多线程(single instruction multiple thread，SIMT)的特点；同时计算过程对数据的访问多是随机访问，不能充分利用GPU合并访存特性。但是声子散射项计算函数保证了在GPU显存中声子散射项的更新，使得整个迭代过程能够完全在GPU上实现，消除了CPU与GPU之间数据传输所需的时间，提升了整体的计算效率。
所述的声子散射项，使用稳定双共轭梯度法(BiCGSTAB)求解线性方程组，求解每个网格单元的能量密度分布，具体为：使用CPU-GPU进行稳定双共轭梯度法计算，其中CPU进行求解过程的流程控制，GPU进行稀疏矩阵向量乘、向量点积、向量数乘和向量加运算，通过迭代计算求得当前收敛的能量密度分布。流程如图2所示。其中关于稀疏矩阵向量乘，向量点击以及向量的乘加等均在GPU内核中计算，而计算过程中流程的判断，如判断迭代是否结束，判断中间标量结果是否为0，是否需要重新初始化并重新进行迭代，这些均由CPU进行。整个求解过程所需数据均保存在GPU显存中，计算过程主机端与设备端的数据传输仅为几个需要判断流程的标量值。
所述的稀疏矩阵向量乘(SpMV)中，每个CUDA线程计算矩阵的一行，每个线程的迭代次数为系数矩阵每行非零数值的个数的最大值m。由于列索引值不是连续的，因此对向量的访问是随机的。对每行所有值计算完成后，计算结果存入结果向量。对结果向量的访问是按照行索引进行，因此warp访问的是一块连续内存，访问过程为合并访问。
所述的向量数乘和向量加(AXPY)计算中，每个CUDA线程计算向量中的一个值，整个计算过程的访存都是合并访问，每个线程完全并行计算。
所述的向量点积，由于涉及不同block之间的数据同步，采用两阶段规约实现，具体包括：
1)每一个线程计算两个数值相乘之积，并将乘积放入共享内存中，之后进行块内规约。将规约过程的循环直接展开，每步计算需进行线程同步。当活动的线程数量小于32时，进行warp内的规约计算。在warp内的规约计算不需要线程间的同步，计算效率更高。块内规约完成后，将每个块的结果按照块号分别保存至临时数组中。
2)对上一阶段的中间结果进行规约，规约过程同第一阶段类似，首先将规约循环展开，每一步的规约计算都需要进行线程同步。当活动线程小于32时，就可进行warp内的规约计算，省略线程同步过程。
优选地，为了优化稀疏矩阵的存储空间，提高GPU上稀疏矩阵的计算效率，所述的稀疏矩阵采用ELLPACK(以下简称ELL)格式存储，该ELL格式使用两个和原始矩阵相同行数的矩阵来存储稀疏矩阵：一个矩阵存储的是列号，一个矩阵存储的是数值，行号用自身所在行表示，该格式的存储方式如图3所示。实际使用过程中，为了简化格式，将二维矩阵转化为一维向量。使用两个一维向量分别存储矩阵的值以及矩阵的列索引，按列方向进行存储。假定稀疏矩阵每行最多m个数值非零，矩阵大小为N*N，则ELL稀疏矩阵中数值向量和列索引向量大小均为N*m。GPU上涉及到稀疏矩阵的计算，使CUDA线程(thread)计算稀疏矩阵的一行，则对稀疏矩阵的数值和列索引访问来说，每一个warp中的线程都是访问连续的一个内存块，实现合并访问，提高访存效率。该方式需确定原矩阵每行的非零元素个数。在本方法中，每行非零元素个数最多为3。实际计算中，若有一行的非零元素个数小于3，在ELL稀疏矩阵对应索引以及数值位置分别置为0，可保证不对计算结果有影响。
优选地，为了避免双倍的数据加载，减少内核数量，节约内核调用的开销，将多个内核融合为一个较大内核。例如在求解线性方程组中，会有形如pi＝ri+βti，yi＝K-1pi等相邻计算。其中K为jacobi预处理矩阵，即为原矩阵的对角矩阵，实际计算过程为一维向量，并行过程与向量数乘向量加(AXPY)一致，因此将两步操作融合成一个GPU内核。这样避免了对pi的多余访问，提高了计算效率。通过将多个操作融合为一个内核函数，减少对数据的额外加载，提高计算效率。
步骤5)根据步骤4得到的每个网格单元的能量密度分布对温度分布和平衡态分布进行更新，具体包括：声子模式温度更新和平衡态分布更新，其中：声子模式温度更新是指：对每一个声子模式来说，将每个网格单元的所有角度进行积分计算得到该声子模式对应的声子模式温度；平衡态分布更新是指：对所有声子模式温度求和得到每个网格单元对应的晶格温度，通过晶格温度求得新的平衡态分布。
优选地，所述的对声子模式温度分布和平衡态分布函数的更新均通过GPU计算，每个CUDA线程(thread)计算一个网格单元的温度值、晶格温度以及平衡态分布的值，每个CUDA块(BLOCK)分配256个线程，总的CUDA网格(Grid)大小为(cell的数量+255)/256。
步骤6)通过比较每个网格单元的能量密度分布更新前后的变化，并当满足收敛条件时停止计算并输出结果，否则返回步骤4，其中变化的计算通过GPU实现，每个CUDA线程(thread)计算一个网格单元差值的平方和，之后通过规约算法求得所有差值的平方和；收敛条件的判断通过CPU执行。
所述的变化，具体为两次更新之间能量密度分布的误差值，即两次更新间能量密度分布差值的平方和的开方，具体采用两步规约的方式，即首先在单个块内的循环展开式的并行规约，每步规约计算都要进行线程同步，当活动线程小于32，即在一个warp内后，可直接进行规约，无需进行线程同步，然后将块内规约结果存入GPU临时向量中并对临时向量进行规约计算，最后将最终规约的标量结果输出至CPU。
所述的收敛条件具体为：两次更新之间能量密度分布的误差值小于网格单元个数与误差系数的乘积即为收敛。
采用以上方法在GPU型号为:NVIDIA Tesla K80，CPU型号为:E5 2680 v3的平台上进行测试。测试将对比采用CPU-GPU并行计算和仅采用CPU进行计算的加速效果。其中GPU计算采用的CUDA框架版本为9.2，CPU中进行线性方程求解采用的是PETSc中库函数，其中PETSc版本为3.10。
CPU与GPU参数信息如下表：
在以上平台上对三个测试算例进行计算，其中每个测试算例网格及网格数量不同，边界条件以及弛豫时间、迁移速度等参数均相同，其中声子模式个数为1，离散角度个数为256。测试算例1的网格大小为24912；测试算例2的网格大小为6282；测试算例3的网格大小为2892。分别实现了通过CPU进行计算求解以及通过GPU对迭代部分进行加速求解。如图5所示为GPU并行加速的加速效果，从图中看出，测试算例越大，加速比越高。这是因为算例越大，所需线程越多，就能充分利用GPU的数量庞大的计算核心，且能够使GPU通过线程快速的切换掩盖对数据的访存的耗时。
总之，本方法提供了基于GPU的有限体积法求解声子玻尔兹曼方程的并行加速方法，在GPU上实现迭代求解过程的并行计算，相较于CPU上串行执行，本方法在2.4万个网格单元的算例上获得了75倍的加速比。
上述具体实施可由本领域技术人员在不背离本发明原理和宗旨的前提下以不同的方式对其进行局部调整，本发明的保护范围以权利要求书为准且不由上述具体实施所限，在其范围内的各个实现方案均受本发明之约束。
Claims (11)
1.一种基于GPU的有限体积法求解玻尔兹曼偏微分方程的并行加速方法，其特征在于，包括以下步骤：
步骤1)划分非结构网格，确定边界条件和计算参数；
步骤2)初始化能量密度分布，求得所有声子模式及其所有方向上的系数矩阵；
步骤3)从CPU内存向GPU显存传输每个网格单元之间的影响系数；
步骤4)计算声子散射项，并使用稳定双共轭梯度法求解线性方程组，对每个网格单元的能量密度分布进行更新；
步骤5)通过GPU对声子模式温度分布和平衡态分布函数进行更新；
步骤6)通过比较每个网格单元的能量密度分布更新前后的变化，并当满足收敛条件时停止计算并输出结果，否则返回步骤4；
所述的边界条件包括固定温度边界、漫反射性边界以及镜面反射型边界；求解计算域为二维矩形；
所述的初始化，即离散，具体包括：①角度的离散，对声子迁移方向的极角和方位角进行离散，确定每个网格单元对应的角度方向和该方向上的控制角的个数；②声子模式的离散，即非灰模型下对色散曲线谱的离散，该离散个数取决于输入文件计算参数中的声子模式个数。
2.根据权利要求1所述的并行加速方法，其特征是，步骤2中所述的所述的初始化是指：将声子玻尔兹曼方程进行角度的离散，声子色散曲线的离散以及计算域的离散：其中：eω，p为能量密度分布，为平衡态能量分布，p为声子模式，vg为声子群速度，为角度的方向，为网格单元法向量，τ为弛豫时间，ΔΩ为每个方向的控制角，A为计算域；因此对于每一个声子模式的每一个角方向及其每一个空间上的网格单元来说，以方向为例，声子玻尔兹曼方程的离散方程为：其中e为能量密度分布，l为长度，V为控制体体积，方程离散之后对所有角度、声子模式下的能量密度分布赋予初值，所有能量密度初值被赋为零。
3.根据权利要求1所述的并行加速方法，其特征是，步骤2中所述的声子模式及其所有方向上的系数矩阵，利用声子玻尔兹曼方程模拟导热问题得到，具体为：在假定当前的平衡态分布函数后，该方程会离散为一个线性方程，而对每个角方向以及整体的区域来说，有线性方程组：该系数矩阵中Kij表示第j个网格单元对第i个网格单元的影响系数。
4.根据权利要求1所述的并行加速方法，其特征是，步骤4中所述的声子散射项是指：通过结合声子散射项和基于能量的声子玻尔兹曼方程，得到在弛豫时间近似下的瞬态项为零的基于能量的声子玻尔兹曼方程：eω，p为能量密度分布，为平衡态能量分布，vg为声子群速度，τω，p为弛豫时间，该等式右边项即为声子散射项。
5.根据权利要求1所述的并行加速方法，其特征是，步骤5中所述的更新是指：通过该能量密度获得给定频率和分支的声子的总的能量Eω，p＝∫4πeω，pdΩ和热流量 Ω为控制体；由此定义整体能量的平均度量值，与声子模式相关的温度Tω，p满足：Eω，p＝Cω，p(Tω，p-Tref)，其中：Cω，p是声子比热容，Tref为决定能量基准线的参考温度，Tω，p为声子模式温度，依此可得整体的平衡温度晶格温度TL满足：由于能量守恒，则散射项在所有声子模式的积包括0，则可知晶格温度TL与Tω，p满足：其中：TL和Tω，p为相关晶格温度和声子模式的温度，通过上述公式完成对温度以及平衡态分布的更新。
6.根据权利要求1所述的并行加速方法，其特征是，为保证数据结构读取和计算的一致性，在GPU上实现了同CPU类似的数据访存函数，即使用一维向量存储多维数据，并且使用多维数据坐标进行访问，使得在GPU计算过程中对数据的访存与CPU一致。
7.根据权利要求1所述的并行加速方法，其特征是，将更新过程中固定不变的系数矩阵绑定到GPU纹理内存，以提高访存效率，进一步提升计算速度。
8.根据权利要求1所述的并行加速方法，其特征是，步骤6中所述的变化，即具体为两次更新之间能量密度分布的误差值，具体为两次更新间能量密度分布差值的平方和的开方，其采用两步规约的方式，即首先在单个块内的循环展开式的并行规约，每步规约计算都要进行线程同步，然后将块内规约结果存入GPU临时向量中并对临时向量进行规约计算，最后将最终规约的标量结果输出至CPU。
9.根据上述任一权利要求所述的并行加速方法，其特征是，为了避免双倍的数据加载，减少内核数量，节约内核调用的开销，将多个GPU内核融合为一个较大内核，即将两步操作融合成一个GPU内核进行处理。
10.一种实现上述任一权利要求所述方法的系统，其特征在于，包括：数据读取模块、初始化模块、主机端设备端数据传输模块、声子散射项计算模块、线性方程组求解模块、声子模式温度求解模块、晶格温度及平衡态分布求解模块以及后处理模块，其中：数据读取模块根据参数文件解析出网格、速度以及弛豫时间信息并输出至初始化模块，初始化模块与主机端设备端数据传输模块相连并传输稀疏矩阵信息，声子散射项计算模块与线性方程求解模块相连并传输声子散射项信息，线性方程组求解模块分别与声子模式温度求解模块、晶格温度及平衡态分布求解模块相连并传输能量密度分布信息，晶格温度及平衡态分布求解模块与后处理模块相连并传递温度及平衡密度分布信息。
11.根据权利要求10所述的系统，其特征是，所述的稀疏矩阵信息采用ELL格式存储；在ELL稀疏矩阵对应索引以及数值位置分别置为0。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910564493.4A CN110275733B (zh) | 2019-06-27 | 2019-06-27 | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910564493.4A CN110275733B (zh) | 2019-06-27 | 2019-06-27 | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110275733A true CN110275733A (zh) | 2019-09-24 |
| CN110275733B CN110275733B (zh) | 2022-11-22 |
Family
ID=67963444
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910564493.4A Active CN110275733B (zh) | 2019-06-27 | 2019-06-27 | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110275733B (zh) |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111858066A (zh) * | 2020-07-30 | 2020-10-30 | 中国空气动力研究与发展中心超高速空气动力研究所 | 气体动理论统一算法中的cpu+gpu异构并行优化方法 |
| CN112257313A (zh) * | 2020-10-21 | 2021-01-22 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
| CN112380793A (zh) * | 2020-11-18 | 2021-02-19 | 上海交通大学 | 基于gpu的湍流燃烧数值模拟并行加速实现方法 |
| CN113987841A (zh) * | 2021-12-24 | 2022-01-28 | 苏州浪潮智能科技有限公司 | 一种求解界面处声子热输运的方法和存储介质 |
| CN114155135A (zh) * | 2021-12-02 | 2022-03-08 | 西北核技术研究所 | 基于gpu集群的相对论性波尔兹曼方程计算方法、存储介质及设备 |
| CN114201887A (zh) * | 2021-12-29 | 2022-03-18 | 苏州浪潮智能科技有限公司 | 声子热输运的稳态模拟方法和装置 |
| CN114329977A (zh) * | 2021-12-29 | 2022-04-12 | 苏州浪潮智能科技有限公司 | 声子热输运的瞬态模拟方法和装置 |
| CN115495919A (zh) * | 2022-09-30 | 2022-12-20 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN115618498A (zh) * | 2022-11-08 | 2023-01-17 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器跨流域流场的预测方法、装置、设备及介质 |
| CN116933553A (zh) * | 2023-08-02 | 2023-10-24 | 上海交通大学 | 数值反应堆中子学的非结构网格体积修正方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101760183A (zh) * | 2009-12-30 | 2010-06-30 | 哈尔滨工业大学 | 计算硅锗超晶格材料界面热阻的方法 |
| CN102681972A (zh) * | 2012-04-28 | 2012-09-19 | 浪潮电子信息产业股份有限公司 | 一种利用GPU加速格子-Boltzmann的方法 |
| CN102944574A (zh) * | 2012-11-14 | 2013-02-27 | 北京科技大学 | 一种介孔异质复合材料热物性的计算方法 |
| CN102945295A (zh) * | 2012-10-15 | 2013-02-27 | 浪潮(北京)电子信息产业有限公司 | 一种格子玻尔兹曼方法的并行加速方法及系统 |
-
2019
- 2019-06-27 CN CN201910564493.4A patent/CN110275733B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101760183A (zh) * | 2009-12-30 | 2010-06-30 | 哈尔滨工业大学 | 计算硅锗超晶格材料界面热阻的方法 |
| CN102681972A (zh) * | 2012-04-28 | 2012-09-19 | 浪潮电子信息产业股份有限公司 | 一种利用GPU加速格子-Boltzmann的方法 |
| CN102945295A (zh) * | 2012-10-15 | 2013-02-27 | 浪潮(北京)电子信息产业有限公司 | 一种格子玻尔兹曼方法的并行加速方法及系统 |
| CN102944574A (zh) * | 2012-11-14 | 2013-02-27 | 北京科技大学 | 一种介孔异质复合材料热物性的计算方法 |
Non-Patent Citations (2)
| Title |
|---|
| ZAHIRI ET AL.: "Numerical investigation of ballistic-diffusive heat transfer through a constriction with the Boltzmann transport equation", 《WEB OF SCIENCE》 * |
| 张立民等: "基于GPU的受限玻尔兹曼机并行加速", 《电子设计工程》 * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111858066B (zh) * | 2020-07-30 | 2022-07-15 | 中国空气动力研究与发展中心超高速空气动力研究所 | 气体动理论统一算法中的cpu+gpu异构并行优化方法 |
| CN111858066A (zh) * | 2020-07-30 | 2020-10-30 | 中国空气动力研究与发展中心超高速空气动力研究所 | 气体动理论统一算法中的cpu+gpu异构并行优化方法 |
| CN112257313A (zh) * | 2020-10-21 | 2021-01-22 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
| CN112257313B (zh) * | 2020-10-21 | 2024-05-14 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
| CN112380793A (zh) * | 2020-11-18 | 2021-02-19 | 上海交通大学 | 基于gpu的湍流燃烧数值模拟并行加速实现方法 |
| CN112380793B (zh) * | 2020-11-18 | 2024-02-13 | 上海交通大学 | 基于gpu的湍流燃烧数值模拟并行加速实现方法 |
| CN114155135A (zh) * | 2021-12-02 | 2022-03-08 | 西北核技术研究所 | 基于gpu集群的相对论性波尔兹曼方程计算方法、存储介质及设备 |
| CN113987841B (zh) * | 2021-12-24 | 2022-04-19 | 苏州浪潮智能科技有限公司 | 一种求解界面处声子热输运的方法和存储介质 |
| CN113987841A (zh) * | 2021-12-24 | 2022-01-28 | 苏州浪潮智能科技有限公司 | 一种求解界面处声子热输运的方法和存储介质 |
| CN114329977A (zh) * | 2021-12-29 | 2022-04-12 | 苏州浪潮智能科技有限公司 | 声子热输运的瞬态模拟方法和装置 |
| CN114329977B (zh) * | 2021-12-29 | 2024-01-23 | 苏州浪潮智能科技有限公司 | 声子热输运的瞬态模拟方法和装置 |
| CN114201887B (zh) * | 2021-12-29 | 2024-01-23 | 苏州浪潮智能科技有限公司 | 声子热输运的稳态模拟方法和装置 |
| CN114201887A (zh) * | 2021-12-29 | 2022-03-18 | 苏州浪潮智能科技有限公司 | 声子热输运的稳态模拟方法和装置 |
| CN115495919A (zh) * | 2022-09-30 | 2022-12-20 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN115495919B (zh) * | 2022-09-30 | 2023-05-26 | 中国船舶科学研究中心 | 一种基于格子玻尔兹曼的时域缓坡方程的数值求解方法 |
| CN115618498A (zh) * | 2022-11-08 | 2023-01-17 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器跨流域流场的预测方法、装置、设备及介质 |
| CN116933553A (zh) * | 2023-08-02 | 2023-10-24 | 上海交通大学 | 数值反应堆中子学的非结构网格体积修正方法 |
| CN116933553B (zh) * | 2023-08-02 | 2024-02-13 | 上海交通大学 | 数值反应堆中子学的非结构网格体积修正方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110275733B (zh) | 2022-11-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110275733A (zh) | 基于有限体积法求解声子玻尔兹曼方程的gpu并行加速方法 | |
| Li et al. | GPU-accelerated preconditioned iterative linear solvers | |
| Komatitsch et al. | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA | |
| Ingram et al. | Glimmer: Multilevel MDS on the GPU | |
| Yamazaki et al. | Improving the performance of CA-GMRES on multicores with multiple GPUs | |
| Rustico et al. | Advances in multi-GPU smoothed particle hydrodynamics simulations | |
| Gaburov et al. | Gravitational tree-code on graphics processing units: implementation in CUDA | |
| Dimond et al. | Accelerating large-scale HPC Applications using FPGAs | |
| US20130194016A1 (en) | System and method for generating a clock gating network for logic circuits | |
| Owens et al. | Discontinuous isogeometric analysis methods for the first-order form of the neutron transport equation with discrete ordinate (SN) angular discretisation | |
| JP2011514587A (ja) | 原子炉炉心をモデル化するためのコンピュータで実行する方法及び対応するコンピュータプログラム製品 | |
| CN110516316A (zh) | 一种间断伽辽金法求解欧拉方程的gpu加速方法 | |
| US20150160371A1 (en) | Gpu accelerated deflation in geomechanics simulator | |
| Dimarco et al. | Towards an ultra efficient kinetic scheme. Part III: High-performance-computing | |
| Bernaschi et al. | A factored sparse approximate inverse preconditioned conjugate gradient solver on graphics processing units | |
| Hassanein et al. | Parallel hardware implementation of the brain storm optimization algorithm using FPGAs | |
| Emelyanov et al. | Analysis of impact of general-purpose graphics processor units in supersonic flow modeling | |
| Chandar et al. | A GPU-based incompressible Navier–Stokes solver on moving overset grids | |
| Bramas | Optimization and parallelization of the boundary element method for the wave equation in time domain | |
| Michels | Sparse-matrix-CG-solver in CUDA | |
| Vasileska et al. | Particle-in-cell code for gpu systems | |
| CN115758492A (zh) | 一种基于三维cad的fdtd共形网格自动生成方法及装置 | |
| Naher et al. | Using Machine Learning to Estimate Utilization and Throughput for OpenCL-Based SpMV Implementation on an FPGA | |
| Cwik et al. | The application of scalable distributed memory computers to the finite element modeling of electromagnetic scattering | |
| Lindahl | MULTI-DIMENSIONAL RESPONSE MATRIX METHOD. |
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 |