!---------------------------------------------------------------------- ! Copyright (c) 2018--2020 Shin-ichi Takehiro. All rights reserved. !---------------------------------------------------------------------- ! !表題 ut_mpi_module テストプログラム ! ! トロイダルポロイダル関係微分関数のテスト ! !履歴 2018/08/09 竹広真一 ! 2020/11/11 竹広真一 セクター計算オプション導入 ! 2021/11/15 竹広真一 ut_RadRot_RadRotRot テスト追加 ! program ut_mpi_module_derivative_mint_test2 use dc_message, only : MessageNotify use dc_test, only : AssertEqual use ut_mpi_module_mint use mpi implicit none !!$ integer,parameter :: im=32, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) !!$ integer,parameter :: mint=1 ! 経度方向対称性 integer,parameter :: im=16, jm=16, km=16 ! 格子点の設定(経度, 緯度, 動径) integer,parameter :: mint=2 ! 経度方向対称性 integer,parameter :: nm=10, lm=16 ! 切断波数の設定(水平, 動径) !!$ integer,parameter :: im=16, jm=16, km=32 ! 格子点の設定(経度, 緯度, 動径) !!$ integer,parameter :: mint=2 ! 経度方向対称性 !!$ integer,parameter :: nm=10, lm=32 ! 切断波数の設定(水平, 動径) real(8),parameter :: ri=0.5, ro=1.5 ! 内外半径 integer :: npv=2 ! 判定誤差設定 integer, parameter :: check_digits = 7 integer, parameter :: ignore = -8 real(8), allocatable :: ut_Psi(:,:) real(8), allocatable :: ut_Phi(:,:) real(8), allocatable :: xvh_DData(:,:,:) real(8), allocatable :: xvh_Psi(:,:,:) real(8), allocatable :: xvh_VLon(:,:,:) real(8), allocatable :: xvh_VLat(:,:,:) real(8), allocatable :: xvh_VRad(:,:,:) real(8), allocatable :: xvh_RadRot(:,:,:) real(8), allocatable :: xvh_RadRotRot(:,:,:) integer, parameter :: n=2 integer :: iproc, np, ierr !---------------- MPI スタート --------------------- !call get_ENV_INT('NPV',npv) call MPI_INIT(IERR) call MPI_COMM_RANK(MPI_COMM_WORLD,IPROC,IERR) call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR) call MessageNotify('M','ut_mpi_module_derivative_mint_test2', & 'ut_mpi_module_mint derivative function test #2') call ut_mpi_initial(im,jm,km,nm,lm,ri,ro,npv=npv,mint=mint) allocate(ut_Psi(nmc,0:lm)) allocate(ut_Phi(nmc,0:lm)) allocate(xvh_DData(0:im-1,1:jc,kc)) allocate(xvh_Psi(0:im-1,1:jc,kc)) allocate(xvh_VLon(0:im-1,1:jc,kc)) allocate(xvh_VLat(0:im-1,1:jc,kc)) allocate(xvh_VRad(0:im-1,1:jc,kc)) allocate(xvh_RadRot(0:im-1,1:jc,kc)) allocate(xvh_RadRotRot(0:im-1,1:jc,kc)) !========== ut_RadRot_xvh_xvh, ut_RadRotRot_xvh_xvh_xvh ========== !---------- 単純流れ v_r=r ---------- xvh_VLon = 0 ; xvh_VLat = 0 ; xvh_VRad = xvh_Rad xvh_RadRot = 0 ; xvh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xvh_xvh with v_r=r', & answer = xvh_RadRot, & check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xvh_xvh with v_r=r', & answer = xvh_RadRotRot, & check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call ut_RadRot_RadRotRot( & xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRot) with v_r=r', & answer = xvh_RadRot, & check = xvh_ut(ut_Psi), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRotRot) with v_r=r', & answer = xvh_RadRotRot, & check = xvh_ut(ut_Phi), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 剛体回転解非線形項(東西流れ) ---------- xvh_VLon = 0 xvh_VLat = xvh_Rad*sin(xvh_Lat)*cos(xvh_Lat) xvh_VRad = -xvh_Rad*cos(xvh_Lat)**2 xvh_RadRot = 0 ; xvh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xvh_xvh with rigid rotation(E-W)', & answer = xvh_RadRot, & check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xvh_xvh with rigid rotation(E-W)', & answer = xvh_RadRotRot, & check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call ut_RadRot_RadRotRot( & xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRot) with rigid rotation(E-W)', & answer = xvh_RadRot, & check = xvh_ut(ut_Psi), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRotRot) with rigid rotation(E-W)', & answer = xvh_RadRotRot, & check = xvh_ut(ut_Phi), & significant_digits = check_digits, ignore_digits = ignore & ) !!$ !---------- 剛体回転解非線形項(南北流れ) ---------- !!$ xvh_VLon = xvh_Rad*cos(xvh_Lat)*sin(xvh_Lon)*cos(xvh_Lon) !!$ xvh_VLat = -xvh_Rad*sin(xvh_Lat)*cos(xvh_Lat)*sin(xvh_Lon)**2 !!$ xvh_VRad = -xvh_Rad*(sin(xvh_Lat)**2*sin(xvh_Lon)**2 + cos(xvh_Lon)**2) !!$ xvh_RadRot = 0 ; xvh_RadRotRot = 0 !!$ !!$ call AssertEqual(& !!$ message='ut_RadRot_xvh_xvh with rigid rotation(N-S)', & !!$ answer = xvh_RadRot, & !!$ check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ call AssertEqual(& !!$ message='ut_RadRotRot_xvh_xvh with rigid rotation(N-S)', & !!$ answer = xvh_RadRotRot, & !!$ check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !--------- 鉛直渦度を伴うベクトル場 ---------- xvh_Psi = xvh_Rad**3 * cos(xvh_Lat)**2*sin(2*xvh_Lon) ! r**2 P_2^-2 xvh_VLon = xvh_gRadLat_ut(ut_xvh(xvh_Psi*xvh_Rad)) xvh_VLat = - xvh_gRadLon_ut(ut_xvh(xvh_Psi*xvh_Rad)) xvh_VRad = 0 xvh_RadRot = 6 * xvh_Psi ! r・▽×▽×(ψr) = L_2ψ xvh_RadRotRot = 0 call AssertEqual(& message='ut_RadRot_xvh_xvh with r^2 Y_2^-2 toroidal field', & answer = xvh_RadRot, & check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xvh_xvh with r^2 Y_2^-2 toroidal field',& answer = xvh_RadRotRot, & check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call ut_RadRot_RadRotRot( & xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRot) with r^2 Y_2^-2 toroidal field',& answer = xvh_RadRot, & check = xvh_ut(ut_Psi), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRotRot) with r^2 Y_2^-2 toroidal field',& answer = xvh_RadRotRot, & check = xvh_ut(ut_Phi), & significant_digits = check_digits, ignore_digits = ignore & ) !---------- 鉛直速度を伴うベクトル場 ----------- ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ xvh_VRad = xvh_ut(ut_l2_ut(ut_xvh(xvh_Psi/xvh_Rad))) xvh_VLat = xvh_GradLat_ut(ut_DRad_ut(ut_xvh(xvh_Psi*xvh_Rad))) xvh_VLon = xvh_GradLon_ut(ut_DRad_ut(ut_xvh(xvh_Psi*xvh_Rad))) xvh_RadRot = 0 xvh_RadRotRot = -xvh_ut(ut_l2_ut(ut_lapla_ut(ut_xvh(xvh_Psi)))) ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ call AssertEqual(& message='ut_RadRot_xvh_xvh with non-vortical field', & answer = xvh_RadRot, & check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRotRot_xvh_xvh with r^2 Y_2^-2 Poloidal field',& answer = xvh_RadRotRot, & check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), & significant_digits = check_digits, ignore_digits = ignore & ) call ut_RadRot_RadRotRot( & xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRot) with r^2 Y_2^-2 Poloidal field',& answer = xvh_RadRot, & check = xvh_ut(ut_Psi), & significant_digits = check_digits, ignore_digits = ignore & ) call AssertEqual(& message='ut_RadRot_RadRotRot_xvh (RadRotRot) with r^2 Y_2^-2 Poloidal field',& answer = xvh_RadRotRot, & check = xvh_ut(ut_Phi), & significant_digits = check_digits, ignore_digits = ignore & ) !========== ut_Potential2Vector, ut_Potential2Rotation ========== !----------------- 剛体回転場 -------------------- ! ut_Psi = ut_xvh(xvh_Rad * sin(xvh_Lat)) ut_Phi = 0.0D0 call ut_Potential2VectorMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = xvh_Rad * cos(xvh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', & answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with rigid rotation field(E-W)', & answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with rigid rotation field(E-W)', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with rigid rotation field(E-W)', & answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 2*cos(xvh_LAT) call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with rigid rotation field(E-W)', & answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 2*sin(xvh_LAT) call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with rigid rotation field(E-W)', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) !!$ !----------------- 剛体回転場(南北流) -------------------- !!$ ! !!$ ut_Psi = ut_xvh(xvh_Rad * cos(xvh_Lat) * sin(xvh_Lon)) !!$ ut_Phi = 0.0D0 !!$ !!$ call ut_Potential2VectorMPI(& !!$ xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) !!$ !!$ xvh_DData = -xvh_Rad*sin(xvh_Lat)*sin(xvh_Lon) !!$ call AssertEqual(& !!$ message='ut_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', & !!$ answer = xvh_VLon, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xvh_DData = -xvh_Rad*cos(xvh_Lon) !!$ call AssertEqual(& !!$ message='ut_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', & !!$ answer = xvh_VLat, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xvh_DData = 0.0D0 !!$ call AssertEqual(& !!$ message='ut_Potential2VectorMPI (Rad) with rigid rotation field(N-S)', & !!$ answer = xvh_VRad, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ call ut_Potential2RotationMPI(& !!$ xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) !!$ !!$ xvh_DData = 2*cos(xvh_Lon) !!$ call AssertEqual(& !!$ message='ut_Potential2RotationMPI (Lon) with rigid rotation field(N-S)', & !!$ answer = xvh_VLon, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xvh_DData = -2*sin(xvh_Lon)*sin(xvh_Lat) !!$ call AssertEqual(& !!$ message='ut_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', & !!$ answer = xvh_VLat, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) !!$ !!$ xvh_DData = 2*sin(xvh_Lon)*cos(xvh_Lat) !!$ call AssertEqual(& !!$ message='ut_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', & !!$ answer = xvh_VRad, & !!$ check = xvh_DData, & !!$ significant_digits = check_digits, ignore_digits = ignore & !!$ ) ! ----------------- 渦無し場 -------------------- ut_Psi = 0.0D0 ut_Phi = ut_xvh(xvh_Rad * sin(xvh_Lat)) call ut_Potential2VectorMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = 0 call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with non-vortical field',& answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 2 * cos(xvh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with non-vortical field',& answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 2 * sin(xvh_Lat) call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with non-vortical field', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with non-vortical field', & answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with non-vortical field', & answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with non-vortical field', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) ! ----------------- 例 4 -------------------- ! ポロイダル速度場 ut_Psi = 0.0D0 ut_Phi = ut_xvh(xvh_Rad**3 * cos(xvh_Lat)**2*cos(2*xvh_Lon)) call ut_Potential2VectorMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = 4*xvh_Rad**2 * cos(xvh_Lat) *(-2)*sin(2*xvh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lon) with poloidal field', & answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 4*xvh_Rad**2 * (-2)*cos(xvh_Lat)*sin(xvh_Lat)*cos(2*xvh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Lat) with poloidal field', & answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 6 * xvh_Rad**2 * cos(xvh_Lat)**2 * cos(2*xvh_Lon) call AssertEqual(& message='ut_Potential2VectorMPI (Rad) with poloidal field', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call ut_Potential2RotationMPI(& xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi ) xvh_DData = 6* xvh_Rad * 2*cos(xvh_Lat)*sin(xvh_Lat)*cos(2*xvh_Lon) call AssertEqual(& message='ut_Potential2RotationMPI (Lon) with poloidal field', & answer = xvh_VLon, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 6* xvh_Rad * cos(xvh_Lat)*(-2)*sin(2*xvh_Lon) call AssertEqual(& message='ut_Potential2RotationMPI (Lat) with poloidal field', & answer = xvh_VLat, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) xvh_DData = 0.0D0 call AssertEqual(& message='ut_Potential2RotationMPI (Rad) with poloidal field', & answer = xvh_VRad, & check = xvh_DData, & significant_digits = check_digits, ignore_digits = ignore & ) call MessageNotify('M','ut_mpi_module_derivative_mint_test2', & 'ut_mpi_module_mint derivative function test #2 succeeded!') !------ MPIの終了 ------ call MPI_FINALIZE(IERR) end program ut_mpi_module_derivative_mint_test2