计算物理已经成为物理学的一个重要分支,经常需要对海量的数据进行复杂的计算。在普通的计算机上研究人员很难通过简单的计算公式对复杂运动进行建模,而精确地模拟却由于巨大的计算量而无法实现。计算机业界正处于并行计算的革命中,强大并且低成本的GPU集群计算能力使得研究者们能够进行快速的实验。

  CUDA介绍

  CUDA(Compute Unified Device Architecture,统一计算设备架构)是一种由NVIDIA公司推出的处理和管理GPU并行计算的硬件和软件的架构。它绕过了图形流水线,直接对GPU的硬件核心做了一层多线程封装,根据其提供的多线程并行编程接口可以很有效地进行多线程编程,开发线程级并行性。这种新型的GPU编程模型语法是对C语言的一个极小扩展,添加了一些使用GPU 硬件特性的专有指令,使得GPU能够更有效地用于通用计算。

  SRIM介绍

  SRIM(the stopping and range of ions into matter)系列软件基于Monte Carlo方法,是目前国际上通用的计算带电粒子在靶材料中入射和能量传递的计算机程序。其原名为TRIM(the Transport of Ions in Matter)1985年由J. F. Ziegler J. P. Biersack和U.littmark编写,主要计算入射粒子在物质中的能损,射程等。以后在此基础上不断升级完善,并于1991年更名为SRIM。

  SRIM对入射原子在靶中的运动分析采用连续慢化假设:

  1.入射原子与材料靶原子核的碰撞采用两体弹性碰撞来描述, 这一部分主要导致入射原子运动轨迹的偏折, 能量损失来自弹性碰撞的能量损失部分。

  2.两次碰撞之间的距离以及碰撞后的参数通过随机抽样得到。

  3.在两次碰撞之间,认为入射原子与材料中的电子作用,连续均匀地损失能量。

  粒子的位置、能量损失以及次级粒子的各种参数都在整个跟踪过程中存储下来,最后得到各种所需物理量的期望值和相应的统计误差。

  SRIM并行化中的一些技术

  并行的基本思想一般有两种:功能并行和数据并行。功能并行是不相关的任务对数据集的不同元素进行不同的操作,而数据并行是不相关的任务对数据集的不同元素进行相同的操作。为了实现入射粒子在物质中运动的高性能模拟,根据SRIM的特点我们对其进行了数据并行。下面以SRIM为物理背景主要介绍采用CUDA并行化时遇到的一些技术问题。  

 共享数据的并发修改

  因为SRIM主要模拟跟踪一大批入射粒子的运动,每个粒子在靶物质中的运动是相互独立的,所以把每个粒子在靶材料中的作用作为线程并行化,在统计能损过程会遇到如此问题:

  MTOT[i]=MTOT[i]+E_loss,若把靶在深度方向上的长度均匀划分成N段,则i为第i段靶材料,MTOT[i]用来存储所有粒子在第i段上的能损总量,E_loss是每个粒子在第i段上的能损。此类问题即为浮点运算共享数据的并发修改问题,是并行开发中最常见的问题。对于整数型变量的并发修改我们可利用CUDA提供的原子加法操作,但由于浮点数加法运算的中间结果可能被截断,并非是可结合的,所以在这里不能应用。我们也可以自己书写高级原子操作通过锁定数据结构来实现,但实践证明其在线程数比较多的情况下很耗时。还有一种可行的方法就是,对变量MTOT分配Yn*blockdim*threaddim大小,其中blockdim=(ns_p+threaddim-1)/threaddim即把每个粒子的对应在靶材料中的能损分别存储,然后再并行相加,这种方法在粒子数较小的情况下是可行的,但等到粒子数很大时改变量会占用很大内存。本文提出一种折中的方法,同样给MTOT分配Yn*blockdim*threaddim大小,其中blockdim,threaddim给于较小的定值,如二者都为16。

  每个线程的起始索引按照如下公式计算int tid=threadIdx.x+blockDim.x*blockIdx.x;

  当粒子数大于线程数时,在每个线程计算完当前索引任务后接着需要对索引进行递增,其中递增的步长为线程格中正在运行的线程数量。即i=i+blockDim.x*gridDim.x;

  那么代码框架可写为如下:

  …
  #define yn (100)
  …
  __global__ void mcc(参数){
  int tid=threadIdx.x+blockDim.x*blockIdx.x;
  int i=tid;
  …
  while(i<particleNum&&tid<blockDim.x*gridDim.x){
  …
  temp=(i%blockDim.x*gridDim.x)*yn+m-1;//m>=0&&m<yn
  MTOT[temp]=MTOT[temp]+E_loss;
  …
  i+= blockDim.x*gridDim.x;
  }
  }
  int main(){
  int blockdim=16;
  int threaddim=16;
  int n_HN_total= yn* blockdim* threaddim;
  …
  mcc<<<blockdim,threaddim>>>( d_MTOT,…参数);
  CUDA_SAFE_CALL(cudaThreadSynchronize());
  CUDA_ SAFE_CAL L (cudaMemcp y( h_MTOT,d_MTOT,
  n_HN_total*sizeof(float),cudaMemcpyDeviceToHost));
  for(int u=0;u<thread_total;u++){
  for(int v=0;v<yn;v++){
  MTOT[v]=MTOT[v]+h_MTOT[u*yn+v];
  …
  }
  }
  …
  }

  这种实现方法的好处是:1.无论线程数增加到多大,该代码都是可行的;2.其占用的内存是固定的不会随着线程数的增加而增大。  

