【Fortranコード】死ぬほどべた書きで書いた「チャンネル流れ」のプログラム
最近出たこちらの参考書をきっかけにFortranを勉強しています。
思いっきりべた書きでチャンネル流れのコードをFortranで書いてみました。
あまりきれいで汎用性があるコードではありませんが追々書き換えていくつもりです。
gfortran 9.4.0
paraview 5.9.0(可視化用）
Fortranコード
program main
implicit none
!----- パラメータ設定 -----
real(kind=8), Parameter :: Lx = 10.0d0, Ly = 1.0d0 !空間領域サイズ
integer, Parameter :: Nx = 150, Ny = 50 !空間分割数
integer, Parameter :: Nt = 1000 !ステップ数
real(kind=8), Parameter :: dx = Lx/real(Nx), dy = Ly/real(Ny) !空間刻み
real(kind=8), Parameter :: dt = 0.01d0 !時間刻み
real(kind=8), Parameter :: Re = 200.0d0 !レイノルズ数
integer :: i, j, t, l
!----- 空間座標 -----
real(kind=8), allocatable, dimension(:) :: x, y
!----- 変数 -----
real(kind=8), allocatable, dimension(:,:) :: u, v, p
!---- 微分 -----
real(kind=8), allocatable, dimension(:,:) :: du, dv, ux, uy, vx, vy, uxx, uyy, vxx, vyy
! ポアソン方程式の収束定義
real(kind=8) :: p_prev,MaxErr,CurErr,f(0:Nx,0:Ny)
real(kind=8) :: MaxP = 1.0d-10
real(kind=8) :: Conv = 1.0d-3
integer,Parameter :: lmax = 5000
! ファイル保存
character(len=40) :: filename
!======== main program =========
! メモリの割付
allocate(x(0:Nx), y(0:Ny))
allocate(u(0:Nx, 0:Ny), v(0:Nx, 0:Ny), p(0:Nx, 0:Ny))
allocate(du(1:Nx, 0:Ny), dv(1:Nx, 0:Ny))
allocate(ux(1:Nx, 0:Ny), uy(1:Nx, 0:Ny), vx(1:Nx, 0:Ny), vy(1:Nx, 0:Ny))
allocate(uxx(1:Nx, 0:Ny), uyy(1:Nx, 0:Ny), vxx(1:Nx, 0:Ny), vyy(1:Nx, 0:Ny))
!x,y座標
do i = 1, Nx
x(i) = dx*real((i-Nx/2))
y(i) = dy*real((i-Ny/2))
end do
!======== 初期状態の設定 ==========
do j = 0, Ny
do i = 0, Nx
u(i,j) = 0.0d0
v(i,j) = 0.0d0
p(i,j) = 0.0d0
end do
end do
! ===== NS方程式 =====
call system("mkdir data")
do t =1, Nt
print*, t
do j=1,Ny-1
do i=1,Nx-1
!======== 境界条件 ==========
!----- inlet flow -----
if (i == 1) then
u(0,j) = 1.0d0
v(0,j) = 0.0d0
p(0,j) = p(1,j) !gradient 0
end if
!----- outlet flow -----
if (i == Nx-1) then
u(Nx,j) = u(Nx-1,j) !gradient 0
v(Nx,j) = v(Nx-1, j) !gradient 0
p(Nx,j) = p(Nx-1,j) ! 大気圧基準0kPa
end if
! ----- wall -----
if (j == 1) then
u(i,0) = 0.0d0 !noslip
v(i,0) = 0.0d0 !noslip
p(i,0) = P(i,1) !gradient 0
end if
if (j == Ny-1) then
u(i,Ny) = 0.0d0 !noslip
v(i,Ny) = 0.0d0 !noslip
p(i,Ny) = p(i,Ny-1) !gradient 0
end if
! NS方程式
ux (i,j) = ( u(i+1,j )-u(i-1,j ) ) / (2.0d0*dx)
uxx(i,j) = ( u(i+1,j )-2.0d0*u(i,j)+u(i-1,j ) )/(dx*dx)
uy (i,j) = ( u(i ,j+1)-u(i ,j-1) ) / (2.0d0*dy)
vx (i,j) = ( v(i+1,j )-v(i-1,j ) ) / (2.0d0*dx)
vy (i,j) = ( v(i ,j+1)-v(i ,j-1) ) / (2.0d0*dy)
uyy(i,j) = ( u(i ,j+1)-2.0d0*u(i,j)+u(i ,j-1) )/(dy*dy)
vxx(i,j) = ( v(i+1,j )-2.0d0*v(i,j)+v(i-1,j ) )/(dx*dx)
vyy(i,j) = ( v(i ,j+1)-2.0d0*v(i,j)+v(i ,j-1) )/(dy*dy)
du(i,j) = (1.d0/Re)*( uxx(i,j)+uyy(i,j) ) - ( u(i,j)*ux(i,j)+v(i,j)*uy(i,j) )
dv(i,j) = (1.d0/Re)*( vxx(i,j)+vyy(i,j) ) - ( u(i,j)*vx(i,j)+v(i,j)*vy(i,j) )
! 速度の更新
u(i,j) = u(i,j) + du(i,j)*dt
v(i,j) = v(i,j) + dv(i,j)*dt
end do
end do
! ポアソン方程式
do l=1,lmax ! 収束イタレーション
MaxErr = 0.0d0
CurErr = 0.0d0
do j=1,Ny-1
do i=1,Nx-1
f(i,j) = (1.d0/dt)*( ( u(i+1,j )-u(i-1,j ) ) / (2.0d0*dx) &
+ ( v(i ,j+1)-v(i ,j-1) ) / (2.0d0*dy) )
p_prev = p(i,j)
p(i,j) = ( dx*dx*(p(i,j+1)+p(i,j-1)) &
+dy*dy*(P(i+1,j)+p(i-1,j)) &
-dx*dx*dy*dy*f(i,j) ) / (2.0d0*(dx*dx+dy*dy))
if (MaxP.LT.abs(p(i,j))) MaxP = p(i,j)
CurErr = ( abs(p(i,j)-p_prev) ) / MaxP
if (MaxErr.LT.CurErr) MaxErr = CurErr
end do
end do
If (MaxErr.LT.Conv) Exit
end do
do j=1,Ny-1
do i=1,Nx-1
u(i,j) = u(i,j) - dt*( p(i+1,j )-p(i-1,j ) ) / (2.0d0*dx)
v(i,j) = v(i,j) - dt*( p(i ,j+1)-p(i ,j-1) ) / (2.0d0*dy)
end do
end do
!--- output ---
If ( mod(t,200)==0 ) then
Write(filename,"(a,i5.5,a)") "data/output",int(t),".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+1, Ny+1, 1
Write(10,"('POINTS ',i9,' float')") (Nx+1)*(Ny+1)*1
do j=0,Ny ; do i=0,Nx
Write(10,"(3(f9.4,1x))") i*dx, j*dy, 0.0d0
end do ; end do
Write(10,"('POINT_DATA ',i9)") (Nx+1)*(Ny+1)*1
!date input
Write(10,"('SCALARS U float')")
Write(10,"('LOOKUP_TABLE default')")
do j=0,Ny ; do i=0,Nx
Write(10,*) u(i,j)
end do; end do
Close(10)
end if
end do ! time step do文
! メモリ解放
! deallocate(x, y)
! deallocate(u, v, p)
! deallocate(du, dv)
! deallocate(ux, uy, vx, vy)
! deallocate(uxx, uyy, vxx, vyy)
end program mainparaviewで出力した結果はこちらです。
慣れないうちはべた書きでコードを書いた方が流れがわかりやすいと思い、ひとまずべた書きでコードを書いてみました。
追ってsubroutineやmodulleといった機能ごとにプログラムを分けて管理しやすいようにして置こうと思います。
そうすることで初期状態、境界条件の設定のみを変えることで他のプログラム部分は意識せずにプログラムを組むことができるということになります。
さらにはOpenMPIで並列化するようにもしたいところです。
Twitter➡@t_kun_kamakiri
Instagram➡kamakiri1225
youtube➡https://www.youtube.com/channel/UCbG6_Q9ZRqqVT6YZOpcjDlQ
ブログ➡宇宙に入ったカマキリ（物理ブログ）
ココナラ➡物理の質問サポートサービス
コミュニティ➡製造業ブロガー