!---------------------------------------------------------------------- ! Copyright (c) 2012 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wq_module_sjpack テストプログラム ! ! wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq のテスト ! wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr のテスト ! wq_Potential2Vector, wq_Potential2Rotation のテスト ! !履歴 2011/09/08 竹広真一 wq_module_deriv2_test.f90 を MPI 用に改造 ! 2011/09/12 竹広真一 wtq_mpi_module 用に改造 ! 2012/04/02 竹広真一 wtq_mpi_module_sjpack 用に改造 ! 2012/04/03 竹広真一 wtq_module_sjpack 用に改造 ! program wtq_module_sjpack_deriv_wq2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wtq_module_sjpack implicit none integer,parameter :: im=32, jm=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: kmo=16, kmi=21 ! 格子点の設定(球殻動径, 球動径) integer,parameter :: nm=10 ! 切断波数の設定(水平) integer,parameter :: lmo=10, lmi=32 ! 切断波数の設定(球殻動径, 球動径) real(8),parameter :: ri=0.5, ro=1.5 ! 内外半径 real(8), dimension(0:im-1,1:jm,kmi) :: xyr_Psi real(8), dimension(0:im-1,1:jm,kmi) :: xyr_KxRGrad real(8), dimension(0:im-1,1:jm,kmi) :: xyr_KGrad real(8), dimension(0:im-1,1:jm,kmi) :: xyr_QOperator real(8), dimension(0:im-1,1:jm,kmi) :: xyr_VLon real(8), dimension(0:im-1,1:jm,kmi) :: xyr_VLat real(8), dimension(0:im-1,1:jm,kmi) :: xyr_VRad real(8), dimension(0:im-1,1:jm,kmi) :: xyr_V0Lon ! 速度(正解, 経度) real(8), dimension(0:im-1,1:jm,kmi) :: xyr_V0Lat ! 速度(正解, 緯度) real(8), dimension(0:im-1,1:jm,kmi) :: xyr_V0rad ! 速度(正解,動径) real(8), dimension(0:im-1,1:jm,kmi) :: xyr_RadRot real(8), dimension(0:im-1,1:jm,kmi) :: xyr_RadRotRot real(8), dimension(0:im-1,1:jm,kmi) :: xyr_Torvel real(8), dimension(0:im-1,1:jm,kmi) :: xyr_Polvel real(8), dimension((nm+1)*(nm+1),0:lmi) :: wq_Torvel ! トロイダルポテンシャル real(8), dimension((nm+1)*(nm+1),0:lmi) :: wq_Polvel ! ポロイダルポテンシャル integer, parameter :: n=3 ! 判定誤差設定 integer, parameter :: check_digits = 8 integer, parameter :: ignore = -7 call MessageNotify('M','wtq_module_sjpack_deriv_wq_test2', & 'wtq_module_sjpack wq-derivative function test #2') call wtq_Initial(im,jm,kmi,kmo,nm,lmi,lmo,ri,ro) !========== wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq ========= call MessageNotify('M','wtq_module_sjpack_derivative_test2', & 'Test for wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq') xyr_Psi = xyr_rad**n * cos(xyr_lat)*sin(xyr_lon) ! r**n Y_1^1 xyr_KxRGrad = xyr_Rad**n * cos(xyr_lat)*cos(xyr_lon) ! r**n Y_1^-1 ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1 xyr_KGrad = (n-1)*xyr_rad**(n-1)* cos(xyr_lat)*sin(xyr_lat)*sin(xyr_lon) ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1 xyr_QOperator = - 3*(n-1)*xyr_rad**(n-1)* & cos(xyr_lat)*sin(xyr_lat)*sin(xyr_lon) call AssertEqual(& message='k x r grad with r**n Y_1^1', & answer = xyr_KxRgrad, & check = xyr_wq(wq_KxRGrad_wq(wq_xyr(xyr_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='k grad with r**n Y_1^1', & answer = xyr_KGrad, & check = xyr_KGrad_wq(wq_xyr(xyr_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Q operator with r**n Y_1^1', & answer = xyr_QOperator, & check = xyr_wq(wq_QOperator_wq(wq_xyr(xyr_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) xyr_Psi = 3*cos(xyr_lat)*sin(xyr_lat) * sin(xyr_lon) * xyr_Rad**2 ! Y_2^1 xyr_KxRGrad = 3*cos(xyr_lat)*sin(xyr_lat) * cos(xyr_lon) * xyr_Rad**2 ! Y_2^-1 ! k・▽r^2 Y_2^1 = 3r Y_1^1 xyr_KGrad = 3*xyr_Rad*cos(2*xyr_Lat)*cos(xyr_Lat)*sin(xyr_Lon) & +3*xyr_Rad*sin(2*xyr_Lat)*sin(xyr_Lat)*sin(xyr_Lon) xyr_QOperator = -9*cos(xyr_lat)*sin(xyr_lon)*xyr_rad call AssertEqual(& message='k x r grad with r**2 Y_2^1', & answer = xyr_KxRgrad, & check = xyr_wq(wq_KxRGrad_wq(wq_xyr(xyr_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='k grad with r**2 Y_2^1', & answer = xyr_KGrad, & check = xyr_KGrad_wq(wq_xyr(xyr_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Q operator with r**2 Y_2^1', & answer = xyr_QOperator, & check = xyr_wq(wq_QOperator_wq(wq_xyr(xyr_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) !=========== wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr ============ call MessageNotify('M','wtq_module_sjpack_derivative_wq_test2', & 'Test for wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr') xyr_VLon = 0.0D0 ; xyr_VLat = 0.0D0 ; xyr_VRad = xyr_Rad xyr_RadRot = 0.0D0 ; xyr_RadRotRot = 0.0D0 call AssertEqual(& message='r Rot v with simple example (v_r=r)', & answer = xyr_RadRot, & check = xyr_wr(wr_RadRot_xyr_xyr(xyr_VLon,xyr_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='r Rot Rot v with simple example (v_r=r)', & answer = xyr_RadRotRot, & check = xyr_wr(wr_RadRotRot_xyr_xyr_xyr(xyr_VLon,xyr_VLat,xyr_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) ! 剛体回転解非線形項(東西流れ) xyr_VLon = 0.0D0 xyr_VLat = xyr_Rad*sin(xyr_Lat)*cos(xyr_Lat) xyr_VRad = -xyr_Rad*cos(xyr_Lat)**2 xyr_RadRot = 0.0D0 ; xyr_RadRotRot = 0.0D0 call AssertEqual(& message='r Rot v with rigid rotation(E-W flow)', & answer = xyr_RadRot, & check = xyr_wr(wr_RadRot_xyr_xyr(xyr_VLon,xyr_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='r Rot Rot v with rigid rotation(E-W flow)', & answer = xyr_RadRotRot, & check = xyr_wr(wr_RadRotRot_xyr_xyr_xyr(xyr_VLon,xyr_VLat,xyr_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) ! 剛体回転解非線形項(南北流れ) xyr_VLon = xyr_Rad*cos(xyr_Lat)*sin(xyr_Lon)*cos(xyr_Lon) xyr_VLat = -xyr_Rad*sin(xyr_Lat)*cos(xyr_Lat)*sin(xyr_Lon)**2 xyr_VRad = -xyr_Rad*(sin(xyr_Lat)**2*sin(xyr_Lon)**2 + cos(xyr_Lon)**2) xyr_RadRot = 0.0D0 ; xyr_RadRotRot = 0.0D0 call AssertEqual(& message='r Rot v with rigid rotation(N-S flow)', & answer = xyr_RadRot, & check = xyr_wr(wr_RadRot_xyr_xyr(xyr_VLon,xyr_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='r Rot Rot v with rigid rotation(N-S flow)', & answer = xyr_RadRotRot, & check = xyr_wr(wr_RadRotRot_xyr_xyr_xyr(xyr_VLon,xyr_VLat,xyr_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) ! 鉛直渦度を伴うベクトル場 xyr_Polvel = xyr_Rad**2 * cos(xyr_Lat)*sin(xyr_Lon) ! r**2 P_1^1 !xyr_Psi = xyr_Rad**3 * cos(xyr_Lat)*sin(xyr_Lat)*sin(xyr_Lon) ! r**3 P_2^1 xyr_VLon = xyr_GradLat_wq(wq_xyr(xyr_Polvel*xyr_Rad)) xyr_VLat = - xyr_GradLon_wq(wq_xyr(xyr_Polvel*xyr_Rad)) xyr_VRad = 0.0D0 xyr_RadRot = 2 * xyr_Polvel ! r・▽×▽×(ψr) = L_2ψ !xyr_RadRot = 6 * xyr_Polvel ! r・▽×▽×(ψr) = L_2ψ xyr_RadRotRot = 0.0D0 call AssertEqual(& message='r Rot v with vortical field', & answer = xyr_RadRot, & check = xyr_wr(wr_RadRot_xyr_xyr(xyr_VLon,xyr_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='r Rot Rot v with vortical field', & answer = xyr_RadRotRot, & check = xyr_wr(wr_RadRotRot_xyr_xyr_xyr(xyr_VLon,xyr_VLat,xyr_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) ! 鉛直速度を伴うベクトル場 xyr_Torvel = 0.0D0 xyr_Polvel = xyr_Rad**3 * cos(xyr_Lat)*sin(xyr_Lon) ! r**2 P_1^1 call wq_Potential2Vector(& xyr_VLon,xyr_VLat,xyr_VRad, wq_xyr(xyr_Torvel),wq_xyr(xyr_Polvel) ) xyr_RadRot = 0.0D0 xyr_RadRotRot = -xyr_wq(wq_L2_wq(wq_Lapla_wq(wq_xyr(xyr_Polvel)))) ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ call AssertEqual(& message='r Rot v with non-vortical field', & answer = xyr_RadRot, & check = xyr_wr(wr_RadRot_xyr_xyr(xyr_VLon,xyr_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='r Rot Rot v with non-vortical field', & answer = xyr_RadRotRot, & check = xyr_wr(wr_RadRotRot_xyr_xyr_xyr(xyr_VLon,xyr_VLat,xyr_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) !==================== wq_Potential2Vector ==================== call MessageNotify('M','wtq_module_sjpack_derivative_wq_test2', & 'Test for wq_Potential2Vector, wq_Potential2VectorMPI') ! ----------------- 例 1 -------------------- ! 剛体回転場 wq_Torvel = wq_xyr(xyr_Rad * sin(xyr_Lat)) wq_Polvel = 0.0D0 xyr_V0lon = xyr_Rad * cos(xyr_Lat) xyr_V0lat = 0.0D0 xyr_V0Rad = 0.0D0 call wq_Potential2Vector(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Vector(VLon) with rigid rotation(EW)',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VLat) with rigid rotation(EW)',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VRad) with rigid rotation(EW)',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- ! 剛体回転場(南北流) wq_Torvel = wq_xyr(xyr_Rad * cos(xyr_Lat) * sin(xyr_Lon)) wq_Polvel = 0.0D0 xyr_V0Lon = -xyr_Rad*sin(xyr_Lat)*sin(xyr_Lon) xyr_V0Lat = -xyr_Rad*cos(xyr_Lon) xyr_V0Rad = 0.0D0 call wq_Potential2Vector(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Vector(VLon) with rigid rotation(NS)',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VLat) with rigid rotation(NS)',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VRad) with rigid rotation(NS)',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 3 -------------------- ! 渦無し場 wq_Torvel = 0.0D0 wq_Polvel = wq_xyr(xyr_Rad**3 * sin(xyr_Lat)) xyr_V0lon = 0.0D0 xyr_V0lat = 4 * xyr_Rad**2 * cos(xyr_Lat) xyr_V0rad = 2 * xyr_Rad**2 * sin(xyr_Lat) call wq_Potential2Vector(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Vector(VLon) with no rotation', & answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VLat) with no rotation', & answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VRad) with no rotation', & answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wq_Torvel = 0.0D0 wq_Polvel = wq_xyr(xyr_Rad**5 * cos(xyr_Lat)*sin(xyr_Lon)) xyr_V0lon = 6 * xyr_Rad**4 * cos(xyr_Lon) xyr_V0lat = - 6 * xyr_Rad**4 * sin(xyr_Lat) * sin(xyr_Lon) xyr_V0rad = 2 * xyr_Rad**4 * cos(xyr_Lat) * sin(xyr_Lon) call wq_Potential2Vector(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Vector(VLat) with poloidal field',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VLat) with poloidal field',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Vector(VRad) with poloidal field',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) !==================== wq_Potential2Rotation ==================== call MessageNotify('M','wtq_module_sjpack_derivative_test2', & 'Test for wq_Potential2Rotation, wq_Potential2RotationMPI') ! ----------------- 例 1 -------------------- ! 剛体回転場 wq_Torvel = wq_xyr(xyr_Rad * sin(xyr_Lat)) wq_Polvel = 0.0D0 ! xyr_Vlon = xyr_Rad * cos(xyr_Lat)' ! xyr_Vlat = 0.0' ! xyr_VRad = 0.0' xyr_V0lon = 0.0D0 xyr_V0lat = 2*cos(xyr_LAT) xyr_V0rad = 2*sin(xyr_LAT) call wq_Potential2Rotation(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Rotation(VLon) with rigid rotation(EW)',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VLat) with rigid rotation(EW)',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VRad) with rigid rotation(EW)',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- ! 剛体回転場(南北流) wq_Torvel = wq_xyr(xyr_Rad * cos(xyr_Lat) * sin(xyr_Lon)) wq_Polvel = 0.0D0 ! xyr_VLon = -xyr_Rad*sin(xyr_Lat)*sin(xyr_Lon) ! xyr_VLat = -xyr_Rad*cos(xyr_Lon) ! xyr_VRad = 0.0 xyr_V0Lon = 2*cos(xyr_Lon) xyr_V0Lat = -2*sin(xyr_Lon)*sin(xyr_Lat) xyr_V0Rad = 2*sin(xyr_Lon)*cos(xyr_Lat) call wq_Potential2Rotation(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Rotation(VLon) with rigid rotation(NS)',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VLat) with rigid rotation(NS)',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VRad) with rigid rotation(NS)',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 3 -------------------- ! 渦無し場 wq_Torvel = 0.0D0 wq_Polvel = wq_xyr(xyr_Rad * sin(xyr_Lat)) ! xyr_Vlon = 0 ! xyr_Vlat = 2 * cos(xyr_Lat) ! xyr_Vrad = 2 * sin(xyr_Lat) xyr_V0lon = 0.0D0 xyr_V0lat = 0.0D0 xyr_V0rad = 0.0D0 call wq_Potential2Rotation(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Rotation(VLon) with no rotation', & answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VLat) with no rotation', & answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VRad) with no rotation', & answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wq_Torvel = 0.0D0 wq_Polvel = wq_xyr(xyr_Rad **3 * cos(xyr_Lat)*sin(xyr_Lon)) ! xyr_Vlon = 4 * xyr_Rad**2 * cos(xyr_Lon) ! xyr_Vlat = -4 * xyr_Rad**2 * sin(xyr_Lat) * sin(xyr_Lon) ! xyr_Vrad = 2 * xyr_Rad**2 * cos(xyr_Lat) * sin(xyr_Lon) xyr_V0lon = 10*xyr_Rad*sin(xyr_Lat)*sin(xyr_Lon) xyr_V0lat = 10*xyr_Rad*cos(xyr_Lon) xyr_V0rad = 0.0D0 call wq_Potential2Rotation(& xyr_VLon,xyr_VLat,xyr_VRad, wq_Torvel, wq_Polvel ) call AssertEqual(& message='Potential2Rotation(VLon) with poloidal field',& answer = xyr_VLon, & check = xyr_V0Lon, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VLat) with poloidal field',& answer = xyr_VLat, & check = xyr_V0Lat, & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='Potential2Rotation(VRad) with poloidal field',& answer = xyr_VRad, & check = xyr_V0Rad, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wtq_module_sjpack_derivative_wq_test', & 'wtq_module_sjpack wq-derivative function test #2 succeeded!') end program wtq_module_sjpack_deriv_wq2