参数传递过多问题

  CUDA对kernel函数中参数大小要求不超过128kb,但是在SRIM程序的并行化中需要传入和传出很多参数,假如只对这些参数简单地罗列,其远远大于128kb,那么我们可以采用结构体将这些参数传入。用结构体从CPU到GPU进行传值时,对于单个变量是很简单的,但如果对数组传值就需要一定的技巧了,以下举例说明。

  若结构体定义为

  typedef struct
  {float E0kev;
  int *N;
  …
  }Test;Kernel函数定义为
  __global__ void D_test(Test *d_test , int particleNum)
  {
  int i=threadIdx.x+blockDim.x*blockIdx.x;
  if(i<particleNum)
  {
  cuPrintf("%d ",(*d_test).N[i];
  }
  }

  结构体在CPU 端的定义及内存分配为:

  Test *test1=(Test *)malloc(sizeof(Test));//cpu
  (*test1).N=(int *)malloc(maxLayerNum*sizeof(int));

  若在GPU端结构体内存分配写为:

  Test * test2;
  CUDA_SAFE_CALL(cudaMalloc((void**)&test2,sizeof(Test )));
  CUDA_SAFE_CALL(cudaMal loc((void**)&(*t est2).N,
  maxLayerNum*sizeof(int )));

  则系统会报Segmentation fault错误。

  在实现时需要定义中间变量,并在GPU端为数组分配内存,代码如下:

  Test test2;
  CUDA_ SAFE_CALL (cudaMal l oc(( void* *)&t es t2 .N,
  maxLayerNum*sizeof(int )));
  CUDA_SAFE_CALL(cudaMemcpy(test 2.N,(* t est1).N,
  maxLayerNum*sizeof(int),cudaMemcpyHostToDevice));

  然后对GPU上的结构体变量定义,并为其传值:

  Test * d_test ;
  CUDA_SAFE_CALL(cudaMalloc((void**)&d_test,sizeof(Test )));
  CUDA_SAFE_CALL(cudaMemcpy(d_test,&test2,sizeof(Test),
  cudaMemcpyHostToDevice));

  SRIM并行化模拟结果

图1 H—>C在1000Kev时能损曲线对比

  图1为H 打C在1000kev时能损曲线对比图,横轴表示靶材料的深度,纵轴表示能损,是由前面总能损MTOT[i]除以总粒子数所得。左图为运行SRIM软件的能损曲线图,右图为用CUDA书写的GPU并行程序能损曲线图,可以看出来两个图的结果基本相同,证明了并行程序的正确性。

表1 P—>C的CPU与GPU时间对比

  表1列出了质子打碳时入射能分别在10kev、100kev、1000kev,入射粒子数分别为10000、100000、1000000时只统计能损和射程,运行全部程序的时间对比。可以看出,在粒子数比较少的情况下在CPU上运行要明显快于GPU,在粒子数比较多时GPU体现出了它的优越性,这是因为在粒子数比较少时,计算量很小,GPU在启动内核时占用很大一部分时间。并且采用该并行算法GPU的时间增长比较缓慢,那么我们可以看出粒子数越多该并行程序的优越性越明显。

  本文主要以SRIM为物理背景,讨论了在CUDA架构下的并行编程,通过串行程序与并行程序的时间对比,可以看出GPU并行程序的优越性。CUDA仅是并行计算的基础,为了提高效率我们试图可以考虑MPI+OpenMP+CUDA的三级并行编程模式,即引入节点内的细粒度并行以及节点间的粗粒度并行机制,形成一种更高效的并行策略。

  (作者单位:1为中国科学院近代物理研究所;2为兰州大学)