ポアソン方程式を計算するのに向上するためred-black-SOR法を用いて計算を行おうとしたところ偶数奇数に分け計算を行うことはできましたがスレッディングやプロセスプールを用いたところエラーは出ませんが計算結果が0として出力されます。何か情報等でも構いません、教えていただけると幸いです
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import itertools
import csv
import time
import concurrent.futures
import threading
import random
# 初期条件
delta_L = 1.0
delta = 1.0
n_inter = 0
nx = 100 #xの範囲
ny = 100 #yの範囲
xmin = 0 #xの最小値
xmax = 100*delta #xの最大値
ymin = 0 #yの最小値
ymax = 100*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-10 #前回との差
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)
# 電荷設置
eps0=1
charge= np.zeros((ny, nx))
Q1=1
Q2=Q1/delta**2
charge[50,50] = Q2
#Main
conv_check = []
class SampleThreading(threading.Thread):
def __init__(self, thread_name):
self.thread_name = str(thread_name)
threading.Thread.__init__(self)
def __str__(self):
return self.thread_name
def run(self,delta,n_inter):
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])
for i in range(nx - 1):
for j in range(i % 2, ny - 1, 2):
if i == xmin or i == xmax or j == ymax or j == ymin:
phi[i, j] = 0
else:
phi[i,j] = phi[i,j]+omega_SOR_recta *((phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4-phi[i,j]+ (dx*dy/(4*eps0))*charge[i,j])
delta = np.max(abs(phi - phi_in))
n_inter += 1
sleep_seconds = random.randint(5, 10)
time.sleep(sleep_seconds)
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])
for i in range(nx-1):
for j in range(i % 2 + 1, ny - 1, 2):
if i == xmin or i == xmax or j == ymax or j == ymin:
phi[i, j] = 0
else:
phi[i,j] = phi[i,j]+omega_SOR_recta *((phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4-phi[i,j]+ (dx*dy/(4*eps0))*charge[i,j])
delta = np.max(abs(phi - phi_in))
n_inter += 1
thread_list = []
for i in range(nx - 1):
thread = SampleThreading(thread_name=i)
thread.start()
thread_list.append(thread)
for thread in thread_list:
thread.join()
# CSV出力
print("The number of total iteration =", n_inter)
print("data_points=", nx*ny)
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.gca(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 = 100
ymin = 0
ymax = 100
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[50+i,50]-phi[50+int(1/delta_L),50])
#電場計算
L = 99
Ey = np.zeros([L,L])
for i in range(L):
for j in range(L):
Ey[i,j] = -(phi[i+1,j]-phi[i-1,j])/2/delta_L
print("Keisan E")
for i in range(1,6):
print("x=",i*delta_L,"E=",Ey[50+i,50])