<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Równanie dyfuzji - dekompozycja siatki z zastosowaniem MPI
Problem
Step1:  non-contigous slice
Step2:  N - slices
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os 
print os.getenv("HOME")
wd = os.path.join( os.getenv("HOME"),"mpi_tmpdir")
if not os.path.isdir(wd):
    os.mkdir(wd)
os.chdir(wd)
print "WD is now:",os.getcwd()
%%writefile mpi002.py
from mpi4py import MPI
import numpy as np 
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
def numpy_diff2d(u,dx2,dy2,c):
    A = (1.0-2.0*(c/dx2+c/dy2))
    u[1:-1,1:-1] =A*u[1:-1,1:-1] + c/dy2*(u[2:,1:-1] + u[:-2,1:-1]) + \
                                   c/dx2*(u[1:-1,2:] + u[1:-1,:-2])
        
N=52
Niter=211
dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy
dt = 0.01
D = 0.1
c = D*dt
u = np.zeros([N, N])       
if rank == 0:
    u[-2,u.shape[1]/2] = 1.0/np.sqrt(dx2*dy2)
    print "CLF = ",c/dx2,c/dy2
    
for i in range(Niter):
    
    if rank == 0:
                        
        comm.Send([u[-2,:], MPI.FLOAT], dest=1)
        comm.Recv([u[-1,:], MPI.FLOAT], source=1)     
        
    elif rank == 1:               
        comm.Recv([u[0,:], MPI.FLOAT], source=0)     
        comm.Send([u[1,:], MPI.FLOAT], dest=0)
       
    numpy_diff2d(u,dx2,dy2,c)
        
#np.savez("udata%04d"%rank, u=u)    
U = comm.gather(u[1:-1,1:-1])
if rank==0:
    np.savez("Udata", U=U)    
    
!mpirun -n 2 python mpi002.py
data = np.load("Udata.npz")
plt.imshow(np.vstack(data['U']))
print data['U'].shape
!pwd
Explanation: Równanie dyfuzji - dekompozycja siatki z zastosowaniem MPI
Problem:
Chcemy rozwiązać równanie dyfuzji wykorzystując wiele procesorów. 
Dzielimy siatkę na obszary o w każdym obszarze rozwiązujemy równanie niezależnie. Po każdym kroku czasowym wykorzystujemy komunikację w MPI do wymiany informacji o przylegających  brzegach odpowiednich obszarów.
End of explanation
%%writefile mpi003.py
from mpi4py import MPI
import numpy as np 
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
def numpy_diff2d(u,dx2,dy2,c):
    A = (1.0-2.0*(c/dx2+c/dy2))
    u[1:-1,1:-1] =A*u[1:-1,1:-1] + c/dy2*(u[2:,1:-1] + u[:-2,1:-1]) + \
                                   c/dx2*(u[1:-1,2:] + u[1:-1,:-2])
        
N=52
Niter=211
dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy
dt = 0.01
D = 0.1
c = D*dt
u = np.zeros([N, N])       
if rank == 0:
    u[u.shape[1]/2,-2] = 1.0/np.sqrt(dx2*dy2)
    print "CLF = ",c/dx2,c/dy2
  
for i in range(Niter):
    
    if rank == 0:
        OUT = u[:,-2].copy()  
        IN = np.empty_like(OUT)
        comm.Send([OUT, MPI.FLOAT], dest=1)
        comm.Recv([IN, MPI.FLOAT], source=1)     
        u[:,-1] = IN
    elif rank == 1: 
        OUT = u[:,1].copy()
        IN = np.empty_like(OUT)               
        comm.Recv([IN, MPI.FLOAT], source=0)     
        comm.Send([OUT, MPI.FLOAT], dest=0)
        u[:,0] = IN
        
    numpy_diff2d(u,dx2,dy2,c)
        
np.savez("udata%04d"%rank, u=u)    
!mpirun -n 2 python mpi003.py
u1 = np.load('udata0000.npz')['u']
u2 = np.load('udata0001.npz')['u']
plt.imshow(np.hstack([u1[:,:-1],u2[:,1:]]))
Explanation: non-contigous slice
End of explanation
%%writefile mpi004.py
from mpi4py import MPI
  
import numpy as np 
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
Nproc = comm.size
def numpy_diff2d(u,dx2,dy2,c):
    A = (1.0-2.0*(c/dx2+c/dy2))
    u[1:-1,1:-1] = A*u[1:-1,1:-1] + c/dy2*(u[2:,1:-1] + u[:-2,1:-1]) + \
                                    c/dx2*(u[1:-1,2:] + u[1:-1,:-2])
N = 16*128        
Nx = N
Ny = N/Nproc
Niter=200
dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy
dt = 0.01
D = 0.2
c = D*dt
u = np.zeros([Ny, Nx])       
if rank == 0:
    u[-2,u.shape[1]/2] = 1.0/np.sqrt(dx2*dy2)
    print "CLF = ",c/dx2,c/dy2
t0 = MPI.Wtime()
for i in range(Niter):
    
    if Nproc>1:
        if rank == 0:       
            comm.Send([u[-2,:], MPI.FLOAT], dest=1)
        if rank >0  and rank < Nproc-1:
            comm.Recv([u[0,:], MPI.FLOAT], source=rank-1)
            comm.Send([u[-2,:], MPI.FLOAT], dest=rank+1)    
        if rank == Nproc - 1:               
            comm.Recv([u[0,:], MPI.FLOAT], source=Nproc-2)     
            comm.Send([u[1,:], MPI.FLOAT], dest=Nproc-2)
        if rank >0  and rank < Nproc-1:
            comm.Recv([u[-1,:], MPI.FLOAT], source=rank+1)
            comm.Send([u[1,:], MPI.FLOAT], dest=rank-1)
        if rank == 0:
            comm.Recv([u[-1,:], MPI.FLOAT], source=1)   
    
    #print rank
    
    comm.Barrier()
    numpy_diff2d(u,dx2,dy2,c)
t1 = MPI.Wtime()    
print rank,t1-t0
    
#np.savez("udata%04d"%rank, u=u)    
if Nproc>1:
    U = comm.gather(u[1:-1,1:-1])
    if rank==0:
         np.savez("Udata", U=U)    
    
!mpirun -H gpu2,gpu3 python mpi004.py
!mpirun -n 4 python mpi004.py
data = np.load("Udata.npz")
plt.imshow(np.vstack(data['U']))
print data['U'].shape
a = np.arange(0,16).reshape(4,4)
b = a[:,2]
c = a[2,:]
np.may_share_memory(a,b),np.may_share_memory(a,c)
a.flags
b.flags
c.flags
a=np.array(range(6))
b = a[2:4]
b=666
print a
np.may_share_memory?
Explanation: N - slices
End of explanation
<END_TASK>