概要
OpenFOAMの時系列データとsetsデータ(sampleコマンドで出力する空間分布のデータ)をPythonで読み込むための関数をつくって,matplotlibでのプロットを簡略化しました。
前提として,OpenFOAMの時系列データはリスタートの度に新しいディレクトリに保存されますが,それをcatコマンドでくっつけてひとつづきのファイルになってるとします。
Pythonのコード
時系列データの読み込み関数用関数
- caseはケースファイル名
- positionsには時系列データの場所の情報をリストでいれます。
- variablesが時系列データの物理量(p, Uとか)のリスト
- 各位置の時系列データを一括で読み込みます。
- 「物理量_場所」の文字列をkey,実際のデータをndarrayとしてvalueにもつディクショナリデータを返します。
def readHistory(case, positions, variables):
import numpy as np
"""
case: case name (string data)
positions: position name (list data)
variables: physical quantities (list data)
"""
data ={}
for position in positions:
for variable in variables:
key = variable+"_"+position
value = np.genfromtxt(case+"/"+position+"/"+variable, unpack=True)
data[key] = value
return data
setsデータの読み込み用関数
引数の説明
- caseはケースファイル名
- fileはsetsデータのファイル名
- variablesが出力したい物理量(p, Uとか)のリスト
- setsデータは複数の物理量のデータが1つのファイルに出力されるので,各物理量を個別に読み込みます。
- 「物理量_場所」の文字列をkey,実際のデータをdataframeとしてvalueにもつディクショナリデータを返します。
def readSets(case, file, variables):
import numpy as np
import pandas as pd
import glob
import os
"""
case: case name (string data)
file: sets data file name (string data)
variables: physical quantities (list data)
"""
# list of timing
list1 = glob.glob(case + "/sets/*")
variable_names = file.split("_")
# remove ".xy"
variable_names[-1] = variable_names[-1].split(".")[0]
position, variable_names[0] = variable_names[0], "x"
#print(variable_names)
data = {}
for variable in variables:
key = variable+"_"+position
data[key] = pd.DataFrame()
for i in list1:
path = i+"/" + file
df = pd.read_table(path, header = None, names = variable_names)
data[key] = pd.concat([data[key], df[variable]], axis=1,)
list_index=[]
for i in list1:
list_index.append(str(float(os.path.basename(i))))
data[key].columns = list_index
data[key] = pd.concat([df[variable_names[0]], data[key].sort_index(axis=1)], axis=1)
# return dictionary of dataframe
return data
matplotlibでプロット
matplotlibでデータをプロットします。
データの配置は以下のようになってるとします。
test01がケース名。
test01/position-A/p
test01/position-A/U
test01/position-B/p
test01/position-B/U
test01/sets/0/center_p_k.xy
test01/sets/10/center_p_k.xy
p, Uはposition-A, Bの時系列データ
center_p_k.xyは時刻0,10の分布データ
です。
ケースが複数あって,ケース名だけ変えたプロットを簡単に作ることを想定した作りになってます。
import matplotlib.pyplot as plt
%matplotlib inline
case = "test01"
history_position = ["position-A", "position-B"]
history_variables = ["p", "U"]
sets_variables = ["p", "k"]
sets_position = "center"
dataSet_test01= [readHistory(case, history_position, history_variables),
readSets(case, "center_p_k.xy", sets_variables)]
# 時系列データ一式のプロット
for i in history_position:
for j in history_variables:
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.plot(dataSet_test01[0][j+"_"+i][0], dataSet_test01[0][j+"_"+i][1], label=j)
ax.legend(bbox_to_anchor=(1.01, 1., 0., 0), loc='upper left', borderaxespad=0.,)
# 空間分布のプロット
for sets in sets_variables:
for column in dataSet_test01[1][sets+"_"+sets_position]:
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.plot(dataSet_test01[1][sets+"_"+sets_position]["x"],
dataSet_test01[1][sets+"_"+sets_position][column],
label=sets+"_"+column)
ax.legend(bbox_to_anchor=(1.01, 1., 0., 0), loc='upper left', borderaxespad=0.,)