92299


ポアソン方程式を計算するのに向上するためred-black-
SOR法を用いて計算を行おうとしたところ偶数奇数に分け計算を行うことはできましたがスレッディングやプロセスプールを用いたところエラーは出ませんが計算結果が０として出力されます。何か情報等でも構いません、教えていただけると幸いです

[code]

    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])
    
[/code]


CC BY-SA 4.0
2022-11-14T03:20:28.513
92212
2022-11-24T15:23:18.313
2022-11-14T04:16:55.687
54717
54717
question
python spyder
ポアソン方程式の並列化ができない


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

[code]

    # 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出力
    
[/code]

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

[code]

    # 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出力
    
[/code]

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

[code]

    # 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出力
    
[/code]


CC BY-SA 4.0
2022-11-18T10:39:14.150
92299
2022-11-24T15:23:18.313
2022-11-24T15:23:18.313
54588
54588
92212
answer