w3hello.com logo
Home PHP C# C++ Android Java Javascript Python IOS SQL HTML videos Categories
How to use a wrapper function that i generated using swig in python?

Okay so i was able to accomplish what i wanted to in Cython, thank god. In the process I found Cython much more powerful of a tool !

So the result is a 2D wave PDE solved using a finite difference. The main computation has been exported to a Cythonized function. The function takes in three 2D np.ndarrays and returns one back. Much easier to work with numpy and other data types as well.

Here is the Cythonized function:

from numpy import *
cimport numpy as np

def cwave_prop( np.ndarray[double,ndim=2] u_prev, np.ndarray[double,ndim=2]
u, np.ndarray[double,ndim=2] u_next):

    cdef double t = 0
    cdef double t_old = 0
    cdef double t_end = 100
    cdef int i,j
    cdef double c = 1
    cdef double Lx = 100
    cdef double Ly = 100
    cdef double dx = 1
    cdef double dy = 1
    cdef double dt = (1/(c))*(1/(sqrt(1/dx**2 + 1/dy**2)))

    while t<t_end:
        t_old=t; t +=dt

        #wave steps through time and space

        for i in xrange(1,99):
            for j in xrange(1,99):
                u_next[i,j] = - u_prev[i,j] + 2*u[i,j] + 
                        ((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + 
                ((c*dt/dx)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1])

        #set boundary conditions of grid to 0

        for j in xrange(100): u_next[0,j] = 0
        for i in xrange(100): u_next[i,0] = 0
        for j in xrange(100): u_next[Lx-1,j] = 0
        for i in xrange(100): u_next[i,Ly-1] = 0

        #set prev time step equal to current one
        for i in xrange(100):
            for j in xrange(100):       
                u_prev[i,j] = u[i,j]; 
                u[i,j] = u_next[i,j]; 



    print u_next

And here is my python script that calls it and also returns it back then plots the result. Any suggestions on writing better code are welcome...

'''George Lees Jr.
2D Wave pde '''

from numpy import *
import numpy as np
import matplotlib.pyplot as plt 
import cwave2
np.set_printoptions(threshold=np.nan)

#declare variables
#need 3 arrays u_prev is for previous time step due to time derivative

Lx=Ly = (100)                       #Length of x and y dims of grid
dx=dy = 1                       #derivative of x and y respectively
x=y = np.array(xrange(Lx))              #linspace to set the initial
condition of wave
u_prev=np.ndarray(shape=(Lx,Ly), dtype=np.double)   #u_prev 2D grid for
previous time step needed bc of time derivative
u=np.ndarray(shape=(Lx,Ly), dtype=np.double)        #u 2D grid
u_next=np.ndarray(shape=(Lx,Ly), dtype=np.double)   #u_next for advancing
the time step #also these are all numpy ndarrays
c = 1                           #setting constant velocity of the wave
dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2))       #we have to set dt
specifically to this or numerical approx will fail!
print dt

#set Initial Conditions and Boundary Points
#I(x) is initial shape of the wave

def I(x,y): return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0)

#set up initial wave shape

for i in xrange(100):
    for j in xrange(100):
        u[i,j] = I(x[i],y[j])


#set up previous time step array

for i in xrange(1,99):
    for j in xrange(1,99):
            u_prev[i,j] = u[i,j] + 0.5*((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j]
+ u[i+1,j]) + 
            0.5*((c*dt/dy)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) 

#set boundary conditions to 0

for j in xrange(100): u_prev[0,j] = 0
for i in xrange(100): u_prev[i,0] = 0
for j in xrange(100): u_prev[Lx-1,j] = 0
for i in xrange(100): u_prev[i,Ly-1] = 0

#call C function from Python
cwave2.cwave_prop( u_prev , u , u_next )
#returned u (2D np.ndarray)

from tempfile import TemporaryFile
outfile = TemporaryFile()
np.save(outfile,u)
fig = plt.figure()
plt.imshow(u,cmap=plt.cm.ocean)
plt.colorbar()
plt.show()




© Copyright 2018 w3hello.com Publishing Limited. All rights reserved.