QiitaのJulia記事の新着を見ていると，コロケート格子＋PISO＋TVDの熱対流解析コードという記事がありました。 「改善点があれば教えて下さい」とのことなので，自分でも分かるような改善点があればと思ってソースを見ました。
普通のプログラミング言語なら明らかにパフォーマンスに影響が出そうなのは大きいfor文内のif文の部分です。
(でもJuliaなら最適化するとあまり変わらないのかも)
例えば，solveUA関数中で，次のように大きなループの中にif文があって，もったいないなぁと。
for j=1:ny+2 for i=1:nx+2 if i==nx+2 aeuva[i,j] = 0.0 elseif i==nx+1 aeuva[i,j] = 2nu*dy/dx else aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0)#- 0.25*dy*(ua[i+1, j] + ua[i,j]) end if i==1 awuva[i,j] = 0.0 elseif i==2 awuva[i,j] = 2nu*dy/dx else awuva[i,j] = nu*dy/dx + dy*max(0.5*(ua[i,j] + ua[i-1,j]), 0.0) end if j==ny+2 anuva[i,j] = 0.0 elseif j == ny+1 anuva[i,j] = 2nu*dx/dy else anuva[i,j] = nu*dx/dy + dx*max(-0.5*(va[i,j+1] + va[i,j]), 0.0)#- 0.25dx*(va[i,j+1] + va[i,j]) end if j==1 asuva[i,j] = 0.0 elseif j==2 asuva[i,j] = 2nu*dx/dy else asuva[i,j] = nu*dx/dy + dx*max(0.5*(va[i,j] + va[i,j-1]), 0.0) end apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + dx*dy/dt end end
ここは，
- aeuva
- awuva
- anuva
- asuva
- apuva
と分けて，それぞれでif文を無くせばいいのかなぁと。
まずは次の部分を確認してみます。
for j=1:ny+2 for i=1:nx+2 if i==nx+2 aeuva[i,j] = 0.0 elseif i==nx+1 aeuva[i,j] = 2nu*dy/dx else aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0)#- 0.25*dy*(ua[i+1, j] + ua[i,j]) end end end
iの条件でnx+2，nx+1，1:nxの3つに分かれます。次のような感じでしょうか。
for j=1:ny+2 aeuva[nx+2, j] = 0.0 end for j=1:ny+2 aeuva[nx+1, j] = 2nu*dy/dx end for j=1:ny+2 for i=1:nx aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0) end end
一つ目と二つ目の境界条件の式は全ての値を固定値に設定するので ループ文にするほどではないと個人的に感じるので pythonのnumpyっぽく書いてみます。
三つ目は単純な二重ループなので，Julia独特の書き方で書いてみます。
aeuva[nx+2, 1:ny+2] .= 0.0 aeuva[nx+1, 1:ny+2] .= 2nu*dy/dx for j=1:ny+2, i=1:nx aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0) end
というのをそれぞれ適用すると全体は次のような感じで書けるでしょうか?
aeuva[nx+2, 1:ny+2] .= 0.0 aeuva[nx+1, 1:ny+2] .= 2nu*dy/dx for j=1:ny+2, i=1:nx aeuva[i,j] = nu*dy/dx + dy*max(-0.5*(ua[i+1,j] + ua[i,j]), 0.0) end awuva[1, 1:ny+2] .= 0.0 awuva[2, 1:ny+2] .= 2nu*dy/dx for j=1:ny+2, i=3:nx+2 awuva[i,j] = nu*dy/dx + dy*max(0.5*(ua[i,j] + ua[i-1,j]), 0.0) end anuva[1:nx+2, ny+2] .= 0.0 anuva[1:nx+2, ny+1] .= 2nu*dx/dy for j=1:ny, i=1:nx+2 anuva[i,j] = nu*dx/dy + dx*max(-0.5*(va[i,j+1] + va[i,j]), 0.0) end asuva[1:nx+2, 1] .= 0.0 asuva[1:nx+2, 2] .= 2nu*dx/dy for j=3:ny+2, i=1:nx+2 asuva[i,j] = nu*dx/dy + dx*max(0.5*(va[i,j] + va[i,j-1]), 0.0) end for j=1:ny+2, i=1:nx+2 apuva[i,j] = aeuva[i,j] + awuva[i,j] + anuva[i,j] + asuva[i,j] + dx*dy/dt end
あとは自分の好み(実行速度はやや遅くなる)として書くコードをなるべく少なくしたい人なので，
ua = zeros(Float64, nx+2, ny+2) ua2 = zeros(Float64, nx+2, ny+2) ua3 = zeros(Float64, nx+2, ny+2) utemp = zeros(Float64, nx+2, ny+2) uold = zeros(Float64, nx+2, ny+2) apuva = zeros(Float64, nx+2, ny+2) aeuva = zeros(Float64, nx+2, ny+2) awuva = zeros(Float64, nx+2, ny+2) anuva = zeros(Float64, nx+2, ny+2) asuva = zeros(Float64, nx+2, ny+2) bua = zeros(Float64, nx+2, ny+2)
のところは
(ua, ua2, ua3, utemp, uold, apuva, aeuva, awuva, anuva, asuva, bua) = (zeros(Float64, nx+2, ny+2) for k in 1:11)
とかにしそうな感じです。 (11行 → 3行)
ただ，
# 次のように書くと uaとua2が同じものを指し示すのでバグの元 ua = ua2 = zeros(Float64, nx+2, ny+2)
と書くとバグの元なので注意が必要です。それぞれでzeros関数を呼ぶ必要があります。
0に初期化しないのであれば，
ua = zeros(Float64, nx+2, ny+2) ua2 = similar(ua)
とかで同じ型・サイズの行列を確保できて楽ちんなんですけどね。
歳をとると書くのがめんどくさいというより，見返した時に長い行数を読むのが面倒なのでとにかく見づらくても 短いソースを書こうとする傾向にあります。
「他の人が作ったライブラリを使って自分のソースはなるべく短く」それが年寄りJulia言語ユーザーujimushi流です。
なお，
for j=2:ny+1 for i=2:nx+1 ptemp[i,j] = pc[i,j] end end
とかただ要素のコピーをしている部分があるのですが，自分なら
ptemp[2:nx+1, 2:ny+1] .= pc[2:nx+1, 2:ny+1]
とかしてしまいそうです。ただの要素のコピーに5行もかけるのはもったいないとか思ってしまいます。
もちろん複雑な処理をしているところはfor文を使ってじっくり書いた方がいいと思います。
後，行列各要素の値を0に初期化する時は，例えば次のようなコードの時は =でつなげても問題ありません。
ua[nx+2, j] = 0.0 va[nx+2, j] = 0.0 ua[1, j] = 0.0 va[1, j] = 0.0 ua2[nx+2, j] = 0.0 va2[nx+2, j] = 0.0 ua2[1, j] = 0.0 va2[1, j] = 0.0
たとえば次のようにも書けます。
ua[nx+2, j] = va[nx+2, j] = ua[1, j] = va[1, j] = 0.0 ua2[nx+2, j] = va2[nx+2, j] = ua2[1, j] = va2[1, j] = 0.0
次のような感じにも書けますね。
ua[[nx+2, 1], j] .= 0.0 va[[nx+2, 1], j] .= 0.0 ua2[[nx+2, 1], j] .= 0.0 va2[[nx+2, 1], j] .= 0.0
まさにコードの行数を短くするためには何でもやるというコードの書き方です。
