MapReduceとN体問題

天体や分子動力学（MD）といったN体問題は計算量が多い問題としてよく知られている．近年流行しているHadoopのMapReduceの計算モデルがN体問題にも適用できると気づいたので，記録しておく．

N体問題とは，複数の質点が各々相互に作用して運動する系について解く問題である．この問題がやっかいなのは，位置と運動量が与えられた１個の質点に関して，次の瞬間の位置と運動量を求めるためには，他のN-1個の質点の位置を知る必要があることだ．１個の質点に関してN-1個の質点との作用を計算しなければならないので，全体ではN*(N-1)回の計算を行う必要がある．計算量の多い問題は，複数のマシンで並列に計算をしたいが，１ステップの計算に全ての質点の情報が必要なので，マシン間の通信が問題となる．こういったプログラミングをするのは非常に面倒だ．

最近にわかに注目を集めているHadoopはMapReduceという計算モデルに落とし込める計算に関しては，簡単に並列分散処理を行うプログラムを作ることができる．MapReduceは演算をMapとShuffleとReduceの３段階に分ける．Map，ReduceフェイズにもちいるMap関数とReduce関数はユーザーが独自にプログラミングできる．

Map関数は，KeyとValueのデータの組から，異なるKeyとValueの組のリストを作成する．

Map(<K1, V1>) -> list <K2, V2>

Shuffleのフェイズで同じKeyをもつ組が集められ，<K2, list V2 >の形になる．

Reduce関数で，Valueが一つの値になる．

Reduce(<K2, list V2>) –> <K2, V3>

この一連のMapReduceの流れをN体問題に置き換える．iが質点の番号，xi, piが質点iの位置と運動量である．

Map関数で，質点iの位置を全てのj = 0…N-1に対して割り当てる．

Map(<i, xi>)　-> list <j, xi> (0<j<N-1, j != i)

そうするとShuffleで全ての質点の位置が集まるので，Reduce関数で質点jの次のステップでの位置と運動量を求めることができる．

Reduce(<j, list xi>, xj, pj) –> <j, x’j, p’j>

Reduce関数の入力値には，当然質点jの現在の位置と運動量が必要とされる．

MapReduceジョブ１回が１ステップの計算となる．現実的には，HadoopではMap, ReduceでそれぞれN回JVMを立ち上げることになるので（JVMの再利用を行うならプロセス数分），1ステップで2N回JVMを起動して落とす．また，MapReduceの度にHDFSに書き込みが行われる．インメモリで計算すべきだが，HDFSのキャッシュ機能がどう働くかは分からない．MapReduceジョブはオーバーヘッドが大きいから，現実的ではないと思う．

移動一般曲線座標系での流体の支配方程式

これはもっと車輪の再発明．

PDFにしたのをこちらに置いておく．

なお，教科書であれば梶島 岳夫『乱流の数値シミュレーション』，藤井孝蔵『乱流の数値計算法』が詳しい．

移動一般曲線座標系での流体の支配方程式の導出

正規直交座標空間x_iで定義された領域V_x内での物理量φの保存式を考える．V_xは時間によって移動・変形するとする．f_jを流束，Fを湧き出しとする．φの保存式は次のように表される．

 移動一般曲線座標ξ^iへ座標変換をする．V_ξをV_xの計算空間への写像とすると，

となる．ここに，Jは座標変換のヤコビアンであり，

である．V_ξは時間によって変化しないとすると，

となる．

 (SCL)

を用いると，

 

と変形できる．これは強保存型と呼ばれる形である．

u_iを流速の正規直交座標成分（共変成分），u^G_iを計算格子の移動速度の正規直交座標成分とする．教保存型の保存式に

を代入すると，

となり，質量保存式となる．密度が一定のとき，うは以下のように書き換えられる．

流速の反変成分

を用いた場合には，

となる．

強保存型の保存式に

を代入すると，

となり，運動量保存式となる．σ_{ij}は応力テンソルであり，

である．pは圧力，μは粘性係数である．密度および粘性係数を一定とすると，上式は以下のように書き換えられる．

νは動粘性係数である．流速の反変成分U^k，格子歪みテンソル

，高さポテンシャルを含んだ圧力

を用いると，

勾配型の場合には，

となる．

同様に，スカラ保存式は，cを濃度，Dを拡散係数として，

となる．

なお，移動一般曲線座標での支配方程式はSCLおよび，

 (VCL)

が成立していれば，正規直交座標での質量保存式，運動量保存式，スカラ保存式

と等価となる．離散化された状態では，SCLおよびVCLは必ずしも成立するわけではないので注意が必要である．

水面波の計算における離散化の影響

2015-03-09 英語にしました。

---

車輪の再発明である．混相流の計算をしている人は経験的に知っていることだが，体系的に調べたものを見たことがなかったので，必要だから計算した．TeXで書いて，PDFにしたのがこちらにある．PDFの方が読みやすいと思う．

水面波の計算における離散化の影響

計算の安定性

差分法を用いて時間発展型の方程式を数値的に解くとき，陽解法の場合，安定した計算を行うためには，格子幅Δに対して，時間刻みΔtに上限があることが知られている．移流方程式でのCourant数uΔt/Δ<1（CFL条件）や，拡散方程式での拡散数DΔt/Δ^2<0.5などが有名である．uは流速，Dは拡散係数である．しかし，水面波の計算を行う際には，これらの条件を満たしていても計算が不安定となる場合があり，別の条件が必要とされる．そこで，微小振幅波の場合に，von Neumannの安定性解析法を適用し，計算の安定性について検証する．

