Accelerating thePerformance of Modal
Aerosol Module (MAM)of E3SM Using OpenACC
Hongzhang Shan1), Zhengji Zhao2),and Marcus Wagner3)
1)

Lawrence Berkeley National Laboratory,2)
National Energy Scientific Computing Center
3)
Cray, Inc

WACCPD 19, Denver COOutline● Introduction● MAM Algorithms and Kernels● Offloading Four Representative Kernels in MAM to GPU Using
OpenACC directives● MAM Performance on Summit● Summary and Conclusion

Introduction● Energy Exascale Earth System Model (E3SM) is a state-of-the-artearth system simulation code
○It has a large code base with over a million lines of Fortran code

○Production code are currently optimized for advanced CPU systems● Making effective use of GPUs, however, remains a challenge● In this work, using the modal aerosol module (MAM) of E3SM as adriving example, we investigated how to effectively offload
computational tasks to GPUs

● We chose to work with OpenACC directives

MAM Algorithms and Kernels-4-Modal Aerosol Module (MAM) of E3SM●●

●E3SM was developed to reliably project decade-to-century scale changes that
could critically impact the U.S. energy sector. It combined the atmosphere,ocean, land, river, ice, and other components.
The computation of the atmosphere component is based upon the Spectral
Element (SE) numerical discretization of underlying PDEs for stratified,hydrostatic fluid dynamics on rotating spheres.
MAM is a submodule from the atmosphere component that plays an importantrole in the climate system by influencing the Earth's radiation budgets andmodifying cloud properties. It predicts the mass and mixing ratios of cloudliquid and cloud ice, diagnoses the mass and mixing ratios of rain and snow,and handles complicated conversions between cloud hydrometeors.

MAM in E3SMa)

●●
●
●●b)

E3SM models the Earth with a cubed-sphere grid (6 faces) as shown in Fig. a).
The resolution of the meshes is defined as the number of spectral elements ne along theedge of each cube face. (6ne2 elements total in the mesh)
Each element contains a np*np tensor product of Gauss-Lobatto-Legendre (GLL)
points depicted in Figure b), the number of unique points (physics columns)
There exists another dimension, namely the vertical direction (except the sphere faces).
Computations between physical columns are independent

