!---------------------------------------------------------------------- ! Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 wsc_module_svpack テストプログラム ! ! 微分のテスト ! !履歴 2017/05/21 竹広真一 ! program wsc_module_svpack_derivative1 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use wsc_module_svpack use wa_module_svpack, only : wa_Lapla_wa implicit none integer,parameter :: im=32, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: nm=10, lm=10 ! 切断波数の設定(水平, 動径) real(8),parameter :: ri=0.5, ro=1.5 ! 内外半径 real(8), dimension(0:im-1,1:jm,0:km) :: xyz_Data real(8), dimension(0:im-1,1:jm,0:km) :: xyz_DData real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VLon real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VLat real(8), dimension(0:im-1,1:jm,0:km) :: xyz_VRad real(8), dimension(0:im-1,1:jm,0:km) :: xyz_Div real(8), dimension(0:im-1,1:jm,0:km) :: xyz_Psi real(8), dimension(0:im-1,1:jm,0:km) :: xyz_xi ! 判定誤差設定 integer, parameter :: check_digits = 10 integer, parameter :: ignore = -11 integer, parameter :: nmin=1, nmax=9 integer :: n real(8) :: pi pi = atan(1.0D0)*4.0D0 call MessageNotify('M','wsc_module_svpack_derivative_test1', & 'wsc_module_svpack derivative function test #1') call wsc_Initial(im,jm,km,nm,lm,ri,ro) xyz_xi = (xyz_Rad - ri)/(ro-ri) !==================== wc_DRad_ws, ws_DRad_wc ==================== xyz_Data = sin(PI*xyz_xi) xyz_DData = PI*cos(PI*xyz_xi) ! call AssertEqual(& message='wc_DRad_ws with sin(z)', & answer = xyz_DData, & check = xyz_wc(wc_DRad_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(2*PI*xyz_xi) xyz_DData = -2*PI*sin(2*PI*xyz_xi) ! call AssertEqual(& message='ws_DRad_wc with cos(2z)', & answer = xyz_DData, & check = xyz_ws(ws_DRad_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== wc_DivRad_ws, ws_DivRad_wc ==================== xyz_Data = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * cos(PI*xyz_xi) ! Y_2^-1 cos(z) xyz_DData = -PI* cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) ! call AssertEqual(& message='ws_DivRad_wc with Y_2^-1 cos(z)', & answer = xyz_DData, & check = xyz_ws(ws_DivRad_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(xyz_Lat)*sin(xyz_Lat) * cos(xyz_Lon) * sin(2*PI*xyz_xi) ! Y_2^1 sin(2z) xyz_DData =2*PI*cos(xyz_Lat)*sin(xyz_Lat) * cos(xyz_Lon) * cos(2*PI*xyz_xi) ! call AssertEqual(& message='wc_DivRad_ws with Y_2^1 sin(2z)', & answer = xyz_DData, & check = xyz_wc(wc_DivRad_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== wc_RotRad_ws, ws_RotRad_wc ==================== xyz_Data = sin(xyz_Lat) * cos(3*PI*xyz_xi) ! Y_2^0 cos(3z) xyz_DData = -3*PI* sin(xyz_Lat) * sin(3*PI*xyz_xi) ! call AssertEqual(& message='ws_RotRad_wc with Y_2^0 cos(3z)', & answer = xyz_DData, & check = xyz_ws(ws_RotRad_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = sin(xyz_Lat) * sin(4*PI*xyz_xi) ! Y_2^0 sin(4z) xyz_DData = 4*PI* sin(xyz_Lat) * cos(4*PI*xyz_xi) ! call AssertEqual(& message='wc_RotRad_ws with Y_2^0 sin(4z)', & answer = xyz_DData, & check = xyz_wc(wc_RotRad_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== ws_Lapla_ws, wc_Lapla_wc ==================== xyz_Data = cos(xyz_Lat)*sin(xyz_Lon)*sin(3*PI*xyz_xi) ! Y_1^-1 sin(3z) xyz_DData = -( 9*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='ws_Lapla_ws with Y_1^-1 sin(3z)', & answer = xyz_DData, & check = xyz_ws(ws_Lapla_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(xyz_Lat)*cos(xyz_Lon)*cos(4*PI*xyz_xi) ! Y_1^1 cos(4z) xyz_DData = -( 16*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='wc_Lapla_wc with Y_1^1 cos(4z)', & answer = xyz_DData, & check = xyz_wc(wc_Lapla_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !==================== ws_LaplaInv_ws, wc_LaplaInv_wc ==================== xyz_Data = cos(xyz_Lat)*sin(xyz_Lon)*sin(3*PI*xyz_xi) ! Y_1^-1 sin(3z) xyz_DData = -1.0D0/( 9*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='ws_LaplaInv_ws with Y_1^-1 sin(3z)', & answer = xyz_DData, & check = xyz_ws(ws_LaplaInv_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(xyz_Lat)*sin(xyz_Lon)*sin(2*PI*xyz_xi) ! Y_1^-1 sin(2z) xyz_DData = -1.0D0/( 4*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='ws_LaplaInv_ws with Y_1^-1 sin(2z)', & answer = xyz_DData, & check = xyz_ws(ws_LaplaInv_ws(ws_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(xyz_Lat)*cos(xyz_Lon)*cos(4*PI*xyz_xi) ! Y_1^1 cos(4z) xyz_DData = -1.0D0/( 16*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='wc_LaplaInv_wc with Y_1^1 cos(4z)', & answer = xyz_DData, & check = xyz_wc(wc_LaplaInv_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Data = cos(xyz_Lat)*cos(xyz_Lon)*cos(5*PI*xyz_xi) ! Y_1^1 cos(4z) xyz_DData = -1.0D0/( 25*PI**2 + 1*2/Ro**2) * xyz_Data ! call AssertEqual(& message='wc_LaplaInv_wc with Y_1^1 cos(5z)', & answer = xyz_DData, & check = xyz_wc(wc_LaplaInv_wc(wc_xyz(xyz_Data))), & significant_digits = check_digits, ignore_digits = ignore & ) !========== xyz_GradLon_ws, xyz_GradLat_ws ========== ! ------ Y_1^-1 sin(z) Toridal field ----- xyz_Psi = cos(xyz_Lat)*sin(xyz_Lon) * sin(PI*xyz_xi) ! Y_1^-1 sin(z) xyz_VRad = 0 xyz_VLon = cos(xyz_Lon) * sin(PI*xyz_xi) /Ro call AssertEqual(& message='xyz_GradLon_ws with Y_1^-1 sin(z)', & answer = xyz_VLon, & check = xyz_GradLon_ws(ws_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_VLat = -sin(xyz_Lat)*sin(xyz_Lon) * sin(PI*xyz_xi) /Ro call AssertEqual(& message='xyz_GradLat_ws with Y_1^-1 sin(z)', & answer = xyz_VLat, & check = xyz_GradLat_ws(ws_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) !========== xyz_GradLon_wc, xyz_GradLat_wc, wc_Div_xyz_xyz_xyz, wc_RotRad_xyz_xyz ========== ! ------ Y_1^1 cos(2z) Toridal field ----- xyz_Psi = cos(xyz_Lat)*cos(xyz_Lon) * cos(2*PI*xyz_xi) ! Y_1^1 cos(2z) xyz_VRad = 0 xyz_VLon = - sin(xyz_Lon) * cos(2*PI*xyz_xi) /Ro call AssertEqual(& message='xyz_GradLon_wc with Y_1^1 cos(2z)', & answer = xyz_VLon, & check = xyz_GradLon_wc(wc_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_VLat = -sin(xyz_Lat)*cos(xyz_Lon) * cos(2*PI*xyz_xi) /Ro call AssertEqual(& message='xyz_GradLat_wc with Y_1^1 cos(2z)', & answer = xyz_VLat, & check = xyz_GradLat_wc(wc_xyz(xyz_Psi)), & significant_digits = check_digits, ignore_digits = ignore & ) xyz_Div = - 2 * xyz_Psi/Ro**2 call AssertEqual(& message='wc_Div_xyz_xyz_xyz with Toroiral Y_1^-1 cos(2z)', & answer = xyz_Div, & check = xyz_wc(wc_Div_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)),& significant_digits = check_digits, ignore_digits = ignore & ) xyz_VLat = - xyz_Vlat call AssertEqual(& message='wc_RotRad_xyz_xyz_xyz with Toroiral Y_1^-1 cos(2z)', & answer = xyz_Div, & check = xyz_wc(wc_RotRad_xyz_xyz(xyz_VLon,xyz_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) ! ------ Y_2^-1 sin(z) Poloidal field ----- xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) ! Y_2^-1 sin(z) xyz_VLon = xyz_GradLon_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z) xyz_VLat = xyz_GradLat_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z) xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro ! (2*3/Ro) Psi xyz_DData = 0.0D0 call AssertEqual(& message='wc_Div_xyz_xyz_xyz with Y_2^-1 sin(z)', & answer = xyz_DData, & check = xyz_wc(wc_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), & significant_digits = check_digits, ignore_digits = ignore & ) ! ------ Y_2^-1 cos(z) Poloidal field ----- xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * cos(PI*xyz_xi) ! Y_2^-1 sin(z) xyz_VLon = xyz_GradLon_ws(ws_DRad_wc(wc_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z) xyz_VLat = xyz_GradLat_ws(ws_DRad_wc(wc_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z) xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro ! (2*3/Ro) Psi xyz_DData = 0.0D0 call AssertEqual(& message='ws_Div_xyz_xyz_xyz with Y_2^-1 sin(z)', & answer = xyz_DData, & check = xyz_ws(ws_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), & significant_digits = check_digits, ignore_digits = ignore & ) !!$ !========== xyz_RotLon_ws_wc, xyz_RotLat_wc_ws ============ !!$ !!$ ! ------ Y_2^-1 sin(z) Poloidal field ----- !!$ !!$ xyz_Psi = cos(xyz_Lat)*sin(xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) ! Y_2^-1 sin(z) !!$ !!$ xyz_VLon = xyz_GradLon_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*sin(lat)cos(lon))cos(z) !!$ xyz_VLat = xyz_GradLat_wc(wc_DRad_ws(ws_xyz(xyz_Psi)))*Ro ! PI*cos(2lat)sin(lon)cos(z) !!$ xyz_Vrad = -xyz_wz(wa_Lapla_wa(wz_xyz(xyz_Psi)))/Ro ! (2*3/Ro) Psi !!$ !!$ xyz_DData = 0.0D0 !!$ call AssertEqual(& !!$ message='wc_Div_xyz_xyz_xyz with Y_2^-1 sin(z)', & !!$ answer = xyz_DData, & !!$ check = xyz_wc(wc_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xyz_DData = (2*3/Ro**2 + PI**2) * cos(2*xyz_Lat) * sin(xyz_Lon) * sin(PI*xyz_xi) !!$ call AssertEqual(& !!$ message='xyz_RadLon_wt with Y_2^-1 sin(z)', & !!$ answer = xyz_DData, & !!$ check = xyz_RotLon_ws_wc(ws_xyz(xyz_Vrad),wc_xyz(xyz_Vlat)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xyz_DData = -(2*3/Ro**2 + PI**2) * sin(xyz_Lat)*cos(xyz_Lon)*sin(PI*xyz_xi) !!$ call AssertEqual(& !!$ message='xyz_RadLat_wt with Y_2^-1 sin(z)', & !!$ answer = xyz_DData, & !!$ check = xyz_RotLat_wc_ws(wc_xyz(xyz_Vlon),ws_xyz(xyz_Vrad)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xyz_VLat = xyz_Rad**(n-1)*cos(2*xyz_Lat)*sin(xyz_Lon) !!$ call AssertEqual(& !!$ message='xyz_GradLat_wt with r^2 Y_2^-1', & !!$ answer = xyz_VLat, & !!$ check = xyz_GradLat_wt(wt_xyz(xyz_Psi)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xyz_Div = - 6* xyz_Psi/xyz_Rad**2 !!$ call AssertEqual(& !!$ message='wt_Div_xyz_xyz_xyz_wt with Toroiral r2 Y_1^-1', & !!$ answer = xyz_Div, & !!$ check = xyz_wt(wt_Div_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)),& !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ !==================== wt_DivLon_xyz, wt_DivLat_xyz ==================== !!$ n=2 !!$ write(6,*)'n =',n !!$ !!$ xyz_VLon = xyz_rad**n * cos(xyz_Lat)**2*sin(xyz_Lon) !!$ xyz_DData = xyz_rad**(n-1)*cos(xyz_Lat)*cos(xyz_Lon) !!$ call AssertEqual(& !!$ message='wt_DivLon_xyz with r^n cos^2 φ sinλ', & !!$ answer = xyz_DData, & !!$ check = xyz_wt(wt_DivLon_xyz(xyz_VLon)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xyz_VLat = xyz_rad**n * cos(xyz_Lat)**2*sin(xyz_Lon) !!$ xyz_DData = -3*xyz_rad**(n-1)*cos(xyz_Lat)*sin(xyz_Lat)*sin(xyz_Lon) !!$ call AssertEqual(& !!$ message='wt_DivLat_xyz with r^n cos^2 φ sinλ', & !!$ answer = xyz_DData, & !!$ check = xyz_wt(wt_DivLat_xyz(xyz_VLat)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ call MessageNotify('M','wsc_module_svpack_derivative_test1', & 'wsc_module_svpack derivative function test #1 succeeded!') end program wsc_module_svpack_derivative1