はじめに
定常流れの数値流体計算は反復計算を行い、十分に収束させる必要がある。
この反復計算は、上手く収束する場合もあれば、発散する場合もある。
最初の100stepの残差の振る舞いから収束するか発散するか予測できないかやってみた。
環境
OS: macOS High Sierra v10.13.2
OpenFOAM v4.0
python 3.6.4
Keras 2.0.2
TensorFlow 1.0.1
ファイル構成
case/
├ converge.py (本プログラム)
├ hassan.csv (学習用発散データ)
├ hassan_test.csv (テスト用発散データ)
├ syuusoku.csv (学習用収束データ)
└ syuusoku_test.csv (テスト用収束データ)
↓github
https://github.com/matsubaraDaisuke/my_qiita/tree/master/converge
簡単な説明
機械学習に関する部分はこちらの資料を参考にした。
https://qiita.com/hiroeorz@github/items/33b85529be0829f34973
入力層はノード数：100
中間層はノード数：400、活性化関数:relu
出力層はノード数：2(収束or発散をone-hot-vectorで表現)、活性化関数:softmax
損失関数：binary_corossentropy
性能指標：accuracy
学習アルゴリズム：SGD
学習データ
p,U,kの空間に対する残差を記録し、100Stepごとにデータを分割し学習させる。
とりあえずp,U,kはゴチャ混ぜにした。
収束した場合と発散した場合でラベル付けを行う。
流体計算の対象はOpenFOAMのチュートリアルである、pitzDailyとmotorBikeとした。
発散データはmotorBikeの安定性を下げて発散させた。
なお、残差の値は直接使わずに対数Log10でとった値を用いている。
ソースコードの説明
# -*- coding: utf-8 -*-
from keras.models import Sequential
from keras.layers import Activation, Dense
from keras.utils.np_utils import to_categorical
import numpy as np
import pandas as pd
# X:[収束] => Y:0
# X:[発散] => Y:1
syuusoku_data = pd.read_csv('syuusoku.csv', header=None)
hassan_data = pd.read_csv('hassan.csv', header=None)
syuusoku_array = syuusoku_data.as_matrix()
hassan_array = hassan_data.as_matrix()
syuusoku_array = syuusoku_array.transpose()
hassan_array = hassan_array.transpose()
syuusoku_array = np.array(syuusoku_array)
hassan_array = np.array(hassan_array)
X_list = np.concatenate([syuusoku_array, hassan_array],axis=0)
syuusoku_label = np.zeros(260)
hassan_label = np.ones(20)
Y_list = np.concatenate([syuusoku_label, hassan_label],axis=0)
# kerasのmodelに渡す前にXをnumpyのarrayに変換する。
X = np.array(X_list)
Y = to_categorical(Y_list) #one-hot化
# 学習のためのモデルを作る
model = Sequential()
# 全結合層(100層->400層)
model.add(Dense(input_dim=100, output_dim=400))
# 活性化関数(ReLu関数)
model.add(Activation("relu"))
# 全結合層(400層->2層)
model.add(Dense(output_dim=2))
# 活性化関数(softmax関数)
model.add(Activation("softmax"))
# モデルをコンパイル
#model.compile(loss="categorical_crossentropy", optimizer="sgd", metrics=["accuracy"])
model.compile(loss="binary_crossentropy", optimizer="sgd", metrics=["accuracy"])
# 学習を実行
model.fit(X, Y, nb_epoch=1000, batch_size=100)
# 学習したモデルで予測する。
# [発散]=> 1([0,1]) 1番目のビットが立っている
# [収束]=> 0([1,0]) 0番目のビットが立っている
# という予測になるはず...
# 発散ケースのテスト
hassan_test_data = pd.read_csv('hassan_test.csv', header=None)
hassan_test_array = hassan_test_data.as_matrix()
hassan_test_array = hassan_test_array.transpose()
hassan_test_array = np.array(hassan_test_array)
#print(hassan_test_data)
results = model.predict_proba(hassan_test_array)
print("Predict(発散):\n", results)
# 収束ケースのテスト
syuusoku_test_data = pd.read_csv('syuusoku_test.csv', header=None)
syuusoku_test_array = syuusoku_test_data.as_matrix()
syuusoku_test_array = syuusoku_test_array.transpose()
syuusoku_test_array = np.array(syuusoku_test_array)
results1 = model.predict_proba(syuusoku_test_array)
print("Predict(収束):\n", results1)
結果
...
...
...
Epoch 992/1000
280/280 [==============================] - 0s - loss: 0.1356 - acc: 0.9571
Epoch 993/1000
280/280 [==============================] - 0s - loss: 0.1395 - acc: 0.9571
Epoch 994/1000
280/280 [==============================] - 0s - loss: 0.1336 - acc: 0.9571
Epoch 995/1000
280/280 [==============================] - 0s - loss: 0.1436 - acc: 0.9571
Epoch 996/1000
280/280 [==============================] - 0s - loss: 0.1327 - acc: 0.9571
Epoch 997/1000
280/280 [==============================] - 0s - loss: 0.1351 - acc: 0.9571
Epoch 998/1000
280/280 [==============================] - 0s - loss: 0.1382 - acc: 0.9571
Epoch 999/1000
280/280 [==============================] - 0s - loss: 0.1376 - acc: 0.9571
Epoch 1000/1000
280/280 [==============================] - 0s - loss: 0.1476 - acc: 0.9571
そこそこ良さそう。
別の計算条件に対して、収束ケースと発散ケースの最初の100stepのUx
に対して予測してみた。
Predict(発散):
[[0.04552955 0.95447046]]
Predict(収束):
[[9.9921906e-01 7.8098016e-04]]
Predict(発散)に対して[0 1]
Predict(収束)に対して[1 0]
に近い値となっており、予測ができている。
最後に
今回は最初の100stepを対象にしたが、10stepくらいで予測できたら実用的かも。
時系列データとみることもできるので、RNNの方が適切?
テストケースが少ないのでもっと検証が必要。
多分、収束しそうだけど急に発散的なものには対応してないと思う。
離散化条件や、メッシュ品質も絡めたいけど面倒なのでやってません。
やっぱり、みんなの計算条件と結果を集めるGUIや環境作りの方が大切かも。