!-- !---------------------------------------------------------------------- ! Copyright (c) 2016 SPMODEL Development Group. All rights reserved. !---------------------------------------------------------------------- ! !表題 ee_mpi_module ! ! spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の ! 流体運動をスペクトル法により数値計算するための Fortran90 関数を ! 提供する. ! ! 内部で ISPACK/FTPACKI の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法については... ! !履歴 2016/08/08 竹広真一 eee_mpi_module より改造 ! 2016/09/22 竹広真一 ef_Lapla_ef, ef_Laplainv_ef 初期化 ! !++ module ee_mpi_module ! != ee_mpi_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== 概要 ! ! spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の ! 流体運動をスペクトル法により数値計算するための Fortran90 関数を ! 提供する. ! ! 内部で ISPACK/FTPACKI の Fortran77 サブルーチンを呼んでいる. ! スペクトルデータおよび格子点データの格納方法については... ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (ef_, vx_, x_, v_) は, 返す値の形を示している. ! ef_ :: スペクトルデータ ! (第 1,2 次元がそれぞれ Y,X 方向波数(X について分割配置)) ! vx_ :: XY 方向 2 次元格子点データ(Y について分散配置 ! x_ :: X 方向 1 次元格子点データ ! v_ :: Y 方向 1 次元格子点データ(分散配置) ! ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv)は, ! その関数の作用を表している. ! ! * 関数名の最後 (_ef_ef,_ef,_vx, _x, _v) は, 入力変数の形が ! スペクトルデータおよび格子点データであることを示している. ! _ef :: スペクトルデータ ! _ef_ef :: 2 つのスペクトルデータ ! _vx :: 2 次元格子点データ ! _x :: X 方向 1 次元格子点データ ! _v :: Y 方向 1 次元格子点データ ! !=== 各データの種類の説明 ! ! * vx : 2 次元格子点データ. ! * 変数の種類と次元は real(8), dimension(,js(ip):je(ip),0:im-1). ! * im は X 座標の格子点数であり, サブルーチン ee_mpi_initial にて ! Y 方向の全格子点数 jm とともにあらかじめ設定しておく. ! * ee_mpi_Initial によってプロセス毎に Y 方向の ! 分割配置が設定され, その情報が public 変数 js,je,jc に設定される. ! * js が分割配置の最小格子点, je が最大格子点, jc が格子点数である. ! ! * ef : スペクトルデータ. ! * 変数の種類と次元は real(8), dimension(-lm:lm,2*kc). ! * mm は Y 方向の最大波数であり, サブルーチン ee_mpi_initial にて ! X 方向の最大波数 km とともにあらかじめ設定しておく. ! (X, Y 方向波数の順ではない)ことに注意. ! * ee_mpi_Initial によってプロセス毎に X 方向波数の ! 分割配置が設定され, その情報が public 変数 ls,le,lc に設定される. ! * ks が分割配置の最小波数, ke が最大波数, kc が波数の数である. ! 実際の格納のされ方は ks, ks+1,...,ke,-ke,-ke+1,...,-ks の順である. ! * 分割配置される前のスペクトルデータの格納のされかたは IPSACK の ! P2PACK に基づく. ! ee_S( L, K) : s_{kl} の実部 ! ee_S(-L,-K) : s_{kl} の虚部 ! ee_S(-L, K) : s_{k(-l)} の実部 ! ee_S( L,-K) : s_{k(-l)} の虚部 ! ee_S( L, 0) : s_{0l} の実部 ! ee_S(-L, 0) : s_{0l} の虚部 ! ee_S( 0, K) : s_{k0} の実部 ! ee_S( 0,-K) : s_{k0} の虚部 ! ee_S( 0, 0) : s_{00} (実数) ! * 元の格子点データが実数であるため, s_{-k-l} = s^*_{kl} なる関係がある. ! ただし ^* は複素共役を表す. ! ! * x, v : X, Y 方向 1 次元格子点データ. ! * 変数の種類と次元はそれぞれ real(8), dimension(0:im-1), ! real(8), dimension(js(ip):je(ip)) ! ! * ef_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * vx_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * x_, v_ で始まる関数が返す値は 1 次元格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! ee_mpi_Initial :: スペクトル変換の格子点数, 切断波数の設定 ! kf_k :: X 方向波数の位置情報 ! !==== 座標変数 ! ! x_X, v_Y :: 格子点座標(X,Y座標)を格納した 1 次元配列 ! x_X_Weight, v_Y_Weight :: 重み座標を格納した 1 次元配列 ! vx_X, vx_Y_Z :: 格子点データの XY 座標(X,Y) ! (格子点データ型 3 次元配列) ! !==== 基本変換 ! ! vx_ef :: スペクトルデータから格子データへの変換 ! ef_vx :: 格子データからスペクトルデータへの変換 ! xy_vx :: 分散格子点データの集積 ! ee_ef :: 分散スペクトルデータの集積 ! ef_ee :: 集積スペクトルデータの分散 ! !==== 微分 ! ! ef_Lapla_ef :: スペクトルデータにラプラシアンを作用させる ! ef_LaplaInv_ef :: スペクトルデータにラプラシアンの逆変換を作用させる ! ef_Dx_ef :: スペクトルデータに X 微分を作用させる ! ef_Dy_ef :: スペクトルデータに Y 微分を作用させる ! !==== 非線形計算 ! ! ef_Joacobian_ef_ef :: 2 つのスペクトルデータからヤコビアンを計算する. ! ef_JoacobianZ_ef :: 渦度スペクトルデータから非発散方程式の非線形項を計算する ! !==== 積分・平均 ! ! IntVX_vx, AvrVX_vx :: 2 次元(XY)格子点データの X,Y 方向積分および平均 ! v_IntX_vx, v_AvrX_vx :: 2 次元(XY)格子点データの X 方向積分および平均 ! x_IntV_vx, x_AvrV_vx :: 2 次元(XY)格子点データの Y 方向積分および平均 ! ! IntX_x, AvrX_x :: 1 次元(X)格子点データの X 方向積分および平均 ! IntV_v, AvrV_v :: 1 次元(Y)格子点データの Y 方向積分および平均 ! !==== スペクトル解析 ! ! ef_EnergyFromStreamfunc_ef :: 流れ関数からエネルギースペクトルを計算 ! ef_EnstrophyFromStreamfunc_ef :: 流れ函数からエンストロフィースペクトルを計算 ! !==== 補間計算 ! ! Interpolate_ef :: 任意の点の値をスペクトルデータから計算する ! use dc_message, only : MessageNotify use mpi implicit none private public ks, ke, kc public js, je, jc public kf_k public ee_mpi_Initial ! 初期化ルーチン public vx_ef, ef_vx ! 基本変換 !!$ public xy_vx ! 基本変換 !!$ public ee_ef, ef_ee ! 基本変換 public ef_Dx_ef, ef_Dy_ef ! 微分 public ef_Lapla_ef, ef_LaplaInv_ef ! 微分 public ef_Jacobian_ef_ef ! ヤコビアン public ef_JacobianZ_ef ! ヤコビアン public IntYX_vx, AvrYX_vx ! 積分・平均 public v_IntX_vx, v_AvrX_vx ! 積分・平均 public x_IntY_vx, x_AvrY_vx ! 積分・平均 public IntX_x, AvrX_x ! 積分・平均 public IntY_v, AvrY_v ! 積分・平均 public Interpolate_ef ! 補間 public x_X, v_Y ! 座標変数 public x_X_Weight, v_Y_Weight ! 座標変数 public vx_X, vx_Y ! 座標変数 public ef_EnergyFromStreamfunc_ef ! エネルギー public ef_EnstrophyFromStreamfunc_ef ! エンストロフィー integer :: im=32, jm=32 ! 格子点の設定(X,Y) integer :: km=10, lm=10 ! 切断波数の設定(X,Y) integer, dimension(:), allocatable :: ke ! X 最大波数 integer, dimension(:), allocatable :: ks ! X 最小波数 integer, dimension(:), allocatable :: kc ! X 波数の数 integer, dimension(:), allocatable :: je ! Y 座標最大格子点 integer, dimension(:), allocatable :: js ! Y 座標最小格子点 integer, dimension(:), allocatable :: jc ! Y 座標格子点数 integer :: np ! MPI プロセス数 integer :: ip ! MPI プロセス ID integer, dimension(:), allocatable :: itj ! FFT 用配列(Y) real(8), dimension(:), allocatable :: tj ! FFT 用配列(Y) integer, dimension(:), allocatable :: iti ! FFT 用配列(X) real(8), dimension(:), allocatable :: ti ! FFT 用配列(X) real(8), dimension(:), allocatable :: x_X ! 格子点座標(X) real(8), dimension(:), allocatable :: v_Y ! 格子点座標(Y) real(8), dimension(:), allocatable :: x_X_Weight ! 格子点重み(X) real(8), dimension(:), allocatable :: v_Y_Weight ! 格子点重み(Y) real(8), dimension(:,:), allocatable :: vx_X ! 格子点(X)座標(2 次元) real(8), dimension(:,:), allocatable :: vx_Y ! 格子点(Y)座標(2 次元) real(8) :: xmin, xmax ! X 座標範囲 real(8) :: ymin, ymax ! Y 座標範囲 real(8) :: xl, yl ! 領域の大きさ(X,Y) real(8), parameter :: pi=3.1415926535897932385D0 save im, jm, km, lm, itj, tj, iti, ti save ks, ke, kc, js, je, jc, np, ip save x_X, v_Y, x_X_Weight, v_Y_Weight, vx_X, vx_Y save xmin, xmax, ymin, ymax, xl, yl contains !--------------- 初期化 ----------------- subroutine ee_mpi_Initial(i_in,j_in,k_in,l_in,xmin_in,xmax_in,ymin_in,ymax_in) ! ! スペクトル変換の格子点数, 波数を設定する. ! 領域の範囲は [0,2pi]x[0,2pi] に決め打ち ! ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで ! 初期設定をしなければならない. ! integer,intent(in) :: i_in ! 格子点数(X) integer,intent(in) :: j_in ! 格子点数(Y) integer,intent(in) :: k_in ! 切断波数(X) integer,intent(in) :: l_in ! 切断波数(Y) real(8),intent(in) :: xmin_in, xmax_in ! X 座標範囲 real(8),intent(in) :: ymin_in, ymax_in ! Y 座標範囲 integer :: ii, jj !!$ integer :: kks, kke, kkc, jjs, jje, jjc integer :: iip integer :: kp, jp integer :: kint, kmod integer :: ierr im = i_in ; jm = j_in km = k_in ; lm = l_in xmin = xmin_in ; xmax = xmax_in ymin = ymin_in ; ymax = ymax_in xl = xmax-xmin ; yl = ymax-ymin CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR) CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR) !!$ kp=km/np+1 !!$ kks=kp*ip !!$ kke=min(kp*(ip+1)-1,km) !!$ if(kke.ge.kks) then !!$ kkc=kke-kks+1 !!$ else !!$ kkc=0 ; kks=0 ; kke=0 !!$ end if !!$ !!$ jp=(jm-1)/np+1 !!$ jjs=jp*ip !!$ jje=min(jp*(ip+1)-1,jm-1) !!$ if(jje.ge.jjs) then !!$ jjc=jje-jjs+1 !!$ else !!$ jjc=0 ; jjs=0 ; jje=0 !!$ end if allocate(ks(0:np-1)) allocate(ke(0:np-1)) allocate(kc(0:np-1)) allocate(js(0:np-1)) allocate(je(0:np-1)) allocate(jc(0:np-1)) allocate(itj(5)) allocate(iti(5)) allocate(tj(jm*2)) allocate(ti(im*2)) call fttrui(im,iti,ti) call fttzui(jm,itj,tj) kint = lm/np kmod = mod(lm,np) kc = kint do iip=0,kmod-1 kc(iip) = kc(iip) + 1 end do ks(0)=0 ; ke(0) = ks(0) + kc(0) -1 do iip=1,np-1 ks(iip) = ks(iip-1) + kc(iip-1) ke(iip) = ks(iip) + kc(iip) - 1 end do do iip=0,np-1 !!$ kp=lm/np+1 !!$ ks(iip)=kp*iip !!$ ke(iip)=min(kp*(iip+1)-1,km) !!$ !!$ !!$ if(ke(iip).ge.ks(iip)) then !!$ kc(iip)=ke(iip)-ks(iip)+1 !!$ else !!$ kc(iip)=0 ; ks(iip)=0 ; ke(iip)=0 !!$ end if jp=(jm-1)/np+1 js(iip)=jp*iip je(iip)=min(jp*(iip+1)-1,jm-1) if(je(iip).ge.js(iip)) then jc(iip)=je(iip)-js(iip)+1 else jc(iip)=0 ; js(iip)=0 ; je(iip)=0 end if enddo allocate(x_X(0:im-1)) allocate(x_X_Weight(0:im-1)) allocate(v_Y(jc(ip))) allocate(v_Y_Weight(jc(ip))) allocate(vx_X(jc(ip),0:im-1)) allocate(vx_Y(jc(ip),0:im-1)) do ii=0,im-1 x_X(ii) = xmin + xl/im*ii enddo x_X_Weight = xl/im do jj=1,jc(ip) v_Y(jj) = ymin + yl/jm*(js(ip)+jj-1) enddo v_Y_Weight = yl/jm vx_X = spread(x_X,1,jc(ip)) vx_Y = spread(v_Y,2,im) call MessageNotify('M','ee_mpi_initial', & 'ee_mpi_module (2020/02/13)is initialized') end subroutine ee_mpi_Initial !--------------- 波数位置問い合わせ ----------------- function kf_k(k) ! ! スペクトルデータ eef(-lm,lm,2*kc) での ! X 方向波数 L の位置情報(第 2 添え字)を調べる. ! ! X 方向波数については ls,ls+1,... le, -le,-le+1,...,-ls ! の順に格納されている. ! integer, intent(IN) :: k ! X 方向波数 integer :: kf_k ! スペクトルデータ第 2 添え字 if ( abs(k) < ks(ip) .OR. ke(ip) < abs(k) ) then kf_k = -1 else if ( k >= 0 ) then kf_k = k-ks(ip)+1 else kf_k = ks(ip)+2*kc(ip)-abs(k) endif end function kf_k !--------------- 基本変換 ----------------- function vx_ef(ef) ! ! スペクトルデータから格子データへ変換する. ! real(8), dimension(js(ip):je(ip),0:im-1) :: vx_ef !(out) 格子点データ real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef !(in) スペクトルデータ real(8), dimension(js(ip):je(ip),2,0:im/2-1) :: vax_work real(8), dimension(kc(ip),0:jm-1,2) :: fya_work real(8), dimension(jc(ip)*im) :: y2 real(8), dimension(kc(ip)*jm*2) :: y3 integer :: k, l fya_work = 0.0D0 do k=1,kc(ip) do l=1,lm fya_work(k,l,1) = ef(l,k) fya_work(k,l,2) = ef(-l,2*kc(ip)-k+1) fya_work(k,jm-l,1) = ef(-l,k) fya_work(k,jm-l,2) = ef(l,2*kc(ip)-k+1) enddo enddo do k=1,kc(ip) fya_work(k,0,1) = ef(0,k) fya_work(k,0,2) = ef(0,2*kc(ip)-k+1) enddo if ( ip == 0 ) then do l=1,lm fya_work(1,l,1) = ef(l,1) fya_work(1,l,2) = ef(-l,1) fya_work(1,jm-l,1) = ef(l,1) fya_work(1,jm-l,2) = -ef(-l,1) end do fya_work(1,0,1) = ef(0,1) endif call fttzub(kc(ip),jm,fya_work,y3,itj,tj) vax_work = vax_fya(fya_work) call fttrub(jc(ip),im,vax_work,y2,iti,ti) vx_ef = reshape(vax_work,(/jc(ip),im/)) end function vx_ef function ef_vx(vx) ! ! 格子データからスペクトルデータへ変換する. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_vx !(out) スペクトルデータ real(8), dimension(js(ip):je(ip),0:im-1), intent(in) :: vx !(in) 格子点データ real(8), dimension(js(ip):je(ip),2,0:im/2-1) :: vax_work real(8), dimension(kc(ip),0:jm-1,2) :: fya_work real(8), dimension(jc(ip)*im) :: y2 real(8), dimension(kc(ip)*jm*2) :: y3 integer :: k, l vax_work = reshape(vx,(/jc(ip),2,im/2/)) call fttruf(jc(ip),im,vax_work,y2,iti,ti) fya_work = fya_vax(vax_work) if ( ip == 0 ) fya_work(1,:,2) = 0.0D0 call fttzuf(kc(ip),jm,fya_work,y3,itj,tj) do k=1,kc(ip) do l=1,lm ef_vx(l,k) = fya_work(k,l,1) ef_vx(-l,2*kc(ip)-k+1) = fya_work(k,l,2) ef_vx(-l,k) = fya_work(k,jm-l,1) ef_vx(l,2*kc(ip)-k+1) = fya_work(k,jm-l,2) enddo enddo do k=1,kc(ip) ef_vx(0,k) = fya_work(k,0,1) ef_vx(0,2*kc(ip)-k+1) = fya_work(k,0,2) enddo if ( ip == 0 ) then do l=1,lm ef_vx(l,1) = fya_work(1,l,1) ef_vx(-l,1) = fya_work(1,l,2) ef_vx(l,2*kc(0)) = 0.0D0 ef_vx(-l,2*kc(0)) = 0.0D0 end do ef_vx(0,1) = fya_work(1,0,1) ef_vx(0,2*kc(0)) = 0.0D0 endif end function ef_vx function vax_fya(fya) real(8), dimension(kc(ip),0:jm-1,2),intent(IN) :: fya real(8), dimension(jc(ip),2,0:im/2-1) :: vax_fya real(8), dimension(0:jm-1,2,kc(ip)) :: yaf_work real(8), dimension(:,:,:), allocatable :: vaf_work integer :: j, k, ips, ipr, ierr, istatus(MPI_STATUS_SIZE) vax_fya = 0.0D0 do k=1,kc(ip) do j=0,jm-1 yaf_work(j,:,k) = fya(k,j,:) enddo enddo do ips=0,np-1 do ipr=0,np-1 if ( ipr /= ips ) then if ( ip == ips ) then allocate(vaf_work(jc(ipr),2,kc(ip))) vaf_work = yaf_work(js(ipr):je(ipr),:,:) call MPI_Send(vaf_work,jc(ipr)*2*kc(ips),MPI_REAL8,& ipr,ips,MPI_COMM_WORLD,IERR) deallocate(vaf_work) endif if ( ip == ipr ) then call MPI_Recv(vax_fya(1,1,ks(ips)),jc(ipr)*2*kc(ips),MPI_REAL8,& ips,ips,MPI_COMM_WORLD,ISTATUS,IERR) endif else if ( ip == ipr ) then vax_fya(:,:,ks(ip):ke(ip)) = yaf_work(js(ip):je(ip),:,1:kc(ip)) endif enddo enddo end function vax_fya function fya_vax(vax) real(8), dimension(jc(ip),2,0:im/2-1),intent(IN) :: vax real(8), dimension(kc(ip),0:jm-1,2) :: fya_vax real(8), dimension(0:jm-1,2,kc(ip)) :: yaf_work real(8), dimension(:,:,:), allocatable :: vaf_work integer :: j, k, ips, ipr, ierr, istatus(MPI_STATUS_SIZE) do ips=0,np-1 do ipr=0,np-1 if ( ipr /= ips ) then if ( ip == ips ) & call MPI_Send(vax(1,1,ks(ipr)),jc(ips)*2*kc(ipr),MPI_REAL8,& ipr,ips,MPI_COMM_WORLD,IERR) if ( ip == ipr ) then allocate(vaf_work(jc(ips),2,kc(ip))) call MPI_Recv(vaf_work(1,1,1),jc(ips)*2*kc(ipr),MPI_REAL8,& ips,ips,MPI_COMM_WORLD,ISTATUS,IERR) yaf_work(js(ips):je(ips),:,:) = vaf_work deallocate(vaf_work) endif else if ( ip == ipr ) then yaf_work(js(ip):je(ip),:,1:kc(ip)) = vax(:,:,ks(ip):ke(ip)) endif enddo enddo do k=1,kc(ip) do j=0,jm-1 fya_vax(k,j,:) = yaf_work(j,:,k) enddo enddo end function fya_vax !--------------- 微分計算 ----------------- function ef_Lapla_ef(ef) ! ! 入力スペクトルデータにラプラシアン(∂xx+∂yy)を作用する. ! ! スペクトルデータのラプラシアンとは, 対応する格子点データに ! ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) をかける ! 計算を行っている. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_Lapla_ef !(out) スペクトルデータのラプラシアン real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef !(in) 入力スペクトルデータ integer kf, k,l ! 作業変数 ef_Lapla_ef = 0.0D0 do kf=1,kc(ip) k = ks(ip) + kf - 1 do l=-lm,lm ef_Lapla_ef(l,kf) = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ef(l,kf) ef_Lapla_ef(l,2*Kc(ip)-kf+1) & = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ef(l,2*kc(ip)-kf+1) enddo enddo end function ef_Lapla_ef function ef_LaplaInv_ef(ef) ! ! 入力スペクトルデータに逆ラプラシアン(∂xx+∂yy)**(-1)を作用する. ! ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに全波数 (k**2 + l**2) で割る ! 計算を行っている. k=l=0 成分には 0 を代入している. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_LaplaInv_ef !(out) スペクトルデータの逆ラプラシアン real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef !(in) スペクトルデータ integer kf, k, l ef_LaplaInv_ef = 0.0D0 do kf=1,kc(ip) k = ks(ip) + kf - 1 do l=-lm,lm if ( k == 0 .AND. l == 0 ) then ef_LaplaInv_ef(l,kf) = 0.0 else ef_LaplaInv_ef(l,kf) = -ef(l,kf)/((2*pi*k/xl)**2+(2*pi*l/yl)**2) ef_LaplaInv_ef(l,2*kc(ip)-kf+1) & = -ef(l,2*kc(ip)-kf+1)/((2*pi*k/xl)**2+(2*pi*l/yl)**2) endif enddo enddo end function ef_LaplaInv_ef function ef_Dx_ef(ef) ! ! 入力スペクトルデータに X 微分(∂x)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 k をかけて ! sin(kx) <-> cos(kx) 成分に入れ換える計算を行っている. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_Dx_ef !(out) スペクトルデータの X 微分 real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef !(in) 入力スペクトルデータ integer kf, k, l ! 作業変数 do kf=1,kc(ip) k = ks(ip) + kf - 1 do l=-lm,lm ef_Dx_ef(l,kf) = -(2*pi*k/xl)*ef(-l,2*kc(ip)-kf+1) ef_Dx_ef(l,2*kc(ip)-kf+1) = -(2*pi*(-k)/xl)*ef(-l,kf) enddo enddo end function ef_Dx_ef function ef_Dy_ef(ef) ! ! 入力スペクトルデータに Y 微分(∂y)を作用する. ! ! スペクトルデータの X 微分とは, 対応する格子点データに Y 微分を ! 作用させたデータのスペクトル変換のことである. ! ! 実際にはスペクトルデータに X 方向波数 l をかけて ! sin(ky) <-> cos(ky) 成分に入れ換える計算を行っている. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_Dy_ef !(out) スペクトルデータの Y 微分 real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef !(in) 入力スペクトルデータ integer kf,l ! 作業変数 do kf=1,kc(ip) if ( ip == 0 .AND. kf == 1 ) then do l=-lm,lm ef_Dy_ef(l,kf) = -(2*pi*l/yl)*ef(-l,kf) ef_Dy_ef(l,2*kc(ip)-kf+1) = 0.0D0 enddo else do l=-lm,lm ef_Dy_ef(l,kf) = -(2*pi*l/yl)*ef(-l,2*kc(ip)-kf+1) ef_Dy_ef(l,2*kc(ip)-kf+1) = -(2*pi*l/yl)*ef(-l,kf) enddo endif enddo end function ef_Dy_ef function ef_Jacobian_ef_ef(ef_a,ef_b) ! ! 2 つのスペクトルデータからヤコビアン ! ! J(A,B)=(∂xA)(∂yB)-(∂yA)(∂xB) ! ! を計算する. ! ! 2 つのスペクトルデータのヤコビアンとは, 対応する 2 つの ! 格子点データのヤコビアンのスペクトル変換のことである. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_Jacobian_ef_ef !(out) 2 つのスペクトルデータのヤコビアン real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_a !(in) 1つ目の入力スペクトルデータ real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_b !(in) 2つ目の入力スペクトルデータ ef_Jacobian_ef_ef = ef_vx( vx_ef(ef_Dx_ef(ef_a))*vx_ef(ef_Dy_ef(ef_b)) & - vx_ef(ef_Dy_ef(ef_a))*vx_ef(ef_Dx_ef(ef_b)) ) end function ef_Jacobian_ef_ef function ef_JacobianZ_ef(ef_zeta) ! ! 渦度スペクトルデータ ζ から流線関数と渦度のヤコビアン ! ! -J(ψ,ζ)=-[(∂xψ)(∂yζ)-(∂yψ)(∂xζ)] ! ! を計算する. ただしψ は (∂xx+∂yy)ψ=ζ を満たす流線関数である. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_JacobianZ_ef !(out) 流線関数と渦度のヤコビアン real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_Zeta !(in) 渦度スペクトルデータ real(8), dimension(jc(ip),0:im-1) :: vx_U ! 速度 X 成分 real(8), dimension(jc(ip),0:im-1) :: vx_V ! 速度 Y 成分 real(8), dimension(-lm:lm,2*kc(ip)) :: ef_UV ! vx_U = - vx_ef(ef_Dy_ef(ef_LaplaInv_ef(ef_Zeta))) vx_V = vx_ef(ef_Dx_ef(ef_LaplaInv_ef(ef_Zeta))) ef_UV = ef_vx(vx_U*vx_V) ef_JacobianZ_ef = - ( ef_Dx_ef(ef_Dx_ef(ef_UV))-ef_Dy_ef(ef_Dy_ef(ef_UV)) )& - ef_Dx_ef(ef_Dy_ef(ef_vx(vx_V**2 - vx_U**2))) end function ef_JacobianZ_ef !--------------- 積分計算 ----------------- function IntYX_vx(vx) ! ! 2 次元格子点データの全領域積分および平均. ! ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた ! 総和を計算している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8) :: IntYX_vx !(out) 積分値 real(8) :: IntYXTmp integer :: i, j, ierr ! 作業変数 IntYXTmp = 0.0d0 do i=0,im-1 do j=1,jc(ip) IntYXTmp = IntYXTmp + vx(j,i) * v_Y_Weight(j) * x_X_Weight(i) enddo enddo CALL MPI_ALLREDUCE(IntYXTmp,IntYX_vx,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntYX_vx function v_IntX_vx(vx) ! ! 2 次元格子点データの X 方向積分 ! ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8), dimension(jc(ip)) :: v_IntX_vx !(out) 積分された 1 次元(Y)格子点データ integer :: i ! 作業変数 v_IntX_vx = 0.0d0 do i=0,im-1 v_IntX_vx(:) = v_IntX_vx(:) + vx(:,i) * x_X_Weight(i) enddo end function v_IntX_vx function x_IntY_vx(vx) ! ! 2 次元格子点データの Y 方向積分 ! ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8), dimension(0:im-1) :: x_IntY_vx !(out) 積分された 1 次元(X)格子点データ real(8), dimension(0:im-1) :: x_IntYTmp integer :: j, ierr ! 作業変数 x_IntYTmp = 0.0d0 do j=1,jc(ip) x_IntYTmp(:) = x_IntYTmp(:) + vx(j,:) * v_Y_Weight(j) enddo CALL MPI_ALLREDUCE(x_IntYTmp,x_IntY_vx,im,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function x_IntY_vx function IntX_x(x) ! ! 1 次元(X)格子点データの X 方向積分 ! ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している. ! real(8), dimension(0:im-1) :: x !(in) 1 次元格子点データ real(8) :: IntX_x !(out) 積分値 IntX_x = sum(x*x_X_Weight) end function IntX_x function IntY_v(v) ! Y 方向積分 ! ! 1 次元(Y)格子点データの Y 方向積分 ! ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している. ! real(8), dimension(jc(ip)) :: v !(in) 1 次元格子点データ real(8) :: IntY_v !(out) 積分値 real(8) :: IntYTmp !(out) 積分値 integer :: ierr IntYTmp = sum(v*v_Y_Weight) CALL MPI_ALLREDUCE(IntYTmp,IntY_v,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntY_v !--------------- 平均計算 ----------------- function AvrYX_vx(vx) ! ! 2 次元格子点データの全領域平均 ! ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8) :: AvrYX_vx !(out) 平均値 AvrYX_vx = IntYX_vx(vx)/(xl*yl) end function AvrYX_vx function v_AvrX_vx(vx) ! ! 2 次元格子点データの X 方向平均 ! ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し, ! x_X_Weight の総和で割ることで平均している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8), dimension(jc(ip)) :: v_AvrX_vx !(out) 平均された 1 次元(Y)格子点データ v_AvrX_vx = v_IntX_vx(vx)/xl end function v_AvrX_vx function x_AvrY_vx(vx) ! ! 2 次元格子点データの Y 方向平均 ! ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し, ! y_Y_Weight の総和で割ることで平均している. ! real(8), dimension(jc(ip),0:im-1) :: vx !(in) 2 次元格子点データ real(8), dimension(0:im-1) :: x_AvrY_vx !(out) 平均された 1 次元(X)格子点データ x_AvrY_vx = x_IntY_vx(vx)/yl end function x_AvrY_vx function AvrX_x(x) ! ! 1 次元(X)格子点データの X 方向平均 ! ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し, ! x_X_Weight の総和で割ることで平均している. ! real(8), dimension(0:im-1) :: x !(in) 1 次元格子点データ real(8) :: AvrX_x !(out) 平均値 AvrX_x = IntX_x(x)/xl end function AvrX_x function AvrY_v(v) ! ! 1 次元(Y)格子点データの Y 方向平均 ! ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し, ! y_Y_Weight の総和で割ることで平均している. ! real(8), dimension(jc(ip)) :: v !(in) 1 次元格子点データ real(8) :: AvrY_v !(out) 平均値 AvrY_v = IntY_v(v)/yl end function AvrY_v !--------------- スペクトル計算 ----------------- function ef_EnergyFromStreamfunc_ef(ef_StrFunc) ! ! 流線関数からエネルギースペクトルを計算する. ! ! E_kl = (1/2)(k^2+l^2)|\psi_kl|^2 ! ! * E_kl の総和が場の平均運動エネルギーとなる. ! * それに領域の面積をかけると全運動エネルギーとなる. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_EnergyFromStreamfunc_ef ! エネルギースペクトル real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_StrFunc ! 流線関数 integer kf, k,l ! 作業変数 do kf=1,kc(ip) k = ks(ip) + kf -1 do l=-lm,lm ef_EnergyFromStreamfunc_ef(l,kf) & = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 ) & * (ef_StrFunc(l,kf)**2 + ef_StrFunc(-l,2*kc(ip)-kf+1)**2) ef_EnergyFromStreamfunc_ef(l,2*kc(ip)-kf+1) & = ef_EnergyFromStreamfunc_ef(l,kf) enddo enddo end function ef_EnergyFromStreamfunc_ef function ef_EnstrophyFromStreamfunc_ef(ef_StrFunc) ! ! 流線関数からエンストロフィースペクトルを計算する. ! ! Q_kl = (1/2)(k^2+l^2)^2|\psi_kl|^2 ! ! * Q_kl の総和が場の平均エンストロフィーとなる. ! * それに領域の面積をかけると全エンストロフィーとなる. ! real(8), dimension(-lm:lm,2*kc(ip)) :: ef_EnstrophyFromStreamfunc_ef ! エンストロフィーースペクトル real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_StrFunc ! 流線関数 integer kf, k,l ! 作業変数 do kf=1,kc(ip) k = ks(ip) + kf -1 do l=-lm,lm ef_EnstrophyFromStreamfunc_ef(l,kf) & = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 )**2 & * (ef_StrFunc(l,kf)**2 + ef_StrFunc(-l,2*kc(ip)-kf+1)**2) ef_EnstrophyFromStreamfunc_ef(l,2*kc(ip)-kf+1) & = ef_EnstrophyFromStreamfunc_ef(l,kf) enddo enddo end function ef_EnstrophyFromStreamfunc_ef !--------------- 補間計算 ----------------- function Interpolate_ef( ef_Data, x, y ) real(8), intent(IN) :: ef_data(-lm:lm,2*kc(ip)) ! スペクトルデータ real(8), intent(IN) :: x ! 補間する点の x 座標 real(8), intent(IN) :: y ! 補間する点の y 座標 real(8) :: Interpolate_ef ! 補間した値 real(8) :: InterpolateTmp ! 作業変数 integer :: kf, k, l, ierr real(8) :: xx, yy xx =(2*PI/xl)*(x - xmin) yy =(2*PI/yl)*(y - ymin) if ( ip == 0 ) then InterpolateTmp = ef_Data(0,1) else InterpolateTmp = 0.0D0 endif ! l=0 if ( ip == 0 ) then do kf=2,kc(ip) k = ks(ip) + kf -1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(0,kf)*cos(k*xx) - ef_Data(0,kf_k(-k))*sin(k*xx) ) end do else do kf=1,kc(ip) k = ks(ip) + kf -1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(0,kf)*cos(k*xx) - ef_Data(0,kf_k(-k))*sin(k*xx) ) end do endif ! k=0 if ( ip == 0 ) then do l=1,lm InterpolateTmp = InterpolateTmp & + 2*( ef_Data(l,kf_k(0))*cos(l*yy) - ef_Data(-l,kf_k(0))*sin(l*yy) ) end do end if ! k*l > 0 do l=1,lm if ( ip == 0 ) then do kf=2,kc(ip) k = ks(ip) + kf -1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(l,kf)*( cos(k*xx)*cos(l*yy) & - sin(k*xx)*sin(l*yy) ) & -ef_Data(-l,kf_k(-k))*( sin(k*xx)*cos(l*yy) & + cos(k*xx)*sin(l*yy) ) ) end do else do kf=1,kc(ip) k = ks(ip) + kf - 1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(l,kf)*( cos(k*xx)*cos(l*yy) & - sin(k*xx)*sin(l*yy) ) & -ef_Data(-l,kf_k(-k))*( sin(k*xx)*cos(l*yy) & + cos(k*xx)*sin(l*yy) ) ) end do endif end do ! k*l < 0 do l=1,lm if ( ip == 0 ) then do kf=2,kc(ip) k = ks(ip) + kf - 1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(-l,kf)*( cos(k*xx)*cos(l*yy) & + sin(k*xx)*sin(l*yy) ) & -ef_Data(l,kf_k(-k))*( sin(k*xx)*cos(l*yy) & - cos(k*xx)*sin(l*yy) ) ) end do else do kf=1,kc(ip) k = ks(ip) + kf - 1 InterpolateTmp = InterpolateTmp & + 2*( ef_Data(-l,kf)*( cos(k*xx)*cos(l*yy) & + sin(k*xx)*sin(l*yy) ) & -ef_Data(l,kf_k(-k))*( sin(k*xx)*cos(l*yy) & - cos(k*xx)*sin(l*yy) ) ) end do endif end do CALL MPI_ALLREDUCE(InterpolateTmp,Interpolate_ef,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function Interpolate_ef end module ee_mpi_module