Zenn
🗂

流体シミュレーション結果にグラフ信号処理を適用してみる

2024/06/09に公開
グラフフーリエ変換
tech

はじめに

グラフ信号処理とはグラフフーリエ変換を利用して、グラフ上の各頂点上に配置された離散信号を分析することです。グラフ構造としては、SNSにおける人物関係や、分子構造、画像（各ピクセルを頂点とみなしグラフとする）が代表的ですが、FEM（有限要素法）やFVM（有限体積法）などの物理シミュレーションで使用されるデータ（メッシュデータ）もグラフ構造とみることができます。（下図参照）

今回は簡易的な流体の数値シミュレーションデータにグラフ信号処理を適用してみて、その結果を確認してみたいと思います。


メッシュデータの例。図の青線がエッジに対応し各交点が頂点に対応するグラフ構造とみなせる

グラフ信号処理概要

理論的な枠組みについて最低限のみ以下に示します。より詳細は、参考文献[1] [2] [3] を参照してください。

計算方法

グラフフーリエ変換

記号を以下のように定義します。

  • N : 節点数
  • \bm L : グラフラプラシアン
  • \bm u_{\lambda_i} : グラフラプラシアンの i 番目の固有値 \lambda_iに対応する固有ベクトル（グラフフーリエ基底）
  • \bm U = {\bm u_{\lambda_0}, \bm h_{\lambda_1}, ..., \bm h_{\lambda_N}} : グラフフーリエ基底をまとめた行列
  • \bm f : 頂点領域におけるグラフ信号
  • \bm H : スペクトル領域におけるフィルタカーネル

グラフフーリエ変換、逆グラフフーリエ変換は以下のように定義されます。

  • グラフフーリエ変換
\bm F = \bm U^{T} \bm f
  • 逆グラフフーリエ変換
\bm f = \bm U \bm F

フィルタリング

スペクトル領域でのフィルタリングを考えます。スペクトル領域におけるフィルタカーネル \bm H は以下のように表されます。

\begin{pmatrix} H(\lambda_0) \\ & H(\lambda_1) \\ & & H(\lambda_2) & \\ & & & \ddots \\ & & & & H(\lambda_{N - 2}) \\ & & & & & H(\lambda_{N - 1}) \end{pmatrix}

フィルタリングされた頂点領域での信号を \bm f_{\text{out}} とすると、以下のように表されます。

\bm f_{\text{out}} = \bm U \bm H \bm U^{T} \bm f_{\text{in}}

図化してまとめると下図のようになります。

フーリエ変換による信号処理との差異

グラフフーリエ変換と（通常の）フーリエ変換との違いを整理します。

最大のポイントは、信号の観測位置同士の「繋がり方」の扱いにあります。k-インデックスの頂点における信号 f(k) が与えられたとき、フーリエ変換では、f(k - 1) は、f(k) の時間的に１サンプル前の信号であることを意味します。一方で、グラフフーリエ変換では、時間的な「前後」関係を仮定しません。グラフ構造では、頂点のインデックスが連続しているとしても直接つながっているとは限らないことや、グラフ構造が空間的な広がりを持っているため、「前後」関係を定義できません。

以下に、グラフフーリエ変換と、通常のフーリエ変換との差異をまとめます。

グラフフーリエ変換 通常のフーリエ変換
信号の位置 グラフ構造をもとに判断する １インデックスが時間的に１サンプル前であることを前提とする
|固有値| << 1 であることの意味 直接接続している頂点の中で距離が長い頂点同士であっても、似たような値を取る。 前後の信号が似た値をとる（長周期の振動）
|固有値| >> 1 であることの意味 直接接続している頂点の中で距離が短い頂点同士であっても、大きく異なる値を取る。 前後の信号が大きく異なる値をとる（短周期の振動）

計算量について

上記の通り、グラフフーリエ変換を実施するにはグラフラプラシアンの固有値、固有ベクトルを求める必要があり、この計算量がボトルネックとなります。固有値、固有ベクトルに関するアルゴリズムは多く考案されていますが、密行列の全固有値を求める代表的なアルゴリズムであるQRアルゴリズムを例にとると、一般に \Omicron(N^3) であるとされます。[4]

ただし、疎行列の場合や固有値の一部のみを計算する場合には計算量が緩和されることもあります。詳細は、scipyの粗行列固有値計算モジュールとして使用されているARPACKのマニュアルなどが参考になります。[5] [6]

グラフ信号処理の適用

シミュレーション例 1（Cavity流れ）

シミュレーション概要

簡単な流れ場として Cavity 流れから見てみたいと思います。本シミュレーションは、箱の上面が一定速度で動くことで、内部の流体が動く様子をシミュレーションしたものです。

項目
節点数 882
ソルバー icoFoam
ソフトウェア OpenFOAM-v2012
時間発展 あり（0.5秒間の非定常シミュレーション）

