analytical and numerical solutions match exactly for alpha zero

moving_domain/moving_heat_equation_analytical.py

+import dolfin as df
 import numpy as np
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 from matplotlib.widgets import Slider
 import progressbar
 import progressbar
-import dolfin as df
+
 df.set_log_level(df.LogLevel.ERROR)
 df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
 df.parameters['form_compiler']['optimize'] = True
 
 
-# parameters
-alpha = 1.0
-T = 1
-dt = 0.001
-L0 = 1
-D = 1.0
-Nx = 2000
-Nt = 1000
-t = np.linspace(0,T,Nt)
-
-# diffusion and advection
-def diffusion(c, tc):
-        return (D * df.inner(c.dx(0), tc.dx(0)))
-
-def advection(c, tc, v):
-        u = df.interpolate(v, c.function_space())
-        return (df.inner((u*c).dx(0),tc))
-
-# create mesh
-mesh = df.IntervalMesh(Nx, 0, L0)
-x = mesh.coordinates()
-# v = df.Constant(1.0)
-v = df.Expression('alpha*x[0]',alpha = alpha,degree=1)
-
-# create function space
-conc_element = df.FiniteElement('P', mesh.ufl_cell(), 1)
-function_space = df.FunctionSpace(mesh,conc_element)
-
-# initial condition
-c0 = df.interpolate(df.Expression('1 + 0.2*cos(pi*x[0]/L0)', pi = np.pi, L0 = L0,degree=1),function_space)
-c0_array = c0.compute_vertex_values(mesh)
-
-# define variational problem
-c = df.Function(function_space)
-tc = df.TestFunction(function_space)
-form = (df.inner((c - c0)/dt, tc)               
-+ diffusion(c,tc)
-+ advection(c, tc, v)
-)
-form = form * df.dx
-
-# time stepping
-ctot = np.zeros_like(t)
-x_array = np.zeros((len(t), mesh.num_vertices()))
-
-x_array[0] = mesh.coordinates()[:,0]
-c_array = np.zeros((len(t), len(c0_array)))
-c_array[0] = c0_array
-ctot[0] = df.assemble(c0*df.dx(mesh))
-
-for n in progressbar.progressbar(range(1, len(t))):
-    df.solve(form == 0, c)
-    c_array[n] = c.compute_vertex_values(mesh)
+def advection_diffusion(Nx, L, Nt, tmax, D, alpha):
+    # mesh, function space, function, test function
+    mesh = df.IntervalMesh(Nx, 0, L)
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    c = df.Function(SFS)
+    tc = df.TestFunction(SFS)
+
+    # x and t arrays
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+
+    # initial condition
+    c0 = df.Function(SFS)
+    c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=L, degree=1))
+
+    # arrays
+    c_array = np.zeros((Nt+1, Nx+1))
+    x_array = np.zeros((Nt+1, Nx+1))
+    x_array[0] = mesh.coordinates()[:, 0]
+    c_array[0] = c0.compute_vertex_values(mesh)
+
+    # velocity
+    v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
+    u = df.interpolate(v, SFS)
+
+    # form
+    cform = (df.inner((c - c0)/dt, tc) 
+            + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+            + df.inner((u*c).dx(0), tc) )* df.dx
+
+    # solve
+    for i in progressbar.progressbar(range(1, Nt+1)):
+        df.solve(cform == 0, c)
+        c_array[i] = c.compute_vertex_values(mesh)
         c0.assign(c)
         c0.assign(c)
-    ctot[n] = df.assemble(c0*df.dx(mesh))
+        df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
+        x_array[i] = mesh.coordinates()[:,0]
 
 
-    df.ALE.move(mesh, df.Expression('v*dt',v=v,dt=dt,degree=1))
-    x_array[n] = mesh.coordinates()[:,0]
+    return c_array, x_array
 
 
-# analytical solution
-c_exact= np.zeros((len(t), len(x)))
+# plot c(x,t) numerical and analytical for given dt
+Nx, L, Nt, tmax, D, alpha = 64, 1, 100, 1, 1, 1
+x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
 
 
-for i in range(len(t)):
-        xprime = x_array[0] / (L0 * np.exp(alpha * t[i]))
-        tprime = ( D / (2 * alpha * L0**2)) * ( 1 - np.exp(-2 * alpha * t[i]))
-        int = np.exp(-alpha * t[i])
-        c_exact[i] = int * (1 + 0.2 * np.cos(np.pi * xprime) * np.exp(-np.pi**2 * tprime))
+# exact solution
+c_exact = np.zeros((Nt+1, Nx+1))
+times = np.linspace(0, tmax, Nt+1)
+for j in range(Nt+1):
+    if alpha == 0:
+        c_exact[j] = 1 + 0.2 * np.cos(np.pi * x[j]/L) * np.exp(-np.pi**2 * D * times[j]/L**2)
+    else:
+        c_exact[j] = 1 + 0.2 * np.cos(np.pi*x[j]*np.exp(-alpha*times[j])/L) * np.exp(-np.pi**2 * D * (1 - np.exp(-2 * alpha * times[j]))/(2 * alpha * L**2)) * np.exp(-alpha * times[j])
 
 
-# plot c(x,t) computed numerically
-fig, ax_comp = plt.subplots(1,1,figsize=(8,6))
-ax_comp.set_xlabel(r'$x$')
-ax_comp.set_ylabel(r'$c(x,t)$')
-ax_comp.set_xlim(np.min(x_array)-1, np.max(x_array)+1)
-ax_comp.set_ylim(np.min(c_array)-1, np.max(c_array)+1)
-cplot, = ax_comp.plot(x_array[0], c0_array)
-c_exactplot, = ax_comp.plot(x_array[0], c_exact[0], 'ro', markersize=3, markevery=50)
+c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
+times = np.linspace(0, tmax, Nt+1)
 
 
-def update(value):
-        ti = np.abs(t-value).argmin()
-        cplot.set_xdata(x_array[ti])
-        cplot.set_ydata(c_array[ti])
-        c_exactplot.set_xdata(x_array[ti])
-        c_exactplot.set_ydata(c_exact[ti])
+fig, ax = plt.subplots(1,1,figsize=(8,6))
+ax.set_xlabel(r'$x$')
+ax.set_ylabel(r'$c(x,t)$')
+ax.set_xlim(np.min(x)-2, np.max(x)+2)
+ax.set_ylim(np.min(c)-2, np.max(c)+2)
 
 
+cplot, = ax.plot(x[0], c[0], 'go', ms=1)
+cexactplot, = ax.plot(x[0], c_exact[0]) 
+
+def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_xdata(x[ti])
+        cplot.set_ydata(c[ti])
+        cexactplot.set_xdata(x[ti])
+        cexactplot.set_ydata(c_exact[ti])
         plt.draw()
         plt.draw()
+        
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
-slider = Slider(sax, r'$t/\tau$', min(t), max(t),
-        valinit=min(t), valfmt='%3.1f',
+slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+        valinit=min(times), valfmt='%3.1f',
         fc='#999999')
         fc='#999999')
 slider.drawon = False
 slider.drawon = False
 slider.on_changed(update)
 slider.on_changed(update)