流体計算の実行速度比較を再訪してみた: C++, Julia

[流体計算の実行速度比較: Fortran, C++, Rust, Python, Julia](https://qiita.com/shigunodo/items/d693dc03323f9a205bb9)の記事は実アプリケーション的なコードにおける速度比較を行っており、とても有用でした。
この記事の中で、JuliaがC++やFortranよりも遅いという結論が出ていました。自分の肌感覚としてはこのようなシンプルなコードの場合CとJuliaは速度差があまり出ないと思っていたので、意外でした。
そこで、こちらの記事にある計算を、1からJuliaで書いてみました。

こちらの記事では様々な言語での速度比較をしていますが、コードはC++のもののみが公開されています。
そこで、このC++のコードをそのままJuliaに移植してみることで、どちらがどのくらいの速さかを見てみようと思います。

記事では400x400というメッシュの計算結果を載せていましたが、これは計算時間が非常にかかるようですので、githubのC++のコードにあった100x100のバージョンを実行してみようと思います。

# 環境
- MacBook Pro 13-inch, M1 2020
- C++のバージョン

```
g++ --version
Configured with: --prefix=/Library/Developer/CommandLineTools/usr --with-gxx-include-dir=/Library/Developer/CommandLineTools/SDKs/MacOSX12.1.sdk/usr/include/c++/4.2.1
Apple clang version 12.0.5 (clang-1205.0.22.11)
Target: arm64-apple-darwin20.5.0
Thread model: posix
InstalledDir: /Library/Developer/CommandLineTools/usr/bin
```
です。

- Juliaのバージョン
Juliaは

```
Julia Version 1.8.5
Commit 17cfb8e65ea (2023-01-08 06:45 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.5.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores
Environment:
  JULIA_SSL_NO_VERIFY_HOSTS = github.com
  LD_LIBRARY_PATH = :/usr/local/lib/
  DYLD_LIBRARY_PATH = /opt/intel_org/oneapi/mkl/2021.3.0/lib
```
です。


# C++
まず、同じ結果が返ってきていることを確認するために、計算で生成された最後のファイル```b0000100.dat```の一番最後の値をみてみます。

```
g++ main.cpp -std=c++17   
```
でコンパイルした場合

最後から10番目は

```
2.016334709006551762e+01
2.016484612166407686e+01
2.016605065164424460e+01
2.016670718736106593e+01
2.016298471334668108e+01
2.017677248704771387e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
```
です。

結果の出力は

```
elapsed:   0 h 01 m 51 s | tstep =    97 | t =      2.911 | iter =      28 | rest:   0 h 00 m 03 s 
elapsed:   0 h 01 m 52 s | tstep =    98 | t =      2.941 | iter =      28 | rest:   0 h 00 m 02 s 
elapsed:   0 h 01 m 54 s | tstep =    99 | t =      2.971 | iter =      28 | rest:   0 h 00 m 01 s 
elapsed:   0 h 01 m 55 s | tstep =   100 | t =      3.001 | iter =      28 | rest:   0 h 00 m 00 s 
```
となりました。1分55秒ほどかかっているようです。

次に、
```
g++ main.cpp -std=c++17 -O3 
```
にすると、

```
elapsed:   0 h 00 m 18 s | tstep =    95 | t =       2.85 | iter =      27 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 18 s | tstep =    96 | t =       2.88 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 19 s | tstep =    97 | t =      2.911 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 19 s | tstep =    98 | t =      2.941 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 19 s | tstep =    99 | t =      2.971 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 19 s | tstep =   100 | t =      3.001 | iter =      28 | rest:   0 h 00 m 00 s 
```
19秒で終わりました。


アウトプットファイルの末尾は
```
2.016334709006551762e+01
2.016484612166407686e+01
2.016605065164424460e+01
2.016670718736106593e+01
2.016298471334668108e+01
2.017677248704771387e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
```
結果は変わっていません。

```
g++ main.cpp -std=c++17 -O2
```
の場合は、

```
elapsed:   0 h 00 m 19 s | tstep =    96 | t =       2.88 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 20 s | tstep =    97 | t =      2.911 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 20 s | tstep =    98 | t =      2.941 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 20 s | tstep =    99 | t =      2.971 | iter =      28 | rest:   0 h 00 m 00 s 
elapsed:   0 h 00 m 20 s | tstep =   100 | t =      3.001 | iter =      28 | rest:   0 h 00 m 00 s 
```
で20秒ほどかかっています。
出力は

```
2.016334709006551762e+01
2.016484612166407686e+01
2.016605065164424460e+01
2.016670718736106593e+01
2.016298471334668108e+01
2.017677248704771387e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
2.009126984126984539e+01
```
で同じです。

最適化オプションによってずいぶん速度が変わるようです。

# Juliaの場合
C++のコードを見ながら、写経するようにJuliaに変換してみました。クラスの構造などはJuliaの型を使うような形でなるべく維持するようにしました。コードは[こちら](https://github.com/cometscome/Juliafulidtest/tree/master/juliafluid2d)にあります。Juliaっぽい書き方をしておらず、C++と同じ結果が出るようにコードを書いたためにC++に似ている形になっています。


出力結果は
```
tstep = 99, iter = 28
tstep = 100, iter = 28
 40.251492 seconds (15.11 M allocations: 2.073 GiB, 0.42% gc time, 1.24% compilation time)
```
で、

アウトプットファイルの末尾は
```
20.163347090065518
20.164846121664077
20.166050651644245
20.166707187361066
20.16298471334668
20.176772487047714
20.091269841269845
20.091269841269845
20.091269841269845
20.091269841269845
```
です。答えがC++と完全に一致していることがわかりますから、コード自体の移植には成功しています。
速度は40秒ほどでした。C++のデフォルトは約120秒ですからそれよりは速いです。一方、O2の最適化のC++コードと比べると、2倍遅い、という結果でした。2倍なら悪くないように思います。


## さらなる改良
次に、inboundsをつけてみました。
その結果は

```
tstep = 96, iter = 28
tstep = 97, iter = 28
tstep = 98, iter = 28
tstep = 99, iter = 28
tstep = 100, iter = 28
 30.240109 seconds (15.23 M allocations: 2.077 GiB, 0.53% gc time, 0.97% compilation time)
```

```
20.163347090065518
20.164846121664077
20.166050651644245
20.166707187361066
20.16298471334668
20.176772487047714
20.091269841269845
20.091269841269845
20.091269841269845
20.091269841269845
```
です。src_inboundsというディレクトリに入っている方がinboundsが入っている方です。これで、c++のO2と比べて1.5倍遅い、くらいになりました。

あとは何をすれば速くなりますかね。もう少しいじれば1倍くらいまで行きそうな気はしているのですが。