1

わからない点
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()
2
  • 質問文の 書式の指定は Markdown や HTML を使用して投稿を書式設定するには?にあるので参考にどうぞ。行が溢れて(?)コードの一部が次の行に回り込んでたようなのである程度は対処してみたけど, 一度質問者さんから内容確認してみてください
    – oriri
    11月8日 6:41
  • 当サイトはあくまでプログラミングのQAサイトなので、「プラズマ生成」が主軸な質問文だと伝わらない可能性があります。
    – cubick
    11月8日 9:55

0

この質問に回答するには、ログインする必要があります。

のタグが付いた他の質問を参照する。