以降では、シミュレーションの最終時刻（0.5秒）での結果を用います。


Cavity流れのシミュレーション結果
（コンターは流速の大きさを表し、矢印が方向を表しています。）

モードを見る

固有モードを確認します。アップロードできるファイルサイズの関係上、下のgifでは低周波数のモードを小さい順に20個表示しています。


固有モード（20個分）

スペクトル分析

圧力および流速ノルムをグラフ信号だとみなして、スペクトルを分析した結果を下図に示します。以下のことがわかります。

  • 低固有値成分が支配的である。
  • 速度ノルムに関するスペクトルは周期性のようなパターンが確認できる。


圧力のスペクトル分析結果


流速ノルムのスペクトル分析結果

フィルタリング

スペクトル分析の結果をもとに、固有値が1.08 ( 400番目の固有値 ) 以下の成分のみを取り出す Low-Passフィルタと、1.08より大きい成分のみを取り出す High-Pass フィルタを、速度ベクトルにかけてみます。

以下のことがわかります。

  • スペクトル分析の結果の通り、低周波数成分で概ねもとの結果を再現できる。
  • 高周波数成分は矩形形状上部の両隅部に現れている。


フィルタリングの結果（左：Low-Pass フィルタ、右: High-Pass フィルタ）

シミュレーション例 2（ mixing-elbow ）

シミュレーション概要

Cavity流れよりも節点数が大きく、全固有値を計算しきれないケースとしてmixing-elbowモデルを用意しました。本シミュレーションは、パイプのような形状をした構造にZ方向下部の流入口から流体を流し入れたときの様子をシミュレーションしたものです。[7]

項目
節点数 15246
ソルバー simpleFoam
ソフトウェア OpenFOAM-v2012
時間発展 なし（定常）


mixing_elbowモデル （左：3Dモデル外観、右：Y=0の断面図。コンターが流速の大きさ、矢印が流速の方向を表す。）

固有モードの確認

固有モードを確認します。アップロードできるファイルサイズの関係上、下のgifでは低周波数のモードを小さい順に20個表示しています。


モード図（同じ固有モードを異なる方向から見ている。）

フィルタリング

節点数が大きいため、すべての固有値を求めるには時間がかかります。そこで、固有値（周波数）の小さい方もしくは大きい方から数個求めることが可能であれば、それぞれLow-Passフィルタ、High-Passフィルタとすることができることを利用し、フィルタリングをかけてみます。固有値をそれぞれ30個ずつもとめて速度場にフィルタリングをかけた結果を下図に示します。

  • 低周波数に対応する信号は、パイプ上部で相対的に高い値をとっている。
  • 高周波数に対応する信号は、流入口付近で相対的に高い値をとっている。


フィルタリングの結果（左：Low-Pass フィルタ、右: High-Pass フィルタ）

さいごに

以上の議論では、OpneFOAMによる流体数値シミュレーションデータに対して、グラフフーリエ変換を適用し、その結果を確認しました。

今回は単純にグラフフーリエ変換を適用してみただけでしたが、以下のようなことが分かると数値計算結果の検証や、ニューラルネットワークにおける損失関数の設定などに役立てる可能性があります。

  • Low-Passフィルタを通した信号の強度と、境界条件の影響力との関係（低周波数成分が大きい箇所は境界条件の影響を強くうけているということなのか）
  • スペクトル領域の高周波数成分と、乱流との関係

参考資料

脚注

  1. David I Shuman, Sunil K. Narang, Pascal Frossard, Antonio Ortega and Pierre Vandergheynst, "The Emerging Field of Signal Processing on Graphs - Extending High-Dimensional Data Analysis to Networks and Other Irregular Domain", https://arxiv.org/pdf/1211.0053, 2012, IEEE Signal Processing Magazine ↩︎

  2. 田中雄一, "グラフ信号処理のすゝめ", https://www.jstage.jst.go.jp/article/essfr/8/1/8_15/_pdf/-char/ja, アクセス日2024.6.9 ↩︎

  3. 佐藤竜馬, "グラフニューラルネットワーク", 講談社, 2024, p146-175 ↩︎

  4. "Eigenvalue algorithms: The QR algorithm", https://www.sci.utah.edu/~akil/docs/courses/2020fall/math6610/lec20.pdf, アクセス日2024.6.9 ↩︎

  5. "ARPACK Users' Guide", http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf, アクセス日2024.6.9 ↩︎

  6. "固有値計算法", RIKEN AICS HPC Spring School, 今村俊幸, https://aics.riken.jp/aicssite/wp-content/uploads/2014/03/sps14_enshu2-1.pdf, アクセス日2024.6.9 ↩︎

  7. OpenFOAM-Primer-Support, https://gitlab.com/OpenCAE/OpenFOAM-Primer-Support, アクセス日2024.6.9 ↩︎

Discussion