２次元熱拡散シミュレーション結果をPython VTKで出力する

はじめに

熱拡散シミュレーションで得たndarray形式の2次元格子の温度データを、VTKデータで出力する方法を記載します。

方法

Python VTKライブラリを利用して、StructuredGrid形式で出力します。

環境　

• Windows 10 home（FEniCSの計算ではWSLを利用）
• Anaconda(python 3.7.6)
• FEniCS 2018

シミュレーション結果

温度拡散シミュレーションは、次の熱拡散方程式をPythonインターフェースのFEniCSを利用して、有限要素法で計算しました。
$$
\frac{\partial}{\partial t} u = \alpha\nabla ^2 u
$$
$u$:温度スカラー $\alpha$:熱拡散係数です。計算結果は、次の図のようになります。

図：計算結果

FEniCSは、以下の記事を参考にしてください。
https://qiita.com/matsxxx/items/b696d2fe730be1683d8d
今回利用したシミュレーション結果を作るコードは、以下になります。

import fenics as fe
import numpy as np

tol = 1e-14 #
Lx = 80 #x方向長さ
Ly = 100 #y方向長さ
Nx = 80 #x軸セル数
Ny = 100 #y軸セル数
T = 10.0 #計算する時間
num_steps = 10 #時間ステップ
dt = T/num_steps #時間刻み
alpha = 100.#熱拡散係数

# メッシュ作成
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(Lx, Ly), Nx, Ny)
V = fe.FunctionSpace(mesh, "P", 1)#有限要素空間

# fenicsの配列を並べ替えるのに利用する。
vertex_to_dof_map = fe.vertex_to_dof_map(V)

# 境界条件の位置設定
def boundary_right(x, on_boundary):
    #x[0]:x座標、x[1]:y座標
    return on_boundary and fe.near(x[0], Lx, tol)


def boundary_top(x, on_boundary):
    #x[0]:x座標、x[1]:y座標
    return on_boundary and fe.near(x[1], Ly, tol)

# 第一種境界条件の設定
bc_right = fe.DirichletBC(V, fe.Constant(100), boundary_right)#Lx側を100℃
bc_top = fe.DirichletBC(V, fe.Constant(80), boundary_top)#Ly側を80℃

# 試行/試験関数の設定
u = fe.TrialFunction(V) #解きたい変数
v = fe.TestFunction(V) #重み関数

# 初期条件
u_0 = fe.Expression('0',degree=2)#初期温度0
u_n = fe.interpolate(u_0, V)

# 変分法の熱拡散方程式
F= u * v * fe.dx + alpha * dt * fe.dot(fe.grad(u), fe.grad(v)) * fe.dx - (u_n) * v * fe.dx 
a, L = fe.lhs(F), fe.rhs(F)

# 変数
u = fe.Function(V)
t = 0#時間

# 非定常計算の実行
for n in range(num_steps):
    #時間更新
    t += dt
    
    #方程式を解く
    fe.solve(a == L, u, [bc_right, bc_top])
    
    #解を更新
    u_n.assign(u)

# 最終結果をndarrayで取り出す。
u_np_temp = u.vector().get_local()
# fenics内の変数の並びからT[0,0], T[0,1]・・・T[Nx-1,Ny-2],T[Nx-1,Ny-1]の並びに変換
u_np = np.zeros_like(u_np_temp)
for i,v in enumerate(vertex_to_dof_map):
    u_np[i] = u_np_temp[v]
    
# npyで保存
np.save("solution", u_np)

# fenicsにはVTKデータで出力することが可能なので、通常はこちらを使う。
# vtk = fe.File("output.pvd")
# vtk << u

VTKライブラリを用いた出力

ndarray形式のシミュレーション結果を、VTKライブラリのStructuredGrid形式で出力します。

import numpy as np
import vtk
Nx = 80
Ny = 100

# シミュレーション結果読み込み
u = np.load("solution.npy")

# 点データの作成
points = vtk.vtkPoints()

for iy in range(Ny+1):
    for ix in range(Nx+1):
        x = ix 
        y = iy
        points.InsertNextPoint(x,y,0)#2次元なので、z座標は0にする
    
# 点データをstructuredGridにセットする
structured_mesh = vtk.vtkStructuredGrid()
structured_mesh.SetExtent(0, Nx, 0, Ny, 0, 0)
structured_mesh.SetPoints(points)

# 温度データをvtk arrayに入力する
temp_value = vtk.vtkDoubleArray()
temp_value.SetName("Temperature")
for t in u:
    temp_value.InsertNextValue(t)
# StructuredGridのPointDataに温度データを追加
structured_mesh.GetPointData().AddArray(temp_value)

# データ出力
writer = vtk.vtkXMLDataObjectWriter()
writer.SetFileName("temperature.vts")
writer.SetInputData(structured_mesh)
writer.Update()

おまけ：PyEVTKを利用して出力

VTKデータで出力するだけなら、PyEVTKを利用したほうが簡単です。

import numpy as np 
import evtk

Nx = 80
Ny = 100
u = np.load("solution.npy")
u = u.reshape(Ny+1, -1)#x,yの並び注意
u = u[:,:,np.newaxis]

# meshgridで座標作成
x = np.arange(0, Nx+1)
y = np.arange(0, Ny+1)
X, Y = np.meshgrid(x,y)
X = X[:,:,np.newaxis]
Y = Y[:,:,np.newaxis]

Z = np.zeros_like(X)#２次元のため、Z座標は0

evtk.hl.gridToVTK("temperature_evtk", X, Y, Z, pointData={"temperature":u})