CAMx平流模块GPU CUDA移植优化
本文是对Cao et al.(2023)论文的进一步解读。
GPU移植方案
移植路线：Fortran代码先改成C代码，然后进一步使用CUDA C扩展改写。
Fortran到C转换
Fortran到C手动重写（可以借助一些转换工具）需要注意以下几点：
实际上就是C/Fortran的符号链接、传参方式、数据类型对应、等价库函数以及存储方式。
标准C到CUDA C转换
由于Fortran不能直接调用CUDA C kernels，需要一个C接口程序用于传递参数和数据。
优化手段
减少CPU-GPU通信频率
原来PPM平流方案Fortran实现方式不适合GPU计算。如下图所示，PPM水平平流以一维方式实现，先做x方向，再做y方向。比如x方向，核心程序hadvppm是针对x方向的一行数据进行计算，在四重嵌套循环（垂直层/平流时间子步/y方向/平流物种）内调用该程序。
如果简单改写为C和CUDA C，每次只传输一行数据到GPU计算，CPU-GPU通信频繁且传输数据量小。因此，重构PPM实现，将hadvppm改写为对3D场一次性做x方向的PPM计算，这就要求将传递数据保存到三维或四维数组中，一次性传给GPU。
GPU线程组织
原先CUDA版本对三维空间点(i,j,k)的循环替换为只用线程索引的索引计算，同时计算x或者y方向的格点。
i=threadIdx.x + blockIdx.x*blockDim.xCUDA线程组织方式如下，每个线程块256个线程，根据水平格点数向上取整线程块数，Grid/Block都是一维组织方式
新方案使用thread和block索引同时计算所有水平格点(i,j)，每个线程块由1024个线程组成（一维），而Grid是二维组织方式。
i=threadIdx.x + blockIdx.x*blockDim.x
j=blockIdx.y
对应CUDA kernel中循环使用方式如下，其中i/j循环替换为if判断，确保不越界。
GPU大规模并行在编程上将循环改成了IF判断（确保不越界）。
优化GPU数据访存效率
循环顺序与数据在内存中布局（Fortran列主元）相同。
MPI + CUDA
- 利用CUDA库函数查询每个节点的GPU卡数
- 使用MPI进程号和取余操作确定每个节点启动的GPU的ID
- 使用CUDA库cudaSetDevice函数绑定CPU核心和GPU卡
绑定GPU代码：
其它
Fortran多维数组索引和C一维索引转换可以使用宏定义切换
#define I3(i,j,k) ((k) * (jdim) * (idim) + (j) * (idim) + (i))
#define I4(i,j,k,ispc) ((ispc)*(kdim)*(jdim)*(idim)+(k)*(jdim)*(idim)+(j)*(idim)+(i))使用起来符合Fortran多维索引习惯，
d_fm[I4(i,j,k,ispc)] = 0.0f;文献
Cao, K., Wu, Q., Wang, Lingling, Wang, N., Cheng, H., Tang, X., Li, D., Wang, Lanning, 2023. GPU-HADVPPM V1.0: a high-efficiency parallel GPU design of the piecewise parabolic method (PPM) for horizontal advection in an air quality model (CAMx V6.10). Geoscientific Model Development 16, 4367–4383. https://doi.org/10.5194/gmd-16-4367-2023