MAM in E3SM
Loop structure for computations
do j = 1, nchunks
do k = 1, nlev
do i = 1, ncols(j)!number of chunks!vertical levels, may have data dependency!number of columns in chunk j!sum(ncols(j) = total physical columns
computation_kernels() !many different kernels
enddo
enddo
enddo

•
•
•

In the parallel implementation physical columns are distributed among theprocesses based on a set of load balancing strategies.
To get better caching effects, all the columns assigned to a process will begrouped in a data structure called a chunk.
In each chunk, a maximum number of columns PCOL is specified atcompilation time.

Experiment Configuration-8-Experiment System – Summit at ORNLo Theoretical peak ~200 PF (dp),o Each Summit node has two IBM Power9processors with six Nvidia V100 GPUs.

o Power9 CPUs are connected with GPUs
through dual NVLINKo 512 GB of DDR4 memory for Power9CPUS and 96 GB of HBM for GPUs.

o Each Nvidia V100 has 80 SMs, 16 GB ofHBM, and a 6 MB L2 cache.

o Each SM contains 64 FP32, 64 INT32, 32
FP64 cores; partitioned to four processing
blocks, each with a warp scheduler.

Experiment Setup●E3SM: https://github.com/E3SM-Project/E3SM.git, branch shz0116/cam/cam_openacc● We used the PGI compiler version 19.4, Spectrum MPI version
10.3.0.1, and CUDA 10.1.168● Other libraries used in the E3SM code included NETCDF 4.6.1,
NETCDF-FORTRAN 4.4.4, ESSL 6.1.0 , Parallel NETCDF 1.8.1, andHDF5 1.10.3● The data set for E3SM is SMS_PS_Ld5.ne16_ne16.FC5AV1C-L,
which stresses the atmosphere physics. Here, ne16 ne16 defines thecubed sphere grid resolution

Offloading Four Representative Kernels in
MAM to GPU Using OpenACC directives- 11 -

Offloading MAM Kernels to GPUs
● Data transfer was not trivial in MAM○

MAM has a large code base with tens of thousands of lines of source code○

MAM does checkpointing with various I/O operations scattered all over the code○an excessive number of temporary subroutines or function variables need to bepromoted and explicitly allocated on the GPU memory as well

● MAM has a flat profile, its run time is distributed across many functions,meaning we could not focus on just a couple of loops● The programming effort needed to optimize different kernels also varies
significantly by kernel. Some required a significant code refactoringKernel: subgrid_mean_updraft
1 !$acc parallel loop collapse(2) copyin(wsig,w0) copyout(ww) private(zz,wa)
2 do k = 1, pver3
do i = 1, ncol
4sigma = max(0.001_r8, wsig(i,k))5
wlarge = w0(i,k)
6xx = 6._r8 * sigma / nbin7
do ibin = 1, nbin!constant nbin=50
8yy = wlarge - 3._r8*sigma + 0.5*xx
9yy = yy + (ibin-1)*xx
10zz(ibin) = yy * exp(-1.*(yy-wlarge)**2/(2*sigma**2))/(sigma*sqrt(2*pi))*xx
11
if (zz(ibin) .gt. 0._r8) then12wa(ibin) = zz(ibin)
13
else
14
wa(ibin) = 0._r8
15
endif
16
end do17
sum_wa = sum( wa(:))
18
if (sum_wa .gt. 0._r8) then19
ww(i,k) = sum_wa
20
else
21
ww(i,k) = 0.001_r8
22
end if23
enddo
24 enddo-calculates the mean updraft velocity

Kernel Performance on a Summit NodeKernel: hetfrz_classnuc_cam_calc
1 !$acc declare create(ncnst, nmodes, ...)!create module variables on device
2 !$acc update device(ncnst, nmodes, ...)
3 ...
4 !$acc enter data create(total_aer_num, ...)
5 ...
6 !$acc parallel loop collapse(2) private(fn), copyin(t,pmid) &7 !$acc& copyout(frzbccnt,...) default(present)8
do k = top_lev, pver9
do i = 1, ncol
10if (t(i,k) .gt. 235.15_r8 .and. t(i,k) .lt. 269.15_r8) thenLines 11-34 in next slide35
end if36
end do37
end do38 ...
39 !$acc exit data delete(total_aer_num, ...)

-calculates the heterogeneous freezing
rates from classical nucleation theory

Kernel: hetfrz_classnuc_cam_calc
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

qcic = min(qc(i,k)/lcldm(i,k), 5.e-3_r8)
ncic = max(nc(i,k)/lcldm(i,k), 0._r8)
con1 = 1._r8/(1.333_r8*pi)**0.333_r8
r3lx = con1*(rho(i,k)*qcic/(rhoh2o*max(ncic*rho(i,k), 1.0e6_r8)))**0.333_r8
r3lx = max(4.e-6_r8, r3lx)
supersatice = svp_water(t(i,k))/svp_ice(t(i,k))!svp_water and svp_ice are two subroutines
fn(1) = factnum(i,k,mode_accum_idx)
if (nmodes == MAM3_nmodes .or. nmodes == MAM4_nmodes) thenfn(2) = factnum(i,k,mode_accum_idx)
fn(3) = factnum(i,k,mode_coarse_idx)
else if (nmodes == MAM7_nmodes) thenfn(2) = factnum(i,k,mode_finedust_idx)
fn(3) = factnum(i,k,mode_coardust_idx)
end ifcall hetfrz_classnuc_calc( &deltatin, t(i,k), pmid(i,k), supersatice,&fn, r3lx, ncic*rho(i,k)*1.0e-6_r8, frzbcimm(i,k), frzduimm(i,k),&frzbccnt(i,k), frzducnt(i,k), frzbcdep(i,k), frzdudep(i,k), hetraer(:,i,k), &awcam(:,i,k), awfacm(:,i,k), dstcoat(:,i,k), total_aer_num(:,i,k), &coated_aer_num(:,i,k), uncoated_aer_num(:,i,k), &total_interstitial_aer_num(:,i,k), &total_cloudborne_aer_num(:,i,k), errstring)!hetfrz_classnuc_calc is a sequential routine with hundreds of lines}

