ただのメモ

他の人に見せても良い方のメモ

山と谷

 

Phase field法なる手法を使えば、2相の境界を上手い具合に図示できるらしい。

それを利用して、A0相とA1相がある状況と、B0相とB1相がある状況を考えて、A0相からA1相に境界近辺で変化しつつ、B0相からB1相でもその境界付近で変化しつつ、B1相はA1相とは共存しがたいため、B0相になっていく。ということを表現できそう。

コードだが、はてなブログの仕様によって、記事に※が入ってしまう問題がある。これの解決策がわからない。編集時点では、そんなことはないのだが、、、。

”””

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import math

nx,ny=32,32;dx,dy=0.5e-6,0.5e-6;
eee=1.0e+6;sigma=1.0;delta=4.*dx;amobi=4.e-14
ram=0.1;bbb=2.*np.log*1/(1.-(1.-2.*ram)))/2. 

aaa=np.sqrt(3.*delta*sigma/bbb);www=6.*sigma*bbb/delta;
pmobi=amobi*math.sqrt(2.*www)/(6.*aaa)

dt=dx*dx/(5.*pmobi*aaa*aaa)/2;nsteps=500

p= np.zeros*2;q=np.zeros*3
r_nuclei=5.*dx
for i in range(0,nx):
    for j in range(0,ny):
        rp = np.sqrt( (i *dx)**2 +(j*dy)**2 ) - r_nuclei
        rq1=np.sqrt( (i *dx-30.*dx)**2 +(j*dy-20.*dy)**2 ) - 0.3*r_nuclei
        rq2=np.sqrt( (i *dx-16.*dx)**2 +(j*dy-30.*dy)**2 ) - 0.3*r_nuclei
        rq3=np.sqrt( (i *dx-10.*dx)**2 +(j*dy-4.*dy)**2 ) - 0.3*r_nuclei
        rq4=np.sqrt( (i *dx-32.*dx)**2 +(j*dy-1.*dy)**2 ) - 0.3*r_nuclei
        p[0,i,j] = 0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*rp))
        q[0,i,j] = max(0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*rq1)),0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*rq2)),0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*rq3)),0.5*(1.-np.tanh(np.sqrt(2.*www)/(2.*aaa)*rq4)))

def do_timestep(p,q):
    for t in range(nsteps-1):
        p1=np.zeros*4;p2=np.zeros*5;p3=np.zeros*6;p4=np.zeros*7;
        q1=np.zeros*8;q2=np.zeros*9;q3=np.zeros*10;q4=np.zeros*11;
        for j in range(ny):
            for i in range(nx):
                ip=i+1;im=i-1;jp=j+1;jm=j-1;
                if ip > nx - 1:
                    ip = nx -1
                if im < 0:
                    im = 0
                if jp > ny - 1:
                    jp = ny -1
                if jm < 0:
                    jm = 0
#                p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+  aaa*aaa**12+  aaa*aaa**13+  aaa*aaa**14*(p[t,i,j]+p1[i,j]/2-0.5+3*(p[t,i,j]+p1[i,j]/2)*(1-(p[t,i,j]+p1[i,j]/2)))+  aaa*aaa**15*(q[t,i,j]+q1[i,j]/2-0.5-6*(p[t,i,j]+p1[i,j]/2-q[t,i,j]+q1[i,j]/2)*(q[t,i,j]+q1[i,j]/2)*(1-(q[t,i,j]+q1[i,j]/2)))+  aaa*aaa**16*(p[t,i,j]+p2[i,j]/2-0.5+3*(p[t,i,j]+p2[i,j]/2)*(1-(p[t,i,j]+p2[i,j]/2)))+  aaa*aaa**17*(q[t,i,j]+q2[i,j]/2-0.5-6*(p[t,i,j]+p2[i,j]/2-q[t,i,j]+q2[i,j]/2)*(q[t,i,j]+q2[i,j]/2)*(1-(q[t,i,j]+q2[i,j]/2)))+  aaa*aaa**18*(p[t,i,j]+p3[i,j]-0.5+3*(p[t,i,j]+p3[i,j])*(1-(p[t,i,j]+p3[i,j])))+  aaa*aaa**19*(q[t,i,j]+q3[i,j]-0.5-6*(p[t,i,j]+p3[i,j]-q[t,i,j]+q3[i,j])*(q[t,i,j]+q3[i,j])*(1-(q[t,i,j]+q3[i,j])))+  aaa*aaa*((q[t,ip,j]+q3[ip,j] - 2*(q[t,i,j]+q3[i,j]) + q[t,im,j]+q3[im,j])/dx/dx + (q[t,i,jp]+q3[i,jp] - 2*(q[t,i,j]+q3[i,j])+ q[t,i,jm]+q3[i,jm])/dy/dy) ) * dt

        for j in range(ny):
            for i in range(nx):
#                p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+  aaa*aaa*((p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                p[t+1,i,j] = p[t,i,j] +(p1[i,j]+p2[i,j]*2+p3[i,j]*2+p4[i,j])/6
                q[t+1,i,j] = q[t,i,j] +(q1[i,j]+q2[i,j]*2+q3[i,j]*2+q4[i,j])/6

do_timestep(p,q)

x = np.linspace(0, nx, nx);y=np.linspace(0, ny, ny)
x, y = np.meshgrid(y, x)

fig = plt.figure()
fig.set_dpi(100)
ax = Axes3D(fig)

