!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wsc_mpi_module_supack テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2017/05/31 竹広真一 wt_module_nmomp_test_derivavit2.f90 より改造 ! program wsc_mpi_module_supack_derivative2_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 ! 内外半径 ! 判定誤差設定 integer, parameter :: check_digits = 8 integer, parameter :: ignore = -9 real(8), allocatable :: wc_Psi(:,:) real(8), allocatable :: ws_Phi(:,:) real(8), allocatable :: xvz_DData(:,:,:) real(8), allocatable :: xvz_Psi(:,:,:) real(8), allocatable :: xvz_VLon(:,:,:) real(8), allocatable :: xvz_VLat(:,:,:) real(8), allocatable :: xvz_VRad(:,:,:) real(8), allocatable :: xvz_RadRot(:,:,:) real(8), allocatable :: xvz_RadRotRot(:,:,:) real(8), allocatable :: xvz_Xi(:,:,:) integer, parameter :: n=2 integer :: iproc, np, ierr real(8), parameter :: pi=3.1415926535897932385D0 !---------------- MPI スタート --------------------- call MPI_INIT(IERR) call MPI_COMM_RANK(MPI_COMM_WORLD,IPROC,IERR) call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR) call MessageNotify('M','wsc_mpi_module_supack_derivative2_test', & 'wsc_mpi_module_supack derivative function test #2') call wsc_mpi_initial(im,jm,km,nm,lm,ri,ro) allocate(wc_Psi(nc,0:lm)) allocate(ws_Phi(nc,lm)) allocate(xvz_DData(0:im-1,1:jc,0:km)) allocate(xvz_Psi(0:im-1,1:jc,0:km)) allocate(xvz_VLon(0:im-1,1:jc,0:km)) allocate(xvz_VLat(0:im-1,1:jc,0:km)) allocate(xvz_VRad(0:im-1,1:jc,0:km)) allocate(xvz_RadRot(0:im-1,1:jc,0:km)) allocate(xvz_RadRotRot(0:im-1,1:jc,0:km)) allocate(xvz_Xi(0:im-1,1:jc,0:km)) xvz_xi = (xvz_Rad - ri)/(ro-ri) !========== wc_KxRGrad_wc, ws_KxRGrad_ws ========== ! ----------------- 例 1 -------------------- xvz_Psi = sin(Pi*xvz_xi) * cos(xvz_lat)*sin(xvz_lon) ! Y_1^-1 sin(z) xvz_DData = sin(Pi*xvz_xi) * cos(xvz_lat)*cos(xvz_lon) call AssertEqual(& message='ws_KxRGrad_ws with Y_1^-1 sin(z)', & answer = xvz_DData, & check = xvz_ws(ws_KxRGrad_ws(ws_xvz(xvz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- xvz_Psi = cos(xvz_lat)*sin(xvz_lat) * sin(xvz_lon) * cos(PI*xvz_xi) ! Y_2^1 cos(z) xvz_DData = cos(xvz_lat)*sin(xvz_lat) * cos(xvz_lon) * cos(PI*xvz_xi) ! Y_2^1 call AssertEqual(& message='wc_KxRGrad_wc with Y_2^1 cos(z)', & answer = xvz_DData, & check = xvz_wc(wc_KxRGrad_wc(wc_xvz(xvz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== wsc_Potential2VectorMPI, wsc_Potential2RotationMPI, wc_RadRot, ws_RadRotRot ========== !----------------- 剛体回転場 -------------------- ! wc_Psi = wc_xvz(sin(xvz_Lat)) ws_Phi = 0.0D0 call wsc_Potential2VectorMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = cos(xvz_Lat) call AssertEqual(& message='wsc_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2VectorMPI (Lat) with rigid rotation field(E-W)', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2VectorMPI (Rad) with rigid rotation field(E-W)', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRot = xvz_wc(wc_L2_wc(wc_Psi)) call AssertEqual(& message='wc_RadRot_xvz_xvz with with rigid rotation field(E-W)', & answer = xvz_RadRot, & check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRotRot = 0.0D0 call AssertEqual(& message='ws_RadRotRot_xvz_xvz with rigid rotation field(E-W)', & answer = xvz_RadRotRot, & check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2RotationMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2RotationMPI (Lon) with rigid rotation field(E-W)', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2RotationMPI (Lat) with rigid rotation field(E-W)', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = xvz_RadRot/Ro call AssertEqual(& message='wsc_Potential2RotationMPI (Rad) with rigid rotation field(E-W)', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) !----------------- 剛体回転場(南北流) -------------------- ! wc_Psi = wc_xvz(cos(xvz_Lat) * sin(xvz_Lon)) ws_Phi = 0.0D0 call wsc_Potential2VectorMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = -sin(xvz_Lat)*sin(xvz_Lon) call AssertEqual(& message='wsc_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = -cos(xvz_Lon) call AssertEqual(& message='wsc_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2VectorMPI (Rad) with rigid rotation field(N-S)', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRot = xvz_wc(wc_L2_wc(wc_Psi)) call AssertEqual(& message='wc_RadRot_xvz_xvz with with rigid rotation field(N-S)', & answer = xvz_RadRot, & check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRotRot = 0.0D0 call AssertEqual(& message='ws_RadRotRot_xvz_xvz with rigid rotation field(N-S)', & answer = xvz_RadRotRot, & check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2RotationMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2RotationMPI (Lon) with rigid rotation field(N-S)', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = xvz_RadRot/Ro call AssertEqual(& message='wsc_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 渦無し場 -------------------- wc_Psi = 0.0D0 ws_Phi = ws_xvz(sin(PI*xvz_xi) * sin(xvz_Lat)) call wsc_Potential2VectorMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = 0 call AssertEqual(& message='wsc_Potential2VectorMPI (Lon) with non-vortical field', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = PI * cos(xvz_Lat) * cos(PI*xvz_xi) call AssertEqual(& message='wsc_Potential2VectorMPI (Lat) with non-vortical field', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 2 * sin(xvz_Lat) * sin(PI*xvz_xi)/Ro call AssertEqual(& message='wsc_Potential2VectorMPI (Rad) with non-vortical field', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRot = 0.0D0 call AssertEqual(& message='wc_RadRot_xvz_xvz with non-vortical field', & answer = xvz_RadRot, & check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRotRot = - xvz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi)))) call AssertEqual(& message='ws_RadRotRot_xvz_xvz with non-vortical field', & answer = xvz_RadRotRot, & check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2RotationMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = (2/Ro**2 + PI**2) * cos(xvz_Lat) * sin(PI*xvz_xi) call AssertEqual(& message='wsc_Potential2RotationMPI (Lon) with non-vortical field', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2RotationMPI (Lat) with non-vortical field', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = xvz_RadRot/Ro call AssertEqual(& message='wsc_Potential2RotationMPI (Rad) with non-vortical field', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wc_Psi = 0.0D0 ws_Phi = ws_xvz(sin(2*PI*xvz_xi) * cos(xvz_Lat)*sin(xvz_Lon)) call wsc_Potential2VectorMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = 2*PI * cos(2*PI*xvz_xi) * cos(xvz_Lon) call AssertEqual(& message='wsc_Potential2VectorMPI (Lon) with poloidal field', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = - 2*PI* cos(2*PI*xvz_xi) * sin(xvz_Lat)*sin(xvz_Lon) call AssertEqual(& message='wsc_Potential2VectorMPI (Lat) with poloidal field', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = 2*sin(2*PI*xvz_xi) * cos(xvz_Lat)*sin(xvz_Lon)/Ro call AssertEqual(& message='wsc_Potential2VectorMPI (Rad) with poloidal field', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRot = 0.0D0 call AssertEqual(& message='wc_RadRot_xvz_xvz with with poloidal field', & answer = xvz_RadRot, & check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_RadRotRot = - xvz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi)))) call AssertEqual(& message='ws_RadRotRot_xvz_xvz with with poloidal field', & answer = xvz_RadRotRot, & check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2RotationMPI(& xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi ) xvz_DData = -(2/Ro**2 + 4*PI**2)*sin(2*PI*xvz_xi)*sin(xvz_Lat)*sin(xvz_Lon) call AssertEqual(& message='wsc_Potential2RotationMPI (Lon) with poloidal field', & answer = xvz_VLon, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = -(2/Ro**2+4*PI**2) * sin(2*PI*xvz_xi) * cos(xvz_Lon) call AssertEqual(& message='wsc_Potential2RotationMPI (Lat) with poloidal field', & answer = xvz_VLat, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvz_DData = xvz_RadRot/Ro call AssertEqual(& message='wsc_Potential2RotationMPI (Rad) with poloidal field', & answer = xvz_VRad, & check = xvz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wsc_mpi_module_supack_derivative_test2', & 'wsc_mpi_module_supack derivative function test #2 succeeded!') !------ MPIの終了 ------ call MPI_FINALIZE(IERR) end program wsc_mpi_module_supack_derivative2_test