Kernel Performance on a Summit Node

Kernel: ccncalc
1 do k=top_lev,pver2
do i=1,ncol
3a(i)=surften_coef/tair(i,k)
4smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
5
enddo
6
do m=1,ntot_amode7
phase=38
call loadaer(state, pbuf, 1, ncol, k, &9m, cs, phase, naerosol, vaerosol, hygro)
10!get data from pbuf to naerosol, vaerosol, and hygro
11where(naerosol(:ncol) .gt. 1.e-3_r8)
12amcube(:ncol)=amcubecoef(m)*vaerosol(:ncol)/naerosol(:ncol)
13
sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol)*amcube(:ncol))14
elsewhere15
sm(:ncol)=1._r8
16
endwhere17
do l=1,psat18
do i=1,ncol
19
arg(i)=argfactor(m)*log(sm(i)/super(l))
20
ccn(i,k,l)=ccn(i,k,l)+naerosol(i)*0.5_r8*(1._r8-erf(arg(i)))
21
enddo
22
enddo
23
enddo
24 enddo-calculate the number of concentrations ofaerosols activated when cloud condensation
nuclei are at supersaturation

1 do k=top_lev,pver
2 do m = 1, ntot_amode3
call loadaer(state, pbuf, 1, ncol, k, &4
m, cs, phase, naerosol(:,m,k), vaerosol(:,m,k), hygro(:,m,k))
5!define naerosol, vaeroosol, hygro as 3D instead of 1D
6 enddo
7 enddo
9 $acc data copy(ccn) copyin(vaerosol, naerosol, hydro) &8
10 $acc& copyin(super,amcubecoef,argfactor,tair,smccoefcoef,surften_coef)

GPU version

11 $acc parallel loop private(a,smcoef,arg,sm,amcube,m,i,l) default(present)
12 do k=top_lev,pver13
do i=1,ncol
14a(i)=surften_coef/tair(i,k)
15smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
16
enddo
17
do m=1,ntot_amode18
phase=3
19
where(naerosol(:ncol) .gt. 1.e-3_r8)
20
amcube(:ncol)=amcubecoef(m)*vaerosol(:ncol.m.k)/naerosol(:ncol,m,k)21sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol,m,k)*amcube(:ncol))22
elsewhere23
sm(:ncol)=1._r8
24
endwhere25
do l=1,psat26
do i=1,ncol27
arg(i)=argfactor(m)*log(sm(i)/super(l))
28
ccn(i,k,l)=ccn(i,k,l)+naerosol(i,m,k)*0.5_r8*(1._r8-erf(arg(i)))
29
enddo
30
enddo
31
enddo
32 enddo
33 ...
34 $acc end data

Kernel Performance on a Summit Node1 do i = 1, ncol
2 ... !more than 400 lines of code3
do n = 1, nsubmix
4
qncld(:) = qcld(:)
5nnew <--> nsav !nnew = 1, nnsav=06srcn(:) = 07
do m = 1, ntot_amode
8mm = mam_idx(m,0)
9srcn(top_lev:pver-1) = srcn(top_lev:pver-1) + &10nact(top_lev:pver-1,m)*raercol(top_lev+1:pver,mm,nsav)
11
tmpa = raercol(pver,mm,nsav)*nact(pver,m) + &12
raercol_cw(pver, mm, nsav) * (...)
13
srcn(pver) = srcn(pver) + max(0.0_r8, tmpa)
14
enddo
15
call explmix(qcld, srcn, ..., qncld) !compute qcld from qncld
1617do m = 1, ntot_amode
18mm = mam_idx(m,0)
19source(top_lev:pver-1) = &20
nact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
21
tmpa = ... !same as line 922source(pver) = max(0.0_r8, tmpa)
23
call explmix(raercol_cw(:, mm, nnew), source, ..., raercol_cw(:, mm, nsav), ...)
24!compute raercol_cw(,,nnew) from raercol_cw(,,nsav)
25
call explmix(raercol(:, mm, nnew), source, ..., raercol(:, mm, nsav), &26
raercol_cw(:, mm, nsav))
27!compute raercol(,,nnew) from raercol(,,nsav) and raercol_cw(,,nsav)
28
do l = 1, nspec_amode(m)29mm = mam_idx(m, l)30source(top_lev:pver-1) = !same as line 17 except using mact instead nact31
tmpa = !same as line 19 except using mact instead nact variable
32
source(pver) = max(0.0_r8, tmpa)
33
call explmix !same as line 21
34
call explmix !same as line 23
35
enddo
36
enddo
37
enddo
38 ...

