はじめてJuliaさわってみた！　Juliaで流体シミュレーション　①移流方程式の写経

## はじめに
数値流体力学やってみたいけど、ハードル高いしPythonおそいしなーとか思っている中で以下のような記事を見かけました。

- 【Python】流体シミュレーション：移流方程式を実装する
https://qiita.com/KQTS/items/354c85adb7d46e28e8da

- 1次元シュレーディンガー方程式を有限要素法で解いてみる in Julia
https://qiita.com/cometscome_phys/items/fbeaa92fcb4b6d8c6240

Juliaで書いたらそれっぽいコードでもそれなりに早く動くのでは？  
ということで、自分の知識には数値流体力学もPythonもJuliaも、何ならMarkupも微分方程式もまったくわかってないけど、とりあえずやってみようという記事です。


とりあえずの目標はPythonでやってるやつをJuliaで動かしてみて、Pythonのコード（高速化は全然されていない）と速度比較や書きやすさ比較してみようというものです。@KQTSさんのコードそのままうつしてます。ありがとうございます！（勝手に使ってごめんなさい）


Juliaの操作は


- Juliaで最低限やっていくための配列操作まとめ
https://qiita.com/A03ki/items/007be353411d19952ef7

- 【Julia初心者向け】Juliaで内包表記
https://qiita.com/keta-sts136/items/c079aab0b4edf29c0472

- JuliaでPyPlot.jl (matplotlib) を使ったヒストグラムのアニメーションを効率的に作成する
https://qiita.com/skotaro/items/57f426d617d42aa4cbaf

- PyPlot and Plots
https://nbviewer.jupyter.org/gist/genkuroki/3d6dbf52a3e52eb7c664bc88632c81d3#PyPlot-and-Plots

などなどみてます。いろいろツッコミどころ多いと思うのでこうしたほうがいいよってのあったら気軽にコメントもらえるとうれしいです。

なお環境はJulia1.5.2, Windows10 HomeでJupyterでIJuliaで動かしてます。Windows版ダウンロードしてきてインストールしてPkg add IJuliaだけしたらAnacondaのJupyterから動いた！環境変数いれてないのでJupyter上でPkgまわり詰まるけど気にしない！

