CFDデータをグラフで表現する方法

2次元の円柱周り流れを対象として、グラフニューラルネットワークでCFDのサロゲートモデルを構築していきます。第一弾として、この記事では、円柱周りのメッシュをグラフで表現するとはどういうことか、Pythonを使って解説していきます。numpyだけでメッシュ作成し、都度可視化して中身を確かめていきます。

メッシュ作成の流れ

グラフに必要なのは計算領域のノード、エッジ、ノード特徴量、エッジ特徴量です。4ステップでメッシュを作成します。

ノードの作成：計算領域内の点を表すノードを作成します。今回は約700点ほど作成します。
エッジの作成：隣接ノードを接続するエッジを作成します。
ノード特徴量の作成：ノードには特徴量を設定します。今回は流速と境界タイプ（interior, inlet, outlet, wall）を割り当てます。
エッジ特徴量の作成：エッジには相対位置、距離を割り当てます。

ノードの作成

まず、計算領域の境界上と内部領域にノード（点）を生成する関数make_points()を定義します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
 
rng = np.random.default_rng(0)

LX, LY = 3.0, 1.0            # 流路幅
CX, CY, R = 0.8, 0.5, 0.15   # 円柱の中心座標・半径

def make_points():
    pts = []
    # 流入・流出・壁面
    for y in np.linspace(0, LY, 21):
        pts.append((0.0, y))          # inlet
        pts.append((LX, y))           # outlet
    for x in np.linspace(0, LX, 61)[1:-1]:
        pts.append((x, 0.0))          # bottom wall
        pts.append((x, LY))           # top wall
    # 円柱表面
    for th in np.linspace(0, 2 * np.pi, 20, endpoint=False):
        pts.append((CX + R * np.cos(th), CY + R * np.sin(th)))
    # 内部点（円柱の内側は除外、円柱近傍をやや密に）
    n_trial = 2000

    # 内部ノードは境界から少し離して生成
    cand = rng.uniform([0.02, 0.02], [LX - 0.02, LY - 0.02], size=(n_trial, 2))

    # 円柱表面に近すぎるノードは除外
    d = np.hypot(cand[:, 0] - CX, cand[:, 1] - CY)
    keep = d > R * 1.15
    cand = cand[keep]
    selected = []

    # ノード同士の距離が近すぎる場合も削除
    min_dist = 0.05
    for p in cand:
        if all(np.hypot(p[0] - q[0], p[1] - q[1]) > min_dist for q in selected):
            selected.append(tuple(p))
    pts.extend(selected)
    return np.array(pts)

ノードを生成し可視化してみます。

pos = make_points()  # [N, 2] ノード座標
N = len(pos)  # 680ぐらい

fig, ax = plt.subplots(figsize=(10, 4))
ax.scatter(pos[:, 0], pos[:, 1])
ax.legend()
ax.set_aspect("equal")
fig.tight_layout()
fig.show()

生成したノード pos の可視化

円柱周り流れのノードが作成できました。流路の外周と円柱表面には境界ノードが等間隔に並び、内部にはノード同士の最小距離を保ったままランダムにノードが配置されています。円柱の内側にノードがないことも確認できます。

エッジの作成

続いて、ドロネー分割というアルゴリズムを利用して、ノード同士を繋ぐエッジを作成します。ドロネー分割は、（2次元の場合）ある平面上に定義された複数の点を使って、その平面を小さな三角形群に分割するアルゴリズムです。先ほど作成したposを使って計算領域をドロネー分割し、三角形群の辺をエッジにします。

ドロネー分割はscipy.spatial.Delaunayで簡単に実施できます。

tri = Delaunay(pos)
tri.simplices
# array([[363, 175, 235],
#        [493, 675, 498],
#        [549, 572, 246],
#        ...,
#        [178, 167, 177],
#        [167, 160, 161],
#        [165, 167, 164]], dtype=int32)

分割後の三角形群はtri.simplicesに格納されます。上記だと、例えば363番、175番、235番のノードで三角形がひとつ定義されていることを意味します。

