Diffusion and Advection

Materials for Diffusion and Advection classes on days 7 and 8.

Author: Jordi

day_7/Diffusion/Poisson1D

@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+"""
+Solution of a 1D Poisson equation: -u_xx = f
+Domain: [0,1]
+BC: u(0) = u(1) = 0
+with f = (3*x + x^2)*exp(x)
+
+Analytical solution: -x*(x-1)*exp(x)
+
+Finite differences (FD) discretization: second-order diffusion operator
+"""
+__author__ = 'Jordi Vila-Pérez'
+__email__ = 'jvilap@mit.edu'
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation    #We have to load this
+from math import pi
+%matplotlib qt
+plt.close()
+
+"Number of points"
+N = 8
+Dx = 1/N
+x = np.linspace(0,1,N+1)
+
+"System matrix and RHS term"
+A = (1/Dx**2)*(2*np.diag(np.ones(N-1)) - np.diag(np.ones(N-2),-1) - np.diag(np.ones(N-2),1))
+F = (3*x[1:N] + x[1:N]**2)*np.exp(x[1:N])
+
+"Solution of the linear system AU=F"
+U = np.linalg.solve(A,F)
+u = np.concatenate(([0],U,[0]))
+ua = -x*(x-1)*np.exp(x)
+
+"Plotting solution"
+plt.plot(x,ua,'-r',linewidth=2,label='$u_a$')
+plt.plot(x,u,':ob',linewidth=2,label='$\widehat{u}$')
+plt.legend(fontsize=12,loc='upper left')
+plt.grid()
+plt.axis([0, 1,0, 0.5])
+plt.xlabel("x",fontsize=16)
+plt.ylabel("u",fontsize=16)
+
+"Compute error"
+error = np.max(np.abs(u-ua))
+print("Linf error u: %g\n" % error)
+
+
+
+
+
+
+

day_7/Diffusion/PoissonSpherical1D

@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+"""
+Solution of a 1D Poisson equation in spherical coordinates: -u_rr - (2/r)*u_r = f
+    Domain: [r0,r0+1]   IMPORTANT r>0
+    BC: u(r0) = u0, u'(r0+1) = 0
+    with f = 2*(2*A(rt)/r + B(rt))*erp(rt), being A(rt) = 2*rt^2 + rt - 3, and B(rt) = 2*rt^2 + 5*rt -2
+    and rt = r-r0
+    
+Analytical solution: 2*(rt)*(3-2*(rt))*exp(rt) + u0
+
+Finite differences (FD) discretization: second-order daccuracy
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Introduce first-order derivative with second-order accuracy
+
+"""
+__author__ = 'Jordi Vila-Pérez'
+__email__ = 'jvilap@mit.edu'
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation    #We have to load this
+from math import pi
+%matplotlib qt
+plt.close()
+
+"Number of points"
+N = 16
+r0 = 13
+Dr = 1/N
+r = np.linspace(0,1,N+1) + r0
+rN = np.concatenate((r, [r[N]+Dr]))
+rNt = rN-r0
+
+"Dirichlet boundary condition at r0"
+u0 = 1
+
+"System matrir and RHS term"
+A = (1/Dr**2)*(2*np.diag(np.ones(N+2)) - np.diag(np.ones(N+1),-1) - np.diag(np.ones(N+1),1))
+"First-order term"
+A = A + (1/Dr)*(np.diag(1/rN[1:N+2],-1) - np.diag(1/rN[0:N+1],1))
+
+"New source term"
+Ar = 2*rNt**2 + rNt - 3
+Br = 2*rNt**2 + 5*rNt - 2
+F = 2*(2*Ar/rN + Br)*np.exp(rNt)
+
+"Boundary condition at r=0"
+A[0,:] = np.concatenate(([1], np.zeros(N+1)))
+F[0] = u0
+
+"Boundary condition at r=0"
+A[N+1,:] = (0.5/Dr)*np.concatenate((np.zeros(N-1),[-1, 0, 1]))
+F[N+1] = 0
+
+"Solution of the linear system AU=F"
+u = np.linalg.solve(A,F)
+u = u[0:N+1]
+
+rt = r-r0
+ua = 2*rt*(3-2*rt)*np.exp(rt)+u0
+
+plt.plot(r,ua,'-r',linewidth=2,label='$u_a$')
+plt.plot(r,u,':ob',linewidth=2,label='$\widehat{u}$')
+plt.legend(fontsize=12,loc='upper left')
+plt.grid()
+plt.xlabel("x",fontsize=16)
+plt.ylabel("u",fontsize=16)
+
+"Compute error"
+error = np.max(np.abs(u-ua))
+print("Linf error u: %g\n" % error)
+

