Nugroho's blog.: The Wrong Code Will often Provide Beautiful Result, :)

Thursday, November 26, 2015

The Wrong Code Will often Provide Beautiful Result, :)

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

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')

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()

.