CN103440163A - 使用gpu并行实现的基于pic模型的加速器仿真方法 - Google Patents

Abstract

本发明公开了使用GPU并行实现的基于PIC模型的加速器仿真方法，包括：将初始化信息从主机复制到计算节点的GPU中；根据初始化信息，确定粒子位置与网格的对应关系；根据粒子位置与网格的对应关系，计算每个网格内所有粒子在网格上的电荷密度权重，得到网格的电荷密度分布；根据网格的电荷密度分布计算网格的电势分布，并根据网格的电势分布计算网格电场分布；计算每个粒子在电场作用下的运动变化，并更新每个粒子的运动状态；以及用每个粒子的更新的运动状态代替初始化信息，迭代地执行以上步骤，直到粒子运动状态满足设计需求。根据本发明实施例，能够解决现有PIC模型加速器仿真算法运算速度低、成本高等技术问题。

Description

使用GPU并行实现的基于PIC模型的加速器仿真方法

技术领域

本发明涉及加速器物理模拟仿真技术领域，具体涉及一种使用图形处理单元(GPU)加速卡的计算机实现基于质点网格(PIC)模型下对加速器中粒子加速过程模拟仿真方法。

背景技术

加速器模拟是加速器物理中一个非常重要的组成部分，细致而科学的加速器模拟对加速器设计和加速器调试起着至关重要的作用。由于在加速器中粒子数量巨大，粒子间相互作用紧密，以及加速器组成结构复杂，都对加速器模拟技术提出了很高的要求。选取合适的模拟模型，采用高效的模拟算法，利用高性能的计算硬件，可以更加真实、细微和直观地模拟出真实加速器环境中粒子的运动形态。

目前，加速器模拟算法的实现以CPU为主，这些方法由于CPU计算能力不足导致计算规模不足，在模拟过程中，极大的限制了计算速度和模拟规模，在可以接受的机时内只能计算很小的空间尺寸和时间尺寸。而基于NVIDIA公司开发的CUDA架构的GPU并行编程为高性能运算提供一种方便、直接又经济的解决方案。CUDA架构的编程模式，依靠其高效的并行模式，方便的库函数，为加速器模拟提供了新的契机。

发明内容

本发明的目的在于，提供一种使用GPU加速卡并行实现基于PIC模型的加速器过程模拟的方法，以解决现有PIC模型加速器模拟算法运算速度低、成本高等技术问题。

根据本发明一方面，提出了一种使用图形处理单元(GPU)实现的基于质点网格(PIC)模型的加速器仿真方法，包括：

a.在主机中产生初始化信息，并将初始化信息从主机复制到计算节点的GPU中，GPU包括多个流处理器；

在GPU中，利用多个流处理器并行执行以下步骤：

b.根据初始化信息，确定粒子位置与网格的对应关系；

c.根据粒子位置与网格的对应关系，计算每个网格内所有粒子在网格上的电荷密度权重，得到网格的电荷密度分布；

d.根据网格的电荷密度分布计算网格的电势分布，并根据网格的电势分布计算网格电场分布；

e.计算每个粒子在电场作用下的运动变化，并更新每个粒子的运动状态；以及

f.用每个粒子的更新的运动状态代替初始化信息，迭代地执行步骤b到e，直到粒子运动状态满足设计需求。

在一个实施例中，在步骤b和e中按照一个流处理器对应一个粒子的方式进行GPU并行处理，并且在步骤c和d中按照一个流处理器对应一个网格的方式进行GPU并行处理。

在一个实施例中，初始化信息包括三维仿真空间划分所得的网格数目、粒子的数目、粒子的三维位置和速度。

在一个实施例中，步骤b包括：确定每个粒子的位置所在的网格的编号并存储在数组中；根据确定的编号对数组中的粒子位置进行排序，使在同一个网格内的所有粒子位置连续排列；在排序后的数组中获取每个网络内粒子的开始位置和结束位置。

在一个实施例中，在步骤b中采用包括多个并行线程的线程块，每个线程块处理预定数目的网格，并且线程块共享访问GPU中的共享内存。