円柱内部の三角形は不要なので除外します。

centers = pos[tri.simplices].mean(axis=1)
inside_cyl = np.hypot(centers[:, 0] - CX, centers[:, 1] - CY) < R
simplices = tri.simplices[~inside_cyl]

続いて、円柱内部を除外したあとの三角形群simplicesからエッジを集めます。

edges = set()
for a, b, c in simplices:
    for i, j in ((a, b), (b, c), (c, a)):
        # まず番号の小さいノード → 大きいノードへのエッジだけ集める
        edges.add((min(i, j), max(i, j)))
edges = np.array(sorted(edges))

# 流体の場合、ノード間情報は双方向でやりとりされる方が良さそう
# なので、番号の大きいノード → 小さいノードへのエッジも追加
edge_index = np.concatenate([edges, edges[:, ::-1]]).T   # [2, E]
E = edge_index.shape[1]

集めたエッジはedgesに格納しています。エッジは[始点ノード, 終点ノード]の形で表現されます。

edges
# array([[  0,   2],
#        [  0,  42],
#        [  1,   3],
#        ...,
#        [661, 675],
#        [669, 684],
#        [679, 685]], dtype=int32)

エッジも可視化してみましょう。

生成したエッジの可視化

隣接ノードを結ぶエッジがいい感じに生成できています。円柱内部の三角形を除外したので、円柱の内側にはエッジが張られていません。

ノード特徴量の作成

メッセージパシングの記事で述べているように、ノードには特徴量ベクトル $\mathbf{x}_u$ を定義する必要があります。今回は流速( $u$ , $v$ )と境界タイプ（内部、入口、出口、壁）を特徴量とします。

まず各ノードへの境界タイプの割り当てです。ノード座標が各境界の座標に一致するかどうかで割り当てています。

INTERIOR, INLET, OUTLET, WALL = 0, 1, 2, 3
node_type = np.full(N, INTERIOR)
eps = 1e-6
node_type[np.abs(pos[:, 0] - 0.0) < eps] = INLET
node_type[np.abs(pos[:, 0] - LX) < eps] = OUTLET
node_type[(np.abs(pos[:, 1]) < eps) | (np.abs(pos[:, 1] - LY) < eps)] = WALL
on_cyl = np.abs(np.hypot(pos[:, 0] - CX, pos[:, 1] - CY) - R) < 1e-3
node_type[on_cyl] = WALL

one_hot = np.eye(4)[node_type]
one_hot # one_hot.shape = (N, 4)
# array([[0., 0., 0., 1.],
#        [0., 0., 0., 1.],
#        [0., 1., 0., 0.],
#        ...,
#        [1., 0., 0., 0.],
#        [1., 0., 0., 0.],
#        [1., 0., 0., 0.]])

境界タイプはone-hotエンコーディングで4次元のベクトルとして扱います。

境界タイプの可視化

可視化すると、左端の青いノードがinlet、右端の赤いノードがoutlet、上下の壁面と円柱表面の緑のノードがwall、それ以外の灰色のノードがinteriorに割り当てられていることが確認できます。

続いて流速です。入口流速 $u$ は次の放物線分布で与えます。

u(y) = -\frac{4}{L_y^2}\left(y - \frac{L_y}{2}\right)^2 + 1 \tag{1}

ここで $y$ は流路高さ方向の座標、 $L_y$ は流路高さです。式(1)は下壁 $y=0$ と上壁 $y=L_y$ で $u=0$ 、流路中央 $y=L_y/2$ で最大値 $u=1$ となる分布です。壁面はすべりなし条件で $u=0$ とし、内部と出口のノードには適当な初期値 $u=0.5$ を与えます。

u = np.zeros(N)
v = np.zeros(N)
u[node_type == INTERIOR] = 0.5 # 適当な内部初期値
inlet_mask = node_type == INLET
u[inlet_mask] = - 4 / LY**2 * (pos[inlet_mask, 1] - LY/2)**2 + 1  # 放物線
u[node_type == OUTLET] = 0.5  # 適当な内部初期値
u[node_type == WALL] = 0.0  # no-slip
u[:10]
# array([0.  , 0.  , 0.19, 0.5 , 0.36, 0.5 , 0.51, 0.5 , 0.64, 0.5 ])

