ループ途中のNumpy配列を保存したいが、エラーが出てしまう。

Question

前提・実現したいこと

python3.6を用いて数値計算を行い、その過程でループ途中のNumpy配列を保存したいと考えています。その際、np.loadtxtおよびsavetxtを使っているのですが、保存の際にエラーが出てしまいます。

発生している問題・エラーメッセージ

C:\Users\USER\python\FDTD\シン・ELECTROMAGNETIC SIMULATION USING THE FDTD(202201
10より)\レーザーのシミュレーション\テスト\増幅の周波数特性測定\ローレンツ媒質導
入>python fd2d_sub_500cell.py
  0%|                                       | 2/4000 [00:04<2:44:14,  2.46s/it]
Traceback (most recent call last):
  File "fd2d_sub_500cell.py", line 82, in <module>
    sz2 = np.loadtxt(os.path.join(dir2,'sz_field_{0}.txt').format(before_time))
  File "C:\Users\USER\AppData\Local\Programs\Python\Python36\lib\site-packages\n
umpy\lib\npyio.py", line 961, in loadtxt
    fh = np.lib._datasource.open(fname, 'rt', encoding=encoding)
  File "C:\Users\USER\AppData\Local\Programs\Python\Python36\lib\site-packages\n
umpy\lib\_datasource.py", line 195, in open
    return ds.open(path, mode, encoding=encoding, newline=newline)
  File "C:\Users\USER\AppData\Local\Programs\Python\Python36\lib\site-packages\n
umpy\lib\_datasource.py", line 535, in open
    raise IOError("%s not found." % path)
OSError: sz/sz_field_1.txt not found.

該当のソースコード

以下はソースコードです。
szの初期値は全要素が0の500x500のNumpy配列です。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
fd2d 細々といじる用
"""

import os
import  numpy  as  np 
from  math  import  pi, sin, cos, exp, sqrt
from  matplotlib  import  pyplot  as  plt 
from  mpl_toolkits.mplot3d.axes3d  import  Axes3D, get_test_data
from tqdm import tqdm

import fd2d_init_value as fi

# データ用のファイル作成
os.mkdir('data')
os.mkdir('sz')

# ディレクトリ指定用
dir1 = 'data/'
dir2 = 'sz/'

# 初期時間
tau = 20

# 周波数対強度のリスト
amp = []

# 強度計算用の配列(実部と虚部)
re_ez = []
im_ez = []

# 特定の点を入れておくためのリスト
point_x = range(1, fi.nsteps + 1)
point_y = []

# 補助項の定義
sz = fi.sz
sz2 = np.zeros((fi.ie, fi.je))

# 各パラメータの読み込み
epsr = fi.epsr
delta0 = fi.delta0
f0 = fi.f0
omega0 = fi.omega0
eps1 = fi.eps1

# 対数スケールの配列作成
start = 10^6
stop = 10^10
freq_list = np.logspace(start,stop)

# 補助項の係数の計算
base = 1 + fi.dt * delta0 * omega0
top1 = 2 - (fi.dt**2) * (omega0**2)
top2 = 1 - (fi.dt * delta0 * omega0)
top3 = (fi.dt**2) * (omega0**2) * eps1

# Main  FDTD Loop 
for  time_step  in  tqdm( range(1,  fi.nsteps  + 1) ):
    
    # 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]) 

    # Source1
    if abs(time_step - tau*3) < tau*3:
        gause_pulse =  ( exp( - ( time_step - tau*3 ) **2 / tau ) )
        fi.dz[250][250] = gause_pulse
    
    else:
        continue

    #補助項Sの2ステップ前の読み込み
    if time_step < 3:
        continue
    else:
        before_time = time_step - 2
        sz2 = np.loadtxt(os.path.join(dir2,'sz_field_{0}.txt').format(before_time))

    # 補助項Sの計算
    for i in range(fi.mir_y_bey, fi.mir_y_ove):
        for j in range(fi.mir_x_sta, fi.mir_x_fin):
            sz[i,j] = (sz[i,j] * top1 / base) - (sz2[i,j] * top2 / base) + (fi.ez[i,j] * top3 / base)
    
    # Sの保存
    np.savetxt(os.path.join(dir2,'sz_field_{0}.txt').format(time_step), sz) 

    # Calculate the Ez field
    for j in range(0, fi.je):
        for i in range(0, fi.ie):
            fi.ez[i, j] = (fi.dz[i, j] - sz[i,j]) / epsr

    # 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])

    # 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])

    # 強度分布の計算
    for i, freq in enumerate(freq_list):
        arg = 2 * pi * freq * fi.dt * time_step
        re_ez[i] = re_ez[i] + cos( arg ) * fi.ez[300][300]
        im_ez[i] = im_ez[i] - sin( arg ) * fi.ez[300][300]
        amp[i] = sqrt( (re_ez[i]) **2 + (im_ez[i]) **2 )

    point_y.append(fi.ez[fi.ie - 30, fi.je - 30])

    if time_step % 100 ==0:
        e_field = np.copy(fi.ez)
        np.savetxt(os.path.join(dir1, 'e_result_{0}.txt'. format(time_step)), e_field)
        e_center = np.copy(fi.ez[:,250])
        np.savetxt(os.path.join(dir1, 'e_center_{0}.txt'. format(time_step)), e_center)

    else :
        continue

f = open('a_point_list(raw).txt', 'w')
for x in point_y:
    f.write(str(x) + "\n")
f.close()

f = open('freq_vs_amp.txt', 'w')
for x in amp:
    f.write(str(x) + "\n")
f.close()

# 特定の点の変化を表すグラフの作成
fig =  plt.figure(figsize=(8,  8)) 
ax = fig.add_subplot(1,1,1)
ax.plot(point_x, point_y)

fig.tight_layout()
plt.savefig("point_vary.png")

試したこと

savetxtと.formatの相性が悪いのかと思い、書き方をいくつか変更してみましたが動きませんでした。ディレクトリ指定の方法が間違っているのかとも思ったのですが、何度確認しても間違いが見つけられませんでした。

補足情報（FW/ツールのバージョンなど）

numpy 1.19.5

Answer 1

理由を知りたければ、savetxtとloadtxtの前に、ファイル名を表示するコードを入れておいて動かしてください。

例えば、

    # Sの保存
    np.savetxt(os.path.join(dir2,'sz_field_{0}.txt').format(time_step), sz) 

の部分に以下のprint文を追加するなどです。

    # Sの保存
    print(os.path.join(dir2,'sz_field_{0}.txt').format(time_step))
    np.savetxt(os.path.join(dir2,'sz_field_{0}.txt').format(time_step), sz) 

Comments

[ShibaSamo]

ppaul様の助言通り、print文を追加してみたところ、「print文があるのに文が表示されない」という結果になりました。無学で恐縮なのですが、特定の関数はループ中では異なる挙動をするということなのでしょうか？

[ppaul]

「print文があるのに文が表示されない」というのは、print文が実行されると期待しているのに、実際は実行されていないという場合が多いです。その原因は、if文の条件が間違っているとか、for文で回しているつもりなのに実は回転数がゼロだったりします。if文の前やfor文の前に、そういう情報を出力するようなprint文を入れてみてください。

[ShibaSamo]

ありがとうございます!
基本を忘れてしまっていました。
ループ中のcontinueをpassにすることで問題を解消することができました!