From 51637605b7bd4912d078919febf740b807b65de1 Mon Sep 17 00:00:00 2001
From: Anne de Jong <anne@amdj.nl>
Date: Thu, 29 Jan 2015 07:36:01 +0100
Subject: [PATCH] Parallel plate flow also working

---
 blflow.py       |  5 ++--
 parplateflow.py | 74 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 76 insertions(+), 3 deletions(-)
 create mode 100755 parplateflow.py

diff --git a/blflow.py b/blflow.py
index 8798c65..491f56f 100755
--- a/blflow.py
+++ b/blflow.py
@@ -35,8 +35,7 @@ def u_np1(un,tn,dt):
     Kn=K(tn)
     unp1=un
     unp1[0]=0                             #Velocity zero ver here
-    for i in range(1,un.size-1):
-        unp1[i]=dt*Kn+un[i]+l*(un[i-1]-2*un[i]+un[i+1])
+    unp1[1:-1]=un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:])
     unp1[-1]=unp1[-2]                     #Approximate 'infinity' bc
     return unp1
 
@@ -61,7 +60,7 @@ while(True):
     t+=dt
     uold=un
     un=u_np1(uold,t,dt)
-    if(i%20==0):
+    if(i%100==0):
         linefd.set_xdata(un)
         linee.set_xdata(u_ex(t))
         p.draw()
diff --git a/parplateflow.py b/parplateflow.py
new file mode 100755
index 0000000..5897ae5
--- /dev/null
+++ b/parplateflow.py
@@ -0,0 +1,74 @@
+#!/usr/bin/python
+
+# Boundary layer flow
+from numpy import *
+import time 
+import matplotlib
+matplotlib.use('TkAgg')
+# from matplotlib.pylab import *
+import pylab as p
+
+
+
+# import matplotlib.animation as animation
+def K(t):                                 #Forcing function
+    return (1-exp(-0.1*t))*cos(t)
+
+
+s=10
+#Define domain
+n=50                                      #Number of gridpoints
+y=linspace(-1,1,n)
+
+dy=y[1]-y[0]
+dt=0.0005
+
+l=(dt/(s**2*dy**2))
+
+hnu=cosh(sqrt(1j)*s*y)/cosh(sqrt(1j)*s)
+fnu=tanh(sqrt(1j)*s)/(sqrt(1j)*s)
+
+
+def u_ex(tn):
+    return (((1-hnu)/(1-fnu))*exp(1j*(tn))/1j).real
+    
+
+def u_np1(un,tn,dt):
+    Kn=K(tn)
+    unp1=un
+    unp1[0]=0                                  #Velocity zero ver here
+    unp1[1:-1]=un[1:-1]+dt*Kn+l*(un[0:-2]-2*un[1:-1]+un[2:])
+    unp1[-1]=0                     #Boundary other side
+    return unp1
+
+un0=zeros(n,float)
+t=0
+un=un0
+# un.append(un0)
+
+# Make the plot
+p.ion()
+linefd, = p.plot(un0,y)
+linee, = p.plot(un0,y)
+p.legend(('Finite difference','Periodic exact'))
+p.ylim(-1,1)
+p.xlim(-1.5,1.5)
+p.ylabel('y')
+p.xlabel('u')
+p.grid('on')
+i=0
+uold=un
+while(True):
+    t+=dt
+    uold=un
+    un=u_np1(uold,t,dt)
+    if(i%100==0):
+        linefd.set_xdata(un)
+        linee.set_xdata(u_ex(t))
+        p.draw()
+        # print("Time:",t)
+    i+=1
+
+
+
+