Rust・Gnuplot・ffmpeg で Python・matplotlib よりは速くなった: 1d CFD, バーガース方程式

結論: GitHub への参照

先日この記事で最低のゴミ・産業廃棄物と書いた, 肖鋒, 長崎孝夫, 『数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析』 2020/1/9 の試行錯誤の結果をまとめておく. はじめに書いておくと, GitHub にあげてあるので, コードを見たい人はそちらを見てほしい.

あれほどゴミのような内容で 3600 円などいい商売だと思いつつ, お金がない中で買ったそれなりに値段の張る本なので, 何だかんだで頑張って C のコードを書き直しながら読んでいる.

自分用メモとして, Python を捨てて Rust・Gnuplot・ffmpeg にいたる経緯はあとで書くとして, まずはコードの話をしよう.

png 作成: Gnuplot

次の crate を使った.

これを見てもらえれば全ては終わる. コード全体は GitHub を見てもらうことにして, ここでは該当部分を抜き出しておこう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// cargo-deps: chrono, gnuplot
extern crate gnuplot;
use gnuplot::*;

pub fn write_png(cnf: &Config, cdata: &CalcData) -> Result<(), Box<dyn error::Error>> {
    let file_name: String = format!("{}/img.{:08}.png", &cnf.dir_name, &cdata.n);
    let mut fg = Figure::new();
    fg.axes2d()
        .set_title(&cnf.title, &[])
        .set_legend(Graph(1.0), Graph(1.0), &[], &[])
        .set_y_range(Fix(cnf.graph_ylim_min), Fix(cnf.graph_ylim_max))
        .set_x_label("x", &[])
        .set_y_label("u", &[])
        .lines(
            &cdata.xs,
            &cdata.us,
            &[Caption("numerical"), LineWidth(2.5)],
        )
        .lines(&cdata.xs, &cdata.ues, &[Caption("exact"), LineWidth(2.5)]);
    fg.save_to_png(file_name, 800, 800);

    Ok(())
}

他にもいろいろインポートは必要だろうが省略した. エラー処理もすさまじくひどいが, よくわかっていない. とりあえず最低限を書けるようになるまで読んだ公式のチュートリアルで出てきた記述を適当にコピペで使っている. (これこそ私が批判した上掲書のゴミコードスタイルだが, 全て棚に上げている.)

いま 1 次元のバーガースを解いているので, 2 次元プロットしたいというだけ. 計算法の都合によって途中で最大値がどんどん落ちていき, 動画化したときに y 軸方向にぶれて見づらいので, y_range を設定してある.

動画化: ffmpeg

これも FFI があるようだが, かったるいので直接コマンド実行している.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
use std::process::Command;
fn write_mp4(cnf: &Config, cdata: &CalcData) -> Result<(), Box<dyn error::Error>> {
    let options = ["-r", "10", "-i", &cnf.specify_png, &cnf.movie_name];

    let run_ffmpeg = Command::new("ffmpeg")
        .args(&options)
        .output()
        .expect("failed to start `ffmpeg`");
    println!("{}", String::from_utf8_lossy(&run_ffmpeg.stdout));

    Ok(())
}

本質的なエラー処理をつけていないのに Result を使っていて頭が悪いが, IO 系なのでつけるだけつけた. もちろんコマンド実行でぐぐればすぐに出てくるし, 多少プログラムが書ける人ならすぐわかるレベルの話だろう.

経緯

改めて書いておくと, 世界を汚染する害悪こと肖鋒, 長崎孝夫, 『数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析』 2020/1/9 では, 計算を実行しながら gnuplot で逐次画像出力の形で「動画化」されていて, ファイル出力としては eps だった. postscript は何もわからないし調べるのも面倒, 拡張子 cpp で実際には C++ でさえなく, 本の指示通りに Visual C++ でやるとコンパイルが通らない C 言語のコードのデバッグはやりたくなかった. (この C コードはネットから落とせるのだが, 本の購入者しか開けないようにパスワードがかかっているので, 買わないと地獄は感じ取れない.)

これまでいろいろな都合によって基本的にちょろいスクリプトには Python を使っていて, 数値計算系のコンテンツ作成も Python で慣れていたので, はじめは Python を使っていた. ちなみにその勉強の成果として Python を使った数学・応用数学入門については次のようなミニ通信講座を作っている.

C からの Python 移植で本を読んでいこう, となり, Python で数値計算となると numpy なので numpy でがんばって書き直していた. numpy となると for ではなくスライスを使った記法で書くのがふつうなわけで, それで書き直していたが, for (i = 1; i <= i_max + 1; i++) { }, for (i = i_max / 2 - 10; i <= i_max / 2 + 10; i++) { } みたいなのを適切なスライスで書き直すのがあまりに面倒で, 嫌になった. これはもう速い言語で for を直接移植した方が楽で早いことに気付き, 前から勉強しようと思っていた Rust に移行した. これならほとんど何も考えずに移植できて, 余計なバグ取りに勤しまずに済んだ.

基本的な計算は問題なくなったので問題は可視化の部分になった. 元の本は Gnuplot だったが, 計算と出力を分離していない産業廃棄物である. (工学者が産廃を作ることに抵抗を感じなかったのか不思議でならない.) Rust では CSV を書き出すことにして, 書き慣れていることから可視化は Python・matplotlib で処理して mp4 化することにした. ただし, これがとにかく遅い. 本では「2 秒分」計算しているのだが, これがとにかく遅い. 「0.1 秒分」の計算くらいならもちろん速くなるが, どんどんなまっていったり, 数値拡散を見るためにある程度長時間挙動を見たいので, これでは困る. いい加減嫌になったので, 自分にできる範囲の高速化としてできることを考えた.

