流体力学 on Paraiso が かんせいした!
見よ!流体力学ソルバーのParaiso実装を!
cell2 <- proceedSingle 1 (dt/2) dR cell cell -- まず時間1次精度でdt/2だけ進めておいて、 cell3 <- proceedSingle 2 dt dR cell2 cell -- それを使ってもとの量をdtだけ進める。
ありのままの姿で書かかれた2次のルンゲ・クッタ法を!
interpolate order i cell = do let shifti n = shift $ compose (\j -> if i==j then n else 0) -- i軸方向、nマスの平行移動を定義 a0 <- mapM (bind . shifti ( 2)) cell -- 2マス平行移動を全ての変数に適用、 a1 <- mapM (bind . shifti ( 1)) cell -- 1マス平行移動を全ての変数に適用、 a2 <- mapM (bind . shifti ( 0)) cell -- 0マス平行移動を全ての変数に適用、 a3 <- mapM (bind . shifti (-1)) cell -- -1マス平行移動を全ての変数に適用、 intp <- sequence $ interpolateSingle order <$> a0 <*> a1 <*> a2 <*> a3 -- 各自由度ごとに補完関数を適用。
TraversableとApplicativeを使って一般的に書かれた空間補完を!
let timescale i = dR!i / (soundSpeed cell + abs (velocity cell !i)) -- 各マスごと、各方向ごと、にタイムステップの上限が決まるので、 dts <- bind $ foldl1 min $ compose timescale -- まず各マスの中で方向ごとのminを取る。 dtG <- bind $ cflG * reduce Reduce.Min dts -- さらに全マス上のグローバルなminを取る。 dt <- bind $ broadcast dtG -- 求まった最小のタイムステップをブロードキャスト。
抽象的なreduceとbroadcastで書かれたこれ実装の大変さや計算時間なんてまったく考えてないだろっていうタイムステップ上限の計算を!
lesta <- starState shockStar shockLeft left -- 左の元状態left から衝撃波後の状態lestaを rista <- starState shockStar shockRight right -- 右の元状態rightから衝撃波後の状態ristaを求め let selector a b c d = -- セレクタを定義し select (0 `lt` shockLeft) a $ select (0 `lt` shockStar) b $ select (0 `lt` shockRight) c d mapM bind $ selector <$> left <*> lesta <*> rista <*> right -- 4つのうちのどれかを選ぶ。
流体の全自由度をひとまとめに扱って(ここではleft, lesta, rista, right)記述された計算や分岐の数々を!
いや、正直同じ条件の分岐を複数生成するのはどうにかしたいのですが・・・
次元の数*1、流体の自由度の数、空間の補完、時間積分の次数、リーマンソルバーの実装などの部品がほぼ直交的に書けてしまう。ある部分を実装するときは関係ない部分は意識しないですむ、書いてしまえば他の部分で何を選ぼうと組み合わせられる、という感じなのでとても楽しいです。これに対して出てくるコードはいまのところ同じ計算を何度もするようなぶよぶよのコードなので、大いに改善の余地ありです。
実装について
たぶん以下のようにするとサンプルがうごくはず
git clone https://github.com/nushio3/Paraiso.git cd Paraiso/ cabal configure cabal install cd examples/Hydro/ make lib #環境にもよりますが make #数分かかるかも ./kh.out
Builder
ざっくり説明しますと、まずLanguage.Paraiso.OM.BuilderにいるBuilderモナドというやつが、いわゆるフロントエンド担当です。抽象構文木を構築します。こいつの正体は、構築中のデータフローグラフ+αを持ったStateモナドです。
type Builder vector gauge val = State (BuilderState vector gauge) val
Builderモナドはいくつかの予約語(load, store, imm, reduce, broadcast, shift, loadIndex)の入出力に現れるほか、Additive, Ring, Fieldなど型クラスのインスタンスになって各種の演算子を取り揃えています。Builderモナド同士を演算すると、
- 引数のBuilderモナドをまず実行してデータフローグラフの入り口となるノードを作って返させ、
- そこからエッジを出し、演算を行うノードを新しく追加し、
- 演算結果を表すノードを作って返す
という処理を行うモナドになるので、Buiderモナドで作った複雑な数式をrunするとあら不思議、その数式を計算するデータフローグラフがひとりでに組み上がるというしくみです。
このグラフはOrthotope-Machine(OM)という仮想マシンのレジスタを操作するもので、OMのレジスタは要素の型(Int, Float, Double, ...)と、配列変数なのかただの変数なのか、程度の区別があります。レジスタの表現として、Valueという、レジスタの型や大きさの情報を型として持っているデータ型と、DynValueという、型や大きさの情報を値レベルに落として持っているデータ型があります。Builderモナドが内部で組み立てるグラフの頂点は単一の型を持たないといけないのでDynValueを使います。ただしBuilderモナドが受け付けたり返したりするのは常にValueなので、Builderモナドを使う限り、うっかり型がおかしい演算を組み立てたりはしないようになっています。
Generator
一方、バックエンド担当はBinderモナドです。こいつはexposeすらされてない闇のモナドです。
type Binder v g a = State (BinderState v g) a
Binderも基本構造だけはシンプルで
addBinding cursor = do ... lhs <- leftHandSide cursor ... rhs <- rightHandSide preCursor bindingModify $ Map.insert cursor (lhs ++ " = " ++ rhs ++ ";") leftHandSide = cursorToSymbol LeftHand rightHandSide cur = do ... when (allocStrategy (getA node0) == Alloc.Delayed) $ addBinding cur cursorToSymbol RightHand cur
addBindingはleftHandSideにrightHandSideを代入するコードを生成し、rightHandSideは自分の値を計算してくれるようなaddBindingを呼び出し・・・という仕組みで再帰的にコードを生成します。そこかしこで直接文字列をいじってCのコードを吐いているので気持ち悪いです。文字列じゃなくてCの構文木を組み立てるような形にしたいもんですが、なんか良いライブラリないですかね。
最適化部
そんなものはない
一応必要なとこしかコード生成しないので、dead code eliminationだけは入ってると言えなくもないです。もとのアルゴリズムの並列性をそこなわない形でグラフにもつことには成功したわけで、こっからががんばりどころです!
*1:いま公開してる版は微妙に手を抜いて次元依存の部分がありますがその気になれば完全に消しちゃえるはず