!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wsc_mpi_module_supack テストプログラム :: スペクトル関数のテスト ! !履歴 2017/05/31 竹広真一 ! program wsc_mpi_module_supack_spectrum_test use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wsc_mpi_module_supack use mpi implicit none integer,parameter :: im=32, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: nm=10, lm=16 ! 切断波数の設定(水平, 動径) real(8),parameter :: ri=0.5, ro=1.5 ! 内外半径 real(8), allocatable :: xvz_Torvel(:,:,:) real(8), allocatable :: wc_Torvel(:,:) real(8), allocatable :: xvz_Polvel(:,:,:) real(8), allocatable :: ws_Polvel(:,:) real(8), allocatable :: xvz_Vlon(:,:,:) real(8), allocatable :: xvz_Vlat(:,:,:) real(8), allocatable :: xvz_Vrad(:,:,:) real(8), allocatable :: xvz_Xi(:,:,:) real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_EkTor real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_EkPol real(8), dimension(0:nm,0:km) :: nz_EkTor real(8), dimension(0:nm,0:km) :: nz_EkPol real(8), dimension(0:km) :: z_Data real(8), dimension(0:km) :: z_Xi ! 判定誤差設定 integer, parameter :: check_digits = 10 integer, parameter :: ignore = -11 real(8) :: PI integer :: n,m integer :: iproc, np, ierr PI = 4.0D0*Atan(1.0D0) !---------------- MPI スタート --------------------- call MPI_INIT(IERR) call MPI_COMM_RANK(MPI_COMM_WORLD,IPROC,IERR) call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR) call wsc_mpi_Initial(im,jm,km,nm,lm,ri,ro) call MessageNotify('M','wsc_mpi_module_supack_spectrum_test',& 'wsc_mpi_module_supack spectrum functions tests.') allocate(xvz_Torvel(0:im-1,jc,0:km)) allocate(wc_Torvel(nc,0:lm)) allocate(xvz_Polvel(0:im-1,jc,0:km)) allocate(ws_Polvel(nc,lm)) allocate(xvz_Vlon(0:im-1,jc,0:km)) allocate(xvz_Vlat(0:im-1,jc,0:km)) allocate(xvz_Vrad(0:im-1,jc,0:km)) allocate(xvz_Xi(0:im-1,jc,0:km)) xvz_xi = (xvz_Rad-ri)/(ro-ri) z_xi = xvz_xi(0,1,:) ! (2Y_1^0+3Y_2^1) xvz_Torvel =( 2*sqrt(3.0D0)*sin(xvz_Lat) & + 3*sqrt(5.0D0/6)*3.0*sin(xvz_Lat)*cos(xvz_Lat) * cos(xvz_Lon) ) ! (3*Y_1^0 + 2Y_1^1 + 2Y_2^0 +Y_2^-2) * sin(z) xvz_Polvel =( 3*sqrt(3.0D0)*sin(xvz_Lat) & + 2* sqrt(3.0D0/2)*cos(xvz_Lat)*cos(xvz_Lon) & + 2*sqrt(5.0D0)*(3.0/2*sin(xvz_Lat)**2-1/2.0) & - sqrt(5.0D0/24)*3.0*cos(xvz_Lat)**2*sin(2*xvz_Lon) ) & * sin(PI*xvz_xi) wc_Torvel = wc_xvz(xvz_Torvel) ws_Polvel = ws_xvz(xvz_Polvel) !--------- Energy spectrum --------- nmz_EkTor=wsc_VMiss nmz_EkPol=wsc_VMiss do n=0,nm do m=-n,n nmz_EkTor(n,m,:) = 0.0D0 nmz_EkPol(n,m,:) = 0.0D0 enddo enddo nmz_EkTor(1,0,:) = 4.0D0 * (4*pi)*Ro**2 nmz_EkTor(2,1,:) = 2*3*(3.0D0/2)**2/2 * (4*pi)*Ro**2 nmz_EkTor(2,-1,:)= nmz_EkTor(2,1,:) nmz_EkPol(1,0,:) = 9.0D0 * ((PI*Ro*cos(PI*z_xi))**2 + 2.0D0*sin(PI*z_xi)**2)*(4*pi) nmz_EkPol(1,1,:) = 1*2*(2/2)**2/2 & * ((PI*Ro*cos(PI*z_xi))**2 + 1*2*sin(PI*z_xi)**2) * (4*pi) nmz_EkPol(1,-1,:)= nmz_EkPol(1,1,:) nmz_EkPol(2,0,:) = 12.0D0 *((PI*Ro*cos(PI*z_xi))**2 + 6.0D0*sin(PI*z_xi)**2) * (4*pi) nmz_EkPol(2,-2,:)= 2*3*(1.0D0/2)**2/2 & *((PI*Ro*cos(PI*z_xi))**2 + 2*3*sin(PI*z_xi)**2) * (4*pi) nmz_EkPol(2,2,:) = nmz_EkPol(2,-2,:) call AssertEqual(& message='nmz_ToroidalEnergySpectrum_wc', & answer = nmz_ToroidalEnergySpectrum_wc(wc_Torvel), & check = nmz_EkTor, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='nmz_PoloidalEnergySpectrum_ws', & answer = nmz_PoloidalEnergySpectrum_ws(ws_Polvel), & check = nmz_EkPol, & significant_digits = check_digits, ignore_digits = ignore & ) nz_EkTor = 0.0D0 nz_EkTor(1,:) = 4.0D0 * (4*pi)*Ro**2 nz_EkTor(2,:) = 2 * 3 *(3.0D0/2)**2 * (4*pi)*Ro**2 nz_EkPol = 0.0D0 nz_EkPol(1,:) = (9.0D0+1*2*(2/2)**2/2*2 )*((PI*Ro*cos(PI*z_xi))**2+2*sin(PI*z_xi)**2)*(4*pi) nz_EkPol(2,:) = (12.0D0+2*3*(1.0D0/2)**2)*((PI*Ro*cos(PI*z_xi))**2+2*3*sin(PI*z_xi)**2)*(4*pi) call AssertEqual(& message='nz_ToroidalEnergySpectrum_wc', & answer = nz_ToroidalEnergySpectrum_wc(wc_Torvel), & check = nz_EkTor, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='nz_PoloidalEnergySpectrum_ws', & answer = nz_PoloidalEnergySpectrum_ws(ws_Polvel), & check = nz_EkPol, & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2VectorMPI(xvz_Vlon,xvz_Vlat,xvz_Vrad,wc_Torvel,ws_Polvel) nmz_EkTor = nmz_ToroidalEnergySpectrum_wc(wc_Torvel) nmz_EkPol = nmz_PoloidalEnergySpectrum_ws(ws_Polvel) z_Data = 0.0D0 do n=0,nm do m=-n,n z_Data = z_Data + nmz_EkTor(n,m,:) + nmz_EkPol(n,m,:) enddo enddo call AssertEqual(& message='Total energy by nmz-spectrum', & answer = IntLonLatRad_xvz(xvz_VLon**2+xvz_VLat**2+xvz_VRad**2)/2.0D0, & check = IntRad_z(z_Data/Ro**2), & significant_digits = check_digits, ignore_digits = ignore & ) nz_EkTor = nz_ToroidalEnergySpectrum_wc(wc_Torvel) nz_EkPol = nz_PoloidalEnergySpectrum_ws(ws_Polvel) call AssertEqual(& message='Total energy by nz-spectrum', & answer = IntLonLatRad_xvz(xvz_VLon**2+xvz_VLat**2+xvz_VRad**2)/2.0D0, & check = IntRad_z(sum(nz_EkTor,1)/Ro**2) & + IntRad_z(sum(nz_EkPol,1)/Ro**2), & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wsc_mpi_module_supack_spectrum_test',& 'wsc_mpi_module_supack spectrum functions tests suceeded!') CALL MPI_FINALIZE(IERR) end program wsc_mpi_module_supack_spectrum_test