rhinoでチューリングパターン - single thread

turing_multi_thread.py

#!encoding: utf-8
import rhinoscriptsyntax as rs
import math as ma
import random

# 並列処理用ライブラリインポート
import threading

n=50
dt=0.5
h=0.1
h2=h*h

a=0.024
b=0.078
Du=0.002
Dv=0.001

u=[]
v=[]
u1=[]
v1=[]

def clear(n):
    u = []
    u1 = []
    v = []
    v1 = []
    row=[]
    for i in range(n+3):
        for j in range(n+3):
            row.append(1)
        u.append(row)
        u1.append(row)
        row=[]
    row=[]
    for i in range(n+3):
        for j in range(n+3):
            row.append(0)
        v.append(row)
        v1.append(row)
        row=[]
    return u,u1,v,v1

def init(u, v, n):
    x=random.randrange(1,n+2)
    y=random.randrange(1,n+2)
    for i in range(1,n+2):
        for j in range(1,n+2):
            r = ma.sqrt((x-i)**2+(y-j)**2)
            if r<8:
                u[i][j]=0.6+random.random()*0.06
                v[i][j]=0.2+random.random()*0.06
    return u, v

def boundary():
    for i in range(1,n+2):
        u[i][0]=u[i][n+1]
        u[i][n+2]=u[i][1]
        u[0][i]=u[n+1][i]
        u[n+2][i]=u[1][i]
    for i in range(1,n+2):
        v[i][0]=v[i][n+1]
        v[i][n+2]=v[i][1]
        v[0][i]=v[n+1][i]
        v[n+2][i]=v[1][i]

# 計算処理をcalcとして外に定義
def calcU(i, u, v, u1, n, h2, b, dt):
        for j in range(1,n+2):
            Lu=(u[i+1][j]+u[i][j+1]+u[i-1][j]+u[i][j-1]-4*u[i][j])/h2
            f=-u[i][j]*v[i][j]**2+a*(1-u[i][j])
            g=u[i][j]*v[i][j]**2-b*v[i][j]
            u1[i][j]=u[i][j]+(Du*Lu+f)*dt
            # return u1

def calcV(i, u, v, v1, n, h2, b, dt):
        for j in range(1,n+2):
            Lv=(v[i+1][j]+v[i][j+1]+v[i-1][j]+v[i][j-1]-4*v[i][j])/h2
            f=-u[i][j]*v[i][j]**2+a*(1-u[i][j])
            g=u[i][j]*v[i][j]**2-b*v[i][j]
            v1[i][j]=v[i][j]+(Dv*Lv+g)*dt
            # return v1

            
# updateはu1/v1 -> u/vへのdeep copyだけを行う
def update():
    for i in range(1,n+2):
        for j in range(1,n+2):
            u[i][j]=u1[i][j]
            v[i][j]=v1[i][j]

                
def display(u,xmin,xmax,ymin,ymax,p):
    domain=(xmin,xmax,ymin,ymax)
    xstep=(domain[1]-domain[0])/n
    ystep=(domain[3]-domain[2])/n
    verts=[]
    for i in range(1,n+2):
        x=domain[0]+(i-1)*xstep
        for j in range(1,n+2):
            y=domain[2]+(j-1)*ystep
            z=u[i][j]*p
            verts.append((x,y,z))
    faces=[]
    for i in range(n):
        for j in range(n):
            e=i*(n+1)+j
            faces.append((e,e+1,e+n+2,e+n+1))
    rs.AddMesh(verts,faces)

# u/u1/v/v1はclear内で初期化
u,u1,v,v1 = clear(n)
u, v = init(u, v, n)

for ite in range(200):
    boundary()
    
    # calcを並列化
    threads = []
    for i in range(1, n+2):
        thread = threading.Thread(target=calcU, args=([i, u, v, u1, n, h2, b, dt]),name="thread%d" % i)
        thread.start()
        thread.join()

    #        threads.append(thread)
        
        
    threads = []
    for i in range(1, n+2):
        thread = threading.Thread(target=calcV, args=([i, u, v, v1, n, h2, b, dt]),name="thread%d" % i)
        thread.start()
        thread.join()

    #        threads.append(thread)
    update()
display(u,-100,100,-100,100,10)