在一个实施例中，在步骤d中，对电荷密度分布进行三维傅立叶变换，根据频域电荷密度分布计算网格的频域电势分布，并对频域电势分布进行三维傅立叶逆变换，以得到网格的电势分布。

在一个实施例中，步骤e包括：计算每个粒子在电场作用下的受力和加速度，并更新每个粒子的三维速度和位置。

在一个实施例中，步骤e还包括：更新所述数组中的粒子位置，对数组中所有粒子位置排序，并更新每个网格内粒子的开始位置和结束位置。

在一个实施例中，利用排序后的数组，线程块对连续排列的同一个网格内的所有粒子位置进行合并访问。

在一个实施例中，在步骤d中利用纹理内存绑定的方法，根据网格的电势分布计算网格的电场分布。

在一个实施例中，在步骤e中，改变线程块的大小，并利用改变后的线程块的大小计算每个粒子的运动变化。

本发明提供一种了使用GPU实现基于PIC模型的加速器模拟的方法，利用快速发展的GPU加速卡技术，采用例如NVIDIA公司的CUDA等并行计算架构来并行实现该算法，可以更为真实地模拟加速器中束团运动过程，其实现效率高，成本低，可以更为便捷，快速地完成加速器模拟过程。此外，在算法结构设计中结合GPU硬件，本发明设计出符合GPU并行模式的算法结构，通过对粒子并行和网格并行的把握，充分利用GPU多线程和内存共享的优点，最大限度实现对PIC模型效率的提高。本发明有效地利用了现有基于CUDA架构的CUFFT库(参见“CUDA FFT(CUFFT)Library Documentation”，NVIDIA，2012.)和THRUST库(参见“CUDA TOOLKITDocumentation”，NVIDIA2012.)，采用粒子排序算法，在高效实现加速器模拟的同时，充分提高了程序的可移植性和可扩展性。最后，本发明在实现过程中，利用共享内存提高数据运算速度、优化数据结构，利用符合GPU硬件架构的合并访问技术减少显存读取次数，合理利用纹理内存，提高关键内核函数的运算速度，很大程度提高了整体程序的运算效率。

附图说明

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。附图中：

图1是实施根据本发明实施例的使用GPU的基于PIC模型的加速器仿真方法的系统结构的示意图；以及

图2是根据本发明实施例的使用GPU并行实现的基于PIC模型的加速器仿真方法的示意流程图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

图1是可以用于实施根据本发明实施例的使用GPU的基于PIC模型的加速器仿真方法的系统结构的示意图。如图1所示，该系统包括主机10、和计算节点20。虽然图中未示出，但是可以理解，主机10可以包括用户接口，例如键盘、触摸板等，通过用户接口，用户可以输入仿真所需的初始化信息、条件等。在图1中，主机10和计算节点20之间经由例如缆线等直接连接。然而，本发明实施例也可以采用其他任何适当的连接方式。例如，主机10可以经由局域网或广域网等与计算节点20通信，这使得能够远程进行仿真实验。例如，计算节点中的GPU可以通过PCI-E插口直接与主机的主板相连，完成CPU内存和GPU显存间的数据通信，并由该但GPU完成模拟过程。在仿真过程中，数据一旦从内存中拷贝到GPU显存中，所有计算均在GPU显存中完成，CPU内存可以只负责数据从显存中提取回来和写入磁盘当中，或者进行屏幕打印。可选地，图1的系统还可以包括外部的存储设备(未示出)，可以与主机连接，可以存储例如各个计算节点的计算结果，以防止死机、断电等意外情况发生导致数据丢失。

在一个实施例中，计算节点20可以是有GPU加速卡的高性能集群。在一个实施例中，计算节点20都具有GF110核心以上的NVIDIA通用计算卡。计算节点20可以进行基于NVIDIA公司开发的CUDA架构的GPU并行编程。

在一个实施例中，主机10可以由例如包括中央处理单元(CPU)的通用计算机实现。

以上系统仅仅是本发明基本构思的一种实现。本领域技术人员可以理解，上述各个部件的功能可以进行再分配或组合，以形成其他的系统构架。此外，如果功能足够强大，上述各个部件的功能可以集成到单个计算机或工作站中。

