!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wsc_mpi_module_supack テストプログラム ! ! 微分のテスト ! !履歴 2017/05/31 竹広真一 wt_module_nmomp_test_derivavit1.f90 より改造 ! program wsc_mpi_module_supack_derivative1_test use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wsc_mpi_module_supack use wa_mpi_module_supack, only : wa_Lapla_wa use mpi implicit none integer,parameter :: im=32, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: nmm=10, lm=16 ! 切断波数の設定(水平, 動径) real(8),parameter :: ri=0.5, ro=1.5 ! 内外半径 real(8), allocatable :: xvz_Data(:,:,:) real(8), allocatable :: xvz_DData(:,:,:) real(8), allocatable :: xvz_VLon(:,:,:) real(8), allocatable :: xvz_VLat(:,:,:) real(8), allocatable :: xvz_VRad(:,:,:) real(8), allocatable :: xvz_Div(:,:,:) real(8), allocatable :: xvz_Psi(:,:,:) real(8), allocatable :: xvz_Xi(:,:,:) ! 判定誤差設定 integer, parameter :: check_digits = 10 integer, parameter :: ignore = -11 integer, parameter :: nmin=1, nmax=9 integer :: iproc, np, ierr real(8) :: pi pi = atan(1.0D0)*4.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 MessageNotify('M','wsc_mpi_module_supack_derivative_test1', & 'wsc_mpi_module_supack derivative function test #1') call wsc_mpi_Initial(im,jm,km,nmm,lm,ri,ro) allocate(xvz_Data(0:im-1,1:jc,0:km)) allocate(xvz_DData(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_Div(0:im-1,1:jc,0:km)) allocate(xvz_Psi(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_DRad_ws, ws_DRad_wc ==================== xvz_Data = sin(PI*xvz_xi) xvz_DData = xvz_GradLat_wc(wc_xvz(xvz_Psi)) xvz_DData = PI*cos(PI*xvz_xi) ! call AssertEqual(& message='wc_DRad_ws with sin(z)', & answer = xvz_DData, & check = xvz_wc(wc_DRad_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(2*PI*xvz_xi) xvz_DData = -2*PI*sin(2*PI*xvz_xi) ! call AssertEqual(& message='ws_DRad_wc with cos(2z)', & answer = xvz_DData, & check = xvz_ws(ws_DRad_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== wc_DivRad_ws, ws_DivRad_wc ==================== xvz_Data = cos(xvz_Lat)*sin(xvz_Lat) * sin(xvz_Lon) * cos(PI*xvz_xi) ! Y_2^-1 cos(z) xvz_DData = -PI* cos(xvz_Lat)*sin(xvz_Lat) * sin(xvz_Lon) * sin(PI*xvz_xi) ! call AssertEqual(& message='ws_DivRad_wc with Y_2^-1 cos(z)', & answer = xvz_DData, & check = xvz_ws(ws_DivRad_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(xvz_Lat)*sin(xvz_Lat) * cos(xvz_Lon) * sin(2*PI*xvz_xi) ! Y_2^1 sin(2z) xvz_DData =2*PI*cos(xvz_Lat)*sin(xvz_Lat) * cos(xvz_Lon) * cos(2*PI*xvz_xi) ! call AssertEqual(& message='wc_DivRad_ws with Y_2^1 sin(2z)', & answer = xvz_DData, & check = xvz_wc(wc_DivRad_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== wc_RotRad_ws, ws_RotRad_wc ==================== xvz_Data = sin(xvz_Lat) * cos(3*PI*xvz_xi) ! Y_2^0 cos(3z) xvz_DData = -3*PI* sin(xvz_Lat) * sin(3*PI*xvz_xi) ! call AssertEqual(& message='ws_RotRad_wc with Y_2^0 cos(3z)', & answer = xvz_DData, & check = xvz_ws(ws_RotRad_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = sin(xvz_Lat) * sin(4*PI*xvz_xi) ! Y_2^0 sin(4z) xvz_DData = 4*PI* sin(xvz_Lat) * cos(4*PI*xvz_xi) ! call AssertEqual(& message='wc_RotRad_ws with Y_2^0 sin(4z)', & answer = xvz_DData, & check = xvz_wc(wc_RotRad_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== ws_Lapla_ws, wc_Lapla_wc ==================== xvz_Data = cos(xvz_Lat)*sin(xvz_Lon)*sin(3*PI*xvz_xi) ! Y_1^-1 sin(3z) xvz_DData = -( 9*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='ws_Lapla_ws with Y_1^-1 sin(3z)', & answer = xvz_DData, & check = xvz_ws(ws_Lapla_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(xvz_Lat)*cos(xvz_Lon)*cos(4*PI*xvz_xi) ! Y_1^1 cos(4z) xvz_DData = -( 16*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='wc_Lapla_wc with Y_1^1 cos(4z)', & answer = xvz_DData, & check = xvz_wc(wc_Lapla_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== ws_LaplaInv_ws, wc_LaplaInv_wc ==================== xvz_Data = cos(xvz_Lat)*sin(xvz_Lon)*sin(3*PI*xvz_xi) ! Y_1^-1 sin(3z) xvz_DData = -1.0D0/( 9*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='ws_LaplaInv_ws with Y_1^-1 sin(3z)', & answer = xvz_DData, & check = xvz_ws(ws_LaplaInv_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(xvz_Lat)*sin(xvz_Lon)*sin(2*PI*xvz_xi) ! Y_1^-1 sin(2z) xvz_DData = -1.0D0/( 4*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='ws_LaplaInv_ws with Y_1^-1 sin(2z)', & answer = xvz_DData, & check = xvz_ws(ws_LaplaInv_ws(ws_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(xvz_Lat)*cos(xvz_Lon)*cos(4*PI*xvz_xi) ! Y_1^1 cos(4z) xvz_DData = -1.0D0/( 16*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='wc_LaplaInv_wc with Y_1^1 cos(4z)', & answer = xvz_DData, & check = xvz_wc(wc_LaplaInv_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Data = cos(xvz_Lat)*cos(xvz_Lon)*cos(5*PI*xvz_xi) ! Y_1^1 cos(4z) xvz_DData = -1.0D0/( 25*PI**2 + 1*2/Ro**2) * xvz_Data ! call AssertEqual(& message='wc_LaplaInv_wc with Y_1^1 cos(5z)', & answer = xvz_DData, & check = xvz_wc(wc_LaplaInv_wc(wc_xvz(xvz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== xvz_GradLon_ws, xvz_GradLat_ws ========== ! ------ Y_1^-1 sin(z) Toridal field ----- xvz_Psi = cos(xvz_Lat)*sin(xvz_Lon) * sin(PI*xvz_xi) ! Y_1^-1 sin(z) xvz_VRad = 0 xvz_VLon = cos(xvz_Lon) * sin(PI*xvz_xi) /Ro call AssertEqual(& message='xvz_GradLon_ws with Y_1^-1 sin(z)', & answer = xvz_VLon, & check = xvz_GradLon_ws(ws_xvz(xvz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_VLat = -sin(xvz_Lat)*sin(xvz_Lon) * sin(PI*xvz_xi) /Ro call AssertEqual(& message='xvz_GradLat_ws with Y_1^-1 sin(z)', & answer = xvz_VLat, & check = xvz_GradLat_ws(ws_xvz(xvz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) !========== xvz_GradLon_wc, xvz_GradLat_wc, wc_Div_xvz_xvz_xvz, wc_RotRad_xvz_xvz ========== ! ------ Y_1^1 cos(2z) Toridal field ----- xvz_Psi = cos(xvz_Lat)*cos(xvz_Lon) * cos(2*PI*xvz_xi) ! Y_1^1 cos(2z) xvz_VRad = 0 xvz_VLon = - sin(xvz_Lon) * cos(2*PI*xvz_xi) /Ro call AssertEqual(& message='xvz_GradLon_wc with Y_1^1 cos(2z)', & answer = xvz_VLon, & check = xvz_GradLon_wc(wc_xvz(xvz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_VLat = -sin(xvz_Lat)*cos(xvz_Lon) * cos(2*PI*xvz_xi) /Ro call AssertEqual(& message='xvz_GradLat_wc with Y_1^1 cos(2z)', & answer = xvz_VLat, & check = xvz_GradLat_wc(wc_xvz(xvz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xvz_Div = - 2 * xvz_Psi/Ro**2 call AssertEqual(& message='wc_Div_xvz_xvz_xvz with Toroiral Y_1^-1 cos(2z)', & answer = xvz_Div, & check = xvz_wc(wc_Div_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) xvz_VLat = - xvz_Vlat call AssertEqual(& message='wc_RotRad_xvz_xvz with Toroiral Y_1^-1 cos(2z)', & answer = xvz_Div, & check = xvz_wc(wc_RotRad_xvz_xvz(xvz_VLon,xvz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) ! ------ Y_2^-1 sin(z) Poloidal field ----- xvz_Psi = cos(xvz_Lat)*sin(xvz_Lat) * sin(xvz_Lon) * sin(PI*xvz_xi) ! Y_2^-1 sin(z) xvz_VLon = xvz_GradLon_wc(wc_DRad_ws(ws_xvz(xvz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z) xvz_VLat = xvz_GradLat_wc(wc_DRad_ws(ws_xvz(xvz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z) xvz_Vrad = -xvz_wz(wa_Lapla_wa(wz_xvz(xvz_Psi)))/Ro ! (2*3/Ro) Psi xvz_DData = 0.0D0 call AssertEqual(& message='wc_Div_xvz_xvz_xvz with Y_2^-1 sin(z)', & answer = xvz_DData, & check = xvz_wc(wc_Div_xvz_xvz_xvz(xvz_Vlon,xvz_Vlat,xvz_Vrad)), & significant_digits = check_digits, ignore_digits = ignore & ) ! ------ Y_2^-1 cos(z) Poloidal field ----- xvz_Psi = cos(xvz_Lat)*sin(xvz_Lat) * sin(xvz_Lon) * cos(PI*xvz_xi) ! Y_2^-1 sin(z) xvz_VLon = xvz_GradLon_ws(ws_DRad_wc(wc_xvz(xvz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z) xvz_VLat = xvz_GradLat_ws(ws_DRad_wc(wc_xvz(xvz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z) xvz_Vrad = -xvz_wz(wa_Lapla_wa(wz_xvz(xvz_Psi)))/Ro ! (2*3/Ro) Psi xvz_DData = 0.0D0 call AssertEqual(& message='ws_Div_xvz_xvz_xvz with Y_2^-1 sin(z)', & answer = xvz_DData, & check = xvz_ws(ws_Div_xvz_xvz_xvz(xvz_Vlon,xvz_Vlat,xvz_Vrad)), & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wsc_mpi_module_supack_derivative1_test', & 'wsc_mpi_module_supack derivative function test #1 succeeded!') !------ MPIの終了 ------ call MPI_FINALIZE(IERR) end program wsc_mpi_module_supack_derivative1_test