!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wt_module_svpack テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2017/05/24 竹広真一 ! program wt_module_svpack_derivative2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wt_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) :: wt_Psi real(8), dimension((nm+1)**2,0:lm) :: wt_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 integer, parameter :: n=2 call MessageNotify('M','wt_module_svpack_derivative_test2', & 'wt_module_svpack derivative function test #2') call wt_initial(im,jm,km,nm,lm,ri,ro) !========== wt_KxRGrad_wt, xyz_KGrad_wt, wt_QOperator_wt ========== ! ----------------- 例 1 -------------------- xyz_Psi = xyz_rad**n * cos(xyz_lat)*sin(xyz_lon) ! r**2 Y_1^-1 xyz_DData = xyz_rad**n * cos(xyz_lat)*cos(xyz_lon) call AssertEqual(& message='wt_KxRGrad_wt with r^2 Y_1^-1', & answer = xyz_DData, & check = xyz_wt(wt_KxRGrad_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1 xyz_DData = (n-1)*xyz_rad**(n-1)* cos(xyz_lat)*sin(xyz_lat)*sin(xyz_lon) call AssertEqual(& message='xyz_KGrad_wt with r^2 Y_1^-1', & answer = xyz_DData, & check = xyz_KGrad_wt(wt_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1 xyz_DData = - 3*(n-1)*xyz_rad**(n-1)* cos(xyz_lat)*sin(xyz_lat)*sin(xyz_lon) call AssertEqual(& message='wt_QOperator_wt r^2 Y_1^-1', & answer = xyz_DData, & check = xyz_wt(wt_QOperator_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- xyz_Psi = cos(xyz_lat)*sin(xyz_lat) * sin(xyz_lon) ! Y_2^1 xyz_DData = cos(xyz_lat)*sin(xyz_lat) * cos(xyz_lon) ! Y_2^1 call AssertEqual(& message='wt_KxRGrad_wt with Y_2^1', & answer = xyz_DData, & check = xyz_wt(wt_KxRGrad_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! k・▽ Y_2^1 = (- 4 Y_3^1 / 15 - Y_1^1 /5) xyz_DData = cos(2*xyz_lat)*cos(xyz_lat)*sin(xyz_lon)/xyz_rad call AssertEqual(& message='xyz_KGrad_wt with Y_2^1', & answer = xyz_DData, & check = xyz_KGrad_wt(wt_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = (16*sin(xyz_lat)**2 - 5)*cos(xyz_lat)*sin(xyz_lon)/xyz_rad call AssertEqual(& message='wt_QOperator_wt Y_2^1', & answer = xyz_DData, & check = xyz_wt(wt_QOperator_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== wt_RadRot_xyz_xyz, wt_RadRotRot_xyz_xyz_xyz ========== !---------- 単純流れ v_r=r ---------- xyz_VLon = 0 ; xyz_VLat = 0 ; xyz_VRad = xyz_Rad xyz_RadRot = 0 ; xyz_RadRotRot = 0 call AssertEqual(& message='wt_RadRot_xyz_xyz with v_r=r', & answer = xyz_RadRot, & check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='wt_RadRotRot_xyz_xyz with v_r=r', & answer = xyz_RadRotRot, & check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 剛体回転解非線形項(東西流れ) ---------- xyz_VLon = 0 xyz_VLat = xyz_Rad*sin(xyz_Lat)*cos(xyz_Lat) xyz_VRad = -xyz_Rad*cos(xyz_Lat)**2 xyz_RadRot = 0 ; xyz_RadRotRot = 0 call AssertEqual(& message='wt_RadRot_xyz_xyz with rigid rotation(E-W)', & answer = xyz_RadRot, & check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='wt_RadRotRot_xyz_xyz with rigid rotation(E-W)', & answer = xyz_RadRotRot, & check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 剛体回転解非線形項(南北流れ) ---------- xyz_VLon = xyz_Rad*cos(xyz_Lat)*sin(xyz_Lon)*cos(xyz_Lon) xyz_VLat = -xyz_Rad*sin(xyz_Lat)*cos(xyz_Lat)*sin(xyz_Lon)**2 xyz_VRad = -xyz_Rad*(sin(xyz_Lat)**2*sin(xyz_Lon)**2 + cos(xyz_Lon)**2) xyz_RadRot = 0 ; xyz_RadRotRot = 0 call AssertEqual(& message='wt_RadRot_xyz_xyz with rigid rotation(N-S)', & answer = xyz_RadRot, & check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='wt_RadRotRot_xyz_xyz with rigid rotation(N-S)', & answer = xyz_RadRotRot, & check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !--------- 鉛直渦度を伴うベクトル場 ---------- xyz_Psi = xyz_Rad**2 * cos(xyz_Lat)*sin(xyz_Lon) ! r**2 P_1^1 !xyz_Psi = xyz_Rad**2 * cos(xyz_Lat)*sin(xyz_Lat)*sin(xyz_Lon) ! r**2 P_2^1 xyz_VLon = xyz_gRadLat_wt(wt_xyz(xyz_Psi*xyz_Rad)) xyz_VLat = - xyz_gRadLon_wt(wt_xyz(xyz_Psi*xyz_Rad)) xyz_VRad = 0 xyz_RadRot = 2 * xyz_Psi ! r・▽×▽×(ψr) = L_2ψ !xyz_RadRot = 6 * xyz_Psi ! r・▽×▽×(ψr) = L_2ψ xyz_RadRotRot = 0 call AssertEqual(& message='wt_RadRot_xyz_xyz with r^2 Y_1^-1 toroidal field', & answer = xyz_RadRot, & check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='wt_RadRotRot_xyz_xyz with r^2 Y_1^-1 toroidal field',& answer = xyz_RadRotRot, & check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 鉛直速度を伴うベクトル場 ----------- xyz_VRad = xyz_wt(wt_l2_wt(wt_xyz(xyz_Psi/xyz_Rad))) xyz_VLat = xyz_gRadLat_wt(wt_dRad_wt(wt_xyz(xyz_Psi*xyz_Rad))) xyz_VLon = xyz_gRadLon_wt(wt_dRad_wt(wt_xyz(xyz_Psi*xyz_Rad))) xyz_RadRot = 0 xyz_RadRotRot = -xyz_wt(wt_l2_wt(wt_lapla_wt(wt_xyz(xyz_Psi)))) ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ call AssertEqual(& message='wt_RadRot_xyz_xyz with non-vortical field', & answer = xyz_RadRot, & check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='wt_RadRotRot_xyz_xyz with r^2 Y_1^-1 toroidal field',& answer = xyz_RadRotRot, & check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) !========== wt_Potential2Vector, wt_Potential2Rotation ========== !----------------- 剛体回転場 -------------------- ! wt_Psi = wt_xyz(xyz_Rad * sin(xyz_Lat)) wt_Phi = 0.0D0 call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = xyz_Rad * cos(xyz_Lat) call AssertEqual(& message='wt_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='wt_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='wt_Potential2Vector (Rad) with rigid rotation field(E-W)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call wt_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 0.0D0 call AssertEqual(& message='wt_Potential2Rotation (Lon) with rigid rotation field(E-W)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2*cos(xyz_LAT) call AssertEqual(& message='wt_Potential2Rotation (Lat) with rigid rotation field(E-W)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2*sin(xyz_LAT) call AssertEqual(& message='wt_Potential2Rotation (Rad) with rigid rotation field(E-W)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) !----------------- 剛体回転場(南北流) -------------------- ! wt_Psi = wt_xyz(xyz_Rad * cos(xyz_Lat) * sin(xyz_Lon)) wt_Phi = 0.0D0 call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = -xyz_Rad*sin(xyz_Lat)*sin(xyz_Lon) call AssertEqual(& message='wt_Potential2Vector (Lon) with rigid rotation field(N-S)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = -xyz_Rad*cos(xyz_Lon) call AssertEqual(& message='wt_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='wt_Potential2Vector (Rad) with rigid rotation field(N-S)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call wt_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 2*cos(xyz_Lon) call AssertEqual(& message='wt_Potential2Rotation (Lon) with rigid rotation field(N-S)', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = -2*sin(xyz_Lon)*sin(xyz_Lat) call AssertEqual(& message='wt_Potential2Rotation (Lat) with rigid rotation field(N-S)', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2*sin(xyz_Lon)*cos(xyz_Lat) call AssertEqual(& message='wt_Potential2Rotation (Rad) with rigid rotation field(N-S)', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 渦無し場 -------------------- wt_Psi = 0.0D0 wt_Phi = wt_xyz(xyz_Rad * sin(xyz_Lat)) call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 0 call AssertEqual(& message='wt_Potential2Vector (Lon) with non-vortical field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2 * cos(xyz_Lat) call AssertEqual(& message='wt_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) call AssertEqual(& message='wt_Potential2Vector (Rad) with non-vortical field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call wt_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 0.0D0 call AssertEqual(& message='wt_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='wt_Potential2Rotation (Lat) with non-vortical field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wt_Potential2Rotation (Rad) with non-vortical field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wt_Psi = 0.0D0 wt_Phi = wt_xyz(xyz_Rad**2 * cos(xyz_Lat)*sin(xyz_Lon)) call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 3 * xyz_Rad * cos(xyz_Lon) call AssertEqual(& message='wt_Potential2Vector (Lon) with poloidal field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = - 3 * xyz_Rad * sin(xyz_Lat) * sin(xyz_Lon) call AssertEqual(& message='wt_Potential2Vector (Lat) with poloidal field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 2 * xyz_Rad * cos(xyz_Lat) * sin(xyz_Lon) call AssertEqual(& message='wt_Potential2Vector (Rad) with poloidal field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call wt_Potential2Rotation(& xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi ) xyz_DData = 4*sin(xyz_Lat)*sin(xyz_Lon) call AssertEqual(& message='wt_Potential2Rotation (Lon) with poloidal field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 4*cos(xyz_Lon) call AssertEqual(& message='wt_Potential2Rotation (Lat) with poloidal field', & answer = xyz_VLat, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 0.0D0 call AssertEqual(& message='wt_Potential2Rotation (Rad) with poloidal field', & answer = xyz_VRad, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','wt_module_svpack_derivative_test2', & 'wt_module_svpack derivative function test #2 succeeded!') end program wt_module_svpack_derivative2