Update files

Author: Jordi

…onvection/ConvectionDiffusion1D_limiters …ection/ConvectionDiffusion1D_limiters.pyday_8/Convection/ConvectionDiffusion1D_limiters renamed to day_8/Convection/ConvectionDiffusion1D_limiters.py day_8/Convection/ConvectionDiffusion1D_limiters renamed to day_8/Convection/ConvectionDiffusion1D_limiters.py

@@ -31,7 +31,7 @@
 beta = 1
 
 "Number of points"
-N = 256
+N = 32
 Dx = 1/N
 x = np.linspace(0,1,N+1)
 

day_8/Convection/PollutantTransport1D day_8/Convection/PollutantTransport1D.pyday_8/Convection/PollutantTransport1D renamed to day_8/Convection/PollutantTransport1D.py

@@ -50,22 +50,22 @@
 nt = np.size(time)
 
 "Initialize solution variable"
-U = np.zeros((N+2,nt))
+U = np.zeros((N+1,nt))
 
 
 for it in range(nt-1):
 
     "System matrix and RHS term"
     "Diffusion term"
-    Diff = nu*(1/Dx**2)*(2*np.diag(np.ones(N+2)) - np.diag(np.ones(N+1),-1) - np.diag(np.ones(N+1),1))
+    Diff = nu*(1/Dx**2)*(2*np.diag(np.ones(N+1)) - np.diag(np.ones(N),-1) - np.diag(np.ones(N),1))
 
     "Advection term:"
 
     "Sensor"
     U0 = U[:,it]
-    uaux = np.concatenate(([U0[0]],U0,[U0[N+1]]))
-    Du = uaux[1:N+4] - uaux[0:N+3] + 1e-8
-    r = Du[0:N+2]/Du[1:N+3]
+    uaux = np.concatenate(([U0[0]], U0,[U0[N]]))
+    Du = uaux[1:N+3] - uaux[0:N+2] + 1e-8
+    r = Du[0:N+1]/Du[1:N+2]
 
     "Limiter"
     if beta>0:
@@ -74,8 +74,8 @@
     else:
         phi = 2*r/(r**2 + 1)
 
-    phim = phi[0:N+1]
-    phip = phi[1:N+2]
+    phim = phi[0:N]
+    phip = phi[1:N+1]
 
 
     "Upwind scheme"
@@ -97,25 +97,25 @@
     "Source term"
     sine = np.sin(2*pi*time[it+1])
     sineplus = 0.5*(sine + np.abs(sine))
-    F = 100*np.exp(-((xN-0.8)/0.01)**2)*sineplus
+    F = 100*np.exp(-((x-0.8)/0.01)**2)*sineplus
 
     "Temporal terms"
-    A = A + (1/dt)*np.diag(np.ones(N+2))
+    A = A + (1/dt)*np.diag(np.ones(N+1))
     F = F + U0/dt
 
     "Boundary condition at x=0"
-    A[0,:] = (1/Dx)*np.concatenate(([-1, 0, 1],np.zeros(N-1)))
+    A[0,:] = (1/Dx)*np.concatenate(([1.5, -2, 0.5],np.zeros(N-2)))
     F[0] = 0
 
     "Boundary condition at x=1"
-    A[N+1,:] = np.concatenate((np.zeros(N+1),[1]))
-    F[N+1] = u0
+    A[N,:] = np.concatenate((np.zeros(N),[1]))
+    F[N] = u0
 
 
     "Solution of the linear system AU=F"
     u = np.linalg.solve(A,F)
     U[:,it+1] = u
-    u = u[1:N+2]
+    u = u[0:N+1]
 
 
 "Animation of the results"
@@ -128,17 +128,18 @@
 
 def animate(i):
 
-    u = U[1:N+2,i]
+    u = U[0:N+1,i]
     plt.plot(x,u)
     myAnimation.set_data(x, u)
     return myAnimation,
 
 anim = animation.FuncAnimation(fig,animate,frames=range(1,nt),blit=True,repeat=False)
 
 
-"Peclet number"
-P = np.abs(c*Dx/nu)
-print("Pe number Pe=%g\n" % P);
+if nu>0:
+    "Peclet number"
+    P = np.abs(c*Dx/nu)
+    print("Pe number Pe=%g\n" % P);
 
 "CFL number"
 CFL = np.abs(c*dt/Dx)