はじめは Rust では CSV を作り, 別途 Gnuplot スクリプトで png 化して, これを ffmpeg で mp4 にしようかと思った. png さえ作れれば動画化は ffmpeg でかなり楽になる. 以前使ったときもここ自体はかなり速かった肌感覚があるので, png 化を楽に速くすることを考えた.

しかし, Gnuplot のスクリプト自体の情報が調べにくく, Gnuplot スクリプトとして動的に大量のファイルを読み込んで処理, みたいなこともしにくいようで, シェルスクリプトをかませないと駄目そうなど, プログラミング弱者の市民には本当に鬱陶しい. はじめは mp4 化と同じく直接のコマンド実行を検討したが, 少し調べたら gnuplot ラッパーがあったので, それを使うことにした.

ここまで来たらあとは ffmpeg をやればいい. 別途実行するのも面倒なので, Rust スクリプト内で直接実行することにした.

もう少し経緯, そして今後の予定

実ははじめは Paraview を使おうと思っていた. OpenFOAM で使っていて, 多少なりとも操作への感覚があったからだ. ただ Rust からの出力として vtk なり vtu を作らなければならず, それで二の足を踏んでいた. しかし Paraview は csv 読み込みもできると聞き, それなら, と思って試してみたが, 1 次元だと使い物にならなかった. いま勉強している純粋な 1 次元の数値計算の可視化に 3 次元用のビューアーを使おうとすること自体に無理があるので, その線は諦めた.

3 次元の計算をするときは多少は慣れているし Paraview も考えたい. ただ, csv からだと読み込みが少し面倒そうなので, やはり vtu を作るのがいいように思う. vtu の仕様の読み込みからしないといけなくて面倒で嫌なのだが ライブラリがあると嬉しい.

それはそれとして, 上掲書は後半で 2 次元を議論する. 2 次元もいまの 2d plot を使うか, 見たいモノによって 3d plot にすればいいので, しばらくは Paraview はいらないだろう.

今の本を読み終わったら粒子法または格子ボルツマン法を勉強する予定で, 既に本も買っている. ここの可視化をどうするといいのかが気になっている. 粒子法については SPHysicsPySPH があるので, これの可視化を参考にしようと思っている. まだ全く見ていないが, 何で可視化しているのだろうか?

格子ボルツマン法はまだ何もわかっていないのだが, 何で粒子法をやりたいかというと, 有限要素法などでメッシュを切りたくないからだ. ちょっと複雑な形状になるとメッシュ生成があまりにも面倒で, メッシュデータの読み込み処理なども書きたくない. 粒子法は粒子法でいろいろ鬱陶しい点があるはずだが, CG 界隈への展開もあるようなので, 別の発展方向が転がっているような気がしていて, それも 1 つのポイントだ.

いま, 特に理工系を志す中高生向けの数学・物理・プログラミングコンテンツを作っているし, たまたま地元の広義知人が区議になったこともあり, その人に地元の理工系教育推進のための相談・提案的なこともしている. 子どもの頃にほしかった・やってみたかったことでもあり, (現状専門外もはなはだしいが) まだまだこの辺の技術・研究ニーズもあるはずで, 中高生へのキャリア教育・技能のマネタイズ教育みたいな視点からも大事そうだと思っている. Pixar の展示会などもあり, ゲームやアニメ・映画など中高生にもある程度身近で需要も見やすいネタなので, いろいろ考えている. もちろん理工系大学生にも一定のニーズがある話題だろう.

私は出自的に, 数学と物理に関して, 能力はどもかく人並以上の尋常ならざる耐性はあるが, 数値計算はまた別種の厳しい世界がある. 何より工学者はゴミのようなコードしか書かないという話しか聞かないし, 数学系だと数値関数解析のような収束系の理論が中心のようでコードがあまり出てこないし, 物理系だと「具体的なコードを示すと時代遅れになりがちだし, 言語の選択については人の趣味もある」といってこれまた具体的なコードが出て来にくいし, いま流行りの統計学・機械学習まわりだとその筋の人達が興味を持つような話ばかりで, 私のような数学・物理系市民には不満なことこの上ないと, とにかく勉強がしづらい. 自ら血を流しつつ, きちんと将来的なマネタイズまで考えられるような, 生徒・学生向けコンテンツの整備を進めていきたい.

また, Python は Google Colab なり何なりで, これまたとにかく厳しい作業である環境構築なしでプログラミングに集中しやすい状況も出ていて, Python 利用の検討はこれはこれで外せないとは思っている. ただ, いまの私の力量ではとにかく遅く書くことしかできず, つらい. 数値計算, ちょっと何かやろうと思うとあっという間に重くなるので, そういう意味でもつらい. 地元的な事情からするとマシンを子ども達 (または家庭) 自身で準備するとなると, スペックも厳しいので, その意味でもいろいろ考えることはある.

理想郷は遥か遠い.

ぜひ登録またはツッコミを: GitHub・YouTube

最近の勉強の結果は次の GitHub リポジトリと YouTube リストにまとめている.

励みになるのでぜひチャンネル登録してほしい.

あと Rust 自体まだはじめて 1 ヶ月すら経っておらず, もとがゴミ以外表現しようのない工学者による C コードを基礎にした Rust コードなので, いまのコードもゴミだろう. 必要ならコード自体にもツッコミを入れてほしい. 泣きながら勉強して直していきたい.