Fortranで解いた問題をParaviewを使ってvtkでアニメーション作成

前回はFortranで2次元の熱伝導方程式を解いて、gnuplotを使って温度場の分布のアニメーションを作成したのです。

前回の記事

前回の記事で示したアニメーションはこんな感じ

カマキリ

前回はgnuplotはスクリプトを作ってアニメーションした

今回は、2次元の熱伝導方程式をFortranで解いて、温度場の分布のアニメーションはparaviewで可視化したいと思います。

本記事の内容

2次元の熱伝導方程式をFortranで解いて温度場の分布のアニメーションをparaviewで可視化する。

カマキリ

今回はparaviewを使います

↓paraviewで可視化した温度場分布のアニメーションはこちら。

※Fortranのコードとスクリプトは本記事にアップされています。
※FortranのコードはFortran90/95の書き方にしているつもりですが正確に書けているかは不明です。
※FortranはWSLを使ってインストールしています。
※paraviewはこちらのサイトよりダウンロードしてお使いください。

僕の使用環境はWindows10でWindows Subsystem for Linux（WSL）を使って、UbuntuからFortranをインストールして使用しています。

カマキリ

早く試してみたいよーった人は「まとめ」へGo!

解くべき方程式
Fortranで2次元熱伝導方程式をFortranで解く
1. 2次元熱伝導方程式の離散化
2. Fortranで解く
paraviewで温度分布のアニメーションを確認
1. vtkファイル形式として出力する部分の記述
paraviewでvtkファイルを読み込んでみよう
まとめ

解くべき方程式

まずは解くべき方程式を明らかにしておきましょう。

2次元熱伝導方程式

\begin{align*}\frac{\partial T}{\partial t}=\kappa\bigg(\frac{\partial^2T}{\partial x^2}+\frac{\partial^2T}{\partial y^2}\bigg)\tag{1}\end{align*}

熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。

今回は初期状態を以下のようにしておきます。

中央だけ中央だけ温度が高く、\(x\)方向の両端は温度０(K)の条件となっています。

また、\(x\),\(y\)の端は初期状態から変化させないような境界条件にしているので、時間発展させると中央の高い温度の部分から低い温度へ熱移動が生じて、下記のような変化を示します。

今回は、これをparaviewを使ってアニメーションを作ろうということになります。

Fortranで2次元熱伝導方程式をFortranで解く

偏微分方程式の離散化の方法は様々あるかと思いますが、ここでは離散化の方法にはこだわっていません。

主目的は「vtkファイル」をFortranで出力してparaviewで読み込むことです。

2次元熱伝導方程式の離散化

2次元熱伝導方程式

\begin{align*}\frac{\partial T}{\partial t}=\kappa\bigg(\frac{\partial^2T}{\partial x^2}+\frac{\partial^2T}{\partial y^2}\bigg)\tag{1}\end{align*}

この偏微分方程式を離散化するにあたって、

時間微分の項（左辺）は陽解法
空間微分の項（右辺）は2次精度

で離散化で進めていきましょう。

離散かされた式が↓こちらです。

\begin{align*}T^{n+1}(x_{i},y_{j})=T^{n}(x_{i},y_{j})+\kappa dt\frac{T^{n}(x_{i}+\Delta x,y_{j})-2T^{n}(x_{i},y_{j})-T^{n}(x_{i}-\Delta x,y_{j})}{\Delta x^2}\\
\kappa dt\frac{T^{n}(x_{i},y_{j}+\Delta y)-2T^{n}(x_{i},y_{j})-T^{n}(x_{i},y_{j}-\Delta y)}{\Delta y^2}\end{align*}

では、これをFortranで解きましょう(^^)/

Fortranで解く

離散化された方程式をFortranで解くためのプログラムを書いていきましょう。

このプログラムの中にvtkファイル形式で温度分布データを出力させる記述があります。

以下がプログラム全体です。

module constants
  implicit none
  integer::i,j,nt
  integer, parameter :: imax = 100,jmax=imax,ntmax = 25
  real(8):: T(0:imax,0:jmax),xmax,ymax,x,y,dx,dy,dt,kp
  integer,parameter :: SP = kind(1.0)
  integer,parameter :: DP = selected_real_kind(2*precision(1.0_SP))
end module constants

module print
  use constants
  implicit none
 
 contains
  subroutine print__profile(step,f)
    real(DP), dimension(0:imax,0:jmax), intent(in) :: f
    integer,intent(in) :: step
    character(len=40) filename
  
   write(filename,"(a,i5.5,a)") "data/temp",int(step),".vtk"
   open(10,file=filename)
   write(10,"('# vtk DataFile Version 3.0')")
   write(10,"('test')")
   write(10,"('ASCII ')")

   write(10,"('DATASET STRUCTURED_GRID')")
   write(10,"('DIMENSIONS ',3(1x,i3))") imax+1, jmax+1, 1

   write(10,"('POINTS ',i9,' float')") (imax+1)*(jmax+1)*1
   do j=0,jmax
     do i=0,imax
      write(10,"(3(f9.4,1x))") i*dx, j*dy, 0.d0
     end do
   end do

   write(10,"('POINT_DATA ',i9)") (imax+1)*(jmax+1)*1

!date input
   write(10,"('SCALARS temp float')")
   write(10,"('LOOKUP_TABLE default')")
   do j=0,jmax
     do i=0,imax
      write(10,*)  T(i,j)
     end do
   end do

   close(10)
    
  end subroutine print__profile
end module print