下面结合图1的系统结构，参照图2来描述根据本发明实施例的使用GPU实现的基于PIC模型的加速器仿真方法。图2示出了该加速器仿真方法的示意流程图。

如图2所示，在步骤202，在主机10中产生初始化信息，并将初始化信息从主机复制到计算节点20的GPU中。具体地，初始化信息可以包括关于网格和粒子的信息，例如三维仿真空间划分所得的网格数目、粒子的数目、粒子的三维位置和速度。在一个实施例中，初始化信息从主机内存复制到计算节点中的GPU设备内存中，粒子初始信息可以包括粒子数目ns、粒子三维方向位置信息pos_x，pos_y，pos_z和三维方向速度信息v_z，v_y，v_z。网格初始信息可以包括将仿真空间划分为Nx*Ny*Nz个网格，Nx，Ny，Nz分别表示在x，y，z方向空间划分的网格数目。

在步骤204，在GPU中，根据初始化信息，确定粒子位置与网格的对应关系。计算节点的GPU对每个空间网格中的所有粒子信息进行处理，其中计算节点GPU的每个流处理器负责处理一个对应网格内所有的粒子。需要读取粒子的三维方向位置数组，即pos_x[]，pos_y[]，pos[]，将粒子位置权重与空间网格格点进行对应。结合以上示例，步骤204具体包括如下过程：

(1)首先在GPU中对粒子的三维方向位置pos_x，pos_y，pos_z做预处理，获得粒子所在的网格信息。开辟大小为粒子数目ns的cell_count_list[]数组，cell_count_list[]数组用来存放各个粒子所处在的空间网格的编号：

Pos_to_cell＜＜＜ns/256，256＞＞＞(cell_count_list，pos_x，pos_y,pos_z)；

在一个实施例中，进行cuda相关函数调用时，每一线程块处理256个粒子，共有ns/256个线程块。

(2)对步骤(1)中获得的粒子cell_count_list信息进行连续排序，使在同一个网格内的所有粒子在GPU设备端内存空间连续排列。这一过程，采用TRUST库中的sort函数并行完成。THRUST库函数实现排序可以有效利用GPU设备的流处理器并行计算，不需要对CPU与GPU端的数据进行交换，避免了数据传输过程的时间开销。

同时，在GPU中利用THRUST库中的sort_by_key函数并行初始化排序索引cell_count_index[]，即将数组的第i个值赋值为i，i大于等于0小于等于ns。cell_count_index的意义是可以在排序后找回排序前粒子分布的原始信息：

Thrust：：sort_by_key(cell_count_list,cell_count_list+ns，cell_count_index)。

(3)在步骤(2)中得到cell_count_list排序的信息后，为了方便对同一网格内的所有粒子统计，需要用search_start函数并行获取每个网格内粒子的开始位置cell_count_start[]与结束位置cell_count_end[]，此时GPU使用ns/256个由256个线程组成的线程块执行这一过程，可以表示为：

Search_start＜＜＜ns/256，256＞＞＞(cell_count_list，cell_count_start，cell_count_end)；

此外，在上述获取开始和结束位置开始前，使用cudaMemset函数重置cell_count_start[]和cell_count_end[]全部为0。

在Search_start的kernel中，利用GPU中共享内存(sharedmemory)shared_cell_count来提高计算效率，Shared_cell_count的大小定义为257个int型。

在GPU硬件架构中，共享内存(参见“NVIDIA CUDA CProgramming Guide”，NVIDIA，2012.)位于芯片上，所以共享内存要比本地和全局内存空间快得多。实际上，对于warp的所有线程，只要在线程之间没有任何库冲突，访问共享内存与访问寄存器一样快。所以，线程块中可以共享访问，访问速度也远远高于全局显存，可以极大地提高线程块内共享临时数据的访问速度，这里将每个线程中的cell_count_list对应到shared_cell_count的共享显存中，利用共享内存的优点，高效实现search_start操作。

Search_start的具体实现过程为：

每一个线程负责将第cell_count_list[blockIdx.x*256+threadIdx.x]的数值记录到shared_cell_count[threadIdx.x]上，而第0个线程还需同时负责将shared_cell_count[(blockIdx.x+1)*256]记录到shared_cell_count[256]的位置上。

