20220423_オープンCAE関西_vtkFortran

Transcript

1. Fortranで美しいフラクタル構造を描く 2022年4月24日
2. 目次 1. モチベーション 2. フラクタルとは 3. 代表的なフタクタル構造 4. Fortranコードで実装 5. 可視化(VTK) 6. 今後 Fortranへの書き直し VTKFortranライブラリの使用
3. モチベーション 学生時代にFortranで数値計算をしていた。 最近「C＆Fortran 演習で学ぶ数値計算」が出てプログラミング勉強に火が付いた。
4. フラクタルとは 自然科学・工学は数式で現象を理解してきた。 しかし、従来のアプローチでは解明できない現象を複雑系と呼ばれる。 流体現象 生体系の現象 経済の現象 •現象自体が理解されていない ⇒数理モデルにすることができない。 •多くの現象の複合 ⇒多くの微分方程式を解く必要があ るため複雑で解けない。 （解いても理解が難しい） 乱流遷移 DANの配列 株価の変動
5. フラクタル構造 そうは言っても何かパターンがある。 べき乗 𝑓(𝑥) = 𝑎𝑥−𝑡 ｘ 𝑓 𝑎𝑥−𝑡 𝑥を𝑏倍する→𝑏𝑥 ⇔ 𝑓(𝑏𝑥) = 𝑎𝑏−𝑡𝑥−𝑡 = 𝑐𝑥−𝑡 パターン構造 𝑏ｘ 𝑓 𝑐𝑥−𝑡 x軸を拡大 元と同じ形になった。 元と同じ形
6. 実現象でのパターン構造 乱流現象 地震規模と頻度 山崩れ規模と頻度 単語の頻度 動物の斑点模様 リアス式海岸 葉のフラクタル構造 人の介在無しに組織が自己形成をしている。 自己組織化という。 ミクロな視点では、ある決まったルールに従ってマクロな視点ではパターンを形成している。 決まったルールを解明できれば、現象を数理モデルで再現できる。 大小さまざまな渦構造 コルモゴロフ則 http://www.pricewalk.net/concept.html
7. 代表的なフタクタル構造 Fortran Paraview ソルバ ポスト マンデルブロー集合 ジュリア集合 バーンズリーのシダ
8. 「マンデルブロー集合」のアルゴリズム マンデルブロー集合 𝑧𝑛+1 = 𝑧𝑛 2 + 𝐶 初期値𝑧0 で定義される漸化式において𝑛 → ∞の極限で 𝑧𝑛 が発散しない複素数𝐶 = 𝑎 + 𝑏𝑖全体の集合 【例】 𝑧1 = 𝑧0 2 + 𝐶 = 𝐶 𝑧2 = 𝑧1 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 𝑧3 = 𝑧2 2 + 𝐶 = 𝑧1 2 + 𝐶 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 2 + 𝐶 ・・・ 𝑎 𝑏 複素数平面上 発散=0 収束=1 発散=0 収束=1 漸化式
9. 「ジュリア集合」のアルゴリズム ジュリア集合 𝑧𝑛+1 = 𝑧𝑛 2 + 𝐶 𝐶を固定。漸化式において𝑛 → ∞の極限で 𝑧𝑛 が発散 しない複素数初期値𝑧0 = 𝑎 + 𝑏𝑖全体の集合 【例】 𝑧1 = 𝑧0 2 + 𝐶 = 𝐶 𝑧2 = 𝑧1 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 𝑧3 = 𝑧2 2 + 𝐶 = 𝑧1 2 + 𝐶 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 2 + 𝐶 ・・・ 𝑎 𝑏 複素数平面上 発散=0 収束=1 漸化式
10. 「バーンズリーのシダ」のアルゴリズム バーンズリーのシダ 新たな(𝒙, 𝒚)座標 確率1% 確率85% 確率7% 確率7% 𝑥 𝑦 𝑥, 𝑦 = (0,0) 1 2 3 𝒇𝒊 𝒙𝒊 , 𝒚𝒊 = 𝟏. 𝟎 プロットする
11. Fortranコードで実装 マンデルブロー集合 ジュリア集合 バーンズリーのシダ https://takun- physics.net/13663 https://note.com/kamakiriphysic s/n/n11859282fa79 https://takun- physics.net/13730/
12. 可視化(VTK) !--------------------------------------------------------------------------------- -- !output of movie file if ( mod(m,1000)==0 ) then write(filename,"(a,i5.5,a)") "data/barsley_",int(m),".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))") Nx, Ny, 1 write(10,"('POINTS ',i9,' float')") Nx*Ny*1 do j=0,Ny-1 do i=0,Nx-1 write(10,"(3(f9.4,1x))") x_(i), y_(j), 0.d0 end do end do write(10,"('POINT_DATA ',i9)") Nx*Ny*1 !date input write(10,"('SCALARS barsley float')") write(10,"('LOOKUP_TABLE default')") do j=0,Ny-1 do i=0,Nx-1 write(10,*) f(i,j) end do end do close(10) end if end do # vtk DataFile Version 3.0 test ASCII DATASET STRUCTURED_GRID DIMENSIONS 512 512 1 POINTS 262144 float -3.0000 -3.0000 0.0000 -2.9883 -3.0000 0.0000 -2.9766 -3.0000 0.0000 -2.9648 -3.0000 0.0000 -2.9531 -3.0000 0.0000 -2.9414 -3.0000 0.0000 -2.9297 -3.0000 0.0000 -2.9180 -3.0000 0.0000 -2.9062 -3.0000 0.0000 POINT_DATA 262144 0.0000000000000000 0.0000000000000000 0.0000000000000000 ③ ② ① ① ③ ② ④ ④ ⑥ ⑤ ⑤ ⑥ • https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf • https://itandcfd.com/vtk_file_format/241/ data/ファイル名_0001.vtk
13. 今後  Fortranへの書き直し integer, parameter :: Nmax=500000, Nx=512, Ny=512 double precision, parameter :: Lx=6.d0, Ly=6.d0, dx=Lx/dble(Nx), dy=Ly/dble(Ny) complex(kind(0d0)), parameter :: uso=(0.d0,1.d0) integer :: i, j, n, m, ix, iy real(kind=8)::random_value double precision:: x, y double precision, allocatable, dimension(:) :: x_, y_ double precision, allocatable, dimension(:,:) :: f character(len=40) filename ・・・ do i=0, Nmax call random_seed !乱数の初期値 call random_number(random_value) ! 座標変換 if (0<= random_value .and. random_value<0.01) then x = 0.d0 y = 0.16d0*y ix = int(Nx/2 + x*Nx/10) iy = int(y * Ny /12) f(ix,iy) = 1.0d0  VTKFortranライブラリの使用 module subtoutine, function https://github.com/szaghi/VTKFortran いまさら Fortran？