def animate(i):
    ax.clear();
    plt.ylim([0,ny]);plt.xlim([0,nx])
    plt.xlabel('x direction');plt.ylabel('y direction')
    ax.plot_surface(x, y, p[i,:,:], rstride=1, cstride=1, cmap=plt.cm.Greens,alpha=0.5,vmax=1,vmin=-1)
    ax.plot_surface(x, y, q[i,:,:], rstride=1, cstride=1, cmap=plt.cm.Reds,alpha=0.5,vmax=1,vmin=-1)
    ax.set_zlim(0, 1)

anim = animation.FuncAnimation(fig,animate,frames=nsteps-1,interval=10,repeat=True)
plt.show()

”””

 

*1:1.+(1.-2.*ram

*2:nsteps,nx,ny

*3:nsteps,nx,ny

*4:nx,ny

*5:nx,ny

*6:nx,ny

*7:nx,ny

*8:nx,ny

*9:nx,ny

*10:nx,ny

*11:nx,ny

*12:p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                p1[i,j] = pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3*p[t,i,j]*(1-p[t,i,j]

*13:p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                q1[i,j] = pmobi * ( 4.*www*q[t,i,j]*(1.-q[t,i,j])*(q[t,i,j]-0.5-6*(p[t,i,j]-q[t,i,j])*q[t,i,j]*(1-q[t,i,j]

*14:q[t,ip,j] - 2*q[t,i,j] + q[t,im,j])/dx/dx + (q[t,i,jp] - 2*q[t,i,j] + q[t,i,jm])/dy/dy) ) * dt

        for j in range(ny):
            for i in range(nx):
                ip=i+1;im=i-1;jp=j+1;jm=j-1;
                if ip > nx - 1:
                    ip = nx -1
                if im < 0:
                    im = 0
                if jp > ny - 1:
                    jp = ny -1
                if jm < 0:
                    jm = 0
#                p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+  aaa*aaa*((p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                p2[i,j] = pmobi * ( 4.*www*(p[t,i,j]+p1[i,j]/2)*(1.-(p[t,i,j]+p1[i,j]/2

*15:p[t,ip,j]+p1[ip,j]/2 - 2*(p[t,i,j]+p1[i,j]/2) + p[t,im,j]+p1[im,j]/2)/dx/dx + (p[t,i,jp]+p1[i,jp]/2 - 2*(p[t,i,j]+p1[i,j]/2)+ p[t,i,jm]+p1[i,jm]/2)/dy/dy) ) * dt
                q2[i,j] = pmobi * ( 4.*www*(q[t,i,j]+q1[i,j]/2)*(1.-(q[t,i,j]+q1[i,j]/2

*16:q[t,ip,j]+q1[ip,j]/2 - 2*(q[t,i,j]+q1[i,j]/2) + q[t,im,j]+q1[im,j]/2)/dx/dx + (q[t,i,jp]+q1[i,jp]/2 - 2*(q[t,i,j]+q1[i,j]/2)+ q[t,i,jm]+q1[i,jm]/2)/dy/dy) ) * dt


        for j in range(ny):
            for i in range(nx):
                ip=i+1;im=i-1;jp=j+1;jm=j-1
                if ip > nx - 1:
                    ip = nx -1
                if im < 0:
                    im = 0
                if jp > ny - 1:
                    jp = ny -1
                if jm < 0:
                    jm = 0
#                p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+  aaa*aaa*((p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                p3[i,j] = pmobi * ( 4.*www*(p[t,i,j]+p2[i,j]/2)*(1.-(p[t,i,j]+p2[i,j]/2

*17:p[t,ip,j]+p2[ip,j]/2 - 2*(p[t,i,j]+p2[i,j]/2) + p[t,im,j]+p2[im,j]/2)/dx/dx + (p[t,i,jp]+p2[i,jp]/2 - 2*(p[t,i,j]+p2[i,j]/2)+ p[t,i,jm]+p2[i,jm]/2)/dy/dy) ) * dt
                q3[i,j] = pmobi * ( 4.*www*(q[t,i,j]+q2[i,j]/2)*(1.-(q[t,i,j]+q2[i,j]/2

*18:q[t,ip,j]+q2[ip,j]/2 - 2*(q[t,i,j]+q2[i,j]/2) + q[t,im,j]+q2[im,j]/2)/dx/dx + (q[t,i,jp]+q2[i,jp]/2 - 2*(q[t,i,j]+q2[i,j]/2)+ q[t,i,jm]+q2[i,jm]/2)/dy/dy) ) * dt


        for j in range(ny):
            for i in range(nx):
                ip=i+1;im=i-1;jp=j+1;jm=j-1;
                if ip > nx - 1:
                    ip = nx -1
                if im < 0:
                    im = 0
                if jp > ny - 1:
                    jp = ny -1
                if jm < 0:
                    jm = 0
#                p[t+1,i,j] = p[t,i,j] + pmobi * ( 4.*www*p[t,i,j]*(1.-p[t,i,j])*(p[t,i,j]-0.5+3./(2.*www)*eee)+  aaa*aaa*((p[t,ip,j] - 2*p[t,i,j] + p[t,im,j])/dx/dx + (p[t,i,jp] - 2*p[t,i,j] + p[t,i,jm])/dy/dy) ) * dt
                p4[i,j]=pmobi*(4.*www*(p[t,i,j]+p3[i,j])*(1.-(p[t,i,j]+p3[i,j]

*19:p[t,ip,j]+p3[ip,j] - 2*(p[t,i,j]+p3[i,j]) + p[t,im,j]+p3[im,j])/dx/dx + (p[t,i,jp]+p3[i,jp] - 2*(p[t,i,j]+p3[i,j])+ p[t,i,jm]+p3[i,jm])/dy/dy) ) * dt
                q4[i,j]=pmobi*(4.*www*(q[t,i,j]+q3[i,j])*(1.-(q[t,i,j]+q3[i,j]