【スマホでCFD】OpenFOAM@Google Colaboratoryで可視化まで行う

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

以下の記事でColabでOpenFOAMを実行する方法が紹介されているが、本記事では可視化までColabで行う。
>[【最速】５分で実施OpenFOAM@Google Colaboratory](https://qiita.com/TatsuyaKatayama/items/8c5b906320db24d73949#2-openfoam%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB)

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

# OpenFOAM@Colabの実行
[Colab](https://colab.research.google.com/)にアクセスし、ノートブックを新規作成。

# 1. OpenFOAMのインストール
下記をセルに入力し、ColabのサーバにOpenFOAMをインストール。
```bash
%%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
%%bash
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity . #チュートリアルケースをコピー
```
```
%cd ./cavity
```

# 3. 計算実行
下記セルを実行し、計算を実行する。
```bash
%%bash
. /opt/openfoam9/etc/bashrc #環境変数の読み込み
blockMesh >log.blockMesh 2>&1 #メッシュ作成
icoFoam >log.icoFoam 2>&1 #ソルバー実行
ls #計算結果が出力されていることを確認
```
![CavityResult.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/d1f87d39-02d5-74c3-a0ce-2e862f76765b.png)

# 4. 可視化
以下セルで、計算結果をvtkファイルに変換する。
```bash
%%bash
. /opt/openfoam9/etc/bashrc
foamToVTK #計算結果をVTK形式に変換
```

以下セルで、[こちら](https://qiita.com/hiro_kuramoto/items/f1bb4af055dfe7e9fd7e)を参考に、matplotlibでvtkファイルを可視化する関数を定義する。
```python
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時点を可視化する。
```python
plotVtk(t=100)
```
![Cavity100p.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/2613956e-b092-40c9-a288-783157416a1f.png)
![Cavity100u.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/a6ab1278-8e1e-972a-bd9c-630f003f690b.png)
![Cavity100stream.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/638b8499-c106-b9e1-a7e2-ce705919763b.png)

ノートブックなので、インタラクティブな表示もできる。
```python
interact(plotVtk, t=(0,100,20) )
```
![CavityInteractive.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/9de3584c-9316-adb6-c9df-9c49c3bc06d5.png)

実際のノートブックは[こちら](https://colab.research.google.com/github/nagamachia/OpenFOAM/blob/main/OpenFOAM_Cavity_Colab.ipynb)。

# 拡張性
簡単な計算条件の変更はsedコマンド等で可能。
```bash
!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のノートブック](https://colab.research.google.com/github/nagamachia/OpenFOAM/blob/main/OpenFOAM_DamBreak_Colab.ipynb)）
![damBreak.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/423862/c61f7830-2176-874d-d8f2-1cd20d8fb865.gif)

# 参考
- [【最速】５分で実施OpenFOAM@Google Colaboratory - Qiita](https://qiita.com/TatsuyaKatayama/items/8c5b906320db24d73949#2-openfoam%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB)
- [OpenFOAMの結果をVTKで出力してmatplotlibで表示する - Qiita](https://qiita.com/hiro_kuramoto/items/f1bb4af055dfe7e9fd7e)
- [Jupyter LabでParaViewを扱う - ただし空気抵抗は無視しないものとする](https://inabower.hateblo.jp/entry/2019/01/05/225403)
- [Jupyter NotebookでOpenFOAMの結果をインタラクティブに可視化した - Qiita](https://qiita.com/TatsuyaKatayama/items/20063fecc7be74758735)
- [PyQt5とvtkを使ってOpenFOAMの結果を可視化してみた - Qiita](https://qiita.com/matsubara_d/items/0c1eeaf587ba7676a906)
- [K3D and Google Colab · Issue #220 · K3D-tools/K3D-jupyter · GitHub](https://github.com/K3D-tools/K3D-jupyter/issues/220)
- [matplotlib pyplotのquiver()でベクトル場を書くとクッソ見づらくなる問題の2つの解決法 - Qiita](https://qiita.com/sotatakahashi/items/fa8cbc5b3ea027ee754a)
- [ipywidgetsでグラフをインタラクティブにする - Qiita](https://qiita.com/ground0state/items/58a576cf09a56a0dd425)
- [VTK: vtkFieldData Class Reference](https://vtk.org/doc/nightly/html/classvtkFieldData.html#afb0e37337df123b5b81adb5857f63ad1)
- [キャビティ流れの解析、paraFoamの実習](http://www.opencae.or.jp/wp-content/uploads/2015/06/Course20090513.pdf)
- [Google Colaboratory上での可視化 - Mesh Wiki](rccm.co.jp/icem/pukiwiki/index.php?Google%20Colaboratory上での可視化)
- [02.OpenFOAMの自動化.ipynb - Colaboratory](https://colab.research.google.com/drive/11eRhnfUW2hPqJHtqzRS2aR3Pr91xEpsm#scrollTo=sGA6syJtNz5n)