test for diffusion on fixed domain

euler_error/diffusion_fixed_domain_solution.py

+import dolfin as df
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+import progressbar
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+Nx, L, D, tmax, dt = 32, 1, 0.1, 0.25, 0.25/400
+Nt = int(tmax/dt)
+
+mesh = df.IntervalMesh(Nx, 0, L)
+x = mesh.coordinates()[:, 0]
+times = np.linspace(0, tmax, Nt+1)
+
+SFS = df.FunctionSpace(mesh, 'P', 1)
+c = df.Function(SFS)
+tc = df.TestFunction(SFS)
+c0 = df.Function(SFS)
+
+c0.interpolate(df.Expression('1 + 0.1 * cos(2*pi*x[0]/L) + 0.2 * cos(3*pi*x[0]/L)', pi=np.pi,L=L, degree=1))
+
+c_exact = np.zeros((Nt+1, Nx+1))
+for i in range(Nt+1):
+    c_exact[i] = 1 + 0.1 * np.cos(2*np.pi*x/L) * np.exp(-4*np.pi**2*D*times[i]/L**2) + 0.2 * np.cos(3*np.pi*x/L) * np.exp(-9*np.pi**2*D*times[i]/L**2)
+
+c_array = np.zeros((Nt+1, Nx+1))
+c_array[0] = c0.compute_vertex_values(mesh)
+
+cform = (df.inner((c - c0)/dt, tc) 
+         + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))) * df.dx
+
+for i in progressbar.progressbar(range(1, Nt+1)):
+    df.solve(cform == 0, c)
+    c0.assign(c)
+    c_array[i] = c0.compute_vertex_values(mesh)
+
+fig, (ax, ax1) = plt.subplots(2, 1)
+fig.suptitle(r'$N_x = %d, \Delta t = %4.3f$' % (Nx, dt))
+cplot, = ax.plot(x, c_array[0], 'ro', mfc='none', ms=6, label='Numerics')
+ceplot, = ax.plot(x, c_exact[0], label='Exact')
+ax.set_xlabel(r'$x$')
+ax.set_ylabel(r'$c(x,t)$')
+ax.legend(loc=1)
+
+error = np.abs(c_array - c_exact)
+print('Max error at t=0.25 is',np.max(error[-1]))
+
+err_plot, = ax1.plot(x, error[1], 'bo', mfc='none', ms=6)
+ax1.set_ylim(np.min(error), np.max(error))
+ax1.set_xlabel(r'$x$')
+ax1.set_ylabel(r'$c(x,t)-c_{\rm exact}(x,t)$')
+
+def update(val):
+    ti = (abs(times-val)).argmin()
+    cplot.set_ydata(c_array[ti])
+    ceplot.set_ydata(c_exact[ti])
+    
+    err_plot.set_ydata(error[ti])
+    plt.draw()
+    
+sax = plt.axes([0.1, 0.92, 0.7, 0.025])
+sl = Slider(sax, 't', times.min(), times.max(), valinit=times.min())
+sl.on_changed(update)
+
+plt.show()