非保存型IDO法の説明
IDO-NCF
the no conservative form of the Interpolated Differential Operator (IDO- CF)
the no conservative form of the Interpolated Differential Operator (IDO- CF)
今日(2007.3.1),東京工業大学を訪問し、スパコンツバメなどを見学。
ツバメは一般学生もIDを作れて特別な許可なしで32CPUまで使えるとのこと…すげー
また、SUNの講習会があったとかでいろいろ話を聞けました
ツバメは一般学生もIDを作れて特別な許可なしで32CPUまで使えるとのこと…すげー
また、SUNの講習会があったとかでいろいろ話を聞けました
- ツバメは各CPU毎にメモリが割り与えられているためそれがボトルネックになったりする
- キャッシュをうまく使うことによって8time(8倍)の性能差が現れたりする
- f(i,j)等をループする場合下記の違いでさえかなりの性能差が出る
Do 100 i=min, max Do 100 j=min, max 100 continue
or
Do 100 j=min, max Do 100 i=min, max 100 continue
こんなの全部コンパイラがチューニングしてくれるのが一番いいんだろうけど、やっぱ奥が深いw
Solaris使ってみてはどうか?とのこと、Solarisには前々から興味があったのだが、コンパイラもSUN純正のがついてくるらしいし、iccとかも学生はただらしい(ぜんぜん知りませんでした)。それにクラスターのジョブ管理は管理用のアプリ(ウェブブラウザ上で動く)のがあるらしくて、各ノードがどれぐらい稼動してるか簡単に見れたのがすごかった。
Solaris使ってみてはどうか?とのこと、Solarisには前々から興味があったのだが、コンパイラもSUN純正のがついてくるらしいし、iccとかも学生はただらしい(ぜんぜん知りませんでした)。それにクラスターのジョブ管理は管理用のアプリ(ウェブブラウザ上で動く)のがあるらしくて、各ノードがどれぐらい稼動してるか簡単に見れたのがすごかった。
- スイッチについてもGigaEntherじゃなくてもっといいのがあると言ってたな(今度調べとけ)
ここからが本題Vlasov方程式を解くならIDO使ってみてはとのこと、(位相空間での)1次元(xとvでの二次元)ではIDOで解けない条件はないだろうとまで言われてた、今度IDOについても調べとけ
テストコード
ido_rk4.f90
Program main implicit none double precision :: dt, dx, o_dx, o_dx2, v_v, f_d, dx_isign, xx double precision :: k1, k2, k3, k4, f1, f2, f3, f4 double precision :: kx1, kx2, kx3, kx4, fx1, fx2, fx3, fx4 double precision :: dft, dftx double precision, dimension(:), allocatable :: f, fn, g, gn, ft, ftx, initial double precision, dimension(:), allocatable :: f1_tmp, f1x_tmp, f2_tmp, f2x_tmp double precision, dimension(:), allocatable :: f1t_tmp, f1tx_tmp, f2t_tmp, f2tx_tmp double precision, dimension(:), allocatable :: f3t_tmp, f3tx_tmp, f4t_tmp, f4tx_tmp integer :: i, j, t, x_up, n_x = 20
allocate( initial(0:n_x+1), f(0:n_x+1), fn(0:n_x+1), g(0:n_x+1), gn(0:n_x+1), ft(0:n_x+1), ftx(0:n_x+1) ) allocate( f1_tmp(0:n_x+1), f1x_tmp(0:n_x+1), f2_tmp(0:n_x+1), f2x_tmp(0:n_x+1) ) allocate( f1t_tmp(0:n_x+1), f1tx_tmp(0:n_x+1), f2t_tmp(0:n_x+1), f2tx_tmp(0:n_x+1) ) allocate( f3t_tmp(0:n_x+1), f3tx_tmp(0:n_x+1), f4t_tmp(0:n_x+1), f4tx_tmp(0:n_x+1) )
!=======================================================================
! Initial Setting
!=======================================================================
! Initial Setting
!=======================================================================
dx = 0.01d0 o_dx = 1.d0/dx dt = 0.01d0 dx_isign = 1.d0 v_v = 0.1d0 xx = -v_v*dt dx_isign = -dsign(1.d0, v_v)
do i = 0, n_x+1 initial(i) = 0.d0 f(i) = 0.d0 fn(i) = 0.d0 ft(i) = 0.d0 ftx(i) = 0.d0 f1_tmp(i) = 0.d0 f1x_tmp(i) = 0.d0 f2_tmp(i) = 0.d0 f2x_tmp(i) = 0.d0 gn(i) = 0.d0 if( 5<=i .and. i<=10 ) f(i) = 1.d0 if( 5<=i .and. i<=10 ) initial(i) = 1.d0 enddo do i = 0, n_x g(i) = (f(i) + f(i+1))/dx enddo g(n_x+1) = 0.d0
!=======================================================================
! Main Loop
!=======================================================================
! Main Loop
!=======================================================================
do t = 1, 10
!########################################################################
do i = 1, n_x x_up = i + int(dx_isign) call cip(f(i), f(x_up), g(i), g(x_up), v_v, ft(i), ftx(i)) f1t_tmp(i) = ft(i) f1tx_tmp(i) = ftx(i)
! f1_tmp(i) = dt*ft(i) !k1 = h*f(x_n,y_n)
! f1x_tmp(i) = dt*ftx(i)
! f1x_tmp(i) = dt*ftx(i)
enddo
!########################################################################
do i = 1, n_x x_up = i + int(dx_isign) call cip(f(i)+f1t_tmp(i)*dt/2.d0, f(x_up)+f1t_tmp(x_up)*dt/2.d0, & g(i)+f1tx_tmp(i)*dt/2.d0, g(x_up)+f1tx_tmp(x_up)*dt/2.d0, & v_v, ft(i), ftx(i)) f2t_tmp(i) = ft(i) !k2 = h*f(x_n+dx,y_n+k1*dt) f2tx_tmp(i) = ftx(i) enddo
!########################################################################
do i = 1, n_x x_up = i + int(dx_isign) call cip(f(i)+f2t_tmp(i)*dt/2.d0, f(x_up)+f2t_tmp(x_up)*dt/2.d0, & g(i)+f2tx_tmp(i)*dt/2.d0, g(x_up)+f2tx_tmp(x_up)*dt/2.d0, & v_v, ft(i), ftx(i)) f3t_tmp(i) = ft(i) !k3 = h*f(x_n+dx,y_n+k1*dt) f3tx_tmp(i) = ftx(i) enddo
!########################################################################
do i = 1, n_x x_up = i + int(dx_isign) call cip(f(i)+f3t_tmp(i)*dt, f(x_up)+f3t_tmp(x_up)*dt, & g(i)+f3tx_tmp(i)*dt, g(x_up)+f3tx_tmp(x_up)*dt, & v_v, ft(i), ftx(i)) f4t_tmp(i) = ft(i) f4tx_tmp(i) = ftx(i) enddo
!########################################################################
do i = 1, n_x f(i) = f(i) + (f1t_tmp(i) + 2.d0*f2t_tmp(i) + 2.d0*f3t_tmp(i) + f4t_tmp(i))*dt/6.d0 g(i) = g(i) + (f1tx_tmp(i) + 2.d0*f2tx_tmp(i) + 2.d0*f3tx_tmp(i) + f4tx_tmp(i))*dt/6.d0 enddo enddo
open(99, file="result_ido_rk4.txt") do i = 1, n_x write(99,*) i, initial(i), f(i) write(*,*) i, initial(i), f(i) enddo close(99) End program main
!-----------------------------------------------------------------------
subroutine cip(f_i, f_up, g_i, g_up, v_v, ft, ftx) implicit none double precision :: f_i, f_up, g_i, g_up, v_v, ft, ftx, dx, o_dx, dx_isign, b, f_d dx = 0.01d0 o_dx = 1.d0/dx dx_isign = -dsign(1.d0, v_v) f_d = (f_up - f_i)*o_dx*dx_isign b = (3.d0*f_d - 2.d0*g_i - g_up)*o_dx*dx_isign ft = -v_v*g_i ftx = -2.d0*v_v*b return end subroutine cip
!-----------------------------------------------------------------------