!---------------------------------------------------------------------- ! Copyright (c) 2013 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wt_wave_module_sjpack テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2013/07/03 竹広真一 wt_module_sjpack_derivative_test2.f90 より改造 ! program wt_wave_module_sjpack_derivative_test2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wt_wave_module_sjpack implicit none integer,parameter :: im=2, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: nm=10, lm=16 ! 切断波数の設定(水平, 動径) integer,parameter :: m=1 ! 東西波数の設定 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(2*(nm-m+1),0:lm) :: wt_Psi real(8), dimension(2*(nm-m+1),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 integer :: j call MessageNotify('M','wt_wave_module_sjpack_derivative_test2', & 'wt_wave_module_sjpack derivative function test #2') call wt_initial(im,jm,km,nm,m,lm,ri,ro) !========== wt_KxRGrad_wt, xyz_KGrad_wt, wt_QOperator_wt ========== ! ----------------- 例 1 -------------------- xyz_Psi = xyz_rad**n * cos(xyz_Lat) ! r**2 P_1^1 do j=1,jm xyz_DData(0,j,:) = -m*z_Rad**n * cos(y_Lat(j)) xyz_DData(1,j,:) = m*z_Rad**n * cos(y_Lat(j)) enddo call AssertEqual(& message='wt_KxRGrad_wt with r^2 P_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) call AssertEqual(& message='xyz_KGrad_wt with r^2 P_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) call AssertEqual(& message='wt_QOperator_wt r^2 P_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) ! P_2^1 do j=1,jm xyz_DData(0,j,:) = -m*cos(y_Lat(j))*sin(y_Lat(j)) ! P_2^1 xyz_DData(1,j,:) = m*cos(y_Lat(j))*sin(y_Lat(j)) ! P_2^1 enddo call AssertEqual(& message='wt_KxRGrad_wt with P_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)/xyz_rad call AssertEqual(& message='xyz_KGrad_wt with P_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)/xyz_rad call AssertEqual(& message='wt_QOperator_wt P_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 ========== !--------- 鉛直渦度を伴うベクトル場 ---------- xyz_Psi = xyz_Rad**2 * cos(xyz_Lat) ! r**2 P_1^1 !xyz_Psi = xyz_Rad**2 * cos(xyz_Lat)*sin(xyz_Lat) ! 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 P_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 P_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_xyz with non-vortical 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 ========== ! ----------------- 例 4 -------------------- ! ポロイダル速度場 wt_Psi = 0.0D0 wt_Phi = wt_xyz(xyz_Rad**2 * cos(xyz_Lat)) call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad,wt_Psi, wt_Phi ) do j=1,jm xyz_DData(0,j,:) = -m * 3 * z_Rad xyz_DData(1,j,:) = m * 3 * z_Rad enddo 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) 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) 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) call AssertEqual(& message='wt_Potential2Rotation (Lon) with poloidal field', & answer = xyz_VLon, & check = xyz_DData, & significant_digits = check_digits, ignore_digits = ignore & ) do j=1,jm xyz_DData(0,j,:) = -m*4 xyz_DData(1,j,:) = m*4 enddo 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_sjpack_derivative_test2', & 'wt_module_sjpack derivative function test #2 succeeded!') end program wt_wave_module_sjpack_derivative_test2