可視化（アニメ）¶

惑星座標を使ってアニメーションを作ってみよう．アニメーションであれば視覚的に動きを捉えることができる．

ParaViewの使用¶

ParaViewはオープンソース，マルチプラットフォームの可視化ソフトウェアである．大規模なデータにも対応できるようである．

>>ParaView のDownloadからインストーラがダウンロードできる． Windows,MacOS,Linuxで使える．

ParaView向けのデータ出力¶

ParaViewは様々な形式のデータを扱うことができるが，VTK XMLフォーマットを用いる．中でもvtkUnstructuredGrid形式（拡張子は「.vtu」）を用いる．以下のような形式である．タグで囲んだ間に粒子の座標や速度などの情報を書き入れる．

<VTKFile type=”UnstructuredGrid” ...>
<UnstructuredGrid>
<Piece NumberOfPoints=”#” NumberOfCells=”#”>
<PointData>...</PointData>
<CellData>...</CellData>
<Points>...</Points>
<Cells>...</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>

Note

〔粒子法入門，越塚誠一ら（丸善出版，2014）〕

ということで，プログラミングコードの中にvtkUnstructuredGrid形式を出力できるサブルーチンを加える．

module output

use globals

implicit none

contains
    subroutine write_vtk(step,c,v,m)
        integer,intent(in)::                step
        double precision,intent(in)::       c(ndim,nbody)
        double precision,intent(in)::       v(ndim,nbody)
        double precision,intent(in)::       m(nbody)

        integer::                   unit_vtu=11
        integer                     i,j
        double precision            v2
        double precision            ct(3)
        character(len=27)           file_vtu
        character(len=20)           step_str

        write(step_str,'(i20)')step
        file_vtu = step_str // '.vtu'
        file_vtu = 'pv_' // adjustl(file_vtu)
        file_vtu = trim(adjustl(file_vtu))
        open(unit=unit_vtu,file=file_vtu,status='replace')
        write(unit_vtu,'(a)')"<?xml version='1.0' encoding='UTF-8'?>"
        write(unit_vtu,'(a)')"<VTKFile xmlns='VTK' byte_order='LittleEndian' version='0.1' type='UnstructuredGrid'>"
        write(unit_vtu,'(a)')"<UnstructuredGrid>"
        write(unit_vtu,'(a,i0,a,i0,a)')"<Piece NumberOfCells='",nbody,"' NumberOfPoints='",nbody,"'>"
        write(unit_vtu,'(a)')" "

        write(unit_vtu,'(a)')"<Points>"
        write(unit_vtu,'(a,i0,a)')"<DataArray NumberOfComponents='",3,"' type='Float32' Name='Position' format='ascii'>"
        do i = 1,nbody
            ct(1:ndim) = c(1:ndim,i)
            do j = 1,3
                if(ndim .eq. 2) then
                    ct(3) = 0.0d0
                end if
                write(unit_vtu,'(f8.3)',advance='no')ct(j)
            end do
            write(unit_vtu,*)
        end do
        write(unit_vtu,'(a)')"</DataArray>"

        write(unit_vtu,'(a)')"</Points>"
        write(unit_vtu,'(a)')"<PointData>"
        write(unit_vtu,'(a)')"<DataArray NumberOfComponents='1' type='Float32' Name='Velocity' format='ascii'>"
        do i = 1,nbody
            v2 = 0.0d0
            do j = 1,ndim
                v2 = v2 + v(j,i)**2
            end do
            write(unit_vtu,'(f8.3)')sqrt(v2)
        end do
        write(unit_vtu,'(a)')"</DataArray>"

        write(unit_vtu,'(a)')"<DataArray NumberOfComponents='1' type='Float32' Name='Mass' format='ascii'>"
        do i = 1,nbody
            write(unit_vtu,'(f8.3)')log10(m(i)*10.0d8)
        end do
        write(unit_vtu,'(a)')"</DataArray>"
        write(unit_vtu,'(a)')"</PointData>"
        write(unit_vtu,'(a)')"<Cells>"
        write(unit_vtu,'(a)')"<DataArray type='Int32' Name='connectivity' format='ascii'>"
        do i = 1,nbody
            write(unit_vtu,'(i0)')i-1
        end do
        write(unit_vtu,'(a)')"</DataArray>"

        write(unit_vtu,'(a)')"<DataArray type='Int32' Name='offsets' format='ascii'>"
        do i = 1,nbody
            write(unit_vtu,'(i0)')i
        end do
        write(unit_vtu,'(a)')"</DataArray>"

        write(unit_vtu,'(a)')"<DataArray type='Int32' Name='types' format='ascii'>"
        do i = 1,nbody
            write(unit_vtu,'(i0)')1
        end do
        write(unit_vtu,'(a)')"</DataArray>"

        write(unit_vtu,'(a)')"</Cells>"
        write(unit_vtu,'(a)')"</Piece>"
        write(unit_vtu,'(a)')"</UnstructuredGrid>"
        write(unit_vtu,'(a)')"</VTKFile>"
        close(unit_vtu)
    end subroutine write_vtk
