!---------------------------------------------------------------------- ! Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 ut_mpi_module_supack テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2016/08/04 竹広真一 wt_mpi_module_supack_test_derivavit2.f90 より改造 ! program ut_mpi_module_svpack_derivative2_test use dc_message, only : MessageNotify use dc_test, only : AssertEqual use ut_mpi_module_svpack 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 :: ut_Psi(:,:) real(8), allocatable :: ut_Phi(:,:) real(8), allocatable :: xyh_DData(:,:,:) real(8), allocatable :: xyh_Psi(:,:,:) real(8), allocatable :: xyh_VLon(:,:,:) real(8), allocatable :: xyh_VLat(:,:,:) real(8), allocatable :: xyh_VRad(:,:,:) real(8), allocatable :: xyh_RadRot(:,:,:) real(8), allocatable :: xyh_RadRotRot(:,:,:) integer, parameter :: n=2 integer :: iproc, np, ierr !---------------- 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','ut_mpi_module_svpack_derivative2_test', & 'ut_mpi_module_svpack derivative function test #2') call ut_mpi_initial(im,jm,km,nm,lm,ri,ro) allocate(ut_Psi(nc,0:lm)) allocate(ut_Phi(nc,0:lm)) allocate(xyh_DData(0:im-1,1:jm,kc)) allocate(xyh_Psi(0:im-1,1:jm,kc)) allocate(xyh_VLon(0:im-1,1:jm,kc)) allocate(xyh_VLat(0:im-1,1:jm,kc)) allocate(xyh_VRad(0:im-1,1:jm,kc)) allocate(xyh_RadRot(0:im-1,1:jm,kc)) allocate(xyh_RadRotRot(0:im-1,1:jm,kc)) !========== ut_KxRGrad_ut, xyh_KGrad_ut, ut_QOperator_ut ========== ! ----------------- 例 1 -------------------- xyh_Psi = xyh_rad**n * cos(xyh_lat)*sin(xyh_lon) ! r**2 Y_1^-1 xyh_DData = xyh_rad**n * cos(xyh_lat)*cos(xyh_lon) call AssertEqual(& message='ut_KxRGrad_ut with r^2 Y_1^-1', & answer = xyh_DData, & check = xyh_ut(ut_KxRGrad_ut(ut_xyh(xyh_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1 xyh_DData = (n-1)*xyh_rad**(n-1)* cos(xyh_lat)*sin(xyh_lat)*sin(xyh_lon) call AssertEqual(& message='xyh_KGrad_ut with r^2 Y_1^-1', & answer = xyh_DData, & check = xyh_KGrad_ut(ut_xyh(xyh_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1 xyh_DData = - 3*(n-1)*xyh_rad**(n-1)* cos(xyh_lat)*sin(xyh_lat)*sin(xyh_lon) call AssertEqual(& message='ut_QOperator_ut r^2 Y_1^-1', & answer = xyh_DData, & check = xyh_ut(ut_QOperatorMPI_ut(ut_xyh(xyh_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- xyh_Psi = cos(xyh_lat)*sin(xyh_lat) * sin(xyh_lon) ! Y_2^1 xyh_DData = cos(xyh_lat)*sin(xyh_lat) * cos(xyh_lon) ! Y_2^1 call AssertEqual(& message='ut_KxRGrad_ut with Y_2^1', & answer = xyh_DData, & check = xyh_ut(ut_KxRGrad_ut(ut_xyh(xyh_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! k・▽ Y_2^1 = (- 4 Y_3^1 / 15 - Y_1^1 /5) xyh_DData = cos(2*xyh_lat)*cos(xyh_lat)*sin(xyh_lon)/xyh_rad call AssertEqual(& message='xyh_KGrad_ut with Y_2^1', & answer = xyh_DData, & check = xyh_KGrad_ut(ut_xyh(xyh_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = (16*sin(xyh_lat)**2 - 5)*cos(xyh_lat)*sin(xyh_lon)/xyh_rad call AssertEqual(& message='ut_QOperator_ut Y_2^1', & answer = xyh_DData, & check = xyh_ut(ut_QOperatorMPI_ut(ut_xyh(xyh_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== ut_RadRot_xyh_xyh, ut_RadRotRot_xyh_xyh_xyh ========== !---------- 単純流れ v_r=r ---------- xyh_VLon = 0 ; xyh_VLat = 0 ; xyh_VRad = xyh_Rad xyh_RadRot = 0 ; xyh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xyh_xyh with v_r=r', & answer = xyh_RadRot, & check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xyh_xyh with v_r=r', & answer = xyh_RadRotRot, & check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 剛体回転解非線形項(東西流れ) ---------- xyh_VLon = 0 xyh_VLat = xyh_Rad*sin(xyh_Lat)*cos(xyh_Lat) xyh_VRad = -xyh_Rad*cos(xyh_Lat)**2 xyh_RadRot = 0 ; xyh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xyh_xyh with rigid rotation(E-W)', & answer = xyh_RadRot, & check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xyh_xyh with rigid rotation(E-W)', & answer = xyh_RadRotRot, & check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 剛体回転解非線形項(南北流れ) ---------- xyh_VLon = xyh_Rad*cos(xyh_Lat)*sin(xyh_Lon)*cos(xyh_Lon) xyh_VLat = -xyh_Rad*sin(xyh_Lat)*cos(xyh_Lat)*sin(xyh_Lon)**2 xyh_VRad = -xyh_Rad*(sin(xyh_Lat)**2*sin(xyh_Lon)**2 + cos(xyh_Lon)**2) xyh_RadRot = 0 ; xyh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xyh_xyh with rigid rotation(N-S)', & answer = xyh_RadRot, & check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xyh_xyh with rigid rotation(N-S)', & answer = xyh_RadRotRot, & check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !--------- 鉛直渦度を伴うベクトル場 ---------- xyh_Psi = xyh_Rad**2 * cos(xyh_Lat)*sin(xyh_Lon) ! r**2 P_1^1 xyh_VLon = xyh_gRadLat_ut(ut_xyh(xyh_Psi*xyh_Rad)) xyh_VLat = - xyh_gRadLon_ut(ut_xyh(xyh_Psi*xyh_Rad)) xyh_VRad = 0 xyh_RadRot = 2 * xyh_Psi ! r・▽×▽×(ψr) = L_2ψ xyh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xyh_xyh with r^2 Y_1^-1 toroidal field', & answer = xyh_RadRot, & check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xyh_xyh with r^2 Y_1^-1 toroidal field',& answer = xyh_RadRotRot, & check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 鉛直速度を伴うベクトル場 ----------- ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ xyh_VRad = xyh_ut(ut_l2_ut(ut_xyh(xyh_Psi/xyh_Rad))) xyh_VLat = xyh_gRadLat_ut(ut_dRad_ut(ut_xyh(xyh_Psi*xyh_Rad))) xyh_VLon = xyh_gRadLon_ut(ut_dRad_ut(ut_xyh(xyh_Psi*xyh_Rad))) xyh_RadRot = 0 xyh_RadRotRot = -xyh_ut(ut_l2_ut(ut_lapla_ut(ut_xyh(xyh_Psi)))) ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ call AssertEqual(& message='ut_RadRot_xyh_xyh with non-vortical field', & answer = xyh_RadRot, & check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xyh_xyh with r^2 Y_1^-1 toroidal field',& answer = xyh_RadRotRot, & check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !========== ut_Potential2Vector, ut_Potential2Rotation ========== !----------------- 剛体回転場 -------------------- ! ut_Psi = ut_xyh(xyh_Rad * sin(xyh_Lat)) ut_Phi = 0.0D0 call ut_Potential2VectorMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = xyh_Rad * cos(xyh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with rigid rotation field(E-W)', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with rigid rotation field(E-W)', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with rigid rotation field(E-W)', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2*cos(xyh_LAT) call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with rigid rotation field(E-W)', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2*sin(xyh_LAT) call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with rigid rotation field(E-W)', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) !----------------- 剛体回転場(南北流) -------------------- ! ut_Psi = ut_xyh(xyh_Rad * cos(xyh_Lat) * sin(xyh_Lon)) ut_Phi = 0.0D0 call ut_Potential2VectorMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = -xyh_Rad*sin(xyh_Lat)*sin(xyh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = -xyh_Rad*cos(xyh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with rigid rotation field(N-S)', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 2*cos(xyh_Lon) call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with rigid rotation field(N-S)', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = -2*sin(xyh_Lon)*sin(xyh_Lat) call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2*sin(xyh_Lon)*cos(xyh_Lat) call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 渦無し場 -------------------- ut_Psi = 0.0D0 ut_Phi = ut_xyh(xyh_Rad * sin(xyh_Lat)) call ut_Potential2VectorMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 0 call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with non-vortical field', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2 * cos(xyh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with non-vortical field', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2 * sin(xyh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with non-vortical field', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with non-vortical field', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with non-vortical field', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with non-vortical field', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 ut_Psi = 0.0D0 ut_Phi = ut_xyh(xyh_Rad**2 * cos(xyh_Lat)*sin(xyh_Lon)) call ut_Potential2VectorMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 3 * xyh_Rad * cos(xyh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with poloidal field', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = - 3 * xyh_Rad * sin(xyh_Lat) * sin(xyh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with poloidal field', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 2 * xyh_Rad * cos(xyh_Lat) * sin(xyh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with poloidal field', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi ) xyh_DData = 4*sin(xyh_Lat)*sin(xyh_Lon) call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with poloidal field', & answer = xyh_VLon, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 4*cos(xyh_Lon) call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with poloidal field', & answer = xyh_VLat, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with poloidal field', & answer = xyh_VRad, & check = xyh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','ut_mpi_module_svpack_derivative_test2', & 'ut_mpi_module_svpack derivative function test #2 succeeded!') !------ MPIの終了 ------ call MPI_FINALIZE(IERR) end program ut_mpi_module_svpack_derivative2_test