Nov
26
It means to compute 2d diffusion equation just like previous post in polar/cylindrical coordinate, and all went to wrong direction, :)
Still trying to understand matplotlib mplot3d behavior
import scipy as sp from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm from matplotlib.ticker import LinearLocator, FormatStrFormatter import matplotlib.pyplot as plt import mpl_toolkits.mplot3d.axes3d as p3 import matplotlib.animation as animation #dr = .1 #dp = .1 #nr = int(1/dr) #np = int(2*sp.pi/dp) nr = 10 np = 10 dr = 1./nr dp = 2*sp.pi/np a = .5 tmax = 100 t = 0. dr2 = dr**2 dp2 = dp**2 dt = dr2 * dp2 / (2 * a * (dr2 + dp2) ) dt /=10. print dt ut = sp.zeros([nr,np]) u0 = sp.zeros([nr,np]) ur = sp.zeros([nr,np]) ur2 = sp.zeros([nr,np]) r = sp.arange(0.,1.,dr) p = sp.arange(0.,2*sp.pi,dp) #initial for i in range(nr): for j in range(np): if ( (i>(2*nr/5.)) & (i<(3.*nr/3.)) ): u0[i,j] = 1. #print u0 def hitung_ut(ut,u0): for i in sp.arange (len(r)): if r[i]!= 0.: ur[i,:] = u0[i,:]/r[i] ur2[i,:] = u0[i,:]/(r[i]**2) ut[1:-1, 1:-1] = u0[1:-1, 1:-1] + a*dt*( (ur[1:-1, 1:-1] - ur[:-2, 1:-1])/dr+ (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2,1:-1])/dr2+ (ur2[1:-1, 2:] - 2*ur2[1:-1, 1:-1] + ur2[1:-1, :-2])/dp2) #hitung_ut(ut,u0) #print ut def data_gen(framenumber, Z ,surf): global ut global u0 hitung_ut(ut,u0) u0[:] = ut[:] Z = u0 ax.clear() plotset() surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False, alpha=0.7) return surf, fig = plt.figure() #ax = fig.gca(projection='3d') ax = fig.add_subplot(111, projection='3d') R = sp.arange(0,1,dr) P = sp.arange(0,2*sp.pi,dp) R,P = sp.meshgrid(R,P) X,Y = R*sp.cos(P),R*sp.sin(P) Z = u0 print len(R), len(P) def plotset(): ax.set_xlim3d(-1., 1.) ax.set_ylim3d(-1., 1.) ax.set_zlim3d(-1.,1.) ax.set_autoscalez_on(False) ax.zaxis.set_major_locator(LinearLocator(10)) ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f')) cset = ax.contour(X, Y, Z, zdir='x', offset=0. , cmap=cm.coolwarm) cset = ax.contour(X, Y, Z, zdir='y', offset=1. , cmap=cm.coolwarm) cset = ax.contour(X, Y, Z, zdir='z', offset=-1., cmap=cm.coolwarm) plotset() surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False, alpha=0.7) fig.colorbar(surf, shrink=0.5, aspect=5) ani = animation.FuncAnimation(fig, data_gen, fargs=(Z, surf),frames=500, interval=30, blit=False) #ani.save('2dDiffusionf500b512.mp4', bitrate=512) plt.show()
.
