- ベストアンサー
FDTD法で遠方界を計算するプログラムの実行結果が表示されない問題
- FDTD法で遠方界を計算するプログラムを作成しましたが、実行結果が正しく表示されません。
- 大きな値を与えるとbus errorが発生する問題もあります。
- 遠方界のプログラムに詳しい方にアドバイスをお願いします。
- 528612
- お礼率36% (33/90)
- その他（プログラミング・開発）
- 回答数3
- ありがとう数24
- みんなの回答 （3）
- 専門家の回答
質問者が選んだベストアンサー
私は大学でFDTD法による遠方界計算とMie散乱の理論値との整合性を調べる研究をしています。 そのため解析領域内に設置した媒質は球形です （媒質の形状によってdirsに格納する値は違ってくると思います） 私の場合、球形媒質の中心点をIc, Jc, Kcとし、dirsには三次元極座標系の値を格納しています。 具体的には、φを0度固定でθを360度回転（XZ面）とφを90度固定でθを360度回転(YZ面)の値です。 教科書中のrfの表現あってるんですかね？ rfにセルサイズをかけるとうまくプログラムが回るのですが・・・
その他の回答 (2)
- kuma4578
- ベストアンサー率100% (1/1)
FDTD法を用いて遠方界を計算しているものです （結果は出力できていますがうまくいってるのかは正直微妙です） これは散乱体を囲む閉局面の、y-z 平面の上部と下部を計算しているプログラムですよね？ 疑問１ 上部と下部ではベクトルが反対方向なので計算方法が同じでいいのでしょうか？ 自分は上部と下部では＋ーを逆にした計算を行っています（宇野さんの教科書にはこの点を詳しく書かれていませんよね） 疑問２ 遠方界計算の範囲に吸収境界が含まれていませんか？5~max-5までとなると,例えばPMLを用いている場合4層ということでしょうか？ 疑問３ dirs,rfにはどのような値を格納していますか？
- f272
- ベストアンサー率46% (8468/18130)
> プログラムに問題があるのか と思っているのなら，そのエラーになるプログラムをすべて見せないと分かりません。 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
補足
ご回答ありがとうございます。 疑問１に関して、確かに上下でベクトルの符号が変わるので、上部と下部では符号が逆転しますね。訂正しました。 疑問２に関してですが、私は吸収境界条件はMurの２次を用いています。PMLは未だ勉強不足でして。 疑問３に関してですがdirsは全く値を格納していませんでした。しかし、どのような値を入れてよいかわかりません。rfには閉曲面（立方体）の対角線の1/2の値（sqrt((nx-10)^2+(ny-10)^2+(nz-10)^2)/2を格納しています。 もう一つ質問なのですが、Ic, Jc, Kcにはどのような値を格納していますか？dirsと同様によいアドバイスがあればお願いします。 よろしくお願いします。