流速 u の可視化

入口で放物線分布（流路中央ほど明るい色）、壁面と円柱表面で0、内部で一様に0.5になっていることが確認できます。なお、 $v$ は全ノードで0のままです。

流速と境界タイプを結合すればノード特徴量xが得られます。ノード特徴量は6次元ベクトルです。

x = np.column_stack([u, v, one_hot])  # [N, 6]
x[:5]
# u, v, interior, inlet, outlet, wall
# array([[0.  , 0.  , 0.  , 0.  , 0.  , 1.  ],
#        [0.  , 0.  , 0.  , 0.  , 0.  , 1.  ],
#        [0.19, 0.  , 0.  , 1.  , 0.  , 0.  ],
#        [0.5 , 0.  , 0.  , 0.  , 1.  , 0.  ],
#        [0.36, 0.  , 0.  , 1.  , 0.  , 0.  ]])

エッジ特徴量の作成

エッジにはノード間の相対位置( $dx, dy$ )とノード間距離を特徴量として割り当てることにします。絶対座標を使うことも可能ですが、今回は敢えて相対位置にします。相対位置であれば、グラフ全体を平行移動しても特徴量が変わりません。位置に依存しない、汎化しやすいモデルの学習につながります。

まず相対位置を計算します。始点ノードから終点ノードへ向かうベクトルです。

src, dst = edge_index
rel = pos[dst] - pos[src] # [E, 2]
rel[:10]
# dx, dy
# array([[ 0.        ,  0.05      ],
#        [ 0.05      ,  0.        ],
#        [ 0.        ,  0.05      ],
#        [-0.05      ,  0.        ],
#        [-0.04318306,  0.03521853],
#        [ 0.        ,  0.05      ],
#        [ 0.05      , -0.05      ],
#        [ 0.05518061,  0.02030679],
#        [ 0.        ,  0.05      ],
#        [-0.04318306, -0.01478147]])

ノード間距離も計算し、エッジ特徴量とします。エッジ特徴量edge_attrは3次元ベクトルです。

dist = np.linalg.norm(rel, axis=1, keepdims=True)  # [E, 1]
edge_attr = np.hstack([rel, dist])
edge_attr  # [E, 3]
# dx, dy, |d|
# array([[ 0.        ,  0.05      ,  0.05      ],
#        [ 0.05      ,  0.        ,  0.05      ],
#        [ 0.        ,  0.05      ,  0.05      ],
#        ...,
#        [-0.08514992, -0.00403709,  0.08524557],
#        [-0.06914038,  0.01746942,  0.0713132 ],
#        [ 0.01847789,  0.11067321,  0.11220513]])

まとめ

CFDデータをグラフとして表現するために、numpyだけで以下の4つを作成しました。

構成要素	変数	形状
ノード座標	`pos`	[N, 2]
エッジ	`edge_index`	[2, E]
ノード特徴量（流速＋境界タイプ）	`x`	[N, 6]
エッジ特徴量（相対位置＋距離）	`edge_attr`	[E, 3]

今回、ノードはランダム配置と境界配置の組み合わせで生成し、エッジはドロネー分割で隣接関係を定義しました。この作業は擬似的にCFDデータを作ったことに相当します。実務ではメッシュや物理場といったCFDデータをそのままノードやエッジに割り当てることになります。

次回は実際にCFDを実施し、MeshGraphNetの学習データセットへ変換していきます。

参考文献

GNNの基礎を学べます。

Graph Representation Learning（William L. Hamilton, Morgan & Claypool, 2020）：GNNの理論を体系的に解説した定番書です。本シリーズの理論編でも参照しており、メッセージパッシングやプーリングの数学的な定式化を丁寧に学べます。

流体力学で考える

日常にある流体を“見える化”するブログ