!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  T2PACK/ISPACK-F90
!      2 次元水路領域問題, Fourier 展開 + Chebyshev 展開法
!
!履歴  2002/01/27  竹広真一
!
module t2pack

  use lumatrix
  use ftrpack
  use ctpack, g_y => g_x, t_dy_t => t_dx_t

  implicit none
  private
  public t2initial
  public gg_x, gg_y
  public gg_et, et_gg
  public et_dx_et, et_dy_et, et_lapla_et, et_jacobian_et_et
  public etboundaries
  public et_laplainv_et, eg_vor2strm_eg
  !public et_vor2strm_et, et_vor2strm1_et

  public g_x, g_y
  public e_g, g_e, e_dx_e
  public t_g, g_t, t_dy_t
  public ctboundaries_DD, ctboundaries_DN, ctboundaries_ND, ctboundaries_NN

  integer            :: im=32, jm=8      ! 格子点の設定(X,Y)
  integer            :: km=10, lm=5      ! 切断波数の設定(X,Y)
  real(8)            :: xl=2.0, yl=1.0   ! 領域の大きさ
  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), dimension(:,:), allocatable :: gg_x, gg_y

  save im, jm, km, lm, xl, yl

  contains
  !--------------- 初期化 -----------------
    subroutine t2initial(i,j,k,l,xmin,xmax,ymin,ymax)

      integer,intent(in) :: i, j           ! 格子点の設定(X,Y)
      integer,intent(in) :: k, l           ! 切断波数の設定(X,Y)

      real(8),intent(in) :: xmin, xmax     ! X 座標範囲
      real(8),intent(in) :: ymin, ymax     ! Y 座標範囲

      im = i       ; jm = j
      km = k       ; lm = l
      xl = xmax-xmin ; yl = ymax-ymin

      call ftrinitial(im,km,xmax-xmin)
      call ctinitial(jm,lm,ymin,ymax)

      allocate(gg_x(0:jm,0:im-1),gg_y(0:jm,0:im-1))
      gg_x = spread(g_x,1,jm+1)
      gg_y = spread(g_y,2,im)

    end subroutine t2initial

  !--------------- 基本変換 -----------------

    function gg_et(et)  ! スペクトル -> 台形格子
      real(8), dimension(0:jm,0:im-1)              :: gg_et
      real(8), dimension(-km:km,0:lm), intent(in)  :: et

      gg_et = g_e(transpose(g_t(et)))

    end function gg_et

    function et_gg(gg)  ! 台形格子 -> スペクトル
      real(8), dimension(-km:km,0:lm)              :: et_gg
      real(8), dimension(0:jm,0:im-1), intent(in)  :: gg

      et_gg = t_g(transpose(e_g(gg)))

    end function et_gg

  !--------------- 微分計算 -----------------

    function et_dx_et(et)   ! スペクトルに作用する x 微分演算子
      real(8), dimension(-km:km,0:lm)                :: et_dx_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: et
      integer k

      do k=-km,km
         et_dx_et(k,:)  =  (-2*pi*k/xl)*et(-k,:)
      enddo
    end function et_dx_et

    function et_dy_et(et)   ! スペクトルに作用する y 微分演算子
      real(8), dimension(-km:km,0:lm)               :: et_dy_et
      real(8), dimension(-km:km,0:lm), intent(in)   :: et

      et_dy_et = t_dy_t(et)

    end function et_dy_et

    function et_lapla_et(et)   ! スペクトルに作用する Laplacian 演算子
      real(8), dimension(-km:km,0:lm)                :: et_lapla_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: et
      integer k

      do k=-km,km
         et_lapla_et(k,:) = -(2*pi*k/xl)**2*et(k,:)
      enddo

      et_lapla_et = et_lapla_et + et_dy_et(et_dy_et(et))

    end function et_lapla_et

    function et_jacobian_et_et(et_a,et_b) ! スペクトルに作用する Jacobian
      real(8), dimension(-km:km,0:lm)                :: et_jacobian_et_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: et_a
      real(8), dimension(-km:km,0:lm), intent(in)    :: et_b
      integer k

      et_jacobian_et_et = et_gg(&
           gg_et(et_dx_et(et_a)) * gg_et(et_dy_et(et_b)) &
           -gg_et(et_dy_et(et_a)) * gg_et(et_dx_et(et_b)) )

    end function et_jacobian_et_et


  !--------------- 境界値問題 -----------------

    subroutine etboundaries(et,values,cond)   ! ディリクレ, ノイマン条件
      ! Chebyshev 空間での境界条件適用

      real(8), dimension(-km:km,0:lm),intent(inout)      :: et
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension(-km:km,2), intent(in), optional :: values
              ! 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              ! 省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              ! 境界条件. 省略時は 'DD'
              !   DD : 両端ディリクレ
              !   DN,ND : ディリクレ/ノイマン条件
              !   NN : 両端ノイマン

      if (.not. present(cond)) then
         if (present(values)) then
            call ctboundaries_DD(et,values)
         else
            call ctboundaries_DD(et)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call ctboundaries_NN(et,values)
         else
            call ctboundaries_NN(et)
         endif
      case ('DN')
         if (present(values)) then
            call ctboundaries_DN(et,values)
         else
            call ctboundaries_DN(et)
         endif
      case ('ND')
         if (present(values)) then
            call ctboundaries_ND(et,values)
         else
            call ctboundaries_ND(et)
         endif
      case ('DD')
         if (present(values)) then
            call ctboundaries_DD(et,values)
         else
            call ctboundaries_DD(et)
         endif
      case default
         call msgdmp('E','etboundaries','B.C. not supported')
      end select

    end subroutine etboundaries

    function et_laplainv_et(et,values) ! 逆 Laplacian, Dirichlet 型境界条件
      ! Chebyshev-tau 法による計算

      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      real(8), dimension(-km:km,0:lm)             :: et_laplainv_et
      real(8), dimension(-km:km,2), intent(in), optional :: values   ! 境界値

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(0:lm,0:lm)           :: cheb
      real(8), dimension(0:lm,0:jm)           :: g
      real(8), dimension(-km:km)              :: value1, value2   ! 境界値

      logical :: first = .true.
      integer :: k,l
      save    :: alu, kp, first

      if (.not. present(values)) then
         value1=0 ; value2=0
      else
         value1 = values(:,1) ; value2 = values(:,2)
      endif

      if ( first ) then
         first = .false.

         allocate(alu(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

         cheb=0
         do l=0,lm
            cheb(l,l)=1
         enddo
         g=g_t(cheb)

         do k=-km,km
            alu(k,:,:) = transpose(t_dy_t(t_dy_t(cheb)) - (2*pi*k/xl)**2*cheb)
            alu(k,lm-1,:) = g(:,0)
            alu(k,lm,:)   = g(:,jm)
         enddo

         call ludecomp(alu,kp)
      endif

      et_work = et
      et_work(:,lm-1) = value1
      et_work(:,lm)   = value2
      et_laplainv_et = lusolve(alu,kp,et_work)

    end function et_laplainv_et

    function eg_vor2strm_eg(eg,values,cond)    ! 渦度から流線を求める. 
                                               ! デフォルトは粘着条件
      ! Chebyshev-Collocation 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta, 
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries

      ! 注 : 最初に呼ばれるときの境界条件で以後計算される(要仕様変更)

      real(8), dimension(-km:km,0:jm),intent(in)  :: eg
              ! 入力渦度分布

      real(8), dimension(-km:km,0:jm)             :: eg_vor2strm_eg
              ! 出力流線分布

      real(8), dimension(2), intent(in), optional :: values
              ! 流線境界値. 境界で一定なので波数 0 成分のみ

      character(len=2), intent(in), optional  :: cond
              ! 境界条件スイッチ. 省略時は 'RR'
              !   RR : 両端粘着条件
              !   RF,FR : 粘着/応力なし条件
              !   FF : 両端応力なし条件

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(-km:km,0:jm)         :: eg_work
      real(8), dimension(0:jm,0:jm)           :: grid
      real(8), dimension(0:jm,0:jm)           :: g
      real(8)                                 :: value1, value2   ! 境界値
      logical                                 :: rigid1, rigid2   ! 境界条件

      logical :: first = .true.
      integer :: k,j
      save    :: alu, kp, first

      if (.not. present(values)) then
         value1=0 ; value2=0
      else
         value1 = values(1) ; value2 = values(2)
      endif

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call msgdmp('E','eg_vor2strm_eg','B.C. not supported')
         end select
      endif

      if ( first ) then
         first = .false.

         allocate(alu(-km:km,0:jm,0:jm),kp(-km:km,0:jm))

         grid=0
         do j=0,jm
            grid(j,j)=1
         enddo
         do k=-km,km
            alu(k,:,:) = transpose( &
                 g_t(t_dy_t(t_dy_t(t_g(grid)))) &
                 - (2*pi*k/xl)** 2* grid )
         enddo

         ! 運動学的条件. 流線は境界で一定
         g=grid
         do k=-km,km
            alu(k,0,:)   = g(:,0)
            alu(k,jm,:)  = g(:,jm)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            g=g_t(t_dy_t(t_g(grid)))
         else
            g=g_t(t_dy_t(t_dy_t(t_g(grid))))
         endif
         do k=-km,km
            alu(k,1,:) = g(:,0)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid2 ) then
            g=g_t(t_dy_t(t_g(grid)))
         else
            g=g_t(t_dy_t(t_dy_t(t_g(grid))))
         endif
         do k=-km,km
            alu(k,jm-1,:) = g(:,jm)
         enddo

         call ludecomp(alu,kp)
      endif

      eg_work = eg
      eg_work(:,1)    = 0               ! 力学的条件
      eg_work(:,jm-1) = 0               ! 力学的条件

      eg_work(:,0) = 0            ! 運動学的条件. 波数 0 以外は 0 
      eg_work(0,0) = value1*2     ! 運動学的条件. 波数 0 は重み 1/2

      eg_work(:,jm)   = 0            ! 運動学的条件. 波数 0 以外は 0 
      eg_work(0,jm)   = value2*2     ! 運動学的条件. 波数 0 は重み 1/2

      eg_vor2strm_eg = lusolve(alu,kp,eg_work)

    end function eg_vor2strm_eg

   !-----  注意 : 計算うまくいかず. pending ---- 
    function et_vor2strm_et(et,values,rigid) ! 渦度から流線を求める. 
                                             ! デフォルトは粘着条件
      ! Chebyshev-tau 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta, 
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries
      !
      ! l=0,1,lm-1,lm 成分の式の代わりに境界条件を与える. 
      ! 渦度の低次成分を無視することは
      ! \nabla^4 \psi = \zeta^2 を解いていることに相当. 
      ! 4 階の微分方程式にすることで境界条件の数とつじつまがあう. 

      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      real(8), dimension(-km:km,0:lm)             :: et_vor2strm_et

      ! 流線境界値. 境界で一定なので波数 0 成分のみ
      real(8), dimension(2), intent(in), optional :: values
      ! 境界条件スイッチ
      logical, dimension(2), intent(in), optional :: rigid

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(0:lm,0:lm)           :: cheb
      real(8), dimension(0:lm,0:jm)           :: g
      real(8)                                 :: value1, value2   ! 境界値
      logical                                 :: rigid1, rigid2   ! 境界条件

      logical :: first = .true.
      integer :: k,l
      save    :: alu, kp, first

      if (.not. present(values)) then
         value1=0 ; value2=0
      else
         value1 = values(1) ; value2 = values(2)
      endif

      if (.not. present(rigid)) then
         rigid1=.true. ; rigid2=.true.
      else
         rigid1 = rigid(1) ; rigid2 = rigid(2)
      endif

      if ( first ) then
         first = .false.

         allocate(alu(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

         cheb=0
         do l=0,lm
            cheb(l,l)=1
         enddo
         do k=-km,km
            alu(k,:,:) = transpose(t_dy_t(t_dy_t(cheb)) - (2*pi*k/xl)**2*cheb)
         enddo

         ! 運動学的条件. 流線は境界で一定
         g=g_t(cheb)
         do k=-km,km
            alu(k,lm-1,:) = g(:,0)
            alu(k,lm,:)   = g(:,jm)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            g=g_t(t_dy_t(cheb))
         else
            g=g_t(t_dy_t(t_dy_t(cheb)))
         endif
         do k=-km,km
            alu(k,0,:) = g(:,0)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid2 ) then
            g=g_t(t_dy_t(cheb))
         else
            g=g_t(t_dy_t(t_dy_t(cheb)))
         endif
         do k=-km,km
            alu(k,1,:) = g(:,jm)
         enddo

         call ludecomp(alu,kp)
      endif

      et_work = et
      et_work(:,0) = 0               ! 力学的条件
      et_work(:,1) = 0               ! 力学的条件

      et_work(:,lm-1) = 0            ! 運動学的条件. 波数 0 以外は 0 
      et_work(0,lm-1) = value1*2     ! 運動学的条件. 波数 0 は重み 1/2

      et_work(:,lm)   = 0            ! 運動学的条件. 波数 0 以外は 0 
      et_work(0,lm)   = value2*2     ! 運動学的条件. 波数 0 は重み 1/2

      et_vor2strm_et = lusolve(alu,kp,et_work)

    end function et_vor2strm_et

   !-----  注意 : 計算うまくいかず. pending ---- 
    function et_vor2strm1_et(et,values,rigid) ! 渦度から流線を求める. 
                                              ! デフォルトは粘着条件
      ! Chebyshev-tau 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta, 
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries
      !
      ! \nabla^4 \psi = \nabla^2\zeta を解く
      ! 4 階の微分方程式にすることで境界条件の数とつじつまがあう. 

      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      real(8), dimension(-km:km,0:lm)             :: et_vor2strm1_et

      ! 流線境界値. 境界で一定なので波数 0 成分のみ
      real(8), dimension(2), intent(in), optional :: values
      ! 境界条件スイッチ
      logical, dimension(2), intent(in), optional :: rigid

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(0:lm,0:lm)           :: cheb
      real(8), dimension(0:lm,0:jm)           :: g
      real(8)                                 :: value1, value2   ! 境界値
      logical                                 :: rigid1, rigid2   ! 境界条件

      logical :: first = .true.
      integer :: k,l
      save    :: alu, kp, first

      if (.not. present(values)) then
         value1=0 ; value2=0
      else
         value1 = values(1) ; value2 = values(2)
      endif

      if (.not. present(rigid)) then
         rigid1=.true. ; rigid2=.true.
      else
         rigid1 = rigid(1) ; rigid2 = rigid(2)
      endif

      if ( first ) then
         first = .false.

         allocate(alu(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

         cheb=0
         do l=0,lm
            cheb(l,l)=1
         enddo
         do k=-km,km
            alu(k,:,:) = transpose( &
                 t_dy_t(t_dy_t(t_dy_t(t_dy_t(cheb)))) &
                 - 2 * (2*pi*k/xl)**2 * t_dy_t(t_dy_t(cheb)) &
                 + (2*pi*k/xl)**4*cheb &
                 )
         enddo

         ! 運動学的条件. 流線は境界で一定
         g=g_t(cheb)
         do k=-km,km
            alu(k,lm-1,:) = g(:,0)
            alu(k,lm,:)   = g(:,jm)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            g=g_t(t_dy_t(cheb))
         else
            g=g_t(t_dy_t(t_dy_t(cheb)))
         endif
         do k=-km,km
            alu(k,lm-3,:) = g(:,0)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid2 ) then
            g=g_t(t_dy_t(cheb))
         else
            g=g_t(t_dy_t(t_dy_t(cheb)))
         endif
         do k=-km,km
            alu(k,lm-2,:) = g(:,jm)
         enddo

         call ludecomp(alu,kp)
      endif

      et_work = et_lapla_et(et)
      et_work(:,lm-3) = 0               ! 力学的条件
      et_work(:,lm-2) = 0               ! 力学的条件

      et_work(:,lm-1) = 0            ! 運動学的条件. 波数 0 以外は 0 
      et_work(0,lm-1) = value1*2     ! 運動学的条件. 波数 0 は重み 1/2

      et_work(:,lm)   = 0            ! 運動学的条件. 波数 0 以外は 0 
      et_work(0,lm)   = value2*2     ! 運動学的条件. 波数 0 は重み 1/2

      et_vor2strm1_et = lusolve(alu,kp,et_work)

    end function et_vor2strm1_et

  end module t2pack

