FDTD 遠方界のプログラム

質問

FDTD法で遠方界を計算するプログラムを作ったのですがコンパイルはうまくいっても実行結果がちゃんと表示されません。プログラムに問題があるのか、大きな値を与えるとbus errorになってしまいます。どなたか遠方界のプリグラムに詳しい方がいらっしゃったら教えて下さい。よろしくお願いします。
プログラム内容
ubroutine far_field
implicit none

call far_field_yz_plane
call far_field_zx_plane
call far_field_xy_plane

return
end subroutine

!===================================================================================
! y-z 平面に対して
!===================================================================================
subroutine far_field_yz_plane
use consts
use fdtd
implicit none
integer :: i, j, k,ie
real(8) :: plane_yz, fac_yz

plane_yz = dy*dz
fac_yz = plane_yz/(4.0d0*pi*c*dt)
ie = nx-6
i = ie+1
do j = 5, ny-5
do k = 5, nz-5
x = (i-0.5d0-ic)*dx
y = (j-jc)*dy
z = (k-kc)*dz
eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1))
ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k))
hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k))
hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1))
do l = 1, nang
time_e = time-dt
time_h = time-dt/2.0d0
timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c
tc_e = time_e+timesh
tc_h = time_h+timesh
m_e = int(tc_e/dt+0.5d0)
m_h = int(tc_h/dt+0.5d0)
a_e = (0.5d0+tc_e/dt-m_e)*fac_yz
b_e = (0.5d0-tc_e/dt+m_e)*fac_yz
ab_e = a_e-b_e
a_h = (0.5d0+tc_h/dt-m_h)*fac_yz
b_h = (0.5d0-tc_h/dt+m_h)*fac_yz
ab_h = a_h-b_h
wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h
wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h
uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e
uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e
wy(l,m_h) = wy(l,m_h)-hzs*ab_h
wz(l,m_h) = wz(l,m_h)+hys*ab_h
uy(l,m_e) = uy(l,m_e)+ezs*ab_e
uz(l,m_e) = uz(l,m_e)-eys*ab_e
wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h
wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h
uy(l,m_e+1) = uy(l,m_e+1)-ezs*a_e
uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e
enddo
enddo
enddo

ie = 5
i = ie
do j = 5, ny-5
do k = 5, nz-5
x = (i-0.5d0-ic)*dx
y = (j-jc)*dy
z = (k-kc)*dz
eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1))
ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k))
hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k))
hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1))
do l = 1, nang
time_e = time-dt
time_h = time-dt/2.0d0
timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c
tc_e = time_e+timesh
tc_h = time_h+timesh
m_e = int(tc_e/dt+0.5d0)
m_h = int(tc_h/dt+0.5d0)
a_e = (0.5d0+tc_e/dt-m_e)*fac_yz
b_e = (0.5d0-tc_e/dt+m_e)*fac_yz
ab_e = a_e-b_e
a_h = (0.5d0+tc_h/dt-m_h)*fac_yz
b_h = (0.5d0-tc_h/dt+m_h)*fac_yz
ab_h =a_h-b_h
wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h
wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h
uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e
uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e
wy(l,m_h) = wy(l,m_h)-hzs*ab_h
wz(l,m_h) = wz(l,m_h)+hys*ab_h
uy(l,m_e) = uy(l,m_e)+ezs*ab_e
uz(l,m_e) = uz(l,m_e)-eys*ab_e
wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h
wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h
uy(l,m_h+1) = uy(l,m_e+1)-ezs*a_e
uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e
enddo
enddo
enddo

return
end subroutine

zx平面、xy平面についても同様に書きました。

回答No.3

kuma4578 · ベストアンサー

私は大学でFDTD法による遠方界計算とMie散乱の理論値との整合性を調べる研究をしています。
そのため解析領域内に設置した媒質は球形です
（媒質の形状によってdirsに格納する値は違ってくると思います）

私の場合、球形媒質の中心点をIc, Jc, Kcとし、dirsには三次元極座標系の値を格納しています。
具体的には、φを0度固定でθを360度回転（XZ面）とφを90度固定でθを360度回転(YZ面)の値です。

教科書中のrfの表現あってるんですかね？
rfにセルサイズをかけるとうまくプログラムが回るのですが・・・

回答No.2

kuma4578

FDTD法を用いて遠方界を計算しているものです
（結果は出力できていますがうまくいってるのかは正直微妙です）

これは散乱体を囲む閉局面の、y-z 平面の上部と下部を計算しているプログラムですよね？

疑問１
上部と下部ではベクトルが反対方向なので計算方法が同じでいいのでしょうか？
自分は上部と下部では＋ーを逆にした計算を行っています（宇野さんの教科書にはこの点を詳しく書かれていませんよね）