program thermal_explicit
  use constants
  use print

  implicit none

  dt   = 0.004
  kp   = 1
  xmax = 10
  ymax=xmax
  dx   = xmax/imax
  dy   = ymax/jmax

  
  open(10,file="initial.txt")
  open(11,file="final.txt")

   !initial state
   do j=0,jmax
     do i = 0, imax
       x = dx*i
       y = dy*j
       if (i < imax/2) then
         T(i,j) = x
       else
         T(i,j) = xmax - x
       endif

      write(10,*) ,x,y,T(i,j)
      end do
    end do
   call print__profile(nt,T(:,:))

   
   !Time step
   do nt = 1, ntmax
     do j=1,jmax-1
       do i = 1, imax-1
          T(i,j) = T(i,j) + kp*dt/(dx*dx)*(T(i+1,j)- 2*T(i,j) +T(i-1,j)) &
                      + kp*dt/(dy*dx)*(T(i,j+1)- 2*T(i,j) +T(i,j-1))
       end do
      end do
      if (mod(nt,1)==0) call print__profile(nt,T(:,:))
      write(*,*)nt
   end do
   
   !Final state
   do j=0,jmax
     do i = 0, imax
       x = dx*i
       y = dy*j
       write(11,*) x,y,T(i,j)
    end do 
   end do
end program thermal_explicit

100

101

102

103

104

105

106

module constants

implicit none

integer::i,j,nt

integer, parameter :: imax = 100,jmax=imax,ntmax = 25

real(8):: T(0:imax,0:jmax),xmax,ymax,x,y,dx,dy,dt,kp

integer,parameter :: SP = kind(1.0)

integer,parameter :: DP = selected_real_kind(2*precision(1.0_SP))

end module constants

module print

use constants

implicit none

contains

subroutine print__profile(step,f)

real(DP), dimension(0:imax,0:jmax), intent(in) :: f

integer,intent(in) :: step

character(len=40) filename

write(filename,"(a,i5.5,a)") "data/temp",int(step),".vtk"

open(10,file=filename)

write(10,"('# vtk DataFile Version 3.0')")

write(10,"('test')")

write(10,"('ASCII ')")

write(10,"('DATASET STRUCTURED_GRID')")

write(10,"('DIMENSIONS ',3(1x,i3))") imax+1, jmax+1, 1

write(10,"('POINTS ',i9,' float')") (imax+1)*(jmax+1)*1

do j=0,jmax

do i=0,imax

write(10,"(3(f9.4,1x))") i*dx, j*dy, 0.d0

end do

write(10,"('POINT_DATA ',i9)") (imax+1)*(jmax+1)*1

!date input

write(10,"('SCALARS temp float')")

write(10,"('LOOKUP_TABLE default')")

do j=0,jmax

do i=0,imax

write(10,*) T(i,j)

end do

close(10)

end subroutine print__profile

end module print

program thermal_explicit

use constants

use print

implicit none

dt = 0.004

kp = 1

xmax = 10

ymax=xmax

dx = xmax/imax

dy = ymax/jmax

open(10,file="initial.txt")

open(11,file="final.txt")

!initial state

do j=0,jmax

do i = 0, imax

x = dx*i

y = dy*j

if (i < imax/2) then

T(i,j) = x

else

T(i,j) = xmax - x

endif

write(10,*) ,x,y,T(i,j)

end do

call print__profile(nt,T(:,:))

!Time step

do nt = 1, ntmax

do j=1,jmax-1

do i = 1, imax-1

T(i,j) = T(i,j) + kp*dt/(dx*dx)*(T(i+1,j)- 2*T(i,j) +T(i-1,j)) &

+ kp*dt/(dy*dx)*(T(i,j+1)- 2*T(i,j) +T(i,j-1))

end do

if (mod(nt,1)==0) call print__profile(nt,T(:,:))

write(*,*)nt

end do

!Final state

do j=0,jmax

do i = 0, imax

x = dx*i

y = dy*j

write(11,*) x,y,T(i,j)

end do

end program thermal_explicit

２次元の熱伝導方程式をFortranで解く部分の解説は、下記の記事に詳しく乗っているので順にお読みください。

paraviewで温度分布のアニメーションを確認

ではparaviewを起動して温度分布を確認してみます。

vtkファイル形式として出力する部分の記述

vtkファイル形式で出力する部分はモジュールとして以下のように書いています。

Fortran：vtkファイル形式の出力部分

出力されたvtkの中身の基本構造は以下です。

その他、vtkファイル形式でのparaviewの読み込み方法については多数あるのですが、本記事では1例を出したにすぎません。

ご興味ある方は下記の記事を参考に色々と遊んでみてださい(^^)/

paraviewでvtkファイルを読み込んでみよう

ここからはparaviewを起動してFortranで解いた2次元熱伝導方程式による温度分布のアニメーションを見てみましょう。

↓paraviewによる基本的な操作は動画にしています。

こんな感じでFortranで解いた問題をvtkファイル形式に出力してアニメーションを見ることができます(^^)/

まとめ

色々と説明を書きましたが、以下のzipを解凍して自身の環境でお試しください(^^)/

20190429_thermal2D_paraview

以下のコマンドでFortranを実行

gfortran 2D_Thermal_Eq.f90 
./a.out

1 2	gfortran 2D_Thermal_Eq.f90 ./a.out

paraviewの操作は「paraviewでvtkファイルを読み込んでみよう」の章を参考にしてください(^^)/

数値計算のためのFortran90/95プログラミング入門(第2版)

牛島省

3,465円(11/21 14:11時点)

Amazon 楽天市場 Yahoo

Amazonの情報を掲載しています

Fortran ハンドブック

田口俊弘

3,520円(11/20 20:31時点)

Amazon 楽天市場 Yahoo

Amazonの情報を掲載しています

Fortranで解いた問題をParaviewを使ってvtkでアニメーション作成

解くべき方程式