!----------------------------------------------------------------------
!     Copyright (c) 2008 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  SPMODEL 2 次元水路領域ブシネスク対流モデル
!      チェビシェフ法による計算(温度固定条件)
!
!      時間積分 : Crank Nicolson 法
!
!      計算結果比較のための解析解については解説文書
!
!         「2 次元水路領域での拡散型方程式の解析解」
!
!      (et_diffusion.pdf)を参照のこと. 
!
! Solving 2-D Diffusioon equation (DD problem)
!     dφ/dt = \nabla φ
!     φ = 0 at y=-d/2,d/2
!
!履歴  2008/10/12  竹広真一: et_galerkin_DD_cn_diff.f90 より改造
!
program et_DD_cn_diff
  !
  ! SPMODEL 2 次元水路領域ブシネスク対流モデル
  !      チェビシェフガラーキン法による計算(温度固定条件)
  !
  !      時間積分 : Crank Nicolson 法
  !
  ! Solving 2-D Diffusioon equation (DD problem)
  !     dφ/dt = \nabla φ
  !     φ = 0 at y=-d/2,d/2
  !
  !      計算結果比較のための解析解については解説文書
  !
  !         「2 次元水路領域での拡散型方程式の解析解」
  !
  !      (et_diffusion.pdf)を参照のこと. 
  !
  use lumatrix
  use et_module
  use gt4_history
  implicit none

 !---- 空間解像度設定 ----
  integer, parameter  :: im=64, jm=16      ! 格子点の設定(X,Y)
  integer, parameter  :: km=21, lm=16      ! 切断波数の設定(X,Y)

 !---- 変数 ----
  real(8)             :: yx_Phi(0:jm,0:im-1)     ! 格子データ(温度)
  real(8)             :: et_Phi(-km:km,0:lm)     ! チェビシェフデータ(温度)
  
  real(8)             :: yx_PhiInit(0:jm,0:im-1) ! 格子データ(温度, 初期値)
  real(8)             :: yx_PhiSol(0:jm,0:im-1)  ! 格子データ(温度, 解析解)
  real(8)             :: yx_PhiError(0:jm,0:im-1)! 格子データ(温度, 誤差)

 !---- 時間積分パラメター ----
  real(8), parameter :: dt=1e-3                 ! 時間ステップ間隔
  integer, parameter :: nt=100, ndisp=10        ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
  real(8), parameter          :: XD=2.0         ! x 方向領域の大きさ
  real(8), parameter          :: YD=1.0         ! y 方向領域の大きさ
  real(8), parameter          :: Kappa=1.0      ! 拡散係数
  character(len=2)            :: BC='DD'        ! 境界条件(DD/DN/ND/NN)

 !---- 座標変数など ----
  real(8), parameter  :: xmin=0, xmax=XD        ! X 座標範囲
  real(8), parameter  :: ymin=-YD/2, ymax=YD/2  ! Y 座標範囲


 !--- Crank Nicholson 陰的計算用作業変数 ---
  real(8) :: DiffLUMT(-km:km,0:lm,0:lm)     
                                    ! Crank Nicholson 陰的計算用行列
  integer :: kpivot(-km:km,0:lm)    ! Crank Nicholson 用行列ピボット情報

 !---- その他 ----
  real(8), parameter  :: pi=3.1415926535897932385D0
  integer :: it=0
  real(8) :: time=0.0D0
  integer :: k, l

 !---------------- 座標値の設定 ---------------------
  call et_Initial(im,jm,km,lm,xmin,xmax,ymin,ymax)

  call CNDiffusionMatrix( kappa, dt, DiffLUMT, kpivot )

 !------------------- 初期値設定 ----------------------
  write(6,*) 'Positive k:cos(kx), negative k: sin(kx)'
  write(6,*) 'Positvie l:cos(lx), negative l: sin(lx)'
  write(6,*) 'k,l?'
  read(5,*)  k,l

  if ( k .ge. 0 .AND. l .ge. 0) then
     write(6,*) 'Type of the initial value : cos(kx)*cos(ly)'
     yx_Phi = cos(2*k*Pi*yx_X/XD) * cos((2*l+1)*Pi*yx_Y/YD)
  else if ( k .lt. 0 .AND. l .ge. 0) then
     write(6,*) 'Type of the initial value : sin(kx)*cos(ly)'
     yx_Phi = sin(2*k*Pi*yx_X/XD) * cos((2*l+1)*Pi*yx_Y/YD)
  else if ( k .ge. 0 .AND. l .lt. 0) then
     write(6,*) 'Type of the initial value : cos(kx)*cos(ly)'
     yx_Phi = cos(2*k*Pi*yx_X/XD) * sin(2*l*Pi*yx_Y/YD)
  else 
     write(6,*) 'Type of the initial value : sin(kx)*sin(ly)'
     yx_Phi = sin(2*k*Pi*yx_X/XD) * sin(2*l*Pi*yx_Y/YD)
  endif
  yx_PhiInit = yx_Phi
  
  et_Phi = et_yx(yx_Phi)

  call output_gtool4_init                            ! ヒストリー初期化
  call output_gtool4

 !------------------- 時間積分 ----------------------
  do it=1,nt               ! Adams-Bashfort & Crank-Nicolson 法による時間積分
     time =it*dt

     et_Phi = et_Phi + dt/2 * kappa * et_Lapla_et(et_Phi)
     et_Phi = LuSolve(DiffLUMT,kpivot,et_Phi)

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        call output_gtool4
     endif
  enddo

  call output_gtool4_close                           ! ヒストリー後処理
  stop

  contains

   !------------------- 拡散項 ----------------------
    subroutine CNDiffusionMatrix( Diffc, dt, DiffLUMatrix, kpivot )
      ! 
      ! Crank Nicholson 拡散陰的計算用行列の設定
      ! (1-D dt/2▽^2, LU 分解)
      ! 
      real(8), intent(IN)   :: Diffc    ! 拡散係数
      real(8), intent(IN)   :: dt       ! 時間刻

      ! Crank Nicholson 拡散陰的計算用行列(1-D dt/2▽^2, LU 分解)
      real(8), intent(OUT)  :: DiffLUMatrix(-km:km,0:lm,0:lm)
      integer, intent(OUT)  :: kpivot(-km:km,0:lm)            ! ピボット情報

      real(8)               :: et_I(-km:km,0:lm)              ! 作業用変数
      real(8)               :: et_DI(-km:km,0:lm)             ! 作業用変数

      integer :: l

      DiffLUMatrix = 0.0
      do l=0,lm
         et_I = 0.0 ; et_I(:,l) = 1.0             ! 各波数成分独立
         et_DI =  - Diffc * dt/2.0 * et_Lapla_et(et_I)
         call et_Boundaries(et_DI,cond=BC)
         DiffLUMatrix(:,:,l) = et_I + et_DI
      enddo

      call LuDeComp(DiffLUMatrix,kpivot)

    end subroutine CNDiffusionMatrix

   !------------------- 出力 ----------------------
    subroutine output_gtool4_init
      !
      ! ヒストリー出力初期化ルーチン
      !
      !   ファイル作成
      !   変数定義
      !
      call HistoryCreate( &                               ! ヒストリー作成
           file='et_DD_cn_diff.nc', &
           title='Diffusion equation', &
           source='Sample program of gtool_history/gtool4', &
           institution='GFD_Dennou Club davis/spmodel project',&
           dims=(/'x','y','t'/), dimsizes=(/im,jm+1,0/),    &
           longnames=(/'X-coordinate','Y-coordinate','time        '/),&
           units=(/'1','1','1'/),                           &
           origin=0.0, interval=real(ndisp*dt) )

      call HistoryPut('x',x_X)                            ! 変数出力
      call HistoryAddattr('x','topology','circular')      ! 周期属性
      call HistoryAddattr('x','modulo',xmax-xmin)         ! 周期属性
      call HistoryPut('y',y_Y)                            ! 変数出力

      call HistoryAddVariable( &                          ! 変数定義
           varname='phi', dims=(/'x','y','t'/), & 
           longname='diffusive element', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='phisol', dims=(/'x','y','t'/), & 
           longname='diffusive element(analytic)', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='phierror', dims=(/'x','y','t'/), & 
           longname='diffusive element(Error)', units='1', xtype='double')

     !---- 実験パラメターを属性として定義, 出力(Global 属性) ----
      call HistoryAddAttr('x','+kappa', kappa )

    end subroutine output_gtool4_init

    subroutine output_gtool4
      !
      ! ヒストリー出力
      !
      !   各物理量の gtool4 ファイルへの出力を行う
      !
      yx_Phi = yx_et(et_Phi)

      if ( l .ge. 0 ) then
         yx_PhiSol = yx_PhiInit &
              * exp( - kappa * (   (2*PI*k/XD)**2 &
                                 +((2*l+1)*PI/YD)**2   )*time )
      else
         yx_PhiSol = yx_PhiInit &
              * exp( - kappa * (  (2*PI*k/XD)**2 &
                                 +(2*PI*l/YD)**2     )*time )
      endif

      yx_PhiError = yx_Phi-yx_PhiSol

      write(6,*) 'it = ',it
      call HistoryPut('t',real(time))
      call HistoryPut('phi',transpose(yx_Phi))
      call HistoryPut('phisol',transpose(yx_PhiSol))
      call HistoryPut('phierror',transpose(yx_PhiError))
    end subroutine output_gtool4

    subroutine output_gtool4_close
      !
      ! ヒストリー出力終了処理
      !
      call HistoryClose
    end subroutine output_gtool4_close

end program et_DD_cn_diff