## 写経
とりあえず、[この記事](https://qiita.com/KQTS/items/354c85adb7d46e28e8da)のものを書き写しています。式の意味とかは元の記事をみてくださいね。

Pythonの元の式から変更あった部分についてはコメントアウトして書き直してます。


そもそも自分はPythonもイマイチかけないので、np.arrangeってなんだっけからだったという。Vectorでかいたけど、[1:Num_stencil_x]で十分かも？  配列番号は1から始まるので1からにしてみたけど、ここは0からでよかったかも？後でちょっと出ますが、既に混乱してはまってます。np.whereはリスト内包表記で。

```julia
# Creat square wave
Num_stencil_x = 101
#x_array = np.arange(Num_stencil_x)
x_array =Vector(1:Num_stencil_x)
#u_array = np.where(((x_array >= 30) |(x_array < 10)), 0.0, 1.0)
u_array = [(if (i>=30) | (i < 10) 0.0 else 1.0 end) for i in x_array]
```
maximumにたどり着くのに時間が地味にかかった。しかもこれ、1から始めたので実は間違ってるという。

```julia
u_lower_boundary = 0.0
u_upper_boundary = 0.0
Time_step = 200

#Delta_x = max(x_array) / (Num_stencil_x-1)
Delta_x = maximum(x_array) / (Num_stencil_x-1)

C = 1
Delta_t = 0.2
CFL = C * Delta_t / Delta_x
total_movement = C * Delta_t * (Time_step+1)

#exact_u_array = np.where(((x_array >= 30 + total_movement) |(x_array < 10 + total_movement)), 0.0, 1.0)
exact_u_array = [(if (i>=30 + total_movement) | (i < 10 + total_movement) 0.0 else 1.0 end) for i in x_array]

```
PyPlot使えばそのままmatlibplot使える！便利！

Jなお、PyPlotのロードにちょっと時間がかかる。しかも初回pkgにaddしたり、動かしたら色々ダウンロードはじめてだいぶ時間かかった。

```julia
# plt.plot(x_array, u_array, label="Initial condition")
# plt.plot(x_array, exact_u_array, label="Answer")
# plt.legend(loc="upper right")
# plt.xlabel("x")
# plt.ylabel("u")
# plt.xlim(0, max(x_array))
# plt.ylim(-0.5,1.5)
using PyPlot
PyPlot.plot(x_array, u_array, label="Initial condition")
PyPlot.plot(x_array, exact_u_array, label="Answer")
PyPlot.legend(loc="upper right")
PyPlot.xlabel("x")
PyPlot.ylabel("u")
PyPlot.xlim(0, maximum(x_array))
PyPlot.ylim(-0.5,1.5)
```
![移流1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/772969/1e360388-7467-e8fe-2f69-4121617566ac.png)

forの表記はPythonより好き。でもendをよく忘れちゃう 。配列は1個ずらして…っと。

```julia
#u_ftcs = u_array.copy()
u_ftcs = copy(u_array)
# タイムステップを設定
#for n in range(Time_step):
for n in 1:Time_step
    #u_old = u_ftcs.copy()
    u_old = copy(u_ftcs)
    #u_ftcs[0] = u_old[0] - CFL / 2 * (u_old[1] - u_lower_boundary)
    u_ftcs[1] = u_old[1] - CFL / 2 * (u_old[2] - u_lower_boundary)
    #u_ftcs[-1] = u_old[-1] - CFL / 2 * (u_upper_boundary - u_old[-1])
    u_ftcs[end] = u_old[end] - CFL / 2 * (u_upper_boundary - u_old[end])
    
    #for j in range(1,Num_stencil_x-1, 1):
    for j in 2: Num_stencil_x-1
        u_ftcs[j] = u_old[j] - CFL / 2 * (u_old[j+1] - u_old[j-1])
    end
end

# plt.plot(x_array, exact_u_array, label="Answer")
# plt.plot(x_array, u_ftcs, label="FTCS")
# plt.legend(loc="upper right")
# plt.xlabel("x")
# plt.ylabel("u(FTCS)")
# plt.xlim(0, max(x_array))
# plt.ylim(-0.5,1.5)

PyPlot.plot(x_array, exact_u_array, label="Answer")
PyPlot.plot(x_array, u_ftcs, label="FTCS")
PyPlot.legend(loc="upper right")
PyPlot.xlabel("x")
PyPlot.ylabel("u(FTCS)")
PyPlot.xlim(0, maximum(x_array))
PyPlot.ylim(-0.5,1.5)
```
よしよし、元の記事とおんなじ画像。
![移流2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/772969/18d049cd-e644-b4da-105d-6e2b7f1b0315.png)

この調子でいけば次の記事のコードまでさっといけそう！

```julia
#u_upwind = u_array.copy()
u_upwind = copy(u_array)

# タイムステップを設定
#for n in range(Time_step):
for n in 1:Time_step
#     u_old = u_upwind.copy()
    u_old = copy(u_upwind)
#     u_upwind[0:1] = u_old[0] - CFL * (u_old[1] - u_lower_boundary)
    u_upwind[1:2] .= u_old[1] - CFL * (u_old[2] - u_lower_boundary)
#     for j in range(1,Num_stencil_x):
    for j in 2:Num_stencil_x-1
#         u_upwind[j:j+1] = u_old[j] - CFL * (u_old[j] - u_old[j-1])
        u_upwind[j:j+1] .= u_old[j] - CFL * (u_old[j] - u_old[j-1])
    end
end

# plt.plot(x_array, exact_u_array, label="Answer")
# plt.plot(x_array, u_upwind, label="Upwind differencing")
# plt.legend(loc="upper right")
# plt.xlabel("x")
# plt.ylabel("u(Upwind differencing)")
# plt.xlim(0, max(x_array))
# plt.ylim(-0.5,1.5)
PyPlot.plot(x_array, exact_u_array, label="Answer")
PyPlot.plot(x_array, u_upwind, label="Upwind differencing")
PyPlot.legend(loc="upper right")
PyPlot.xlabel("x")
PyPlot.ylabel("u(Upwind differencing)")
PyPlot.xlim(0, maximum(x_array))
PyPlot.ylim(-0.5,1.5)
```
![移流3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/772969/dba2c72a-1142-53f0-a0b8-cc8c956bfdaa.png)

（中略）

よし最後のCIPのだ！

```julia
# u_cip= u_array.copy()
u_cip= copy(u_array)

# partial_u_cip = ((np.append(u_cip[1:], u_upper_boundary) + u_cip)/2 - (np.append(u_lower_boundary, u_cip[:-1]) + u_cip)/2)/ Delta_x
partial_u_cip = ((vcat(u_cip[2:end], [u_upper_boundary]) .+ u_cip)./2 .- (vcat([u_lower_boundary],u_cip[1:end-1]) .+ u_cip)./2)./ Delta_x
##いいかきかたがわからない！

# 時間発展させる
# for n in range(Time_step):
for n in 1:Time_step
#     u_old = u_cip.copy()
#     partial_u_old = partial_u_cip.copy()
    u_old = copy(u_cip)
    partial_u_old = copy(partial_u_cip)
#     u_cip[0] = 0
#     partial_u_cip[0] = 0
    u_cip[1] = 0
    partial_u_cip[1] = 0

#     for j in range(1, Num_stencil_x):
    for j in 2 : Num_stencil_x-1
#         a = (partial_u_old[j] + partial_u_old[j-1]) / Delta_x**2 - 2.0 * (u_old[j] - u_old[j-1]) / Delta_x**3
        a = (partial_u_old[j] + partial_u_old[j-1]) / Delta_x^2 - 2.0 * (u_old[j] - u_old[j-1]) / Delta_x^3
#         b = 3 * (u_old[j-1] - u_cip[j]) / Delta_x**2 + (2.0*partial_u_old[j] + partial_u_old[j-1]) / Delta_x
        b = 3 * (u_old[j-1] - u_cip[j]) / Delta_x^2 + (2.0*partial_u_old[j] + partial_u_old[j-1]) / Delta_x
        c = partial_u_old[j]
        d = u_old[j]
        xi = - C * Delta_t  # C > 0
#         u_cip[j:j+1] = a * xi**3 + b * xi**2 + c * xi + d
        u_cip[j:j+1] .= a * xi^3 + b * xi^2 + c * xi + d
#         partial_u_cip[j:j+1] = 3 * a * xi**2 + 2 * b * xi + c
        partial_u_cip[j:j+1] .= 3 * a * xi^2 + 2 * b * xi + c
    end
end

# plt.plot(x_array, exact_u_array, label="Answer")
# plt.plot(x_array, u_cip, label="CIP method")
# plt.legend(loc="upper right")
# plt.xlabel("x")
# plt.ylabel("u(CIP)")
# plt.xlim(0, max(x_array))
# # plt.ylim(-0.5,1.5)
PyPlot.plot(x_array, exact_u_array, label="Answer")
PyPlot.plot(x_array, u_cip, label="CIP method")
PyPlot.legend(loc="upper right")
PyPlot.xlabel("x")
PyPlot.ylabel("u(CIP)")
PyPlot.xlim(0, maximum(x_array))
PyPlot.ylim(-0.5,1.5)
```
![移流5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/772969/f8d158c3-d523-6917-6236-71d068a1e3ea.png)

あれえナンデ発散！

```
In [94]:Delta_x
Out[94]:1.01
In [95]:maximum(x_array)
Out[95]:101
```
おまえか！？ いや、この差では発散は…ないやろ…。（pythonでは1.0,100です）
forの中で発散してそうなので中身をチェック。


pythonの配列で[j:j+1]って…jしか見てないのでは…。配列の違いは慣れないとなぁ。というわけで元のコードでは思うところあったのかj+1も書いてるけど実際はjしか参照してなかったのが、そのままコピペで律儀にj+1も更新してたのが間違いでした！ちゃんちゃん！無事元と同じ出力ができました！

```julia
#         u_cip[j:j+1] = a * xi**3 + b * xi**2 + c * xi + d
        u_cip[j] = a * xi^3 + b * xi^2 + c * xi + d 
#         partial_u_cip[j:j+1] = 3 * a * xi**2 + 2 * b * xi + c
        partial_u_cip[j] = 3 * a * xi^2 + 2 * b * xi + c 
```
![移流7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/772969/4f5afe7e-dfad-3280-95e1-1bd320abab2e.png)

気になる計測時間は、
Julia ([これ](#https://twitter.com/genkuroki/status/948582344696373249)みてBenchmarkToolsで測定)
- minimum time:     87.012 ms (0.00% GC)
- median time:      104.454 ms (0.00% GC)
- mean time:        136.777 ms (0.34% GC)
- maximum time:     358.496 ms (0.00% GC)
Python(%%timeで測定)
- 132ms

どう考えてもpyplotの読み込み時間しかかかってなさそうなので当然の結果ですな。[バーガース方程式の数値シミュレーション](https://qiita.com/KQTS/items/0c4f6c47a4d56881a178)はちょっと重かったのでそこでのお楽しみ。

コードは[jupyternotebookのファイルをそのままあげてます](https://github.com/myu65/flow_simulation/blob/main/移流方程式写経Julia版.ipynb)。

今日はここまで！次回はバーガース方程式を写経したいです。