0

ポアソン方程式を計算するのに向上するため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])
3
  • コード中のコメントが一部文字化けしているようです。読めるように修正いただけないでしょうか?
    – merino
    Commented 2022年11月14日 3:25
  • 修正いたしました Commented 2022年11月14日 4:17
  • red-black-SOR法は知らないので細かいロジックはおってないけど if n_inter % 50 == 0: の中の for i in range(nx - 1): のインデントが怪しいような。。。 50回に1回しかこの部分通らないので phi_in = phi.copy() => delta = np.max(abs(phi - phi_in)) は残り49回入力が変化しないのに毎回同じ結果になる無駄な計算してるようにしか見えないので違和感あります。その下の while 文もphi_in = phi.copy()しかやってないのでループ脱出できるように見えないし、これインデントあってますか? 見当違いなこと言ってたらすみません
    – user35303
    Commented 2022年11月14日 6:25

1 件の回答 1

2

現在のコードからは i(x) の値毎にスレッドを作り並行処理させる方針が読み取れたので,その方針で一例を示します。
最初に,核となる計算部分(# main# CSV出力 の間)を赤グループと黒グループが交互に計算されるようにシングルスレッドで記述します。

# main
def calc_checkered(offset):
    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:
                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])

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)  # red group
    calc_checkered(1)  # black group
    delta = np.max(abs(phi - phi_in))
    n_inter += 1

# CSV出力

次に,赤グループと黒グループを交互に計算することでグループ内では計算順序の制約がなくなったことを利用し,i(x) の値毎にスレッドを作り並行処理するように記述を変更します。
一応,シングルスレッド版とマルチスレッド版で収束までの反復回数と最後の結果に差がないことは確認しました。
なお,マルチスレッド版でコメントアウトの3カ所を元に戻すと(スリープでパフォーマンスは落ちますが)マルチスレッド化していることが確認できます。

# main
def calc_alternate(i, offset):
    for j in range((i + offset) % 2, ny, 2):
        if i == 0 or i == nx - 1 or j == 0 or j == ny - 1:
            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])
    # sleep_seconds = random.randint(5, 10) / 100
    # time.sleep(sleep_seconds)

def calc_checkered(offset):
    thread_list = []
    for i in range(nx):
        thread = threading.Thread(target=calc_alternate, args=(i, offset))
        thread.start()
        thread_list.append(thread)
    # print(f'threading.active_count: {threading.active_count()}')
    for thread in thread_list:
        thread.join()

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)  # red group
    calc_checkered(1)  # black group
    delta = np.max(abs(phi - phi_in))
    n_inter += 1

# CSV出力

(追記,御参考)
計算時間の短縮のため numba モジュールを導入(pip install numba)する例を示します。私の環境では,この計算部分が5分の1程度に短縮されました。

# Main
from numba import jit  # Better to write at the beginning of the code

@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

# CSV出力
2
  • 丁寧な回答ありがとうございます。上記のコードを使用しマルチスレッドで計算を行うことができました。しかし、今回は処理速度の向上を目的としてマルチスレッドを取り入れましたがシングルスレッドで計算を行った時とほぼ変わらないかシングルスレッドのほうが処理速度が早くなってしまいます。処理速度を早くするためにマルチスレッドを取り入れるのは間違いなのでしょうか?調べてみるとProcessPoolExecutor()が良いという人もいますが取り入れることは可能でしょうか。以前ProcessPoolExecutor()を使おうとした時は今回の質問と同じく計算結果が0として出力されるので断念しました。。。 Commented 2022年11月24日 12:02
  • 現行の Python では,今回のような(外部とのやりとりによる待ちがなく)ひたすら計算するケースでは「処理速度を速くするためにマルチスレッドを取り入れる」は得策ではないと思われます。理由は「Python GIL」を検索して頂くと詳しい説明が見つけられると思います。次に,マルチプロセス化ですが(プロセス間のメモリ共有の)経験がなく可否判断は難しいです。そこで,現行の計算時間の主要因と思われる「Python の多重 for ループの遅さ」への対策に numba モジュールを導入する例を回答に追記します。
    – Delft View
    Commented 2022年11月24日 15:19

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

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。