day_7/Diffusion/Slides_Diffusion.pdf

@@ -0,0 +1,17 @@
+---
+title: 
+description: Finite Difference Discretization of Elliptic Problems: the Heat Equation
+author: Jordi Vila-Pérez
+keywords: space-weather,finite differences, diffusion, heat equation, Poisson
+---
+
+Jordi Vila-Pérez | jvilap@mit.edu
+
+----------
+
+## Objectives
+
+- Understand basic concepts of discretization of simple PDE models using finite differences.
+- Be able to discretize a 1D elliptic problem (steady and transient).
+- Learn how to model different kinds of boundary conditions and how to approximate them numerically.
+- Get used to write and read code to solve PDEs numerically.

day_7/Diffusion/diffusion.md

@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+"""
+Solution of a 1D Convection-Diffusion equation: -nu*u_xx + c*u_x = f
+Domain: [0,1]
+BC: u(0) = u(1) = 0
+with f = 1
+
+Analytical solution: (1/c)*(x-((1-exp(c*x/nu))/(1-exp(c/nu))))
+
+Finite differences (FD) discretization:
+    - Second-order cntered differences advection scheme
+
+"""
+__author__ = 'Jordi Vila-Pérez'
+__email__ = 'jvilap@mit.edu'
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from math import pi
+%matplotlib qt
+plt.close()
+import matplotlib.animation as animation
+
+"Flow parameters"
+nu = 0.01
+c = 2
+
+"Number of points"
+N = 128
+Dx = 1/N
+x = np.linspace(0,1,N+1)
+
+"System matrix and RHS term"
+"Diffusion term"
+Diff = nu*(1/Dx**2)*(2*np.diag(np.ones(N-1)) - np.diag(np.ones(N-2),-1) - np.diag(np.ones(N-2),1))
+"Advection term: centered differences"
+Advp = -0.5*c*np.diag(np.ones(N-2),-1)
+Advm = -0.5*c*np.diag(np.ones(N-2),1)
+Adv = (1/Dx)*(Advp-Advm)
+A = Diff + Adv
+"Source term"
+F = np.ones(N-1)
+
+
+"Solution of the linear system AU=F"
+U = np.linalg.solve(A,F)
+u = np.concatenate(([0],U,[0]))
+ua = (1/c)*(x-((1-np.exp(c*x/nu))/(1-np.exp(c/nu))))
+
+plt.plot(x,ua,'-r',linewidth=2,label='$u_a$')
+plt.plot(x,u,':ob',linewidth=2,label='$\widehat{u}$')
+plt.legend(fontsize=12,loc='upper left')
+plt.grid()
+plt.axis([0, 1,0, 2/c])
+plt.xlabel("x",fontsize=16)
+plt.ylabel("u",fontsize=16)
+
+
+"Compute error"
+error = np.max(np.abs(u-ua))
+print("Linf error u: %g\n" % error)
+
+"Peclet number"
+P = np.abs(c*Dx/nu)
+print("Pe number Pe=%g\n" % P);
+
+

day_8/Convection/ConvectionDiffusion1D.py