疑問２
遠方界計算の範囲に吸収境界が含まれていませんか？5~max-5までとなると,例えばPMLを用いている場合4層ということでしょうか？

疑問３
dirs,rfにはどのような値を格納していますか？

補足

ご回答ありがとうございます。

疑問１に関して、確かに上下でベクトルの符号が変わるので、上部と下部では符号が逆転しますね。訂正しました。

疑問２に関してですが、私は吸収境界条件はMurの２次を用いています。PMLは未だ勉強不足でして。

疑問３に関してですがdirsは全く値を格納していませんでした。しかし、どのような値を入れてよいかわかりません。rfには閉曲面（立方体）の対角線の1/2の値（sqrt((nx-10)^2+(ny-10)^2+(nz-10)^2)/2を格納しています。

もう一つ質問なのですが、Ic, Jc, Kcにはどのような値を格納していますか？dirsと同様によいアドバイスがあればお願いします。
よろしくお願いします。

回答No.1

f272

> プログラムに問題があるのか

と思っているのなら，そのエラーになるプログラムをすべて見せないと分かりません。
use consts
use fdtd
とあるのに，これらのモジュールはどうなっているのでしょう？
またsubroutine far_field_yz_planeの中身を長々と書いていますが，例えばこれを（計算結果がおかしくなるとしても）適当に短くして，同じエラーを起こす最低の長さのプログラムに書きなおすことができるでしょう。そのようにして，回答者が同じエラーを起こすことのできる完全なプログラムを提示してください。

お礼

すみません。プログラムを載せようとしたら字数制限があり全部入りきりませんでした。また、自分でプログラムをよく見なしてみます。ありがとうございました。

補足

f272さんのおっしゃる通りモジュールの部分も提示しなければわからないですよね。大変失礼いたしました。

module consts
real(8), parameter :: pi = 3.141592653589793d0 !円周率
real(8), parameter :: c = 2.998d8 !光速 c [m/sec]
real(8), parameter :: epsilon0 = 8.854d-12 !真空の誘電率 [F/m]
real(8), parameter :: mu0 = 4.0d-7*pi !真空の透磁率 [H/m]
real(8), parameter :: z0 = 120.0d0*pi !波動インピーダンス [Ω]
end module

module fdtd
! ********** 遠方界 *********
integer, parameter :: nang = 180
integer :: l
real(8) :: dirs(500,3)
integer :: m_e, m_h
integer :: ic, jc, kc
real(8), parameter :: rf = 6075d0
real(8) :: x, y, z
real(8) :: exs, eys, ezs
real(8) :: hxs, hys, hzs
real(8) :: time_e, time_h, timesh
real(8) :: tc_e, tc_h
real(8) :: a_e, b_e, ab_e, a_h, b_h, ab_h
real(8) :: wx(500,500), wy(500,500), wz(500,500)
real(8) :: ux(500,500), uy(500,500), uz(500,500)

end module

subroutine far_field
implicit none

call far_field_yz_plane
call far_field_zx_plane
call far_field_xy_plane

return
end subroutine

!===================================================================================
! y-z 平面に対して
!===================================================================================
subroutine far_field_yz_plane
use consts
use fdtd
implicit none
integer :: i, j, k, ie
real(8) :: plane_yz, fac_yz

plane_yz = dy*dz
fac_yz = plane_yz/(4.0d0*pi*c*dt)
ie = nx-6
i = ie+1
do j = 5, ny-5
do k = 5, nz-5
x = (i-0.5d0-ic)*dx
y = (j-jc)*dy
z = (k-kc)*dz
eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1))
ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k))
hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k))
hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1))
do l = 1, nang
time_e = time-dt
time_h = time-dt/2.0d0
timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c
tc_e = time_e+timesh
tc_h = time_h+timesh
m_e = int(tc_e/dt+0.5d0)
m_h = int(tc_h/dt+0.5d0)
a_e = (0.5d0+tc_e/dt-m_e)*fac_yz
b_e = (0.5d0-tc_e/dt+m_e)*fac_yz
ab_e = a_e-b_e
a_h = (0.5d0+tc_h/dt-m_h)*fac_yz
b_h = (0.5d0-tc_h/dt+m_h)*fac_yz
ab_h = a_h-b_h
wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h
wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h
uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e
uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e
wy(l,m_h) = wy(l,m_h)-hzs*ab_h
wz(l,m_h) = wz(l,m_h)+hys*ab_h
uy(l,m_e) = uy(l,m_e)+ezs*ab_e
uz(l,m_e) = uz(l,m_e)-eys*ab_e
wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h
wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h
uy(l,m_e+1) = uy(l,m_e+1)-ezs*a_e
uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e
enddo
enddo
enddo