この記事は流体解析のオープンソースの一つであるXcompact3d(旧Incompact3d)をインストールして計算を走らせてみるまでの備忘録です．
Xcompact3D
Xcompact3DはImperial CollegeのTurbulence Simulation Groupが開発している高次精度差分法を用いた乱流解析ソルバーである． 2018頃までは非圧縮性用のIncompact3dと低マッハ数近似を用いた圧縮性用のQuasincompact3dが開発されていたが，統合されてXcompact3dとなっている．機能の概要は以下のとおりであり，詳細はXcompact3dやBartholomew P. et al.で知ることができる．
- 直交格子による高次精度差分法を用いており，6次精度コンパクト差分が利用できる．
- DNS及びLESを行うことができる．
- ポアソンソルバーは3次元FFTを用いた高速解法
- MPIによる並列計算
- 物体境界は埋め込み境界法(Immersed Boundary Method, IBM)で表現できる．
- コードはFortran 90-95で記述されている．
- コードはCPUベースのスーパーコンピュータ向けとして開発されており，GPUには対応していない．
今回は一般的なLinux上でIntel Fortran Compilerを用いてコンパイル・インストールして使ってみる．
Xcompact3dのソースダウンロード
Xcompact3dはGithubにて公開されている．今回は非圧縮性用のIncompact3dを使用してみる．レポジトリにREADME.mdがあるので，それに従って作業を行う．
以降，Linux端末でのコマンド操作を順次示す．
まず，Gitレポジトリからソースコードをダウンロードする． SSH接続できる場合は以下のようにする．
$ git clone git@github.com:xcompact3d/Incompact3d.git
SSHが使えない場合はHTTP接続で持ってくる．
$ git clone https://github.com/xcompact3d/Incompact3d.git
レポジトリには以下のディレクトリが含まれている.
-
decomp2d: 2次元分割FFTライブラリ -
docs: ドキュメント．sphinxで作成する． -
examples: いくつかの例題計算 -
post_vtk: 後処理用にVTKファイルを生成するコード -
script: いくつかのpythonスクリプト -
src: 本体ソースコード
コンパイル
Xcompact3dはライブラリとして2DECOMP&FFT(2次元分割FFT)を必要とするが，上記のレポジトリに既に含まれているので，別途ダウンロード＆インストールする必要は無い．また，大規模な並列IOを可能とするライブラリADIOS2を組み入れることが可能であるが，今回はそこまで検証しないことにする．
Makefileが用意されているので，コンパイル・ビルドするためにはmakeを実行さえすればよい．make以外にはCMakeも用意されている．
Makefileはコンパイラとしてgcc(gfortran)またはintel compiler(ifort)に対応しており，インテルコンパイラを利用したい場合は以下のようにMakefile中で指定を修正する．
#=======================================================================
# Makefile for Xcompact3D
#=======================================================================
# Choose pre-processing options
# -DDOUBLE_PREC - use double-precision
# -DSAVE_SINGLE - Save 3D data in single-precision
# -DDEBG - debuggin xcompact3d.f90
# generate a Git version string
GIT_VERSION := $(shell git describe --tag --long --always)
DEFS = -DDOUBLE_PREC -DVERSION=\"$(GIT_VERSION)\"
LCL = local # local,lad,sdu,archer
IVER = 17 # 15,16,17,18
CMP = intel # intel,gcc
FFT = generic # generic,fftw3,mkl
修正・変更が必要な箇所は他に見当たらないので，makeを実行する(コンパイラは既にインストール済みとする)．
$ make
計算機によるが，コンパイルは数分で完了する．エラーが表示されずにコンパイルが終了したら，実行モジュールxcompact3dが生成される．
例題計算(Taylor Green Vortex)の実行
examplesディレクトリには以下のような例題がある．
- ABL: ?
- Channel-Flow: チャネル流
- Cylinder: 円柱周りの流れ
- Dbg-schemes: スキームデバッグ
- Lock-exchange: 密度成層流
- Mixing-layer: 混合層
- Sandbox: Jupyterを利用した解析支援ツール
- Taylor-Green-Vortex: Taylor Green Vortex
- Turbulent-Boundary-Layer: 乱流境界層
- User: ユーザー定義計算
- Wind-Turbine: 風車周りの流れ
README.mdにはTaylor-Green-Vortexが１６コアcpuで5分以下で終わりますよ，と書いてあるのでまずはじめにチュートリアルとして実行してみる．xcompact3dでは計算条件等を記入した入力ファイルとしてinput.i3dが必要であるが，上記の例題ではあらかじめ用意がされている．適当な計算ディレクトリを作成してTaylor Green Vortexの計算を行う．
$ cd ~/Incompact3d
$ mkdir run
$ cd run
$ cp -r ../examples/Taylor-Green-Vortex ./TGV
$ cd TGV
$ mpirun -np 16 ../../xcompact3d < input.i3d > TGV.log &
計算が回っている間に，入力ファイルinput.i3dの中身を調べてみる．
input.i3d
! -*- mode: f90 -*-
!===================
&BasicParam
!===================
! Flow type (1=Lock-exchange, 2=TGV, 3=Channel, 4=Periodic hill, 5=Cylinder, 6=dbg-schemes)
itype = 2
! Domain decomposition
p_row=0 ! Row partition
p_col=0 ! Column partition
! Mesh
nx=65 ! X-direction nodes
ny=65 ! Y-direction nodes
nz=65 ! Z-direction nodes
istret = 0 ! y mesh refinement (0:no, 1:center, 2:both sides, 3:bottom)
beta = 0.259065151 ! Refinement parameter (beta)
! Domain
xlx = 3.14159265358979 ! Lx (Size of the box in x-direction)
yly = 3.14159265358979 ! Ly (Size of the boy in y-direction)
zlz = 3.14159265358979 ! Lz (Size of the boz in z-direction)
! Boundary conditions
nclx1 = 1
nclxn = 1
ncly1 = 1
nclyn = 1
nclz1 = 1
nclzn = 1
! Flow parameters
iin = 1 ! Inflow conditions (1: classic, 2: turbinit)
re = 1600. ! nu=1/re (Kinematic Viscosity)
u1 = 8. ! u1 (max velocity) (for inflow condition)
u2 = 8. ! u2 (min velocity) (for inflow condition)
init_noise = 0.0 ! Turbulence intensity (1=100%) !! Initial condition
inflow_noise = 0.0 ! Turbulence intensity (1=100%) !! Inflow condition
! Time stepping
dt = 0.005 ! Time step
ifirst = 1 ! First iteration
ilast = 4000 ! Last iteration
! Enable modelling tools
ilesmod=0 ! if 0 then DNS
iscalar=0 ! If iscalar=0 (no scalar), if iscalar=1 (scalar)
iibm=0 ! Flag for immersed boundary method
! Enable io
ivisu=1 ! Store snapshots
ipost=1 ! Do online postprocessing
/End
!====================
&NumOptions
!====================
! Spatial derivatives
ifirstder = 4 ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact)
isecondder = 4 ! (1->2nd central, 2->4th central, 3->4th compact, 4-> 6th compact, 5->hyperviscous 6th)
! Time scheme
itimescheme = 5 ! Time integration scheme (1->Euler,2->AB2, 3->AB3, 4->AB4,5->RK3,6->RK4)
/End
!=================
&InOutParam
!=================
! Basic I/O
irestart = 0 ! Read initial flow field ?
icheckpoint = 1000 ! Frequency for writing backup file
ioutput = 1000 ! Frequency for visualization
ilist = 5 ! Frequency for the output
nvisu = 1 ! Size for visualisation collection
/End
!=================
&Statistics
!=================
nstat = 1 ! Size arrays for statistic collection
/End
!#######################
! OPTIONAL PARAMETERS
!#######################
!=================
&LESModel
!=================
jles=1
smagcst=0.1
walecst=0.5
iwall=0
/End
!=================
&CASE
!=================
tgv_twod = .FALSE.
/End
レイノルズ数reを変え，対応して時間刻みdt・計算ステップ数ilast・格子点数(nx,ny,nz)を変えて計算ができそうである．入力ファイル*.i3dは幾つか用意されているので比較できる．
計算が終了したら、標準出力ファイルTGV.logの末尾を確認し，計算が無事に終わっていることを確認する．
TGV.logの末尾
Good job! Xcompact3d finished successfully!
2DECOMP with p_row*p_col= 0 0
nx*ny*nz= 274625
nx,ny,nz= 65 65 65
dx,dy,dz= 4.908738521234047E-002 4.908738521234047E-002
4.908738521234047E-002
Averaged time per step (s): 2.4054682E-02
Total wallclock (s): 96.21873
Total wallclock (m): 1.603645
Total wallclock (h): 2.6727423E-02
計算終了後は以下のようなファイル・ディレクトリが存在する状態になる．
CMakeLists.txt input.i3d probes/ time_evol.dat
TGV.log input_DNS_Re1600.i3d restart0001000.info
adios2_config.xml input_ILES_Re5000.i3d restart0002000.info
checkpoint input_test.i3d restart0003000.info
data/ out/ restart0004000.info
dataディレクトリの中に１０００ステップごとのフィールドデータが生成されるが，可視化用のデータは無いようである．上記のpost_vtkにあるコードを活用して可視化を試みる．
VTKファイルの生成
post_vtkディレクトリの中身を確認する．
$ cd ~/Incompact3d/post_vtk
$ ls
CMakeLists.txt paraview_incompact3d_vtk.f90
paraview_incompact3d_vtk.f90が使えそうなので，コードの中身を見てみる(READMEなどは無いので)．
paraview_incompact3d_vtk.f90
!Copyright (c) 2012-2022, Xcompact3d
!This file is part of Xcompact3d (xcompact3d.com)
!SPDX-License-Identifier: BSD 3-Clause
!##############################################################################
!!
!! PROGRAM: visu_vtk
!! DESCRIPTION: Create vtk/vtr files from XC3D output. Output file is txt
!! always with Float32 precision regardless of the precision used
!! to write Xc3d output
!!
!! AUTHOR:
!! MODIFIED: Stefano Rolfo
!!
!##############################################################################
program visu_vtk
implicit none
integer(4) :: nx,ny,nz
real(4) :: xlx,yly,zlz,dt,dx,dy,dz
real(8),dimension(:,:,:),allocatable :: input_var_64
real(4),dimension(:,:,:),allocatable :: input_var_32
real(4),dimension(:,:,:),allocatable :: output_var_32
integer :: i,j,k,count,nfil, ivar
integer(4) :: istret,nclx, ncly, nclz, nvar
integer(4) :: nfiles, icrfile, file1, filen, ifile
integer(4) :: dig1, dig2, dig3, dig4, dig5, dig6, dig7
integer(4) :: step
real(4),dimension(:),allocatable :: y1
real(4),dimension(:),allocatable :: y2
real(4),dimension(:),allocatable :: y3
real(4) :: pi,xa
character(len=128) :: input_file
character(len=128) :: output_file
integer(4) :: idigits, iprec
integer :: length_input_file
integer(4) :: iparaview
character(len=128), dimension(:), allocatable :: var_name
!IF THE DATA ARE STORED WITH 3 DIGITS, IE UX001,UX002,ETC.
character(3) :: chits3
!IF THE DATA ARE STORED WITH 4 DIGITS, IE UX0001,UX0002,ETC.
character(4) :: chits4
!IF THE DATA ARE STORED WITH 5 DIGITS, IE UX00001,UX00002,ETC.
character(5) :: chits5
!IF THE DATA ARE STORED WITH 5 DIGITS, IE UX00001,UX00002,ETC.
character(6) :: chits6
!IF THE DATA ARE STORED WITH 5 DIGITS, IE UX00001,UX00002,ETC.
character(7) :: chits7
write (*,*) 'nx, ny, nz - Incompact3D'
read (*,*) nx, ny, nz
write (*,*) 'xlx, yly, zlz - Incompact3D'
read (*,*) xlx, yly, zlz
write (*,*) 'nclx, ncly, nclz - Incompact3D'
read (*,*) nclx, ncly, nclz
write (*,*) 'n files, first file, last file'
read (*,*) nfiles,file1, filen
write (*,*) 'Stretching in the y direction (Y=1/N=0)?'
read (*,*) istret
write (*,*) 'Number of digits in the file name (Allowed 3, 4, 5, 6, 7)?'
read (*,*) idigits
if (idigits < 3 .or. idigits > 7) then
write(*,*) 'Numbers of digits has to be between 3 and 7 '
stop
endif
write (*,*) 'Writing precision for the XC3D postprocessing (Single = 4, Double = 8)?'
read (*,*) iprec
print *, 'iprec =', iprec
write (*,*) 'Paraview version target (i.e. 5.8 => 58)'
read (*,*) iparaview
write (*,*) 'How many variables you would like to export'
read (*,*) nvar
allocate(var_name(nvar))
do ivar = 1, nvar
write(*,*) 'XC3D variable name (i.e. ux)'
read(*,*) var_name(ivar)
enddo
do ivar = 1, nvar
write(*,*) 'Var list ', trim(var_name(ivar))
enddo
! Allocate of the in/out variables
allocate(input_var_64(nx,ny,nz))
allocate(input_var_32(nx,ny,nz))
allocate(output_var_32(nx,ny,nz))
if (nclx==0) dx=xlx/real(nx,4)
if (nclx==1 .or. nclx==2) dx=xlx/real(nx-1.,4)
if (ncly==0) dy=yly/real(ny,4)
if (ncly==1.or.ncly==2) dy=yly/real(ny-1.,4)
if (nclz==0) dz=zlz/nz
if (nclz==1.or.nclz==2) dz=zlz/real(nz-1.,4)
allocate(y1(nx))
allocate(y2(ny))
allocate(y3(nz))
do i=1,nx
y1(i)=(i-1)*dx
enddo
if (istret==1) then
print *,'We need to read the yp.dat file'
open(12,file='yp.dat',form='formatted',status='unknown')
do j=1,ny
read(12,*) y2(j)
enddo
close(12)
else
do j=1,ny
y2(j)=(j-1)*dy
enddo
endif
do k=1,nz
y3(k)=(k-1)*dz
enddo
step = int(filen/nfiles)
! Time loop to write the variables
do ifile = file1, filen, step
! Get the digits for the file name
if (idigits == 3) then
dig1 = ifile/100 + 48
dig2 = ( ifile - 100*( ifile/100 ) )/10 + 48
dig3 = ( ifile - 10*( ifile/10 ) )/1 + 48
chits3(1:3) = char(dig1)//char(dig2)//char(dig3)
write(*,*) 'File name has structure var_name'//chits3
output_file = 'XC3D'//chits3
elseif (idigits == 4) then
dig1 = ifile/1000 + 48
dig2 = ( ifile - 1000*( ifile/1000 ) )/100 + 48
dig3 = ( ifile - 100*( ifile/100 ) )/10 + 48
dig4 = ( ifile - 10*( ifile/10 ) )/1 + 48
chits4(1:4) = char(dig1)//char(dig2)//char(dig3)//char(dig4)
write(*,*) 'File name has structure var_name'//chits4
output_file = 'XC3D'//chits4
elseif (idigits == 5) then
dig1 = ifile/10000 + 48
dig2 = ( ifile - 10000*( ifile/10000 ) )/1000 + 48
dig3 = ( ifile - 1000*( ifile/1000 ) )/100 + 48
dig4 = ( ifile - 100*( ifile/100 ) )/10 + 48
dig5 = ( ifile - 10*( ifile/10 ) )/1 + 48
chits5(1:5) = char(dig1)//char(dig2)//char(dig3)//char(dig4)//&
char(dig5)
write(*,*) 'File name has structure var_name'//chits5
output_file = 'XC3D'//chits5
elseif (idigits == 6) then
dig1 = ifile/100000 + 48
dig2 = ( ifile - 100000*( ifile/100000 ) )/10000 + 48
dig3 = ( ifile - 10000*( ifile/10000 ) )/1000 + 48
dig4 = ( ifile - 1000*( ifile/1000 ) )/100 + 48
dig5 = ( ifile - 100*( ifile/100 ) )/10 + 48
dig6 = ( ifile - 10*( ifile/10 ) )/1 + 48
chits6(1:6) = char(dig1)//char(dig2)//char(dig3)//char(dig4)//&
char(dig5)//char(dig6)
write(*,*) 'File name has structure var_name'//chits6
output_file = 'XC3D'//chits6
else
dig1 = ifile/1000000 + 48
dig2 = ( ifile - 1000000*( ifile/1000000 ) )/100000 + 48
dig3 = ( ifile - 100000*( ifile/100000 ) )/10000 + 48
dig4 = ( ifile - 10000*( ifile/10000 ) )/1000 + 48
dig5 = ( ifile - 1000*( ifile/1000 ) )/100 + 48
dig6 = ( ifile - 100*( ifile/100 ) )/10 + 48
dig7 = ( ifile - 10*( ifile/10 ) )/1 + 48
chits7(1:7) = char(dig1)//char(dig2)//char(dig3)//char(dig4)//&
char(dig5)//char(dig6)//char(dig7)
write(*,*) 'File name has structure var_name'//chits7
output_file = 'XC3D'//chits7
endif
! Creation of the file output for the time step
write(*,*) 'output_file ',output_file
nfil=41
if (iparaview >= 58 ) then
output_file=trim(output_file)//'.vtr'
else
output_file=trim(output_file)//'.vtk'
endif
write(*,*) 'output_file ',output_file
open(nfil,file=output_file)
write(nfil,*)'<VTKFile type="RectilinearGrid" version="0.1"',' byte_order="LittleEndian">'
write(nfil,*)' <RectilinearGrid WholeExtent=','"1 ',nx,' 1 ',ny,' 1 ',nz,'">'
write(nfil,*)' <Piece Extent=','"1 ',nx,' 1 ',ny,' 1 ',nz,'">'
write(nfil,*)' <Coordinates>'
write(nfil,*)' <DataArray type="Float32"',' Name="X_COORDINATES"',' NumberOfComponents="1">'
write(nfil,*) (y1(i),i=1,nx)
write(nfil,*)' </DataArray>'
write(nfil,*)' <DataArray type="Float32"',' Name="Y_COORDINATES"',' NumberOfComponents="1">'
write(nfil,*) (y2(j),j=1,ny)
write(nfil,*)' </DataArray>'
write(nfil,*)' <DataArray type="Float32"',' Name="Z_COORDINATES"',' NumberOfComponents="1">'
write(nfil,*) (y3(k),k=1,nz)
write(nfil,*)' </DataArray>'
write(nfil,*)' </Coordinates>'
! Now we can add the variables
write(nfil,*)' <PointData Scalars="scalar">'
do ivar = 1, nvar
if (idigits == 3) then
input_file = trim(var_name(ivar))//'-'//chits3//'.bin'
elseif (idigits == 3) then
input_file = trim(var_name(ivar))//'-'//chits4//'.bin'
else
input_file = trim(var_name(ivar))//'-'//chits5//'.bin'
endif
length_input_file=len_trim(input_file)
write(*,*) 'length_input_file',length_input_file
write(*,*) 'File to open ', input_file
open(10,file=input_file,form='unformatted',&
access='direct', recl=8, status='old')
count = 1
do k=1,nz
do j=1,ny
do i=1,nx
if (iprec == 8) then
read(10,rec=count) input_var_64(i,j,k)
else
read(10,rec=count) input_var_32(i,j,k)
endif
count = count + 1
enddo
enddo
enddo
close(10)
output_var_32=0.
if (iprec == 8) then
output_var_32=real(input_var_64,4)
else
output_var_32=input_var_32
endif
write(nfil,*)' <DataArray Name="',trim(var_name(ivar)),'" type="Float32"', &
' NumberOfComponents="1"',' format="ascii">'
write(nfil,*) (((output_var_32(i,j,k),i=1,nx),j=1,ny),k=1,nz)
write(nfil,*)' </DataArray>'
enddo
! End of writing data
write(nfil,*)' </PointData>'
! Write the end of the file
write(nfil,*)' </Piece>'
write(nfil,*)' </RectilinearGrid>'
write(nfil,*)'</VTKFile>'
close(nfil)
enddo
deallocate(y1)
deallocate(y2)
deallocate(y3)
deallocate(input_var_64)
deallocate(input_var_32)
deallocate(output_var_32)
end program visu_vtk
プログラム実行時に，read文でパラメータをいくつか指定しないとダメなようなので，計算入力ファイルinput.i3dと計算で出力されたデータを参考にして，以下のようなvtk生成プログラム用の入力ファイルinput.pvを作成し，~/Incompact3d/run/TGV/data/に置いておく．
input.pv
65 65 65 ! nx, ny, nz
3.14159265358979 3.14159265358979 3.14159265358979 ! xlx, yly, zlz
1 1 1 ! nclx, ncly, nclz
4 1 4 ! nfiles, file1, filen
0 ! istret
3 ! idigits
8 ! iprec
59 ! iparaview
6 ! nvar
ux
uy
uz
pp
vort
critq
準備ができたので，まずparaview_incompact3d_vtk.f90をコンパイルして実行モジュールを作成する．単一のプログラムなので，Makefileを用意せず単純にコンパイルする．
$ cd ~/Incompact3d/post_vtk/
$ ifort -assume byterecl -o pv_i3d paraview_incompact3d_vtk.f90
コンパイルオプション-assume bytereclは必ず指定すること．このプログラムでは，計算出力ファイルを開く際のopen文でダイレクトアクセスを行いrecl=8としてレコード長を指定しているが，インテルコンパイラのデフォルトの設定ではレコード長の単位が"ワード(4バイト)"であり，このレコード長単位をバイトに指定するコンパイルオプションを付けないと，レコード長が4*8=32となって実行時エラーとなる．
コンパイルが完了して実行モジュールpv_i3dが生成されたら，計算出力データがあるディレクトリでプログラムを実行する．
$ cd ~/Incompact3d/run/TGV/data
$ ~/Incomact3d/post_vtk/pv_i3d < input.pv > pv_i3d.log &
実行が終わると，4ステップ分のvtkファイルXC3D000[1-4].vtrが生成されるので，これらをParaViewで読み込んで以下のような可視化を行うことができる．
Fig.1 1000ステップ(t=5.0)での渦度等値面
Fig.2 4000ステップ(t=20.0)での渦度等値面
Fig.3 1000ステップ(t=5.0)での圧力と速度ベクトル
Fig.4 4000ステップ(t=20.0)での圧力と速度ベクトル
まとめ
今回，Xcompact3d(Incompact3d)のインストールを行い，計算実行・可視化までを行った．Incompact3dのインストールは極めて容易であり，乱流計算・DNSを始めたい人には取りかかりやすく，ソースコードと照らし合わせて学習にもなる．ただし，マニュアル・ユーザーガイドはほとんど無いので，ソースコードの解読，参考論文の参照，解析の実行をまんべんなくできる人向けのオープンソースかもしれない．大規模計算での使い勝手や性能については今後調査してどこかで報告したい．