考察する支配方程式はEuler方程式と連続の式である．

ここに，u=(u,v,w)，P，ρはぞれぞれ流体の速度，重力ポテンシャルを含んだ圧力，流体の密度である．

Eular方程式は非線形項を含むため，von Neumannの安定性解析法を適用できない．そこで速度が充分小さいとして非線形項を消去した以下の式について検討する．

上式より，

が成立する．

境界条件は，水深が十分深く，水位変動ηが微小とすると，

となる．g，γは，それぞれ重力加速度，表面張力を密度で割った値である．

空間x，y方向，時間tに関して離散化し，計算することを考える．z方向には離散化を行わない．x，y方向の格子間隔をそれぞれΔx，Δy，時間刻みをΔtとする．x=lΔx，y=mΔy，t=nΔtの差分点での，u，v，w，P，ηの値をそれぞれ，

とおく．また，

と表すこととする.

空間差分にレギュラ格子系，コロケート格子系，スタガード格子系を用いた場合について調べ，差分の方法による違いについて検討を行う．また，水位変動の時間発展にEuler陽解法，２次精度Adams-Bashforce法を用いた場合について調べ，時間発展の方法による違いについて検討する．

レギュラ格子系の場合

速度の離散化をレギュラ格子系にて行った場合には，離散化された支配方程式は，

となる．離散化された境界条件は，

となる．

u^n_{l,m}，v^n_{l,m}，w^n_{l,m}，P^n_{l,m}，η^n_{l,m}を以下のようにFourier級数展開し，波数ごとの安定性を調べていく．

とする．波数k_x，k_yのFourier係数の方程式は以下のようになる．

上からそれぞれ式(1),(2),(3),(4),(5),(6),(7)とする．

式(1),(2),(3),(4)より，

として，

境界条件(4)(5)より，

である．

式(3)(7)より，

となる．

とおくと，

となる．

α，βを

と定義すると，H^nに関する漸化式は以下のように書き換え
られる．

計算の安定条件はすべての波数k_x，k_yに対して|α|<1かつ
|β|<1となることである．

B^2-1<0の場合，α，βは複素数なので，

となり，計算は安定である．B^2-1>0の場合，B>1のとき，α>1，B<-1のとき，β< -1となるので計算は不安定である．したがって，計算が安定となるのはすべての波数k_x，k_yについてB^2-1<0が成立するときである．

Δを

と定義する．BはK_1=K_2=0のとき，最大値

をとり，K_1 = 2/Δ，K_2 = 1/(2Δ)のとき，最小値

をとる．

したがって，安定条件は

である．整理すると，

となる．

コロケート格子系の場合

コロケート格子系の場合は，連続の式の離散化の方法が変わる．レギュラ位置からスタガード位置に速度を補間したものを用いて連続の式を評価する．２次精度のコロケート格子系で離散化された連続の式は，

となり，等間隔直交格子の場合はレギュラ格子系と全く等価である．

したがって，計算の安定条件はレギュラ格子系の場合と同様，

となる．

スタガード格子系の場合

速度の離散化をスタガード格子系にて行った場合には，離散化された支配方程式は，

となる．境界条件はレギュラ格子系と同様となる．

Fourier係数の方程式は，

となる．

として，H^nに関する漸化式は

となる．

同様の議論により，計算の安定条件は

となる．

スタガード格子系の場合は，レギュラ格子系およびコロケート格子系の場合と比べ
て，計算が安定となる時間刻みΔtの限界が1/2倍となる．

レギュラ格子系およびコロケート格子系では，１回微分を差分近似する際，2Δx離れた点での差分を取る．一方で，スタガード格子系では，１回微分を差分近似する際，Δx離れた点での差分を取る．したがって，実質の解像度はスタガード格子系の方が1/2の細かさであり，時間刻みは1/2の大きさでなければならない．

Adams-Bashforce法の場合

水位変動の時間発展に以下のような２次精度のAdams-Bashforce法を用いた場合について考える．

H^nに関する漸化式は，レギュラ格子系およびコロケート格子系の場合は，

となる．

同様の議論により，安定条件はレギュラ格子系およびコロケート格子系の場合，

，スタガード格子系の場合，

となる．水位変動の時間発展に２次精度のAdams-Bashforce法を用いた場合の条件は，Eular陽解法を用いた場合と同じである．

根号がついた項は波長πΔの水面波の伝播速度を表しており，この条件は波長πΔの水面波に対するCFL条件と解釈できる．最小の波長は2Δであり，最小の波長の伝播を解像できなくとも，計算が安定になることがわかる．この理由は，後述するように，離散化の影響により，水面波の伝播速度が遅れるからと考えられる．

以上の結果は水面波の計算の際に，時間刻みΔtを決定する上で重要な目安となる．

なお，粘性がある場合でも，波の伝播速度はわずかに小さくなるだけであるので，粘性の影響は考慮する必要はない．

位相誤差

この項では，数値計算で水面波を計算した際に，波の位相が差分近似によってどのような影響を受けるか検討する．前節で定義されたαを，

とおくと，

である．計算の安定条件は，すべての波数に対して，θ_αが存在するために条件でもあ
る．θ_αは時間を1ステップ計算を進めた際に波数kの波に加わる位相である．角振動数σは単位時間あたりに変化する位相なので，

となる．βの場合も同様である．

離散化されていない時の角振動数

と比較すると，低波数の領域では角振動数は互いに一致しているが，高波数では離散化によって位相が遅れ，角振動数が小さくなる．高波数の波の挙動が重要となる計算では，位相誤差が計算結果に影響を及ぼす可能性がある．