「IDO Memo」の編集履歴（バックアップ）一覧に戻る

IDO Memo - (2010/03/17 (水) 08:56:08) のソース

*非保存型IDO法の説明

#contents()
&big(){&font(red){IDO-NCF}}
the no conservative form of the Interpolated Differential Operator (IDO- CF)



今日(2007.3.1),東京工業大学を訪問し、スパコンツバメなどを見学。
ツバメは一般学生もIDを作れて特別な許可なしで32CPUまで使えるとのこと…すげー
また、SUNの講習会があったとかでいろいろ話を聞けました
+ツバメは各CPU毎にメモリが割り与えられているためそれがボトルネックになったりする
+キャッシュをうまく使うことによって８time（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
こんなの全部コンパイラがチューニングしてくれるのが一番いいんだろうけど、やっぱ奥が深いｗ
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
!=======================================================================
      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
!=======================================================================
      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)
        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
!-----------------------------------------------------------------------