!---------------------------------------------------------------------- ! Copyright (c) 2024 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wt_zonal_module テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2024/03/19 竹広真一 wt_wave_module_derivative_test2.f90 より改造 ! program wt_zonal_module_derivative_test2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wt_zonal_module implicit none integer,parameter :: im=1, 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,0:lm) :: wt_Psi real(8), dimension(nm+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_zonal_module_derivative_test2', & 'wt_zonal_module 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 * sin(xyz_Lat) ! r**2 P_1^0 xyz_DData = 0.0D0 call AssertEqual(& message='wt_KxRGrad_wt with r^2 P_1^0', & 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^0 xyz_DData = xyz_rad**(n-1)*( cos(xyz_Lat)**2 + n*sin(xyz_Lat)**2) call AssertEqual(& message='xyz_KGrad_wt with r^2 P_1^0', & answer = xyz_DData, & check = xyz_KGrad_wt(wt_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) ! Q r**n Y_1^0 = -k ・▽ r**n Y_1^0 xyz_DData = - 3*xyz_rad**(n-1)*(sin(xyz_Lat)**2-1.0D0/3.0D0) call AssertEqual(& message='wt_QOperator_wt r^2 P_1^0 ', & answer = xyz_DData, & check = xyz_wt(wt_QOperator_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 2 -------------------- xyz_Psi = (3.0D0*sin(xyz_Lat)**2-1.0D0)/2.0D0 ! P_2^0 xyz_DData = 0.0D0 call AssertEqual(& message='wt_KxRGrad_wt with P_2^0', & answer = xyz_DData, & check = xyz_wt(wt_KxRGrad_wt(wt_xyz(xyz_Psi))), & significant_digits = check_digits, ignore_digits = ignore & ) ! k・▽ Y_2^0 xyz_DData = 3*sin(xyz_Lat)*cos(xyz_Lat)**2/xyz_rad call AssertEqual(& message='xyz_KGrad_wt with P_2^0', & answer = xyz_DData, & check = xyz_KGrad_wt(wt_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_DData = 6.0D0/(5.0D0*xyz_Rad)*( & 8*(5*sin(xyz_Lat)**3-3*sin(xyz_Lat))/2.0D0- 3*sin(xyz_Lat) ) call AssertEqual(& message='wt_QOperator_wt P_2^0', & 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 * sin(xyz_Lat) ! r**2 P_1^0 !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^0 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^0 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 * sin(xyz_Lat)) call wt_Potential2Vector(& xyz_VLon,xyz_VLat,xyz_VRad,wt_Psi, wt_Phi ) xyz_DData = 0.0D0 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 * cos(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 * sin(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*cos(xyz_Lat) call AssertEqual(& message='wt_Potential2Rotation (Lon) with poloidal 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 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_derivative_test2', & 'wt_module derivative function test #2 succeeded!') end program wt_zonal_module_derivative_test2