物理学に基づくニューラルネットワーク（PINNs）を使った流体解析
By M.Sato
1. 概要
PINNs（Physics-Informed Neural Networks）は、偏微分方程式や常微分方程式で表される物理法則をニューラルネットワークの損失関数に組み込み、 データと物理法則の両方を考慮して解を導く解析手法である。
💡 PINNsの特徴
✅ メッシュ生成が不要
✅ 複雑な境界条件や不完全なデータに対しても柔軟に対応可能
✅ 逆問題（未知のパラメータを推定）にも応用可能
❌ 損失関数が複雑になりやすい
❌ 計算コストが高い
❌ 精度保証が難しい
ここでは、PINNsを使って蓋駆動正方キャビティの計算を行ってみた。 今回の計算条件では概ね妥当な計算結果が得られた。
なお今回の解析ではここにある資料を参考にし、 プログラムに若干の修正を加えた。
2. PINNsの解析手法
2.1. ニューラルネットワーク
PINNsにおけるニューラルネットワークは基礎となる微分方程式の解を近似するために使われる。
ニューラルネットワークは、入力層、複数の隠れ層、出力層から構成される。
各層には複数のニューロンがあり、ニューロン同士は重みとバイアスを介して接続されている。
隠れ層では、活性化関数（例えばTanh関数）を用いて非線形変換を行う。
ℹ️ 流体解析におけるPINNsの入出力
入力: 空間座標（例：$x$, $y$）
出力: 流速成分（例：$u$, $v$）、圧力（$p$）
2.2. 損失関数
従来のデータ駆動型のニューラルネットワークでは観測データに基づいてモデルを学習するが、 PINNsでは物理法則を損失関数に組み込むことで、 データが不足している場合でも物理的に一貫した解を得ることができる。
💡 PINNsの核心: 物理法則を損失関数に組み込むことで、データ駆動と物理駆動を融合
例えば流体解析で使用するナビエ・ストークス方程式の場合、 以下のような損失関数を定義することができる。
損失の式に含まれる$\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}$ などは、ニューラルネットワークで得られた予測値$u$の空間勾配として計算される。
全損失は大きく２つの部分から構成されている。
支配方程式に従うことを強制する損失関数
微分方程式の残差を最小化する損失関数を定義する。
$$ \begin{align} L_{PDE} &= mean(L_p^2 + L_u^2 + L_v^2) \end{align} $$
ここで、$L_p$、$L_u$、$L_v$はそれぞれ連続式、 $x$方向運動方程式、$y$方向運動方程式の残差である。
$$ \begin{align} L_p &=\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \\ L_u &=u\frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + \frac{1}{\rho} \frac{\partial p}{\partial x} - \nu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\ L_v &=u\frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + \frac{1}{\rho} \frac{\partial p}{\partial y} - \nu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) \end{align} $$
境界条件を満たすことを強制する損失関数
各境界において境界条件を満たすように損失関数を定義する。
例えば以下に示すキャビティモデルでは以下のような境界条件損失を定義できる。
$$ \begin{align} L_{BC} &= L_{top} + L_{bottom} + L_{left} + L_{right} \end{align} $$
ここで、
$$ \begin{align} L_{top} &= mean((u - U_{lid})^2) + mean(v^2) \\ L_{bottom} &= mean(u^2) + mean(v^2) \\ L_{left} &= mean(u^2) + mean(v^2) \\ L_{right} &= mean(u^2) + mean(v^2) \end{align} $$
全体の損失関数
$$ \begin{align} L_{total} &= \lambda L_{PDE} + L_{BC} \end{align} $$
ここで、$\lambda$はPDE損失の重みであり、境界条件損失とPDE損失のバランスを調整するために使用される。
2.3. 最適化
損失関数が定義されたら、最適化アルゴリズムを用いてニューラルネットワークの重みとバイアスを更新し、損失関数を最小化する。
ℹ️ 代表的な最適化手法
✅ Adam法: 確率的勾配降下法の一種。学習初期に有効
✅ L-BFGS法: 準ニュートン法。収束精度を高めるのに有効
2.4. 解析手法の流れ図
解析手順を図で表示すると以下のようになる。
3. 蓋駆動正方キャビティ流れへの適用
3.1. 解析モデル
蓋駆動正方キャビティ問題は、流体力学における古典的なベンチマーク問題であり、2次元の正方形領域内で、 上部の蓋が一定速度で動くことによって流れが駆動される問題である。
⚙️ 計算条件: レイノルズ数 Re = 100
なお参照した実験データは以下の論文中に書かれているものを使用した。
3.2. 解析結果
PINNsを使って得られた流速分布と圧力分布は以下の通りである。
以下の計算ではGPUを使用した。
また、流速の中心線分布を実験データと比較した結果は以下の通りである。
✅ 検証結果: PINNsの計算結果は概ね実験データと一致している
損失関数の履歴は以下の通りである。
今回の計算では、最初にAdam法で1000エポック学習を行い、 その後L-BFGS法で500エポック学習を行った。
3.3. プログラム例
pythonを使ったPINNsの実装例を以下に示す。
PyTorchを用いてニューラルネットワークを構築し、 Adam法とL-BFGS法を組み合わせて学習を行うようにした。
ハイパーパラメータの管理にはHydraを使用し、 設定ファイル(config.yaml)からパラメータを読み込むようにした。
#!/usr/bin/env python
# coding: utf-8
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
from torch.autograd import grad
from torch.optim import Adam, LBFGS
import hydra
from omegaconf import DictConfig, OmegaConf
import os # 出力ディレクトリの管理に必要
#----------------------------------------- Check for GPU
def get_device():
if torch.cuda.is_available():
device = torch.device("cuda")
print("Using GPU:", torch.cuda.get_device_name(0))
else:
device = torch.device("cpu")
print("Using CPU")
return device
#----------------------------------------- define network class
class PINN(nn.Module):
def __init__(self, layer_sizes):
super(PINN, self).__init__()
# ネットワーク構造をリストで動的に定義
layers = []
for i in range(len(layer_sizes) - 1):
layers.append(nn.Linear(layer_sizes[i], layer_sizes[i+1]))
# 最後のレイヤー以外にTanhアクティベーションを追加
if i < len(layer_sizes) - 2:
layers.append(nn.Tanh())
self.net = nn.Sequential(*layers)
def forward(self, x, y):
# Ensure x and y are column vectors for concatenation
inputs = torch.cat((x.view(-1, 1), y.view(-1, 1)), dim=1)
outputs = self.net(inputs)
u = outputs[:, 0].view(-1, 1)
v = outputs[:, 1].view(-1, 1)
p = outputs[:, 2].view(-1, 1)
return u, v, p
#----------------------------------------- 微係数の計算
def compute_derivatives(u, v, p, x, y):
u_x = grad(u, x, torch.ones_like(u), create_graph=True, allow_unused=True)[0]
u_y = grad(u, y, torch.ones_like(u), create_graph=True, allow_unused=True)[0]
v_x = grad(v, x, torch.ones_like(v), create_graph=True, allow_unused=True)[0]
v_y = grad(v, y, torch.ones_like(v), create_graph=True, allow_unused=True)[0]
p_x = grad(p, x, torch.ones_like(p), create_graph=True, allow_unused=True)[0]
p_y = grad(p, y, torch.ones_like(p), create_graph=True, allow_unused=True)[0]
u_x = u_x if u_x is not None else torch.zeros_like(u)
u_y = u_y if u_y is not None else torch.zeros_like(u)
v_x = v_x if v_x is not None else torch.zeros_like(v)
v_y = v_y if v_y is not None else torch.zeros_like(v)
p_x = p_x if p_x is not None else torch.zeros_like(p)
p_y = p_y if p_y is not None else torch.zeros_like(p)
u_xx = grad(u_x, x, torch.ones_like(u_x), create_graph=True, allow_unused=True)[0]
u_yy = grad(u_y, y, torch.ones_like(u_y), create_graph=True, allow_unused=True)[0]
v_xx = grad(v_x, x, torch.ones_like(v_x), create_graph=True, allow_unused=True)[0]
v_yy = grad(v_y, y, torch.ones_like(v_y), create_graph=True, allow_unused=True)[0]
u_xx = u_xx if u_xx is not None else torch.zeros_like(u_x)
u_yy = u_yy if u_yy is not None else torch.zeros_like(u_y)
v_xx = v_xx if v_xx is not None else torch.zeros_like(v_x)
v_yy = v_yy if v_yy is not None else torch.zeros_like(v_y)
return u_x, u_y, v_x, v_y, p_x, p_y, u_xx, u_yy, v_xx, v_yy
#----------------------------------------- pde損失の計算
def pde_loss(u, v, p, x, y, cfg):
# 計算領域内部における微係数を計算
u_x, u_y, v_x, v_y, p_x, p_y, u_xx, u_yy, v_xx, v_yy = compute_derivatives(u, v, p, x, y)
NU = cfg.constants.nu
RHO = cfg.constants.rho
# Navier-Stokes in 2D (Steady, Incompressible)
f_continuity = u_x + v_y
# Momentum X: u*u_x + v*u_y + p_x/rho - nu*(u_xx + u_yy) = 0
f_momentum_x = u * u_x + v * u_y + p_x/RHO - NU * (u_xx + u_yy)
# Momentum Y: u*v_x + v*v_y + p_y/rho - nu*(v_xx + v_yy) = 0
f_momentum_y = u * v_x + v * v_y + p_y/RHO - NU * (v_xx + v_yy)
loss_continuity = torch.mean(f_continuity**2)
loss_momentum_x = torch.mean(f_momentum_x**2)
loss_momentum_y = torch.mean(f_momentum_y**2)
return loss_continuity + loss_momentum_x + loss_momentum_y
#----------------------------------------- 境界条件損失の計算
def boundary_loss(model, num_points_boundary, U_LID_val, device):
# 境界点群の生成は関数内で実行
# Top Boundary (y=1): u=U_LID, v=0
x_top = torch.rand(num_points_boundary // 4, 1, requires_grad=True).to(device)
y_top = torch.ones_like(x_top).to(device)
u_top, v_top, _ = model(x_top, y_top)
loss_top = torch.mean((u_top - U_LID_val)**2) + torch.mean(v_top**2)
# Bottom Boundary (y=0): u=0, v=0
x_bottom = torch.rand(num_points_boundary // 4, 1, requires_grad=True).to(device)
y_bottom = torch.zeros_like(x_bottom).to(device)
u_bottom, v_bottom, _ = model(x_bottom, y_bottom)
loss_bottom = torch.mean(u_bottom**2) + torch.mean(v_bottom**2)
# Left Boundary (x=0): u=0, v=0
y_left = torch.rand(num_points_boundary // 4, 1, requires_grad=True).to(device)
x_left = torch.zeros_like(y_left).to(device)
u_left, v_left, _ = model(x_left, y_left)
loss_left = torch.mean(u_left**2) + torch.mean(v_left**2)
# Right Boundary (x=1): u=0, v=0
y_right = torch.rand(num_points_boundary // 4, 1, requires_grad=True).to(device)
x_right = torch.ones_like(y_right).to(device)
u_right, v_right, _ = model(x_right, y_right)
loss_right = torch.mean(u_right**2) + torch.mean(v_right**2)
return loss_top + loss_bottom + loss_left + loss_right
#----------------------------------------- トレーニングの実行
def train_pinn(model, cfg, device):
# cfgから設定値の抽出
U_LID = cfg.constants.u_lid
NUM_ADAM_EPOCHS = cfg.training.adam.epochs
NUM_LBFGS_EPOCHS = cfg.training.lbfgs.epochs
NUM_COLLOCATION_POINTS = cfg.data.num_collocation_points
NUM_BOUNDARY_POINTS = cfg.data.num_boundary_points
ADAM_LEARNING_RATE = cfg.training.adam.lr
PDE_LOSS_WEIGHT = cfg.training.pde_loss_weight
history = []
# --- Adam Optimizer Phase ---
print("\n--- Starting Adam Optimization ---")
optimizer_adam = Adam(model.parameters(), lr=ADAM_LEARNING_RATE)
for epoch in range(NUM_ADAM_EPOCHS):
optimizer_adam.zero_grad()
# 常に新しい点群を生成 (静的な点群を使用したい場合は、ループ外で生成します)
x_collocation_adam = torch.rand(NUM_COLLOCATION_POINTS, 1, requires_grad=True).to(device)
y_collocation_adam = torch.rand(NUM_COLLOCATION_POINTS, 1, requires_grad=True).to(device)
u_pred_collocation, v_pred_collocation, p_pred_collocation = model(x_collocation_adam, y_collocation_adam)
loss_pde = pde_loss(u_pred_collocation, v_pred_collocation, p_pred_collocation, x_collocation_adam, y_collocation_adam, cfg)
loss_bc = boundary_loss(model, NUM_BOUNDARY_POINTS, U_LID, device)
loss = PDE_LOSS_WEIGHT * loss_pde + loss_bc
loss.backward()
optimizer_adam.step()
if (epoch + 1) % 100 == 0:
print(f"Adam Epoch {epoch+1}/{NUM_ADAM_EPOCHS}, Loss: {loss.item():.6f}, PDE Loss: {loss_pde.item():.6f}, BC Loss: {loss_bc.item():.6f}")
history.append((epoch + 1, loss.item()))
# --- L-BFGS Optimizer Phase ---
print("\n--- Starting L-BFGS Optimization ---")
# L-BFGSは一般的に静的な点群を使用します
x_collocation_lbfgs = torch.rand(NUM_COLLOCATION_POINTS, 1, requires_grad=True).to(device)
y_collocation_lbfgs = torch.rand(NUM_COLLOCATION_POINTS, 1, requires_grad=True).to(device)
# 境界点群もL-BFGSのために静的に生成
x_boundary_lbfgs_top = torch.rand(NUM_BOUNDARY_POINTS // 4, 1, requires_grad=True).to(device)
y_boundary_lbfgs_top = torch.ones_like(x_boundary_lbfgs_top).to(device)
x_boundary_lbfgs_bottom = torch.rand(NUM_BOUNDARY_POINTS // 4, 1, requires_grad=True).to(device)
y_boundary_lbfgs_bottom = torch.zeros_like(x_boundary_lbfgs_bottom).to(device)
y_boundary_lbfgs_left = torch.rand(NUM_BOUNDARY_POINTS // 4, 1, requires_grad=True).to(device)
x_boundary_lbfgs_left = torch.zeros_like(y_boundary_lbfgs_left).to(device)
y_boundary_lbfgs_right = torch.rand(NUM_BOUNDARY_POINTS // 4, 1, requires_grad=True).to(device)
x_boundary_lbfgs_right = torch.ones_like(y_boundary_lbfgs_right).to(device)
optimizer_lbfgs = LBFGS(model.parameters(),
lr=cfg.training.lbfgs.lr,
max_iter=cfg.training.lbfgs.max_iter,
history_size=cfg.training.lbfgs.history_size,
line_search_fn=cfg.training.lbfgs.line_search_fn)
current_global_epoch = NUM_ADAM_EPOCHS
def closure():
nonlocal current_global_epoch # L-BFGSの反復ごとにグローバルエポックを更新
optimizer_lbfgs.zero_grad()
u_pred_collocation, v_pred_collocation, p_pred_collocation = model(x_collocation_lbfgs, y_collocation_lbfgs)
loss_pde = pde_loss(u_pred_collocation, v_pred_collocation, p_pred_collocation, x_collocation_lbfgs, y_collocation_lbfgs, cfg)
# 境界損失はL-BFGSのために静的に生成された点群を使用
u_top, v_top, _ = model(x_boundary_lbfgs_top, y_boundary_lbfgs_top)
loss_top = torch.mean((u_top - U_LID)**2) + torch.mean(v_top**2)
u_bottom, v_bottom, _ = model(x_boundary_lbfgs_bottom, y_boundary_lbfgs_bottom)
loss_bottom = torch.mean(u_bottom**2) + torch.mean(v_bottom**2)
u_left, v_left, _ = model(x_boundary_lbfgs_left, y_boundary_lbfgs_left)
loss_left = torch.mean(u_left**2) + torch.mean(v_left**2)
u_right, v_right, _ = model(x_boundary_lbfgs_right, y_boundary_lbfgs_right)
loss_right = torch.mean(u_right**2) + torch.mean(v_right**2)
loss_bc = loss_top + loss_bottom + loss_left + loss_right
# L-BFGSフェーズではPDE重みは1.0を使用 (cfgで上書き可能)
loss = cfg.training.lbfgs.pde_loss_weight * loss_pde + loss_bc
if loss.requires_grad:
loss.backward()
return loss
for epoch_lbfgs in range(NUM_LBFGS_EPOCHS):
# L-BFGSでは、max_iter回反復するごとにstep()を呼び出す必要があります。
# closure関数内で反復が管理されます。
loss_lbfgs = optimizer_lbfgs.step(closure)
current_global_epoch += 1
if (epoch_lbfgs + 1) % 10 == 0 or epoch_lbfgs == 0:
print(f"L-BFGS Epoch {epoch_lbfgs+1}/{NUM_LBFGS_EPOCHS} (Global: {current_global_epoch}), Loss: {loss_lbfgs.item():.6f}")
history.append((current_global_epoch, loss_lbfgs.item())) # Store (global epoch, loss)
return history
#----------------------------------------- 結果の保存および作図
def save_and_plot_results(model, cfg, training_history_raw, device):
DOMAIN_SIZE = cfg.constants.domain_size
REYNOLDS_NUMBER = cfg.constants.reynolds_number
num_test_points = 100
x_test_coords = torch.linspace(0, DOMAIN_SIZE, num_test_points).view(-1, 1)
y_test_coords = torch.linspace(0, DOMAIN_SIZE, num_test_points).view(-1, 1)
X_grid, Y_grid = torch.meshgrid(x_test_coords.squeeze(), y_test_coords.squeeze(), indexing='xy')
X_flat_tensor = X_grid.flatten().unsqueeze(1).to(device)
Y_flat_tensor = Y_grid.flatten().unsqueeze(1).to(device)
model.eval()
with torch.no_grad():
u_pred_flat, v_pred_flat, p_pred_flat = model(X_flat_tensor, Y_flat_tensor)
U_plot = u_pred_flat.cpu().numpy().reshape(num_test_points, num_test_points)
V_plot = v_pred_flat.cpu().numpy().reshape(num_test_points, num_test_points)
P_plot = p_pred_flat.cpu().numpy().reshape(num_test_points, num_test_points)
# --- SAVE DATA TO .DAT FILE ---
X_np_flat = X_flat_tensor.cpu().numpy()
Y_np_flat = Y_flat_tensor.cpu().numpy()
U_np_flat = u_pred_flat.cpu().numpy()
V_np_flat = v_pred_flat.cpu().numpy()
P_np_flat = p_pred_flat.cpu().numpy()
output_data = np.hstack((X_np_flat, Y_np_flat, U_np_flat, V_np_flat, P_np_flat))
output_filename = f"lid_driven_cavity_re{REYNOLDS_NUMBER}_pinn_results.dat"
np.savetxt(output_filename, output_data, fmt='%.8e', delimiter='\t', header='X\tY\tU\tV\tP', comments='')
print(f"Results saved to {output_filename}")
# --- Plotting 2D contours ---
fig, axs = plt.subplots(1, 3, figsize=(18, 6))
im1 = axs[0].imshow(U_plot, cmap='jet', origin='lower',
extent=[X_grid.min(), X_grid.max(), Y_grid.min(), Y_grid.max()])
axs[0].set_title('U Velocity')
axs[0].set_xlabel('X')
axs[0].set_ylabel('Y')
fig.colorbar(im1, ax=axs[0])
im2 = axs[1].imshow(V_plot, cmap='jet', origin='lower',
extent=[X_grid.min(), X_grid.max(), Y_grid.min(), Y_grid.max()])
axs[1].set_title('V Velocity')
axs[1].set_xlabel('X')
axs[1].set_ylabel('Y')
fig.colorbar(im2, ax=axs[1])
im3 = axs[2].imshow(P_plot, cmap='jet', origin='lower',
extent=[X_grid.min(), X_grid.max(), Y_grid.min(), Y_grid.max()])
axs[2].set_title('Pressure')
axs[2].set_xlabel('X')
axs[2].set_ylabel('Y')
fig.colorbar(im3, ax=axs[2])
skip_quiver = 8
axs[0].quiver(X_grid.cpu().numpy()[::skip_quiver, ::skip_quiver],
Y_grid.cpu().numpy()[::skip_quiver, ::skip_quiver],
U_plot[::skip_quiver, ::skip_quiver],
V_plot[::skip_quiver, ::skip_quiver], color='white', scale=5)
axs[1].quiver(X_grid.cpu().numpy()[::skip_quiver, ::skip_quiver],
Y_grid.cpu().numpy()[::skip_quiver, ::skip_quiver],
U_plot[::skip_quiver, ::skip_quiver],
V_plot[::skip_quiver, ::skip_quiver], color='white', scale=5)
plt.tight_layout()
plt.savefig('pinn_contours.png')
plt.close(fig)
print("Contour plot saved as pinn_contours.png")
# --- Plotting the loss history ---
epochs_for_plot = [item[0] for item in training_history_raw]
losses_for_plot = [item[1] for item in training_history_raw]
loss_data = np.column_stack((epochs_for_plot, losses_for_plot))
np.savetxt("loss.csv", loss_data, delimiter=",", header="epoch,loss", comments="")
plt.figure(figsize=(10, 6))
plt.plot(epochs_for_plot, losses_for_plot, label='Total Loss')
plt.title('PINN Training Loss History (Adam then L-BFGS)')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.yscale('log')
plt.grid(True)
plt.legend()
plt.savefig('loss_history.png')
plt.close()
print("Loss history plot saved as loss_history.png")
# 流速分布の実験結果(Ghia et al.)
ghia_u_y = np.array([1.000000000E+000, 9.765999913E-001, 9.688000083E-001, 9.609000087E-001, 9.531000257E-001, 8.515999913E-001, 7.343999743E-001, 6.172000170E-001, 5.000000000E-001, 4.530999959E-001, 2.813000083E-001, 1.719000041E-001, 1.015999988E-001, 7.029999793E-002, 6.250000000E-002, 5.469999835E-002, 0.000000000E+000])
ghia_u_val = np.array([1.000000000E+000, 8.412299752E-001, 7.887099981E-001, 7.372199893E-001, 6.871700287E-001, 2.315099984E-001, 3.319999902E-003, -1.364099979E-001, -2.058099955E-001, -2.108999938E-001, -1.566199958E-001, -1.014999971E-001, -6.434000283E-002, -4.774999991E-002, -4.191999882E-002, -3.717000037E-002, 0.000000000E+000])
ghia_v_x = np.array([ 1.0000, 0.9688, 0.9609, 0.9531, 0.9453, 0.9063, 0.8594, 0.8047, 0.5000, 0.2344, 0.2266, 0.1563, 0.0938, 0.0781, 0.0703, 0.0625, 0.0000])
ghia_v_val = np.array([ 0.00000, -0.05906, -0.07391, -0.08864, -0.10313, -0.16914, -0.22445, -0.24533, 0.05454, 0.17527, 0.17507, 0.16077, 0.12317, 0.10890, 0.10091, 0.09233, 0.00000])
# Extract centerline data from PINN predictions
x_half_idx = np.argmin(np.abs(X_grid.cpu().numpy()[0, :] - 0.5))
u_centerline_pinn = U_plot[:, x_half_idx]
y_half_idx = np.argmin(np.abs(Y_grid.cpu().numpy()[:, 0] - 0.5))
v_centerline_pinn = V_plot[y_half_idx, :]
# Plotting Comparison
fig, axs = plt.subplots(1, 2, figsize=(12, 6))
axs[0].plot(ghia_u_val, ghia_u_y, 'ro', label='Ghia et al. (1982)')
axs[0].plot(u_centerline_pinn, y_test_coords.cpu().numpy(), 'b-', label='PINN Prediction')
axs[0].set_title(f'U-velocity at X = 0.5 (Re={REYNOLDS_NUMBER})')
axs[0].set_xlabel('U-velocity')
axs[0].set_ylabel('Y-coordinate')
axs[0].grid(True)
axs[0].legend()
axs[0].set_xlim([-0.3, 1.1])
axs[0].set_ylim([0, 1])
axs[1].plot(ghia_v_x, ghia_v_val, 'ro', label='Ghia et al. (1982)')
axs[1].plot(x_test_coords.cpu().numpy(), v_centerline_pinn, 'b-', label='PINN Prediction')
axs[1].set_title(f'V-velocity at Y = 0.5 (Re={REYNOLDS_NUMBER})')
axs[1].set_xlabel('X-coordinate')
axs[1].set_ylabel('V-velocity')
axs[1].grid(True)
axs[1].legend()
axs[1].set_xlim([0, 1])
axs[1].set_ylim([-0.3, 0.3])
plt.tight_layout()
plt.savefig('centerline_comparison.png')
plt.close(fig)
print("Centerline comparison plot saved as centerline_comparison.png")
@hydra.main(config_path=".", config_name="config", version_base="1.3")
def main(cfg: DictConfig):
# Hydraの出力ディレクトリを表示
print(f"Hydra working directory: {os.getcwd()}")
print(OmegaConf.to_yaml(cfg)) # ロードされた設定を出力
# シードの設定
torch.manual_seed(cfg.seed)
np.random.seed(cfg.seed)
device = get_device()
if device.type == "cuda":
torch.cuda.manual_seed_all(cfg.seed)
torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False
# 物理定数の計算
# 運動学的粘度 nu = UL / Re
cfg.constants.nu = (cfg.constants.u_lid * cfg.constants.domain_size) / cfg.constants.reynolds_number
print(f"Calculated Kinematic Viscosity (nu): {cfg.constants.nu:.6e}")
print(f"Starting PINN training for Re={cfg.constants.reynolds_number}...")
# モデルの初期化 (ネットワーク構造はcfgから)
layer_sizes = [2] + [cfg.model.hidden_size] * cfg.model.num_hidden_layers + [3]
model = PINN(layer_sizes).to(device)
# 学習の実行
training_history_raw = train_pinn(model, cfg, device)
print("Training finished.")
# 結果の保存とプロット
save_and_plot_results(model, cfg, training_history_raw, device)
# オプション：モデルの保存
torch.save(model.state_dict(), f"pinn_model_re{cfg.constants.reynolds_number}.pth")
print("Model saved.")
if __name__ == "__main__":
main()
3.4. 設定ファイル(config.yaml)の例
# @package _global_
seed: 42
# 物理定数
constants:
reynolds_number: 100
u_lid: 1.0
domain_size: 1.0
rho: 1.0
nu: ??? # 計算によって上書きされます
# ネットワークモデルの設定
model:
num_hidden_layers: 5
hidden_size: 50
# データポイントの設定
data:
num_collocation_points: 10000
num_boundary_points: 1000
# 学習の設定
training:
# 境界損失に対するPDE損失の重み
pde_loss_weight: 100.0
adam:
epochs: 1000
lr: 0.001
lbfgs:
epochs: 500
lr: 0.5
max_iter: 50
history_size: 100
line_search_fn: "strong_wolfe"
pde_loss_weight: 1.0 # L-BFGSフェーズのPDE損失の重み (通常は1.0)
4. CPU計算とGPU計算の比較
CPUとGPUの計算環境におけるPINNsモデルの性能を比較した。
pytorchモジュールでは、CPUとGPUの両方で計算が可能であり、同一のコードで切り替えができる。
CPUでは64bit浮動小数点数（倍精度）が使用されるが、 GPUでは大部分のテンソル演算は32bit浮動小数点数（単精度）を使って行われる。
今回のテスト計算でCPUとGPUを使った場合の計算時間は以下の通りである。
| 環境 | エポック数(Adam法+LBFGS法) | 計算時間 |
|---|---|---|
| CPU (Intel i7-13700H) | 1000+500 | 41分 |
| GPU (NVIDIA RTX4060) | 1000+500 | 8分 |
⚡ 性能比較結果: GPUを使用することで、CPUの約5倍高速に計算可能
またPINNsを使って得られたキャビティ中心線上の流速分布を比較すると以下のようになった。
$v(x)$分布に若干の違いが見られるが、全体的には両者はほぼ一致していることが分かった。
5. まとめ
PINNsを使って蓋駆動キャビティ問題を解析し、流速分布と圧力分布を得た。
📋 本解析の成果
✅ PINNsにより蓋駆動キャビティ流れの解析に成功
✅ 得られた流速分布が実験データとほぼ一致
✅ GPU使用によりCPUの約5倍高速に計算可能
PINNsは物理法則を損失関数に組み込むことで、データが不足している場合でも物理的に一貫した解を得ることができる。
🚀 今後の展望: より複雑な流体問題への適用や、PINNsの性能向上に向けた研究が期待される