end module output

40行目「write(unit_vtu,’(f8.3)’,advance=’no’)ct(j)」は書式指定付出力である．「’(f8.3)’」はct(j)を文字幅8，小数点以下3桁の実数で印字せよの意味である．つまり「_:_-1.263」「_:_:_0.850」のような数が印字される．これまで数字を出力する際にこの書式のことは気にせず「*」としてきた．この場合は適当な書式（桁数等）がデータの型に応じて自動的に選択される．倍精度実数を使って計算をしていても図のプロットやアニメの作成にはそれほど多くの桁数は必要ない．出力ファイルが大きくなるだけである．

Note

advance=’no’は印字後に改行するなという命令．

Note

二次元の計算しかしていないのだが，3次元の座標データがないとParaViewがうまく使えなかった．34-40行目あたりでz座標0.0を無理矢理書き入れている．

このサブルーチンをメインループの中で時々呼び出す．

if(mod(i-1,paraviewout) .eq. 0) then
    call write_vtk(i,coord,veloc,mass)
end if

paraviewoutは出力間隔である．呼び出した瞬間の位置情報などがvtuファイルに書きこまれる．出力数に応じた複数のvtuファイルが生成される．

ParaViewでアニメ作成¶

全く使いこなせていない．使用方法はwebで検索してなんとかしよう．以下，メモ

Point Spliteで粒子を描きたければ，Tools→Manage Plugins→PointSplite_pluginをLoadしておく必要がある．Auto Loadにしておけば次回起動時はLoad不要．
File→Openで生成したvtuファイルを選択．複数ファイルを一括で選択できる．
Sources→Axesで軸を追加．
Tools→Lock View Size Customで描画領域を設定．Youtubeに投稿するために1280×720に設定したが，Save animationで893×543に調整されてしまう？

Paraviewの使用画面は以下のような感じ．

作成した動画をYoutubeに掲載した． >>Youtube1

粒子の軌跡を追加する場合

File→Openで生成したvtuファイルを選択
PropertiesのRepresentationはPointにする（Point Spliteは使わない）．
PointにしてからApplyを押すとFiltersが使用可能になる．
Filters→Alphabetical→Temporal Particles To Pathlinesを選択する
Pipeline Browserに追加されたPathlinesの目玉マークを押す．
PropertiesのMask Pointsを1に．1粒子ごとに（つまり全ての粒子の）軌跡を描くの意味．10だと10粒子ごとに1本の線になる．
Max Track Lengthは250位に，ながければ長い軌跡になる．
これで再生すると軌跡だけが書けるはず．Particlesはよくわからないので，描画しない（目玉マークを押さない）

きれいな粒子も描きたければ．．．

もう一度データを読み込む．(File→Openで生成したvtuファイルを選択)
今度はPropertiesのRepresentationをPoint Spliteにする．

以下のようなパイプラインになった．

作成した動画をYoutubeに掲載した． >>Youtube2

月が地球の軌道の外側と内側を交互に移動しながら周回するため2番目の軌跡の色が縞状になっている．

可視化（アニメ）¶

ParaViewの使用¶

ParaView向けのデータ出力¶

ParaViewでアニメ作成¶

Table Of Contents

Previous topic

Next topic

This Page

Navigation

可視化（アニメ）¶

ParaViewの使用¶

ParaView向けのデータ出力¶

ParaViewでアニメ作成¶

Table Of Contents

Previous topic

Next topic

This Page

Quick search

Navigation