!---------------------------------------------------------------------- ! Copyright(c) 2002-2014 SPMODEL Development Group All rights reserved. !---------------------------------------------------------------------- ! Solving KdV equation ! ! d\zeta/dt = -\zeta d\zeta/dx + d^3\zeta/dx^3 ! ! Time integration is performed by the _leap frog scheme_. ! program kdv use ae_module use gtool_history use dc_types, only: DP implicit none !---- 空間解像度設定 ---- integer, parameter :: im=128 ! 格子点の設定(X) integer, parameter :: km=42 ! 切断波数の設定(X) !---- 変数 ---- real(DP) :: g_Zeta(0:im-1) ! 格子データ real(DP) :: e_Zeta(-km:km) ! スペクトルデータ real(DP) :: e_Zeta0(-km:km) ! スペクトルデータ real(DP) :: e_Zeta1(-km:km) ! スペクトルデータ !---- 座標変数など ---- real(DP), parameter :: xmin=0.0, xmax=3.0 !---- 時間積分パラメター ---- real(DP), parameter :: dt=1d-6 ! 時間ステップ間隔 integer, parameter :: nt=5000, ndisp=100 ! 時間積分数, 表示ステップ !---- 物理パラメター ---- real(DP), parameter :: X1=(xmax+xmin)*5D-1 ! 初期分布 X 座標 real(DP), parameter :: U1=7.20D2 ! 初期分布の振幅 real(DP), parameter :: X2=(xmax+3.0D0*xmin)*2.5D-1 ! 初期分布 X 座標 real(DP), parameter :: U2=1.44D3 ! 初期分布の振幅 real(DP) :: pi ! 円周率 integer :: it ! DO 変数 pi = 4.0D0*atan(1.0D0) call ae_initial(im,km,xmin,xmax) ! スペクトル初期化 !------------------- 初期値設定 ---------------------- g_Zeta = U1*sech((g_X-X1)/sqrt(12/u1))**2 & & + U2*sech((g_X-X2)/sqrt(12/u2))**2 e_Zeta = e_g(g_Zeta) e_Zeta1 = e_Zeta e_Zeta0 = e_Zeta it = 0 ! 初期化 call output_gtool_init ! ヒストリー初期化 call output_gtool !------------------- 時間積分 ---------------------- do it=1,nt e_Zeta = e_Zeta0 + 2*dt * & ! Leap frog 時間積分 & ( -e_g(g_e(e_Zeta1)*g_e(e_Dx_e(e_Zeta1))) & & - e_Dx_e(e_Dx_e(e_Dx_e(e_Zeta1))) ) e_Zeta0=e_Zeta1 e_Zeta1=e_Zeta if(mod(it,ndisp) .eq. 0)then ! 出力 g_Zeta = g_e(e_Zeta) call output_gtool endif enddo call output_gtool_close ! ヒストリー後処理 contains function sech(x) real(DP) :: x(0:im-1) real(DP) :: sech(0:im-1) sech = 2.0D0/(exp(x)+ exp(-x)) end function sech subroutine output_gtool_init call HistoryCreate( & ! ヒストリー作成 file='kdv1.nc', title='KDV equation model (Leap frog)', & source='Sample program of SPMODEL/KdV equation', & institution='GFD_Dennou Club: SPMODEL Projectt', & dims=(/'x','t'/), dimsizes=(/im,0/), & longnames=(/'X-coordinate','time '/), & units=(/'1','1'/), & origin=0.0, interval=real(ndisp*dt) ) call HistoryPut('x',g_X) ! 変数出力 call HistoryAddattr('x','topology','circular') ! 周期属性 call HistoryAddattr('x','modulo',xmax-xmin) ! 周期属性 call HistoryAddVariable( & ! 変数定義 varname='zeta', dims=(/'x','t'/), & longname='displacement', units='1', xtype='double') end subroutine output_gtool_init subroutine output_gtool write(6,*) 'it = ',it call HistoryPut('t',real(it*dt)) call HistoryPut('zeta',g_Zeta) end subroutine output_gtool subroutine output_gtool_close call HistoryClose end subroutine output_gtool_close end program kdv