@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+"""
+Solution of a 1D Convection-Diffusion equation: -nu*u_xx + c*u_x = f
+Domain: [0,1]
+BC: u(0) = u(1) = 0
+with f = 1
+
+Analytical solution: (1/c)*(x-((1-exp(c*x/nu))/(1-exp(c/nu))))
+
+Finite differences (FD) discretization:
+    - Second-order cntered differences advection scheme
+    - First-order upwind
+
+"""
+__author__ = 'Jordi Vila-Pérez'
+__email__ = 'jvilap@mit.edu'
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from math import pi
+%matplotlib qt
+plt.close()
+import matplotlib.animation as animation
+
+"Flow parameters"
+nu = 0.01
+c = 2
+
+"Scheme parameters"
+beta = 1
+
+"Number of points"
+N = 256
+Dx = 1/N
+x = np.linspace(0,1,N+1)
+
+"Time parameters"
+dt = 0.1
+time = np.arange(0,3+dt,dt)
+nt = np.size(time)
+
+"Initialize solution variable"
+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-1)) - np.diag(np.ones(N-2),-1) - np.diag(np.ones(N-2),1))
+
+    "Advection term:"
+        
+    "Sensor"
+    U0 = U[:,it]
+    uaux = np.concatenate(([0],U0,[0]))
+    Du = uaux[1:N+1] - uaux[0:N] + 1e-8
+    r = Du[0:N-1]/Du[1:N]
+    
+    "Limiter"
+    if beta>0:
+        phi = np.minimum(np.minimum(beta*r,1),np.minimum(r,beta))
+        phi = np.maximum(0,phi)
+    else:
+        phi = 2*r/(r**2 + 1)
+        
+    phim = phi[0:N-2]
+    phip = phi[1:N-1]
+        
+    
+    "Upwind scheme"
+    cp = np.max([c,0])
+    cm = np.min([c,0])
+    
+    Advp = cp*(np.diag(1-phi) - np.diag(1-phip,-1))
+    Advm = cm*(np.diag(1-phi) - np.diag(1-phim,1))
+    Alow = Advp-Advm
+    "Centered differences"
+    Advp = -0.5*c*np.diag(phip,-1)
+    Advm = -0.5*c*np.diag(phim,1)
+    Ahigh = Advp-Advm
+        
+    Adv = (1/Dx)*(Ahigh + Alow)
+    A = Diff + Adv
+    "Source term"
+    F = np.ones(N-1)
+    
+    "Temporal terms"
+    A = A + (1/dt)*np.diag(np.ones(N-1))
+    F = F + U0/dt
+
+
+    "Solution of the linear system AU=F"
+    u = np.linalg.solve(A,F)
+    U[:,it+1] = u
+    u = np.concatenate(([0],u,[0]))
+
+
+ua = (1/c)*(x-((1-np.exp(c*x/nu))/(1-np.exp(c/nu))))
+
+"Animation of the results"
+fig = plt.figure()
+ax = plt.axes(xlim =(0, 1),ylim =(0,1/c)) 
+plt.plot(x,ua,'-r',linewidth=2,label='$u_a$')
+myAnimation, = ax.plot([], [],':ob',linewidth=2)
+plt.grid()
+plt.xlabel("x",fontsize=16)
+plt.ylabel("u",fontsize=16)
+
+def animate(i):
+    
+    u = np.concatenate(([0],U[0:N+1,i],[0]))
+    plt.plot(x,u)
+    myAnimation.set_data(x, u)
+    return myAnimation,
+
+anim = animation.FuncAnimation(fig,animate,frames=range(1,nt),blit=True,repeat=False)
+
+"Compute error"
+error = np.max(np.abs(u-ua))
+print("Linf error u: %g\n" % error)
+
+"Peclet number"
+P = np.abs(c*Dx/nu)
+print("Pe number Pe=%g\n" % P);
+
+"CFL number"
+CFL = np.abs(c*dt/Dx)
+print("CFL number CFL=%g\n" % CFL);
+