此时，编号为blockIdx.x*256+threadIdx.x的流处理器会将shared_cell_count[threadIdx.x]与shared_cell_count[threadIdx.x+1]两个值进行比较，若其值不一样，则分别记录为：

cell_count_start[blockIdx.x*256+threadIdx.x]＝shared_cell_count[threadIdx.x]

和

cell_count_start[blockIdx.x*256+threadIdx.x-1]＝shared_cell_count[threadIdx.x]

(4)、在GPU中利用cell_count_index索引重新更新粒子的位置信息，仍以256个线程为一组线程块，共ns/256个并行线程，编号为I＝blockIdx.x*256+threadIdx.x的流处理器处理第i个粒子更新后的粒子三维位置和速度信息：

Pos_x_new[I]＝pos_x[cell_count_index[i]]；

…

v_x_new[I]＝pos_x[cell_count_index[i]]；

此时每个流处理器的内核函数中包含cell_count_end[i]-cell_count_start[i]次循环。在如下描述的计算网格的电荷密度分布数组rho[]的过程中，利用更新过的粒子位置信息依次将编号从cell_count_start[i]至cell_count_end[i]的粒子的电荷密度权重存储在电荷密度分布数组的相应元素rho[i]中。

在CUDA架构下，由GPU硬件架构设计，全局显存的读入过程，是根据线程块的大小批量读入至内核函数当中的，数据在内存中存储时，如果线程对应数据在内存中连续时，可以减少内存数据的载入次数，这个过程叫做合并访问。在这种算法设计下的，粒子在排序后，存储顺序由粒子的网格编号决定，同一个网格编号的粒子内存连续，在内存读取过程中，达成内存合并访问，可以有效的减少线程块block中的内存读取次数，极大的提高了内存访问速度，有效的提高了GPU并行权重的效率。

在步骤206，根据粒子位置与网格的对应关系，计算每个网格内所有粒子在网格上的电荷密度权重，得到网格的电荷密度分布。具体地，初始化电荷密度数组rho[]，利用cudaMemset函数将初值全部赋值为0。GPU的每个线程块包含256个线程，共分配(Nx*Ny*Nz/256)个线程块，每个流处理器I＝blockIdx.x*256+threadIdx.x处理编号为i的网格内所有粒子的电荷密度权重。

在步骤208，根据网格的电荷密度分布计算网格的电势分布，并根据网格的电势分布计算网格电场分布。在GPU中，根据空间网格的电荷密度分布，使用CUFFT库函数，利用GPU流处理器对应空间网格的并行方式并行地将电荷密度分布转换为频域电荷密度分布，通过频域电荷密度分布求解空间网格格点处的频域电势分布，再使用CUFFT库函数，将频域电势分布转换为网格点上的电势分布，进而并行求解出空间网格的电场分布。步骤208具体包括如下过程。

(1)在GPU设备端，利用cufft库函数，对电荷密度分布作三维傅立叶正变换，变换后，得到频域电荷密度分布。需注意，在变换之前，应将电荷密度分布格式由float*转换为cuffiComplex*，，使之成为可以使用CUFFT库函数的数据格式。

可以采用如下的操作来实现傅立叶变化，

cufftPlan3d(&plan，Nz，Ny，Ny,CUFFT_C2C)；

cufftExecC2C(plan,rho，rho_fft,,CUFFT_FORWARD)；

在这个变换过程中，Nz，Ny，Nx的顺序不得改变，它是由三维网格向一维网格转换过程的顺序决定的；变换过程中，输入数据为rho，输出数据为rho_fft；变换标识为CUFFT_FORWARD，也就是采用傅立叶正变换。

变换后，得到频域的电荷密度分布rho_fft[]；

(2)根据步骤(1)中获得的频域电荷密度，可在GPU设备中并行地求解频域电势分布，即：

Rho_to_phi＜＜＜(Nx*Ny*Nz/256)，256＞＞＞(phi_fit，rho_fft)；

