LoginSignup
2

More than 1 year has passed since last update.

Xcompact3d(Incompact3d)をインストールして使ってみる

Last updated at Posted at 2022-06-10

この記事は流体解析のオープンソースの一つであるXcompact3d(旧Incompact3d)をインストールして計算を走らせてみるまでの備忘録です.

Xcompact3D

Xcompact3DはImperial CollegeのTurbulence Simulation Groupが開発している高次精度差分法を用いた乱流解析ソルバーである. 2018頃までは非圧縮性用のIncompact3dと低マッハ数近似を用いた圧縮性用のQuasincompact3dが開発されていたが,統合されてXcompact3dとなっている.機能の概要は以下のとおりであり,詳細はXcompact3dBartholomew 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が16コア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ディレクトリの中に1000ステップごとのフィールドデータが生成されるが,可視化用のデータは無いようである.上記の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で読み込んで以下のような可視化を行うことができる.

vort001.png

Fig.1 1000ステップ(t=5.0)での渦度等値面

vort004.png

Fig.2 4000ステップ(t=20.0)での渦度等値面

uvec001.png

Fig.3 1000ステップ(t=5.0)での圧力と速度ベクトル

uvec004.png

Fig.4 4000ステップ(t=20.0)での圧力と速度ベクトル

まとめ

今回,Xcompact3d(Incompact3d)のインストールを行い,計算実行・可視化までを行った.Incompact3dのインストールは極めて容易であり,乱流計算・DNSを始めたい人には取りかかりやすく,ソースコードと照らし合わせて学習にもなる.ただし,マニュアル・ユーザーガイドはほとんど無いので,ソースコードの解読,参考論文の参照,解析の実行をまんべんなくできる人向けのオープンソースかもしれない.大規模計算での使い勝手や性能については今後調査してどこかで報告したい.

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2