熱計算によるC++とJuliaの速度比較

はじめに

現在，私は研究で使用するため，熱のシミュレーションをC++を使って開発しています．しかし，C++は高速であるが，それ以外の利点は少ないと感じています．コードは冗長的で，デバッグも煩雑であるため，開発とは別のところで時間が奪われてしまうことが時々あります．またC++は描画ツールが乏しくカバーするためにPythonを使用しており，開発言語が2つあり，手間がかかると感じています．そんな時，先輩がJuliaは速く，書きやすいとJuliaの良さを熱弁していたため，どれくらいの速度がでるのかC++とJuliaで比較してみようと思います．

実行環境

• OS: Ubuntu 24.04 (WSL2)
• CPU: 16 × AMD Ryzen 7 5700X 8-Core Processor
• RAM: 15 GiB

Juliaについて

Juliaについて簡単に説明します．
Juliaは，インタプリタ言語とコンパイル言語の特長を兼ね備えた，高性能なプログラミング言語です．動的型付けをサポートしつつ，コンパイル言語に匹敵する高速な実行速度を実現しています．この高速性を支える仕組みとして，JuliaはJust-In-Time (JIT) コンパイルを採用しています．JITコンパイルでは，コードが実行時にコンパイルされ，最適化されます．その結果，初回の実行には若干の時間がかかるものの，コンパイル後のキャッシュが利用されるため，2回目以降の実行は非常に高速です．この最適化により，JuliaはCやFortranといった従来のコンパイル言語に匹敵する処理速度を持つとされています．

コード

以下のレポジトリにあります．

比較の仕方

比較のために解く方程式は次の2次元の熱方程式とします．

$$
\frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
$$

$u(x, y, t)$は温度分布を表し、$\alpha$は熱伝導率を示します．

以下に偏微分方程式を離散化する方法について説明します．

離散化の方法

時間の離散化:
オイラー法 (Euler Method) を使用する．これは，時間ステップ $\Delta t$ ごとに温度分布を更新する方法です．

$$
u^{n+1} = u^n + \Delta t \cdot f(u^n)
$$

$u^n$ は時間ステップ $n$ における温度分布， $f(u^n)$は熱方程式の右辺です．

空間の離散化:
中心差分法 (Central Difference Method) を使用して空間の2次微分 $\frac{\partial^2 u}{\partial x^2}$ と $\frac{\partial^2 u}{\partial y^2}$を離散化します．

$$
\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \
\frac{\partial^2 u}{\partial y^2} \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\Delta y)^2}
$$

つまり，熱方程式の右辺である$f(u^n)$は次のようになります．

$$
f(u^n) = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\Delta y)^2}
$$

シミュレーションの設定

初期条件の設定
初期条件の温度は中心からの距離に応じて温度が指数関数的に減少するような分布を作成しています．
$r_0$は半径（75），$(x_c, y_c)$は空間の中心の座標，$(x,y)$は任意のセルの中心点を示します．

$$
u(x, y, 0) = 10 \cdot \exp\left(-\frac{\sqrt{(x - x_c)^2 + (y - y_c)^2}}{r_0}\right)
$$

境界条件の設定
境界の温度は0としています

メッシュの設定
空間を150×150と定義し，150×150のメッシュに分割して計算しています．

時間
シミュレーション時間は100秒行い，時間ステップ $\Delta t$は0.1で計算します．

比較方法

計測の仕方

計測はシミュレータを実行した箇所を測定しています．つまり初期化やファイル保存部分は測定外とします．計測は3回測定し，その平均を比較します．ただし，JuliaはJITコンパイルによって実行時間が長くなるため初回実行時ではなく，2回目以降の実行時間を使用します．
またC++の最適化フラグの有無についても調査してみます．最適化フラグは-O3 -march=native -fltoを使用します．（@fujitanozomuさん，コメントありがとうございます．）
以下はC++とJuliaの測定箇所です．

C++の測定箇所

  auto start_time = std::chrono::high_resolution_clock::now();
  for (int t = 0; t < num_time_steps; ++t) {
    solve_heat_equation(temperature, time_step, alpha, deltaX, deltaY);
  //  結果の保存関数
  //  save_to_file(temperature, t);
  }
  auto end_time = std::chrono::high_resolution_clock::now();

  std::chrono::duration<double> execution_time = end_time - start_time;
  std::cout << "Total execution time: " << execution_time.count() << " seconds"
            << std::endl;

Juliaの測定箇所

	start_time = now()
	solution = solve(problem, Euler(), dt = time_step, saveat = time_step)
  end_time = now()  

  println("Total execution time: $(end_time - start_time)")

結果

回数	C++[s] 最適化なし	Julia[s]	C++[s] 最適化あり
1	0.606	0.205	0.033
2	0.680	0.577	0.033
3	0.635	0.176	0.033
平均	0.640	0.319	0.033

コンパイルで最適化なしの場合，C++の実行時間平均は0.640秒で，Juliaと比較して，約2.0倍低速とJuliaのほうが高速でした．しかし，コンパイルで最適化ありだと，C++の実行時間平均は 0.033秒 で，Juliaと比較して約 9.7倍高速 になりました．

結論

C++のコンパイルが最適化される前はJuliaが約2.0倍高速でしたが，コンパイルを最適化すると，Juliaよりも9.7倍高速になり，C++がめちゃくちゃ高速であるとわかりました．ただ，それらしく記述するだけでJuliaも十分高速でしたので，使う機会があれば，使ってみたいと思います．これからもC++からは逃れそうにないです．今後もC++で頑張っていきます．
もしさらなる高速化や改善案があれば、ぜひ教えてください。

追記

記事投稿後，プログラミング言語の速度比較サイトってよくよく考えたら普通にありそうと思い，探してみたらあったので，備忘録のために残しておきます．

• マンデルブロ集合を生成するコードを使用した速度比較
  https://benchmarksgame-team.pages.debian.net/benchmarksgame/performance/mandelbrot.html
• 適応的台形法による数値積分アルゴリズムを用いた速度比較
  https://github.com/astrojhgu/adaptrapezoid_benchmark