其中，phi_fft为频域电势分布，rho_fft为频域电荷密度。在内核函数中，流处理I＝blockIdx.x*256+threadldx.x被用于求解第i个网格内频域电势，通过计算即可得到空间内所有网格格点上的频域电势分布。

该过程需要大量使用GPU中的临时寄存器，为了最高效地发挥GPU的处理能力，采用每个线程块中包含256个线程的线程块划分方式。

(3)根据步骤(2)得到的频域电势分布，继续利用CUFFT库函数，对频域电势分布做傅立叶逆变换：

cufftPlan3d(&plan，Nz，Ny，Ny，CUFFT_C2C)；

cufftExecC2C(plan,phi_fit，phi，CUFFT_INVERSE)；

在这个变换过程中，Nz，Ny，Nx的顺序不得改变，它是由三维网格向一维网格转换过程的顺序决定的；变换过程中，输入数据为phi_fft，输出数据为phi；变换标识为CUFFT_INVERSE，也就是采用傅立叶逆变换。

变换后，需要对网格点上的电势分布进行数据格式的转化，即将cufftComplex*转化为float*，转化后，得到电势分布。

(4)分配(Nx*Ny*Nz/256)个线程块，每块线程块由256个线程组成，将步骤(3)中得到的空间电势分布分别在x，y和z方向上求解空间网格电场分布。

在这个过程中，采用纹理内存(参见“NVIDIA CUDA CProgramming Guide”，NVIDIA，2012)绑定的方法实现。纹理内存空间具有高速缓存，所以纹理拾取仅在高速缓存缺失时，耗费从设备内存中的一个内存读取，否则它仅耗费从纹理高速缓存中的一个读取。纹理高速缓存针对2D空间局部性进行了优化，所以读取紧密相邻的纹理地址的同一WARP的线程将达到最佳性能。此外，它还设计用于流水化具有恒定延迟的拾取。纹理内存绑定，有更高的访问带宽，而且不受访问模式的约束，在更快寻址的同时，可以隐藏计算时间，有时候会改善应用程序之行随机访问数据的性能。结合上述示例，纹理内存绑定的方法可以包括如下步骤：

a)绑定电势分布为纹理内存：

cudaBindTexture(0,rt,phi,Nx*Ny*Nz)；

b)利用电势分布phi[]，求出空间中的电场分布Ex[]，Ey[]，Ez[]：

phi_to_Ex＜＜＜Nx*Ny*Nz/256,256＞＞＞(phi,Ex)；

…

phi_to_Ez＜＜＜Nx*Ny*Nz/256,256＞＞＞(phi,Ez)；

在计算过程中，一共有Nx*Ny*Nz个格点，流处理对应网格，分别在x、y和z方向计算每个格点处的电场分布，GPU内核划分为Nx*Ny*Nz/256个线程块，每个线程块对应256个线程。在三维方向计算电场分布的时候，由于电势phi[]是以一维数组方式存储，电势phi[]在y方向和z方向分布不连续，纹理内存的使用对y方向和z方向的电场分布Ey和Ez的计算会有较大提高。

c)解除纹理内存绑定：

cudaUnBindTexture(rt)；

在使用过纹理内存后，解除phi[]的纹理内存绑定。

在步骤210，计算每个粒子在电场作用下的运动变化，并更新每个粒子的运动状态。具体地，根据步骤208中所求的电场分布，如下并行计算每个粒子在所受电场力作用下的位置和速度变化，并更新粒子位置及速度变化：

Vel_＜＜＜ns/128，128＞＞＞(v_x,v_y,v_z,Ex,Ey,Ez，pos_x，pos_y,pos_z，v_x_new,v_y_new,v_z_new)；

Pos＜＜＜ns/256,256＞＞＞(、pos_x,pos_y,pos_z,v_x,v_y,v_z,pos_x_new,pos_y_new,pos_z_new)；

其中pos_x_new，pos_y_new，pos_z_new，v_x_new，v_y_new，v_z_new为步骤204中更新后得到的粒子位置及粒子速度，pos_x，pos_y，pos_z，v_x，v_y，v_z是此次计算后更新得到的粒子位置及速度。

