はじめに

Google ColabでOpenFOAMチュートリアルを実行し、計算結果の可視化まで行う方法を記載する。

以下の記事でColabでOpenFOAMを実行する方法が紹介されているが、本記事では可視化までColabで行う。

【最速】５分で実施OpenFOAM@Google Colaboratory

OpenFOAM@Colabの利点

クラウドの壊しても問題ない環境で環境構築できる
PCが無い場合でも、スマホでOpenFOAMを体験できる
スマホのブラウザでColabノートブックにアクセスすると、スマホで計算を実行し、その場で可視化まで可能

OpenFOAM@Colabの実行

Colabにアクセスし、ノートブックを新規作成。

1. OpenFOAMのインストール

下記をセルに入力し、ColabのサーバにOpenFOAMをインストール。

%%bash
sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -"
add-apt-repository http://dl.openfoam.org/ubuntu
apt-get update
apt-get -y install openfoam9
pip3 install vtk #可視化のためにvtkをインストール

2. 計算Caseを準備

チュートリアルケースをカレントディレクトリにコピーする。
ここではCavityケースを利用する。

%%bash
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity . #チュートリアルケースをコピー

%cd ./cavity

3. 計算実行

下記セルを実行し、計算を実行する。

%%bash
. /opt/openfoam9/etc/bashrc #環境変数の読み込み
blockMesh >log.blockMesh 2>&1 #メッシュ作成
icoFoam >log.icoFoam 2>&1 #ソルバー実行
ls #計算結果が出力されていることを確認

4. 可視化

以下セルで、計算結果をvtkファイルに変換する。

%%bash
. /opt/openfoam9/etc/bashrc
foamToVTK #計算結果をVTK形式に変換

以下セルで、こちらを参考に、matplotlibでvtkファイルを可視化する関数を定義する。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import vtk
from vtk.util.numpy_support import vtk_to_numpy
from scipy.interpolate import griddata
from ipywidgets import interact

def plotVtk(t):
  filename = f"VTK/cavity_{t}.vtk"

  reader = vtk.vtkUnstructuredGridReader()
  reader.SetFileName(filename)
  reader.Update()
  data = reader.GetOutput()

  # cell data から point data変換
  cell2point = vtk.vtkCellDataToPointData()
  cell2point.SetInputData(reader.GetOutput())
  cell2point.Update()
  # 座標データの配列化
  coord = vtk_to_numpy(cell2point.GetOutput().GetPoints().GetData())

  x = coord[:,0]
  y = coord[:,1]
  z = coord[:,2] 

  # メッシュグリッド用
  xmin, xmax = min(x), max(x)
  ymin, ymax = min(y), max(y)
  xi = np.linspace(xmin, xmax, 100)
  yi = np.linspace(ymin, ymax, 100)

  # GetAbstractArray('p')は圧力、GetAbstractArray('U')は速度ベクトルデータ
  p = vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray('p'))
  u = vtk_to_numpy(cell2point.GetOutput().GetPointData().GetAbstractArray('U'))
  speed = np.sqrt(u[:,0]**2 + u[:,1]**2) # ベクトルをスカラー値に変換
  U = u[:,0] # ベクトルをスカラー値に変換
  V = u[:,1] # ベクトルをスカラー値に変換

  # 圧力のコンター図出力
  plt.axes().set_aspect('equal')
  levels = np.linspace(-5,5,21) # 描画する範囲の設定(最小，最大，色の分割数)
  plt.tricontourf(x,y,p,levels=levels,cmap="jet",vmin=-5, vmax=5) #vmin,vmax で色の範囲を設定
  cbar = plt.colorbar()
  cbar.set_ticks(np.arange(-5,5.1,0.5)) # カラーバーの目盛り
  cbar.set_label("p")
  plt.show()

  # 流速のコンター図出力
  plt.axes().set_aspect('equal')
  levels = np.linspace(-1,1,21) 
  plt.tricontourf(x,y,U,levels=levels,cmap="jet", vmin=-1, vmax=1)
  cbar = plt.colorbar()
  cbar.set_ticks(np.arange(-1.0,1.01,0.1)) 
  cbar.set_label("u")
  # plt.quiver(x,y,U,V,color='red',angles='xy',scale_units='xy', scale=10.0) #長さでベクトル場を可視化
  plt.quiver(x,y,U/speed,V/speed,speed,cmap='jet',scale=15.0) #カラーマップでベクトル場を可視化
  plt.show()

  # 流線出力
  plt.axes().set_aspect('equal')
  levels = np.linspace(0,1,21) 
  velocity = griddata((x, y), u, (xi[None,:], yi[:,None]), method='linear') # 速度ベクトルをグリッドに合わせて線形補間
  speed2 = np.sqrt(velocity[:,:,0]**2 + velocity[:,:,1]**2)
  strm = plt.streamplot(xi, yi, velocity[:,:,0], velocity[:,:,1], linewidth=1, arrowstyle="-", density=1.5)
  plt.tricontourf(x,y,speed,levels=levels,cmap="jet", vmin=0, vmax=1)
  cbar = plt.colorbar()
  cbar.set_ticks(np.arange(0,1.01,0.1)) 
  cbar.set_label("|u|")
  plt.show()

t=100時点を可視化する。

plotVtk(t=100)

ノートブックなので、インタラクティブな表示もできる。

interact(plotVtk, t=(0,100,20) )

実際のノートブックはこちら。

拡張性

簡単な計算条件の変更はsedコマンド等で可能。

!sed -i -e "25c \        value           uniform (2 0 0);" 0/U #上壁面の境界条件を2m/sに変更

まとめ

Google ColabでOpenFOAMを実行し、matplotlibで可視化を行った。
Google ColabではGUIをサポートしていないらしく、直接PawaviewやparaFoamが利用できないため、かなり面倒な方法になった。
Paraviewは偉大。

おまけ

ダムブレークのチュートリアルも同様にOpenFOAM@Colabで可視化することができた。
（DamBreakのノートブック）

参考