5.4 Kernel: nsubmix

Kernel: nsubmix

Kernel: restructured nsubmix
1 do mm = 1, ncnst_tot
2
do k = top_lev, pver3m= mam_idx_1d(1, mm)
4
l= mam_idx_1d(2, mm)
5
kp1 = min(k+1,pver)6
km1 = max(k-1,top_lev)7
if (l == 0) then8tmpa = nact(k,m)*raercol(kp1,mm,nsav)
9
if (k == pver) then10
tmpa = tmpa + raercol_cw(pver,mm,nsav)*(nact(pver,m) - taumix)
11
tmpa = max(0.0_r8, tmpa)
12
endif
13
else
14
tmpa = mact(k,m)*raercol(kp1,mm,nsav)
15
if (k == pver) then16tmpa = tmpa + raercol_cw(pver,mm,nsav)*(mact(pver,m) - taumix)
17
tmpa = max(0.0_r8, tmpa)
18
endif
19
endif
20
call explmix(raercol_cw(k, mm, nnew), source, ...)
21
call explmix(raercol(k, mm, nnew), source, ...)
22end do23 end doKernel Performance on a Summit Node

MAM Kernel Performance on Summit- 24 -Kermel Performance Analysis●●MAM kernels are light, the average run time is within milliseconds.
Secondly, the parallelism mainly comes from the vertical level (pver=72) and thenumber of columns in a block (pcol). Considering that the Nvidia GPUs use a warp sizeof 32 as the scheduling unit, neither pver(= 72) nor ncol(<= 16) is a perfect fit,
resulting in thread resource waste.

●To improve the performance from an application perspective, the size of pcol andpver can therefore be aligned to multiples of warp size.
Kernel Average runtime on GPUs and CPUS

Kernel Performance vs PCOL Values

Kernel Performance vs MPSBoth CPU and CPU+GPU Versions Scale Well across Multiple Nodes

Summary and Conclusion● We investigated if GPUs can be used to accelerate the performanceof MAM, a module of E3SM on Summit using OpenACC.

● We have achieved over a 5X performance speedup by offloadingsome of the kernels to Nvidia Volta GPUs.

● The results revealed that under the current E3SM configuration for
product runs, some parameter settings do not suit offloading, such as
the number of columns per block and the number of vertical levels.
○These settings not only severely limit the degree of parallelism but also fail tomake effective use of GPU thread resources, becoming a performance bottleneck.

Summary and Conclusion (cont.)● Run time scatters across many kernels, and each computational
kernel is relatively light, with average run time in milliseconds or lessper call. Such light computational kernels particularly require
OpenACC implementation to further reduce overhead from kernel
launching, data transfer, etc.

● Performance was primarily limited by the kernel nsubmix, which did
not benefit from GPU offloading.

● We are looking into overlapping computations and data transfer using
async OpenACC directives and possibly merging kernel computations
to improve MAM’s performance.

Acknowledgement● All authors from Lawrence Berkeley National Laboratory were
supported by the Office of Advanced Scientific Computing Researchin the Department of Energy Office of Science under contract number
DE-AC02-05CH11231.

● This research used resources of the Oak Ridge Leadership
Computing Facility at the Oak Ridge National Laboratory, which is
supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC05-00OR22725.

Thank You!