前提・実現したいこと
python3.6を使ってFDTD法の計算結果を描画したいと考えています。
しかし、計算自体は行えたのですが、描画をする際にエラーが出ます。
エラーコードを読むに、Main FDTD Loopの中の計算には問題がなく、その後のmatplotlibの部分に問題があると思うのですが、どうにもできず途方に暮れています。
ソースコードの
gaz_arrange.txt
gbz_arrange.txt
は、200x200の配列(?)を入れたテキストファイルです。
ie,je,ia,ibなどは全て整数値です。
以下がエラーメッセージです。
C:\Users\USER\python\FDTD\ELECTROMAGNETIC SIMULATION USING THE FDTD> python fd2d
_main.py
100%|███████████████████████████████████████| 194/194 [00:01<00:00, 135.79it/s]
100%|████████████████████████████████████████| 300/300 [03:47<00:00, 1.32it/s]
Traceback (most recent call last):
File "fd2d_main.py", line 123, in <module>
plot_e_field(ax, plotting_point[ 'data_to_plot'], plotting_point[ 'num_step
s'], plotting_point['label'])
File "fd2d_main.py", line 105, in plot_e_field
ax.plot_surface(Y, X, data, rstride=1, cstride=1, color='white', edgecolor
='black', linewidth=.25)
File "C:\Users\USER\AppData\Local\Programs\Python\Python36\lib\site-packages\m
pl_toolkits\mplot3d\axes3d.py", line 1555, in plot_surface
if Z.ndim != 2:
AttributeError: 'NoneType' object has no attribute 'ndim'
また、次のものがソースコードです。
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
fd2d メインのプログラム
"""
import numpy as np
import scipy.constants as sc
from math import pi, sin, cos
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D, get_test_data
from tqdm import tqdm
import fd2d_scatter_arrange as fs
import fd2d_init_value as fi
# 誘電体の配置の読み込み
gaz = np.loadtxt('gaz_arrange.txt')
gbz = np.loadtxt('gbz_arrange.txt')
# Dictionary to keep track of desired points for plotting
plotting_points = [{'label': 'a','num_steps':4700, 'data_to_plot':None}, {'label':'b', 'num_steps':4800, 'data_to_plot':None}, { 'label':'c' ,'num_steps':4900, 'data_to_plot':None}, {'label':'d', 'num_steps':5000, 'data_to_plot':None}, ]
# Main FDTD Loop
for time_step in tqdm( range(1, fi.nsteps + 1) ):
# Incident Ez values
for j in range(1, fi.je):
fi.ez_inc[j] = fi.ez_inc[j] + 0.5 * (fi.hx_inc[j - 1] - fi.hx_inc[j])
# Absorbing Boundary Conditions
fi.ez_inc[0] = fi.boundary_low.pop(0)
fi.boundary_low.append(fi.ez_inc[1])
fi.ez_inc[fi.je - 1] = fi.boundary_high.pop(0)
fi.boundary_high.append(fi.ez_inc[fi.je - 2])
# Calculate the Dz field
for j in range(1, fi.je):
for i in range(1, fi.ie):
fi.dz[i, j] = fi.gi3[i] * fi.gj3[j] * fi.dz[i, j] + fi.gi2[i] * fi.gj2[j] * 0.5 * (fi.hy[i, j] - fi.hy[i - 1, j] - fi.hx[i, j] + fi.hx[i, j - 1])
# Source
pulse = sin(2 * pi * fi.freq * fi.dt * (fi.t0 - time_step))
fi.ez_inc[3] = pulse
# Incident Dz values
for i in range(fi.ia, fi.ib):
fi.dz[i, fi.ja] = fi.dz[i, fi.ja] + 0.5 * fi.hx_inc[fi.ja - 1 ]
fi.dz[i, fi.jb] = fi.dz[i, fi.jb] - 0.5 * fi.hx_inc[fi.jb - 1 ]
# Calculate the Ez field
for j in range(0, fi.je):
for i in range(0, fi.ie):
fi.ez[i, j] = gaz[i, j] * (fi.dz[i, j] - fi.iz[i,j])
fi.iz[i, j] = fi.iz[i, j] + gbz[i, j] * fi.ez[i,j]
# Calculate the Incident Hx
for j in range(0, fi.je - 1):
fi.hx_inc[j] = fi.hx_inc[j] + 0.5 * (fi.ez_inc[j] - fi.ez_inc[j + 1])
# Calculate the Hx field
for j in range(0, fi.je - 1):
for i in range(0, fi.ie - 1):
curl_e = fi.ez[i, j] - fi.ez[i, j + 1]
fi.ihx[i, j] = fi.ihx[i, j] + curl_e
fi.hx[i, j] = fi.fj3[j] * fi.hx[i, j] + fi.fj2[j] * (0.5 * curl_e + fi.fil[i] * fi.ihx[i,j])
# Incident Hx values
for i in range(fi.ia, fi.ib):
fi.hx[i, fi.ja - 1] = fi.hx[i, fi.ja - 1] + 0.5 * fi.ez_inc[fi.ja]
fi.hx[i, fi.jb] = fi.hx[i, fi.jb] - 0.5 * fi.ez_inc[fi.jb]
# Calculate the Hy field
for j in range(0, fi.je - 1):
for i in range(0, fi.ie - 1):
curl_e = fi.ez[i, j] - fi.ez[i + 1, j]
fi.ihy[i, j] = fi.ihy[i, j] + curl_e
fi.hy[i, j] = fi.fi3[i] * fi.hy[i, j] - fi.fi2[i] * (0.5 * curl_e + fi.fjl[j] * fi.ihy[i,j])
# Incident Hy values
for j in range(fi.ja, fi.jb):
fi.hy[fi.ia - 1, j] = fi.hy[fi.ia - 1, j] - 0.5 * fi.ez_inc[j]
fi.hy[fi.ib - 1, j] = fi.hy[fi.ib - 1, j] + 0.5 * fi.ez_inc[j]
# Save data at certain points for later plotting
for plotting_point in plotting_points:
if time_step == plotting_point['num_steps']:
plotting_point['data_to_plot'] = np.copy(fi.ez)
# Plot Fig
plt.rcParams['font.size']= 12
plt.rcParams['grid.color']= 'gray'
plt.rcParams['grid.linestyle'] = 'dotted'
fig = plt.figure(figsize=(8, 8))
X, Y = np.meshgrid(range(fi.je), range(fi.ie))
def plot_e_field(ax, data, timestep, label):
"""3d Plot of E field at a single timestep"""
ax.set_zlim(-0.5, 1)
ax.view_init(elev=15, azim=25)
ax.plot_surface(Y, X, np.array(data), rstride=1, cstride=1, color='white', edgecolor='black', linewidth=.25)
ax.zaxis.set_rotate_label(False)
ax.set_zlabel(r' $E_{Z}$', rotation=90, labelpad=10, fontsize=14)
ax.set_zticks([-0.5, 0, 0.5, 1])
ax.set_xlabel('cm')
ax.set_ylabel('cm')
ax.set_xticks(np.arange(0, 200, step=50))
ax.set_yticks(np.arange(0, 200, step=50))
ax.text2D( 0.25, 0.3, "T = {}".format(timestep) , transform=ax.transAxes)
ax.xaxis.pane.fill = ax.yaxis.pane.fill = ax.zaxis.pane.fill = False
plt.gca().patch.set_facecolor( 'white')
ax.text2D(-0.05, 0.8, "( {} ) ".format(label), transform=ax.transAxes)
ax.dist = 11
# Plot the E field at each of the four time steps saved earlier
for subplot_num, plotting_point in enumerate(plotting_points):
ax = fig.add_subplot(2, 2, subplot_num + 1, projection='3d' )
plot_e_field(ax, plotting_point[ 'data_to_plot'], plotting_point[ 'num_steps'], plotting_point['label'])
fig.tight_layout()
plt.savefig("image.png")
plt.show()
試したこと
plot_surfaceの最初の引数3つの内に問題があると思い、dataをNumpy配列にしてみたのですが駄目でした。
また、コードを書いた当初は問題なく実行できていたのですが、zticksとxticks、yticksの値を変更した後にエラーが出るようになりました。
この部分に問題があるのかと値を戻したところ、エラーのみが残ってしまいました。
元々はplt.savefig("image.png")の行で終わっていたのですが、その後にplt.show()が合った方が挙動が安定するという話を聞きましたので、エラーを吐くようになってから文を追加しました。
補足情報(FW/ツールのバージョンなど)
matplotlib 3.3.4
numpy 1.19.5