海洋大循環モデル PFESコードは、POM (Princeton Ocean Model) をベースにした、 大気モデル計算と海洋モデル計算を交互に繰り返すカップリングコードです。
POMは、Princeton大学NOAA-Dynalysis研究所のMellor,G.L.が開発した
3次元潮流解析プログラムで、広域な海洋循環や沖合の近海潮流に適しています。
近年、そのモデルが一般公開されたことから，多くの研究者によって
使用されています。
PFESコードの主構造は、以下のように時間発展ループ中で、
大気モデル処理と海洋モデル処理を交互に繰り返す、というものです。
海洋モデル処理ルーチンからは、さらに、多数の下位計算ルーチンが
呼び出されます。
do iintp = 1, iend ! 海洋モード時間ループ
do iatmos = 1, atm_per_ocn ! 大気モデル時間ループ
call atmos ! 大気モデル処理
enddo
call ocean ! 海洋モデル処理
enddo
PFESコードの主なモジュール構成は以下のとおりです。
定数データ、パラメータの宣言モジュール
海洋モデル用データ宣言モジュール ocean_common
海洋モデルルーチンoceanを含むモジュール ocean_model
海洋モデルの下位計算ルーチン
大気モデルルーチンatmosを含むモジュール atmosphere_model
大気モデルから海洋モデルへの中継ルーチン
カップリングデータ宣言モジュール coupling_data
主プログラム
作業用データの宣言モジュール
まず、主要部分であるatmos、ocean、および oceanの下位計算ルーチン群をみると、 ほとんどのループネストにおいて、2次元配列の2次元目、あるいは3次元配列の2次元目 をアクセスする並列化可能ループが、最外側ループとなっていることが わかります。
例えば、海洋モデル用データ宣言モジュールocean_common 中で宣言されている3次元配列advx,advyは、以下のようなループで 参照されており、
do j = 1, jmx
do k = 1, kbm1
do i = 1, im
adx2d(i,j) = adx2d(i,j) + advx(i,j,k)*dz(k)
ady2d(i,j) = ady2d(i,j) + advy(i,j,k)*dz(k)
drx2d(i,j) = drx2d(i,j) + drhox(i,j,k)*dz(k)
dry2d(i,j) = dry2d(i,j) + drhoy(i,j,k)*dz(k)
aam2d(i,j) = aam2d(i,j) + aam(i,j,k)*dz(k)
enddo
enddo
enddo
2次元配列dx,dyは、以下のようなループで参照されています。
do j = j2, jmx
do i = 2, im
aru(i,j)=.25*(dx(i,j)+dx(i-1,j))*(dy(i,j)+dy(i-1,j))
arv(i,j)=.25*(dx(i,j)+dx(i,j-1))*(dy(i,j)+dy(i,j-1))
enddo
enddo
(ここで、im = 2164、 jmx = 813、 kbm1 = 20です。)
そこで、カップリングデータが宣言されているモジュールcoupling_data 、海洋モデル用データが宣言されているモジュールocean_common 、およびサブルーチンoceanを含むモジュールocean_model の主要配列を、2次元目でマップすることにします。
ここで、注意しなければならないのは、 PFESコードにおける配列の上下限値が、例えば以下のように、 それぞれ異なっている、と言うことです。
integer, parameter :: jbh = 1, jth = 0
integer, parameter :: jbdx = 1, jtdx = 1
integer, parameter :: jbva = 0, jtva = 1
real(kind=rt) :: h(im,1-jbh :jmd+jth ),
& dx(im,1-jbdx :jmd+jtdx ),
& va(im,1-jbva :jmd+jtva ),
このため、DISTRIBUTE指示文でBLOCK分散を指定するだけでは、
各抽象プロセッサへマップされる要素が少しずつずれてしまい、
リモートアクセスが増大してしまいます。
そのため、以下のように、全ての配列の上下限値を包含するような
1次元のテンプレート配列TPLを宣言し、それをBLOCK分散します。
そして、全ての配列をテンプレート配列TPLに対して、整列させることにします。
!HPF$ TEMPLATE,DISTRIBUTE (BLOCK) :: TPL(0:jmd+1)
例えば、カップリングデータ宣言モジュールCPL_DATA.F 中には、 以下のようなマッピング指示文を挿入します。
!HPF$ TEMPLATE,DISTRIBUTE (BLOCK) :: TPL(0:jmd+1) !HPF$ ALIGN (*,i) WITH TPL(i) :: cld,tatm,hum,swrad !HPF$&,workl1,workl2,wusurf,wvsurf,wtsurf,wssurf !HPF$&,egf, egb,Wsw, Wse, Wnw, Wne,alon, alat !HPF$ ALIGN (*,i,*) WITH TPL(i) :: t, tb
ただし、OCEAN.F中の2次元配列sbuf(jmg,4)に関しては、 例えば以下のように、1次元目に対応するループが並列化対象であるため、 一次元目でマップすることにします。
do j = 1, jmx
do k = 1, kbm1
do i = 3, im-2
darea = dx(i,j)*dy(i,j)*fsm(i,j)
dvol = darea*dt(i,j)*dz(k)
sbuf(j,1) = sbuf(j,1) + dvol
sbuf(j,2) = sbuf(j,2) + tb(i,j,k)*dvol
sbuf(j,3) = sbuf(j,3) + sb(i,j,k)*dvol
sbuf(j,4) = sbuf(j,4) + 0.125*( (u(i,j,k)+u(i+1,j,k))**2
& + (v(i,j,k)+v(i,j+1,k))**2)*dvol
enddo
enddo
enddo
さらに、oceanから引用される下位計算ルーチンには 配列の仮引数を持つものがあり、それらの仮引数にも、 対応する実引数と同一のマッピング指示文を指定します。
この段階で、一度コマンドラインオプション-Minfo付で
翻訳してみることにしましょう。
すると、COMMON_REEDFLUX.Fの以下の部分で
call radinp(dodiavg ,clat, alon(i,j), calday, dayspy,
以下のようなエラーメッセージにより翻訳が終了します。
Scalar element of nonsequential array alon cannot be passed to dummy array argument
これは、
非順序的な配列alonの配列要素を実引数として引用する、という
HPF仕様に違反した記述が含まれているためです。
呼出元手続における配列alonの宣言は、以下の通りであり、
real(kind=rt), intent(in) :: alon(im,jmd)
サブルーチンradinpの対応する仮引数は、以下のように宣言されているので、
real(kind=rt), intent(in) :: clon(im)
実引数の形状を仮引数の形状と一致させるため、 以下のように部分配列実引数に書き直します。
call radinp(dodiavg ,clat, alon(:,j), calday, dayspy,
HPFでは、SEQUENCE指示文を指定しない場合、
基本的には実引数と仮引数の形状を一致させる必要がありますので、
他にも、このような不一致がないかどうかチェックしてみるのが良いでしょう。
すると、COUPLER_LIB.F中のサブルーチンmap_atm2ocn
で引用されている、COMMON_REEDFLUX.F中のサブルーチンreedfluxは、
実引数が、3次元配列tb(im,1-jbtb :jmd+jttb ,kb)であるのに対して、
仮引数がtb(im,1-jbtb:jmd+jttb)となっており、形状が
一致していないことが分かります。
そのため、サブルーチンreedfluxの仮引数を、
以下のように、実引数の形状と一致するよう修正し、
! real(kind=rt), intent(in) :: tb(im,1-jbtb:jmd+jttb)
real(kind=rt), intent(in) :: tb(im,1-jbtb:jmd+jttb,kb)
さらに、サブルーチンreedfluxから、配列tbを実引数として 引用されるサブルーチンreedflxdに関しても、 対応する仮引数stempの形状と、 stempの参照の添字を、それぞれ、3次元配列になるよう 以下のように修正しておきます。
! real(kind=rt), intent(in) :: stemp(im,1-jbstmp:jmd+jtstmp)
real(kind=rt), intent(in) :: stemp(im,1-jbstmp:jmd+jtstmp,kb)
! ts = stemp(i,j)+tk + tbias
ts = stemp(i,j,1)+tk + tbias
ここで、改めて、HPF/ESにより、コマンドラインオプション-Minfo付で
翻訳してみましょう。
その結果、翻訳時のメッセージが出力されますが、
まず"array copied"というメッセージを探してみます。
これは、ループを並列化する際、そのままでは
並列化の基準配列とのマッピングが合わなかったため
ブロードキャスト通信が発生したことを意味しています。
すると、例えば、
サブルーチンprofq中の局所配列kn(im,jmd,kb)、boygr(im,jmd,kb)
に対してマッピングを指定しなかったために、以下のdo jのループを
並列化したときに、各抽象プロセッサ上で設定されるboygr
の値を、全プロセッサで共有するためのブロードキャスト通信が
発生していることが分かります。
do j = 1, jmx
do k = 2, kbm1
do i = 1, im
q2b(i,j,k) = abs(q2b(i,j,k))
q2lb(i,j,k) = abs(q2lb(i,j,k))
boygr(i,j,k) = gee*((rho(i,j,k-1)-rho(i,j,k))
& /(dzz(k-1)*tps(i,j)))
& +gee**2*2.*1.025/(dtef(i,j,k-1)**2+dtef(i,j,k)**2)
enddo
enddo
enddo
このような通信は、kn(im,jmd,kb)、boygr(im,jmd,kb) を2次元目でマップすることで抑制できるので、以下のように マッピング指示文を追加します。
!HPF$ ALIGN (*,i,*) WITH TPL(i) :: kn,boygr
同様に、サブルーチンprofu、profv中の局所配列dhに対しても マッピングを指定することで、"array copied"という メッセージを消すことができます。
次に、 "expensive communication"というメッセージが
出力されている個所を探しましょう。
これは、オーバヘッドの高いパタンの通信が生成されたことを意味しています。
このメッセージが出力されているコードとして、
例えば、OCEAN.F中の以下のループがあります。
do j = 1, jmx
do i = 1, im
if(abs(vaf(i,j)).ge.vamax) vamax = abs(vaf(i,j))
enddo
enddo
このループでは、vafの絶対値の最大値をvamaxに求めていますが、 現在のところ HPF/ESは、 自動的にこのような集計演算を検出することができないため、 以下のように、最大値の算出をHPF2.0仕様の集計文の形式に 合うよう書き直して、REDUCTION節付のINDEPENDENT指示文を指定します。
!HPF$ INDEPENDENT,NEW(i),REDUCTION(vamax)
do j = 1, jmx
!HPF$ INDEPENDENT
do i = 1, im
vamax = max(vamax,abs(vaf(i,j)))
enddo
enddo
また、以下のようなループでも、"expensive communication" というメッセージが出力されています。
do j = 1, idwids
if (j .ge. jst(myid) .and. j .le. jend(myid)) then
jx = j - jend(myid-1)
do k = 1, kb-1
do i = 1, im
tsbcv = (1.e0-tday)*tsbc1(i,j,k)
& + tday*tsbc2(i,j,k)
ssbcv = (1.e0-tday)*ssbc1(i,j,k)
& + tday*ssbc2(i,j,k)
uf(i,jx,k) = uf(i,jx,k) +
& gc1*float(idwids-j)
& *(max(tsbcv-tbias,tmin)-t(i,jx,k))*fsm(i,jx)
vf(i,jx,k) = vf(i,jx,k) +
& gc1*float(idwids-j)
& *(ssbcv-sbias-s(i,jx,k))*fsm(i,jx)
enddo
enddo
endif
enddo
このループでは、jend(myid-1)の値が0のため、
例えば、uf(:,j,:)を基準として並列化したときに、全く
通信の必要なくjのループで並列化できることが分かります。
そのため、INDEPENDENT指示文を指定すると共に、
uf(j)を基準として並列化するようON指示文を指定し、
さらに、LOCAL節に、全ての配列を指定してみましょう。
修正後のループは、以下のようになります。
!HPF$ INDEPENDENT,NEW(j,k,i,tnbcv,snbcv)
do j = 1, idwids
!HPFJ ON HOME(uf(:,j,:)), LOCAL(uf,fsm,t,vf,s) BEGIN
if (j .ge. jst(myid) .and. j .le. jend(myid)) then
jx = j - jend(myid-1)
do k = 1, kb-1
do i = 1, im
tsbcv = (1.e0-tday)*tsbc1(i,j,k)
& + tday*tsbc2(i,j,k)
ssbcv = (1.e0-tday)*ssbc1(i,j,k)
& + tday*ssbc2(i,j,k)
uf(i,jx,k) = uf(i,jx,k) +
& gc1*float(idwids-j)
& *(max(tsbcv-tbias,tmin)-t(i,jx,k))*fsm(i,jx)
vf(i,jx,k) = vf(i,jx,k) +
& gc1*float(idwids-j)
& *(ssbcv-sbias-s(i,jx,k))*fsm(i,jx)
enddo
enddo
endif
!HPFJ ENDON
enddo
同様のパタンは、OCEAN.F中に合計3個所見つかります。
次に、"reduction inlined"というメッセージを探しましょう。
これは、集計演算ループとして並列化されたことを意味しますが、
現在のところHPF/ESは、
集計演算を自動的に並列化する場合、
ループ分割を行った後、ループ中の個々の文に関して並列化を行うので、
集計演算を行うDOループには、
以下のように、REDUCTION節付のINDEPENDENT指示文を指定した方が、
並列化の粒度が大きくなります。
このようなループは、OCEAN.F中に3つ見つかります。
!HPF$ INDEPENDENT,REDUCTION(atot,eaver)
do j = 1, jmg
atot = atot +sbuf(j,1)
eaver = eaver+sbuf(j,2)
enddo
最後に、もう1度コマンドラインオプション-Minfo付で 翻訳して、これまでのべたようなメッセージが無くなっていることを 確認します。
オリジナルソースと比較した修正内容は、以下の通りです。
| Fortranソースに対する変更分 | 69行 |
|---|---|
| Fortranソースに対する増加分 | 68行 |
| Fortranソースに対する修正量の割合 | 137 / 8076 = 0.017 |
他の文の修正としては、これまで説明した、実引数と仮引数の形状を 一致させるための修正 などが含まれます。
| 修正内容 | 行数 |
|---|---|
| INCLUDE行の挿入 | 3行 |
| INDEPENDENT指示文の挿入 | 7行 |
| ON指示構文の挿入 | 4行 |
| マッピングのための指示文の挿入 | 27行 |
| 計時ルーチン関連の文の挿入 | 27行 |
| 他の文の修正 | 69行 |