!---------------------------------------------------------------------- ! COPYRIGHT (c) 2006 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 at_aq_galerkin_RRFF2 ! ! チェビシェフ−ガラーキン法 ! 非圧縮流体の流線関数・流れポテンシャル用 ! ! case RR: 両端粘着条件(両端で値と 1 階微分が 0) ! f(:,i=0)=f(:,i=im)=f'(:,i=0)=f'(:,i=im)=0 ! [ f(:,x=xmax)=f(:,x=xmin)=f'(:,x=xmax)=0=f'(:,x=xmin)=0 ] ! ! case FF: 両端自由すべり条件(両端で値と 2 階微分が 0) ! f(:,i=0)=f(:,i=im)=f''(:,i=0)=f''(:,i=im)=0 ! [ f(:,x=xmax)=f(:,x=xmin)=f''(:,x=xmax)=0=f''(:,x=xmin)=0 ] ! ! case RF: 片端粘着条件, 片端自由すべり条件 ! (両端で値が 0, 片側で 1 階微分, もう一方で 2 階微分が 0) ! f(:,i=0)=f(:,i=im)=f'(:,i=0)=f''(:,i=im)=0 ! [ f(:,x=xmax)=f(:,x=xmin)=f''(:,x=xmax)=0=f'(:,x=xmin)=0 ] ! ! case FR: 片端自由すべり条件, 片端粘着条件 ! (両端で値が 0, 片側で 2 階微分, もう一方で 1 階微分が 0) ! f(:,i=0)=f(:,i=im)=f''(:,i=0)=f'(:,i=im)=0 ! [ f(:,x=xmax)=f(:,x=xmin)=f'(:,x=xmax)=0=f''(:,x=xmin)=0 ] ! ! k 次のガラーキン基底が ! \phi_k(x)=T_k(x) + C_3 T_3(x) + C_2 T_2(x) ! + C_1 T_1(x) + C_0 T_0(x) ! の型式(Type2) ! !履歴 2006/01/24 竹広真一 新規作成 ! 2006/01/25 竹広真一 モジュール初期化フラッグ導入 ! 2010/09/21 竹広真一 at_module と境界条件スイッチを揃えた ! module at_aq_galerkin_RRFF use lumatrix use at_module use dc_message implicit none private public :: at_aq_galerkin_RRFF_Initial ! 初期化サブルーチン public :: aq_at, q_t ! チェビシェフ->ガラーキン変換 public :: aq_ag, q_g ! 格子点->ガラーキンデータ変換 public :: at_aq, t_q ! ガラーキン->チェビシェフ変換 public :: ag_aq, g_q ! ガラーキン->格子点データ変換 public :: aq_Dx_aq, q_Dx_q ! X 微分 real(8), allocatable :: TQ(:,:) ! ガラーキン->チェビシェフ変換行列 real(8), allocatable :: QT(:,:) ! チェビシェフ->ガラーキン変換行列 integer, allocatable :: kp(:) ! LU 分解用ピボット格納配列 real(8), allocatable :: alpha(:) ! チェビシェフ<->ガラーキン行列積係数 real(8), allocatable :: beta(:) ! チェビシェフ<->ガラーキン行列和係数 integer :: im ! 格子点数 integer :: km ! チェビシェフ切断波数 integer,parameter :: ks=4 ! ガラーキン基底最低次数 character(LEN=2) :: bc ! 境界条件(RR/FF/RF/FR) logical :: RRFF_Initialized=.false. ! 初期化フラッグ save im, km, bc, TQ, QT, kp, alpha, beta, RRFF_Initialized public TQ, QT, alpha, beta contains ! ! チェビシェフ−ガラーキン法 ! 両端粘着条件用モジュール ! 初期化サブルーチン ! subroutine at_aq_galerkin_RRFF_Initial(i_in,k_in,bc_in) integer, intent(IN) :: i_in ! (in)格子点数 integer, intent(IN) :: k_in ! (in)チェビシェフ切断波数 character(LEN=2), intent(IN) :: bc_in ! (in)境界条件(RR/FF/RF/FR) integer :: k, l, m, n im=i_in ; km=k_in ; bc=bc_in allocate(TQ(0:km,ks:km),QT(ks:km,ks:km),kp(ks:km)) allocate(alpha(0:km),beta(0:km)) select case (bc) !---- 両端粘着条件用変換行列設定 case ('RR','rr') TQ = 0.0D0 do k=ks,km TQ(0,k) = ( (-1)**k + 1 ) * (k**2-4)/8.0D0 TQ(1,k) = -( 1 - (-1)**k ) * (9-k**2)/16.0D0 TQ(2,k) = -( (-1)**k + 1 ) * k**2 /8.0D0 TQ(3,k) = ( 1 - (-1)**k ) * (1-k**2)/16.0D0 TQ(k,k) = 1.0 enddo call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',& 'Conversion matrices for RR-B.C. initialized.') !---- 両端自由すべり条件用変換行列設定 case ('FF','ff','SS','ss') TQ = 0.0D0 do k=ks,km TQ(0,k) = ( (-1)**k + 1 ) * ( k**2*(k**2-1)-12 )/24.0D0 TQ(1,k) = ( 1 - (-1)**k ) * ( k**2*(k**2-1)-72 )/144.0D0 TQ(2,k) = -( (-1)**k + 1 ) * k**2*(k**2-1) /24.0D0 TQ(3,k) = -( 1 - (-1)**k ) * k**2*(k**2-1) /144.0D0 TQ(k,k) = 1.0 enddo call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',& 'Conversion matrices for FF-B.C. initialized.') !---- 片端粘着条件, 片端自由すべり条件変換行列設定 case ('RF','rf','RS','rs') TQ = 0.0D0 do k=ks,km TQ(0,k) = (-1)**k * k**2*(k**2-1)/48.0D0 & + 3.0D0/16.0D0 * k**2 & - 3.0D0/32.0D0 * ( 1 - (-1)**k ) & - 1.0D0/2.0D0 * ( 1 + (-1)**k ) TQ(1,k) = - (-1)**k * k**2*(k**2-1)/96.0D0 & + 1.0D0/32.0D0 * k**2 & - 33.0D0/64.0D0 * ( 1 - (-1)**k ) TQ(2,k) = - (-1)**k * k**2*(k**2-1)/48.0D0 & - 3.0D0/16.0D0 * k**2 & + 3.0D0/32.0D0 * ( 1 - (-1)**k ) TQ(3,k) = 1.0D0/32.0D0 * ( (-1)**k * k**2*(k**2-1)/3.0D0 & - k**2 & + ( 1 - (-1)**k )/2.0D0 ) TQ(k,k) = 1.0 enddo call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',& 'Conversion matrices for RF-B.C. initialized.') !---- 片端自由すべり条件, 片端粘着条件変換行列設定 case ('FR','fr','SR','sr') TQ = 0.0D0 do k=ks,km TQ(0,k) = k**2*(k**2-1)/48.0D0 & - 3.0D0/16.0D0 * (-1)**(k+1) * k**2 & + 3.0D0/32.0D0 * ( 1 - (-1)**k ) & - 1.0D0/2.0D0 * ( 1 + (-1)**k ) TQ(1,k) = k**2*(k**2-1)/96.0D0 & + 1.0D0/32.0D0 * (-1)**(k+1) * k**2 & - 33.0D0/64.0D0 * ( 1 - (-1)**k ) TQ(2,k) = - k**2*(k**2-1)/48.0D0 & + 3.0D0/16.0D0 * (-1)**(k+1) * k**2 & - 3.0D0/32.0D0 * ( 1 - (-1)**k ) TQ(3,k) = -1.0D0/32.0D0 * ( k**2*(k**2-1)/3.0D0 & + (-1)**(k+1) * k**2 & - ( 1 - (-1)**k )/2.0D0 ) TQ(k,k) = 1.0 enddo call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',& 'Conversion matrices for FR-B.C. initialized.') case default call MessageNotify('E','at_aq_galerkin_RRFF_Initial(Type2)',& 'Argument for B.C. not valid. Should be RR/FF/RF/FR.') end select beta=1.0 ; beta(0)=0.5D0 if (im .eq. km ) beta(km)=0.5D0 ! 両端粘着条件用変換逆行列 alpha=1.0 ; alpha(0)=2.0D0 QT = 0.0D0 do m=ks,km do n=ks,km do l=0,km QT(m,n) = QT(m,n) + alpha(l)*TQ(l,m)*TQ(l,n) enddo enddo enddo call LUDecomp(QT,kp) RRFF_Initialized=.true. end subroutine at_aq_galerkin_RRFF_Initial ! ! 両端粘着条件 ! チェビシェフ係数 -> ガラーキン係数変換(2次元データ) ! function aq_at(at_data) real(8), intent(IN) :: at_data(:,0:) !(in) チェビシェフ係数 real(8) :: aq_at(size(at_data,1),ks:km) !(out) ガラーキン係数 real(8) :: aq_work(size(at_data,1),ks:km) ! 作業用配列 integer :: k,m if ( .not. RRFF_Initialized ) & call MessageNotify('E','aq_at',& 'at_aq_galerkin_RRFF_module(Type2) not initialized') aq_work =0.0 do m=ks,km do k=0,km aq_work(:,m) = aq_work(:,m) & + alpha(k) * beta(k) * at_data(:,k) * TQ(k,m) enddo enddo aq_at = LUSolve(QT,kp,aq_work) end function aq_at ! ! 両端粘着条件 ! チェビシェフ係数 -> ガラーキン係数変換(1次元データ) ! function q_t(t_data) real(8), intent(IN) :: t_data(0:km) !(in) チェビシェフ係数 real(8) :: q_t(ks:km) !(out) ガラーキン係数 real(8) :: q_work(ks:km) ! 作業用配列 integer :: k,m if ( .not. RRFF_Initialized ) & call MessageNotify('E','q_t',& 'at_aq_galerkin_RRFF_module(Type2) not initialized') q_work =0.0 do m=ks,km do k=0,km q_work(m) = q_work(m) & + alpha(k) * beta(k) * t_data(k) * TQ(k,m) enddo enddo q_t = LUSolve(QT,kp,q_work) end function q_t ! ! 両端粘着条件 ! ガラーキン係数 -> チェビシェフ係数変換(2次元データ) ! function at_aq(aq_data) real(8), intent(IN) :: aq_data(:,ks:) !(in) ガラーキン係数 real(8) :: at_aq(size(aq_data,1),0:km) !(out) チェビシェフ係数 integer :: m, n if ( .not. RRFF_Initialized ) & call MessageNotify('E','at_aq',& 'at_aq_galerkin_RRFF_module(Type2) not initialized') at_aq = 0.0D0 do m=0,km do n=ks,km at_aq(:,m) = at_aq(:,m) + TQ(m,n)*aq_data(:,n)/beta(m) enddo enddo end function at_aq ! ! 両端粘着条件 ! ガラーキン係数 -> チェビシェフ係数変換(1次元データ) ! function t_q(q_data) real(8), intent(IN) :: q_data(ks:km) !(in) ガラーキン係数 real(8) :: t_q(0:km) !(out) チェビシェフ係数 integer :: m, n if ( .not. RRFF_Initialized ) & call MessageNotify('E','q_t',& 'at_aq_galerkin_RRFF_module(Type2) not initialized') t_q = 0.0D0 do m=0,km do n=ks,km t_q(m) = t_q(m) + TQ(m,n)*q_data(n)/beta(m) enddo enddo end function t_q ! ! 両端粘着条件 ! 格子点データ -> ガラーキン係数変換(2次元データ) ! function aq_ag(ag_data) real(8), intent(IN) :: ag_data(:,0:) !(in) 格子点データ real(8) :: aq_ag(size(ag_data,1),ks:km) !(out) ガラーキン係数 aq_ag = aq_at(at_ag(ag_data)) end function aq_ag ! ! 両端粘着条件 ! 格子点データ -> ガラーキン係数変換(1次元データ) ! function q_g(g_data) real(8), intent(IN) :: g_data(0:im) !(in) 格子点データ real(8) :: q_g(ks:km) !(out) ガラーキン係数 q_g = q_t(t_g(g_data)) end function q_g ! ! 両端粘着条件 ! ガラーキン係数 -> 格子点データ変換(2次元データ) ! function ag_aq(aq_data) real(8), intent(IN) :: aq_data(:,ks:) !(in) ガラーキン係数 real(8) :: ag_aq(size(aq_data,1),0:im) !(out) 格子点データ ag_aq = ag_at(at_aq(aq_data)) end function ag_aq ! ! 両端粘着条件 ! ガラーキン係数 -> 格子点データ変換(1次元データ) ! function g_q(q_data) real(8), intent(IN) :: q_data(ks:km) !(in) ガラーキン係数 real(8) :: g_q(0:im) !(out) 格子点データ g_q = g_t(t_q(q_data)) end function g_q ! ! 両端粘着条件 ! X 微分計算(1 次元) ! function aq_Dx_aq(aq_data) real(8), intent(IN) :: aq_data(:,ks:) !(in) ガラーキン係数 real(8) :: aq_Dx_aq(size(aq_data,1),ks:km) !(out) 微分ガラーキン aq_Dx_aq = aq_at(at_Dx_at(at_aq(aq_data))) end function aq_Dx_aq ! ! 両端粘着条件 ! X 微分計算(1 次元) ! function q_Dx_q(q_data) real(8), intent(IN) :: q_data(ks:km) real(8) :: q_Dx_q(ks:km) q_Dx_q = q_t(t_Dx_t(t_q(q_data))) end function q_Dx_q end module at_aq_galerkin_RRFF