Transcript 高速演算技術の開発
土木分野からの GPGPUへの期待と不安 浅井光輝＠九州大学 1 流体衝撃と防災 山間部での土石流 沿岸域での津波・高浪 最終的に評価したいのは，流体的な物体から受ける衝撃力を見積もり，安全・安心な構造物を． Motivation for water induced impact ② <During a Storm surge(Typhoon Songa 2004) in Hokkaido> wave-dissipating blocks 新たなエネルギーインフラ Wind lens proposed by Prof. Ohya in kyusyu-u. Solar panel Floating structure may be constructed by carbon-fiber reinforced concrete. proposed by Prof. Ohta in kyusyu-u. 4 大規模な数値解析実現に向けて ・３D-CAD ・数値地図・GISとの連携 入出力 ・高速演算技術の開発 ・解析手法の開発 ・３D立体視 ｿﾌﾄ面 対策 ﾊｰﾄﾞ面 対策 ・並列計算（スパコン使用） ・GPGPUの使用 5 大規模な数値解析実現に向けて ・３D-CAD ・数値地図・GISとの連携 入出力 ・高速演算技術の開発 ・解析手法の開発 ・３D立体視 ｿﾌﾄ面 対策 ﾊｰﾄﾞ面 対策 ・並列計算（スパコン使用） ・GPGPUの使用 6 KS-MORによる動的FEMの効率化 Background of research Dynamic analysis in Civil Engineering Bridge Building 3D-CAD Simplify Simplify beam and truss element Overall response Local stress CPU time mass-spring ‘direct’ FEM 橋梁の支承構造の詳細までモデル化 ゴム支承モデル 地震動を上部構造へ逃がすことにより 下部構造の変形を抑制する ゴムと鋼材の剛性比＝10MPa：200GPa ヒンジ支承モデル 橋軸方向の回転を許容、 橋軸直角方向については 変形、回転を拘束する 有限要素メッシュ図 40～50万要素 解析結果 ゴム支承 ヒンジ 1000ステップ2週間程度の計算→2日程度 RESULT: KS-MOR VON MISSES STRESS BASIS=10 BASIS=20 FEM STEPS 10 2.0 ×10-2 1.3 ×10-2 6.8 ×10-3 2.0 ×10-4 MPa Max Stress=3.7021000 STEPS Max Stress= 3.7019999 Max stress= 3.7025001 Max Stress= 3.3640000 Max Stress=3.3636000 Max stress=3.3656001 20 2.0 ×10-2 1.3 ×10-2 6.8 ×10-3 2.0 ×10-4 MPa STEPS 30 2.0 ×10-2 1.3 ×10-2 6.8 ×10-3 2.0 ×10-4 MPa Max Stress= 3.0415999 Max Stress=3.0409998 Max stress= 3.0439999 大規模な数値解析実現に向けて ・３D-CAD ・数値地図・GISとの連携 入出力 ・高速演算技術の開発 ・解析手法の開発 ・３D立体視 ｿﾌﾄ面 対策 ﾊｰﾄﾞ面 対策 ・並列計算（スパコン使用） ・GPGPUの使用 13 Developing tool for a FSI problem SPH FEM 流体問題 構造問題 Fluid (SPH) 粒子法 dynamics Dynamic change Start from 2008 Structure (FEM) Analysis 有限要素法 Slow change Hybrid simulator Nonlinear Structure Analysis due to water induced impact force Strong background including nonlinear material model 粒子法による関数近似の考え方 Physical quantities at target particle Target body Discritize into particles f xi f ( x j )W ( xi x j , h)dx j h: smoothing distance j ：neighbors i ：target particle Distance Weight function （kernel function） Physical quantities at neighbors Cubic-spline ISPH法の主な計算の流れ ①近傍粒子探索 ②粒子ごとに仮の速度更新 ③圧力ポアソン方程式の解法 SMAC法と同様 ④粒子の加速度・速度・位置の更新 16 MPS法（非圧縮性流体用SPH） n 1 i u 1 n1 u t p 2uin n i f ① predictor step ui uin t 2uin f ① predictor ② pressure Predict the current velocity without pressure gradient term ② pressure evaluation p 2 ③ corrector n 1 t u Evaluate “current” pressure Need to solve the Poisson’s eqn. ③ corrector step 17 Continuity cond. 1 u n1 u p n1 t SPH解析例 breaking water column Momochi, Fukuoka breakwater highway baseball stadium The geography is modeled from digital 1.4 million particles Map. Coupler between SPH and FEM Boundary particle on the edge of FEM 1 p① ( p1 p2 p3 ) / 3 2 ① 3 FEM 4 p② ( p4 p5 p6 ) / 3 5 ② 6 7 p③ ( p7 p8 p9 ) / 3 8 ③ 9 19 Time stepping scheme for each tool Courant condition e t c = length of element elastic wave speed 2 cp tsolid ≪ tliquid Explicit integration Explicit integration Implicit integration Explicit integration tsolid = tliquid 20 1 1 2 1 E 21 E Summary of 1way coupling analysis Pressure distribution on the boundary F EM 構造 S PH 流体 Fluid 問題 問題 Structure (FEM) (SPH) Analysis 有限要素法 粒子法 dynamics × Deformation ・I-SPH(using SMAC) ・Ghost particle ・Semi-implicit anal. 21 ・Liner elasticity ・Small deformation ・Implicit analysis Hybrid simulator 大規模な数値解析実現に向けて ・３D-CAD ・数値地図・GISとの連携 入出力 ・高速演算技術の開発 ・解析手法の開発 ・３D立体視 ｿﾌﾄ面 対策 ﾊｰﾄﾞ面 対策 ・並列計算（スパコン使用） ・GPGPUの使用 22 GPGPU計算の高速化と注意点 永井 学志 岐阜大学 工学部 数理デザイン工学科 23 講演内容 1. 自己紹介 2. 計算力学をGPGPUで高速計算する前に ～失敗しないためこれだけは理解して下さい～ • 電子計算機の概要 • CPU bound? もしくは memory bound? 3. GPGPU適用の具体例（研究紹介） ～連立1次方程式の反復求解～ • 一般的な疎行列の場合（実施済み） • 特殊化：イメージベース有限要素法（実施中） 24 講演内容 1. 自己紹介 2. 計算力学をGPGPUで高速計算する前に ～失敗しないためこれだけは理解して下さい～ • 電子計算機の概要 • CPU bound? もしくは memory bound? 3. GPGPU適用の具体例（研究紹介） ～連立1次方程式の反復求解～ • 一般的な疎行列の場合（実施済み） • 特殊化：イメージベース有限要素法（実施中） 25 FLOPS変遷（チャンピオンデータ） semi-log図 kflops Mflops Gflops 103 106 109 • 最新PC: 1010～ 1011 • 数十万円 Tflops 1012 • GPUが一足先に到達！ • 数万～数十万円 Pflops 1015 • 次世代スパコン: 1016 • 1000億円 26 時代は逐次実行からすっかり並列実行へ 電子計算機のクロック数 • 1秒あたりの電圧の切り替え回数 CPUがGHz（109Hz）に留まる理由 • 緑色の波長550nm * 550THz （550 x1012Hz） • PC（CPU，memory）サイズ：0.1mオーダ * 一昔前のケータイ（GHz）のダイポールアンテナ * 0.1mオーダ 縦(クロックup)がダメなら横(コア数)に • メイニーコア by Intel（MPI，OpenMP） • GPU（OpenCL） • 並列化効率 * 1パーセントの逐次実行があれば，x100倍が限界 27 CPUとGPUの比較（昨年度） CPU GPU 数個のコア 数百個のコア メモリ帯域：数十GB/s 百数GB/s 理論値：数十GFLOPS 1TFLOPS↑（単精度） • 100GFLOPS程度（倍精度） 例 メモリ帯域 CPU(DDR3-1066) 約6倍！ GPU(GTX285) 0 0G 20G 40G CPU(Core i7) 60G 80G 100G 120G 140G 160G 180G 倍精度 FLOPS 約20倍！ GPU(GTX285) 0 0G 200G 400G 600G 800G 1000G 単精度 1200G PCを介して使う GPU諸元 数100GFLOPS （倍精度計算） 数100コア 演算器群 演算器＆ 演算器群 演算器群 レジスタ+α 演算器＆ 演算器群 レジスタ キャッシュ キャッシュ (テクスチャメモリ) 160GByte/sec 数10GFLOPS 数コア 25GByte/sec 主メモリ 8GByte/sec （グローバルメモリ） 4万円GPU (NVIDIA GTX285) 10万円PC (Intel Core i7 920) 主メモリ I/O 講演内容 1. 自己紹介 2. 計算力学をGPGPUで高速計算する前に ～失敗しないためこれだけは理解して下さい～ • 電子計算機の概要 • CPU bound? もしくは memory bound? 3. GPGPU適用の具体例（研究紹介） ～連立1次方程式の反復求解～ • 一般的な疎行列の場合 30 疎行列の連立1次方程式の反復求解 Solve Ax=b 科学技術計算で最終的によく出てくる • 汎用的 GPUでPCよりも数倍高速化できる • ∵ Memory boundの典型 • NVIDIAも実装法の論文を発表（2008年末） • 研究としての新規性はない 31 有限要素解析の特徴 連立1次方程式 Ax bの求解がメイン 自由度数 2次元問題 自由度小 3次元問題 自由度大 解法 直接解法 （LU分解等） 反復法 （共役勾配法等） マルチコアCPU or GPUを用いて高速化する 第21話 共役勾配法 前処理演算 …… …………… ………. 33 ベクトルの和 行列とベクトルの積 ベクトルの内積 有限要素法の剛性行列A とその格納法 疎行列 CRS形式で格納 （Compressed Row Storage） 例 1 0 A 2 0 2 0 2 0 0 0 3 4 a (:) 1 2 ・ ・ ・ 5 2 4 2 DOF(:) 1 3 5 2 5 1 ind (1, :) 1 4 6 ・・・ ind (2, :) 3 5 8・・・ 5 4 0 3 3 4 ・・・ 4 ・・・ 第21話 マトリックスの記憶（反復法用） 行列とベクトルの積に必要な情報させ覚えておけばよい。 ＝Fill-inを考えなくてもOK 標準フォーマットなどは以下を参考に！ ＜CRS形式の例＞Matrix Market→http://math.nist.gov/MatrixMarkt/ 10 1 A 2 0 2 2 4 0 0 0 3 0 4 5 0 0 val 10. 2. colind 1 4 1 rowptr 35 2 1 -1. 1 2. 2 4. 3 2. 1 .3 4 -4. 4 -5. 5 3 4 5 6 7 8 9 3 6 8 10 非ゼロの値 列番号 10 行の先頭のアドレス FEMのSMP並列化 格言 Element By Element法は使うな！ Node By Node法を使え！ EBE法 0 = 0 NBN法 アセンブリングした非ゼロ成分のみを圧縮して覚える（CRSなど） 0 0 36 行列ベクトルの代入先が重複 データレース（書き込み漏れ）の可能性が ベンチマーク 3次元弾性体問題 0.1kN 2m Y X Z 10m 2m 隣接行列（1.5k） 自由度数は1.5k，93k，218k 行あたりの平均の非ゼロ成分数は40程度 ベンチマーク結果 3次元弾性体問題 8秒 218k GTX285 ベクトル型(キャッ シュ有) GTX285 ベクトル型(キャッ シュ無) GTX285 スカラ型 93k i7 2.66GHz 4コア 0.2秒 i7 2.66GHz 1コア 1.5k 0G 2G 4G 6G まとめ：疎行列のAx=bの反復求解 •一般的なPCに比べ，3～4倍の計算 速度を得た（∵ memory bound） •自由度が大きいほど計算速度が向 上していた（∵ レイテンシ） GPUを即戦力としたいならば… お勧め：簡単にはなったが専門家と組む • コンピュータに詳しいと自覚できない場合 * 概要を理解したら，詳細は投げるべき * 概要＝memory/CPU boundを考えた上限Gflops • コンサルテーション（各社）の利用 • ソフト（例：MASAMUNE）の利用 • 情報処理系の人と共同戦線 * じゃじゃ馬に「腕が鳴る」 * 「ソフトウェア マネジド キャッシュ」 * 「直線ドッカンのターボ」 40 素人(アサイ)の失敗① ISPH法の高速化に向け、GPUを使ってみよう（外注で） ①近傍粒子探索 ここだけGPUを使ってみよう ②粒子ごとに仮の速度更新 （マルチコア対応） ③圧力ポアソン方程式の解法 (マルチコア対応) ④粒子の加速度・速度・位置の更新 （マルチコア対応） 41 素人(アサイ)の失敗② ISPH法の高速化に向け、GPUを使ってみよう（外注で） なぜ、近傍粒子探索にGPU？ ・サーチ O(N2) ・Treeサーチ O(N) log(N) + マルチコア対応困難 ・バケットサーチ O(N) +CUDAサンプルで使用 Fortranからバケットサーチのカーネル(外部ルーチン)として 呼ぶまでの作業を外注で依頼しよう！（デモコードとして実績あり） 100万粒子がリアルタイムで流れている… (1世代前のもので500GFlops) 42 素人(アサイ)の失敗③ ISPH法の高速化に向け、GPUを使ってみよう（外注で） ①近傍粒子探索 (GPU対応？) 30倍高速化？ ②粒子ごとに仮の速度更新 （マルチコア対応） ③圧力ポアソン方程式の解法 (マルチコア対応) ④粒子の加速度・速度・位置の更新 （マルチコア対応） 43 全体でどれだけ高速化が実現可能か？ GPU⇔CPUのメモリの転送速度は？ 2012年以降、30～50コアへ GPU Xeon, Opteron 200～300コア 6～12コア 演算器 OpenMPで並列 CUDA、OpenCLで並列 C言語の知識がほしい！ キャッシュ （テキスチャメモリ） 演算器 キャッシュ 並立化効率に効く 160GB/sec 25GB/sec 主メモリ（4GBまで） 44 CPU 主メモリ（100GB以上も） 8GB/sec 遅い！ 素人(アサイ)の失敗④ ISPH法の高速化に向け、GPUを使ってみよう（外注で） ①近傍粒子探索 0.3 =>0.01 (GPU対応？) 30倍高速化？ ②粒子ごとに仮の速度更新0.1 （マルチコア対応） ③圧力ポアソン方程式の解法0.5 ホットスポットはここ！ (マルチコア対応) ④粒子の加速度・速度・位置の更新0.1 （マルチコア対応） 45 1.0 ⇒0.71程度？ アサイの反省 1. 2. 3. 4. 冷静にホットスポットを試算 どうせ使うのならフルGPU化 もっとハードを勉強しよう！ Intelの次世代CPU・GPUコアに注目 メモリが統合される！ 46