わからない点
cahrgeのデータ保存場所の名前。外挿した後のcharge(HUT
,HST
,MST
,MUT
)名前はそれぞれ変えたほうがいいのか。各所でオブジェクトを表示する処理を行ったがないほうがいいのか。各格子点に割り振られた電荷を用いて電位計算を行う。計算した電位により発生した電界を内挿により荷電粒子に加わる静電気力として与える。
下記に現時点のプログラムを記載します。
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
import itertools
import csv
import threading
import random
import time
from numba import jit # Better to write at the beginning of the code
import decimal
import math
# 条件設定
delta_L = 1.0
delta = 1.0
n_inter = 0
nx = 1000 #x方向のメッシュの数
ny = 1000 #y方向のメッシュの数
xmin = 0 #xの最小値
xmax = 1000*delta #xの最大値
ymin = 0 #yの最小値
ymax = 1000*delta #yの最大値
n0 = 1e12; #density in [/m^3]
phi0 = 0; #reference potential
Te = 1; #electron temperature in eV
v_drift = 7000; #ion injection velocity, 7[km/s]
phi_p = -7; #plate potential
convergence_criterion = 1.0e-5 #前回との差
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
# 初期状態
phi = np.zeros((ny, nx))
electrode = np.zeros((ny, nx))
x = np.linspace(xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)
# for SOR method
aa_recta=0.5*(np.cos(np.pi/nx)+np.cos(np.pi/ny)) #
omega_SOR_recta = 2/(1+np.sqrt(1-aa_recta**2)) # 長方形領域に対する最適加速パラメータ
print("omega_SOR_rect=",omega_SOR_recta)
T=1 #粒子計算における時間ステップ
# 粒子データを格納するための辞書を初期化
particle_data = {}
# 粒子数
n_particles = 100
# 粒子の座標と名前をランダムに生成して粒子データに追加
for i in range(n_particles):
a = (random.randint(0, 1000000)) / 1000
b = (random.randint(0, 1000000)) / 1000
particle_name = f"Particle_{i}" # 粒子の名前を設定
particle_data[particle_name] = {"a": a, "b": b, "charge": []} # chargeキーを追加し初期化
# 粒子ごとに計算と結果の保存を行う
for name, data in particle_data.items():
a = data["a"]
b = data["b"]
# 粒子ごとの計算
aa = math.floor(a)
bb = math.floor(b)
a = math.floor(a * 1000)
b = math.floor(b * 1000)
aa = aa * 1000
bb = bb * 1000
a1 = (a - aa)
a2 = ((aa + 1000) - a)
b1 = (b - bb) * 1000
b2 = ((bb + 1000) - b)
HU = b2 * a1
MU = a2 * b2
HS = a1 * b1
MS = a2 * b2
MEN = HU + HS + MU + MS
HUT = (HU / MEN * T)
HST = (HS / MEN * T)
MUT = (MU / MEN * T)
MST = (HS / MEN * T)
aa = math.floor(aa / 1000)
bb = math.floor(bb / 1000)
charge = np.zeros((ny, nx))
charge[aa, bb+1] = HUT
charge[aa, bb] = HST
charge[aa+1, bb] = MST
charge[aa+1, bb+1] = MUT
# 粒子の計算結果をデータに保存
data["charge"] = charge
# 計算結果をファイルに保存
with open("{}.csv".format(name), "w") as f:
writer = csv.writer(f)
writer.writerows(charge)
# 計算結果の出力(任意の方法で)
for name, data in particle_data.items():
charge = data["charge"]
print("Particle Name:", name)
print("Charge Matrix:")
print(charge)
# 電荷置
eps0 = 1
# 電位計算 外挿?
@jit(nopython=True)
def calc_checkered(offset, ph):
for i in range(nx):
for j in range((i + offset) % 2, ny, 2):
if i == 0 or i == nx - 1 or j == 0 or j == ny - 1:
ph[i, j] = 0
else:
ph[i, j] = ph[i, j] + omega_SOR_recta * (
(ph[i+1, j] + ph[i-1, j] + ph[i, j+1]
+ ph[i, j-1])/4- ph[i, j]
- (dx*dy/(4*eps0))*charge[i, j])
conv_check = []
while delta > convergence_criterion:
phi_in = phi.copy()
if n_inter % 50 == 0:
print("interaction No =", n_inter, "delta=", delta)
conv_check.append([n_inter, delta])
calc_checkered(0, phi) # red group
calc_checkered(1, phi) # black group
delta = np.max(abs(phi - phi_in))
n_inter += 1
#red-black-SOR法
# 電場計算 内挿?
# ファイルを書き込みモードでオープン
body = (phi)
header = ('Electrical potential of 2D array')
with open('Deni1.csv', 'w') as f:
writer = csv.writer(f) # writerオブジェクトを作成
writer.writerows(header) # ヘッダーを書き込む
writer.writerows(body) # 内容を書き込む
body2 = (x)
header2 = ('Value1 of axis')
with open('jikunoatai1.csv', 'w') as f1:
writer = csv.writer(f1)
writer.writerow(header2)
writer.writerow(body2)
def plot2D(x, y, phi):
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid (x , y)
surf = ax.plot_surface(X, Y, phi[:], rstride=1, cstride=1,
cmap=cm.viridis,linewidth=0, antialiased=False)
ax.view_init(30, 225)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('phi')
plot2D(x, y, phi)
plt.show()
#オブジェクトを生成
fig = plt.figure(figsize=(11,7), dpi=100)
fig
xmin = 0
xmax = 1000
ymin = 0
ymax = 1000
x = np.linspace(xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)
X, Y = np.meshgrid(x, y)
plt.contourf(X, Y, phi, alpha=0.5, cmap=cm.viridis)
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
#電位計算
print("Keisan V")
for i in range(1,6):
print("x=",i*delta_L,"V=",phi[500+i,500]
- phi[500+int(1/delta_L),500])
#電場計算
L = 999
Ey = np.zeros([L,L])
for i in range(1, L-1):
for j in range(1, L-1):
Ey[i,j] = -(phi[i+1,j]-phi[i-1,j])/2/delta_L
print("Electric Field Calculation")
print("Keisan E")
for i in range(1,6):
print("x=",i*delta_L,"E=",Ey[500+i,500])
# 電場の可視化
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(x[1:L-1], y[1:L-1])
surf = ax.plot_surface(X, Y, Ey[1:L-1, 1:L-1], rstride=1, cstride=1,
cmap=cm.viridis, linewidth=0, antialiased=False)
ax.view_init(30, 225)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Electric Field (Ey)')
plt.show()