在CUDA架构中，由于GPU硬件的寄存器数量有限，在内核函数中，如果临时变量过多，则需要减少线程块中的线程数目。若同一个线程块中的线程数目过大，则会造成GPU中临时寄存器数目不够，使得计算速度急速下降甚至计算出错。

在计算速度变化的内核函数中，需要利用GPU中较多的寄存器存储临时变量，需要控制线程块中的线程数目，所以每个线程块大小划分为128个线程，线程块个数为ns/128个。

在步骤212，判断仿真结果是否满足设计需求。如果不满足，则将步骤210中更新后的粒子位置及速度变化值带回到步骤204中继续计算，直到仿真结果满足设计需求。如果满足，则以上过程结束。

可以将GPU设备内存中的结果数据拷贝回主机内存，并释放相应的CPU和GPU端内存。

本发明实施例根据GPU多轻量计算核心的特点，对算法结构及模型结构进行合理的安排与设计，将GPU中流处理器与算法中的网格、粒子有效地结合在一起，使PIC模型更加适合GPU并行模式，同时，充分提高了GPU线程的使用效率，极大地减少了模拟过程的计算时间。

此外，本发明中，通过合理利用高速访问的GPU中的共享内存，优化数据结构，采用符合GPU硬件结构的合并访问，有效的减少显存读取次数，合理利用纹理内存，对关键内核函数提高30％左右的运算效率，很大限度的发挥GPU计算的优点。

本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (11)

1.一种使用图形处理单元GPU实现的基于质点网格PIC模型的加速器仿真方法，包括：

a.在主机中产生初始化信息，并将初始化信息从主机复制到计算节点的GPU中，GPU包括多个流处理器；

在GPU中利用多个流处理器并行执行以下步骤：

b.根据初始化信息，确定粒子位置与网格的对应关系；

c.根据粒子位置与网格的对应关系，计算每个网格内所有粒子在网格上的电荷密度权重，得到网格的电荷密度分布；

d.根据网格的电荷密度分布计算网格的电势分布，并根据网格的电势分布计算网格电场分布；

e.计算每个粒子在电场作用下的运动变化，并更新每个粒子的运动状态；以及

f.用每个粒子的更新的运动状态代替初始化信息，迭代地执行步骤b到e，直到粒子运动状态满足设计需求。

2.根据权利要求1所述的方法，其中在步骤b和e中按照一个流处理器对应一个粒子的方式进行GPU并行处理，并且在步骤c和d中按照一个流处理器对应一个网格的方式进行GPU并行处理。

3.根据权利要求1或2所述的方法，其中初始化信息包括三维仿真空间划分所得的网格数目、粒子的数目、粒子的三维位置和速度。

4.根据权利要求3所述的方法，其中步骤b包括：

确定每个粒子的位置所在的网格的编号并存储在数组中；

根据确定的编号对数组中的粒子位置进行排序，使在同一个网格内的所有粒子位置连续排列；

在排序后的数组中获取每个网络内粒子的开始位置和结束位置。

5.根据权利要求4所述的方法，其中在步骤b中采用包括多个并行线程的线程块，每个线程块处理预定数目的网格，并且线程块共享访问GPU中的共享内存。

6.根据权利要求1或2所述的方法，其中，在步骤d中，对电荷密度分布进行三维傅立叶变换，根据频域电荷密度分布计算网格的频域电势分布，并对频域电势分布进行三维傅立叶逆变换，以得到网格的电势分布。

7.根据权利要求1或2所述的方法，其中步骤e包括：计算每个粒子在电场作用下的受力和加速度，并更新每个粒子的三维速度和位置。

8.根据权利要求4中所述的方法，其中步骤e还包括：更新所述数组中的粒子位置，对数组中所有粒子位置排序，并更新每个网格内粒子的开始位置和结束位置。

9.根据权利要求4所述的方法，其中利用排序后的数组，线程块对连续排列的同一个网格内的所有粒子位置进行合并访问。

10.根据权利要求1或2所述的方法，其中在步骤d中利用纹理内存绑定的方法，根据网格的电势分布计算网格的电场分布。

11.根据权利要求1或2所述的方法，其中在步骤e中，改变线程块的大小，并利用改变后的线程块的大小计算每个粒子的运动变化。

Images