!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wsc_module_svpack テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2017/05/20 竹広真一 ! program wsc_module_svpack_derivative2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wsc_module_svpack 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), dimension(0:im-1,jm,0:km) :: xyz_DData real(8), dimension(0:im-1,jm,0:km) :: xyz_Psi real(8), dimension((nm+1)**2,0:lm) :: wc_Psi real(8), dimension((nm+1)**2,lm) :: ws_Phi real(8), dimension(0:im-1,jm,0:km) :: xyz_VLon real(8), dimension(0:im-1,jm,0:km) :: xyz_VLat real(8), dimension(0:im-1,jm,0:km) :: xyz_VRad real(8), dimension(0:im-1,jm,0:km) :: xyz_RadRot real(8), dimension(0:im-1,jm,0:km) :: xyz_RadRotRot real(8), dimension(0:im-1,jm,0:km) :: xyz_xi real(8), parameter :: pi=3.1415926535897932385D0 call MessageNotify('M','wsc_module_svpack_derivative_test2', & 'wsc_module_svpack derivative function test #2') call wsc_initial(im,jm,km,nm,lm,ri,ro) xyz_xi = (xyz_Rad - ri)/(ro-ri) !========== wc_KxRGrad_wc, ws_KxRGrad_ws ========== ! ----------------- 例 1 -------------------- xyz_Psi = sin(Pi*xyz_xi) * cos(xyz_lat)*sin(xyz_lon) ! Y_1^-1 sin(z) xyz_DData = sin(Pi*xyz_xi) * cos(xyz_lat)*cos(xyz_lon) call AssertEqual(& message='ws_KxRGrad_ws with Y_1^-1 sin(z)', & answer = xyz_DData, & check = xyz_ws(ws_KxRGrad_ws(ws_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- xyz_Psi = cos(xyz_lat)*sin(xyz_lat) * sin(xyz_lon) * cos(PI*xyz_xi) ! Y_2^1 cos(z) xyz_DData = cos(xyz_lat)*sin(xyz_lat) * cos(xyz_lon) * cos(PI*xyz_xi) ! Y_2^1 call AssertEqual(& message='wc_KxRGrad_wc with Y_2^1 cos(z)', & answer = xyz_DData, & check = xyz_wc(wc_KxRGrad_wc(wc_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== wsc_Potential2Vector, wsc_Potential2Rotation, wc_RadRot, ws_RadRotRot ========== !----------------- 剛体回転場 -------------------- ! wc_Psi = wc_xyz(sin(xyz_Lat)) ws_Phi = 0.0D0 call wsc_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = cos(xyz_Lat) call AssertEqual(& message='wsc_Potential2Vector (Lon) with rigid rotation field(E-W)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Vector (Lat) with rigid rotation field(E-W)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Vector (Rad) with rigid rotation field(E-W)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRot = xyz_wc(wc_L2_wc(wc_Psi)) call AssertEqual(& message='wc_RadRot_xyz_xyz with with rigid rotation field(E-W)', & answer = xyz_RadRot, & check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRotRot = 0.0D0 call AssertEqual(& message='ws_RadRotRot_xyz_xyz with rigid rotation field(E-W)', & answer = xyz_RadRotRot, & check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Rotation (Lon) with rigid rotation field(E-W)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Rotation (Lat) with rigid rotation field(E-W)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = xyz_RadRot/Ro call AssertEqual(& message='wsc_Potential2Rotation (Rad) with rigid rotation field(E-W)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) !----------------- 剛体回転場(南北流) -------------------- ! wc_Psi = wc_xyz(cos(xyz_Lat) * sin(xyz_Lon)) ws_Phi = 0.0D0 call wsc_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = -sin(xyz_Lat)*sin(xyz_Lon) call AssertEqual(& message='wsc_Potential2Vector (Lon) with rigid rotation field(N-S)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = -cos(xyz_Lon) call AssertEqual(& message='wsc_Potential2Vector (Lat) with rigid rotation field(N-S)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Vector (Rad) with rigid rotation field(N-S)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRot = xyz_wc(wc_L2_wc(wc_Psi)) call AssertEqual(& message='wc_RadRot_xyz_xyz with with rigid rotation field(N-S)', & answer = xyz_RadRot, & check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRotRot = 0.0D0 call AssertEqual(& message='ws_RadRotRot_xyz_xyz with rigid rotation field(N-S)', & answer = xyz_RadRotRot, & check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Rotation (Lon) with rigid rotation field(N-S)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Rotation (Lat) with rigid rotation field(N-S)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = xyz_RadRot/Ro call AssertEqual(& message='wsc_Potential2Rotation (Rad) with rigid rotation field(N-S)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 渦無し場 -------------------- wc_Psi = 0.0D0 ws_Phi = ws_xyz(sin(PI*xyz_xi) * sin(xyz_Lat)) call wsc_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = 0 call AssertEqual(& message='wsc_Potential2Vector (Lon) with non-vortical field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = PI * cos(xyz_Lat) * cos(PI*xyz_xi) call AssertEqual(& message='wsc_Potential2Vector (Lat) with non-vortical field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2 * sin(xyz_Lat) * sin(PI*xyz_xi)/Ro call AssertEqual(& message='wsc_Potential2Vector (Rad) with non-vortical field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRot = 0.0D0 call AssertEqual(& message='wc_RadRot_xyz_xyz with non-vortical field', & answer = xyz_RadRot, & check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRotRot = - xyz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi)))) call AssertEqual(& message='ws_RadRotRot_xyz_xyz with non-vortical field', & answer = xyz_RadRotRot, & check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = (2/Ro**2 + PI**2) * cos(xyz_Lat) * sin(PI*xyz_xi) call AssertEqual(& message='wsc_Potential2Rotation (Lon) with non-vortical field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wsc_Potential2Rotation (Lat) with non-vortical field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = xyz_RadRot/Ro call AssertEqual(& message='wsc_Potential2Rotation (Rad) with non-vortical field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wc_Psi = 0.0D0 ws_Phi = ws_xyz(sin(2*PI*xyz_xi) * cos(xyz_Lat)*sin(xyz_Lon)) call wsc_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = 2*PI * cos(2*PI*xyz_xi) * cos(xyz_Lon) call AssertEqual(& message='wsc_Potential2Vector (Lon) with poloidal field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = - 2*PI* cos(2*PI*xyz_xi) * sin(xyz_Lat)*sin(xyz_Lon) call AssertEqual(& message='wsc_Potential2Vector (Lat) with poloidal field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2*sin(2*PI*xyz_xi) * cos(xyz_Lat)*sin(xyz_Lon)/Ro call AssertEqual(& message='wsc_Potential2Vector (Rad) with poloidal field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRot = 0.0D0 call AssertEqual(& message='wc_RadRot_xyz_xyz with with poloidal field', & answer = xyz_RadRot, & check = xyz_wc(wc_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_RadRotRot = - xyz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi)))) call AssertEqual(& message='ws_RadRotRot_xyz_xyz with with poloidal field', & answer = xyz_RadRotRot, & check = xyz_ws(ws_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call wsc_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wc_Psi, ws_Phi ) xyz_DData = -(2/Ro**2 + 4*PI**2)*sin(2*PI*xyz_xi)*sin(xyz_Lat)*sin(xyz_Lon) call AssertEqual(& message='wsc_Potential2Rotation (Lon) with poloidal field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = -(2/Ro**2+4*PI**2) * sin(2*PI*xyz_xi) * cos(xyz_Lon) call AssertEqual(& message='wsc_Potential2Rotation (Lat) with poloidal field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = xyz_RadRot/Ro call AssertEqual(& message='wsc_Potential2Rotation (Rad) with poloidal field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wsc_module_svpack_derivative_test2', & 'wsc_module_svpack derivative function test #2 succeeded!') end program wsc_module_svpack_derivative2