max error scaling linearly with dt

diffusion/diffusion_moving_domain.py

@@ -7,61 +7,61 @@ import progressbar
 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
 
 
-Nx, L, D, tmax, dt = 32, 1.0, 0.1, 2.0, 0.05
-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)', 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)
-
-c_array = np.zeros((Nt+1, Nx+1))
-c_array[0] = c0.compute_vertex_values(mesh)
-
-cform = (df.inner((c - c0)/dt, tc) 
+def advection_diffusion(Nx, L, Nt, tmax, D):
+    # 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
+    x = mesh.coordinates()[:, 0]
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+
+    # initial condition
+    c0 = df.Function(SFS)
+    c0.interpolate(df.Expression('1 + 0.1 * cos(2*pi*x[0]/L)', pi=np.pi,L=L, degree=1))
+
+    # arrays
+    c_array = np.zeros((Nt+1, Nx+1))
+    c_array[0] = c0.compute_vertex_values(mesh)
+
+    # form
+    cform = (df.inner((c - c0)/dt, tc) 
             + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))) * df.dx
             + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))) * df.dx
 
 
-for i in progressbar.progressbar(range(1, Nt+1)):
+    # solve
+    for i in progressbar.progressbar(range(1, Nt+1)):
         df.solve(cform == 0, c)
         df.solve(cform == 0, c)
         c0.assign(c)
         c0.assign(c)
         c_array[i] = c0.compute_vertex_values(mesh)
         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 = c_array-c_exact
-print(np.max(error))
+    return c_array
 
 
-err_plot, = ax1.plot(x, error[0], '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)$')
+# parameters
+Nx, L, D, tmax = 64, 1.0, 0.1, 1.0
+nt_array = np.array([50, 100, 200, 400, 600, 800, 1600])
+dt_array = tmax/nt_array
+mesh = df.IntervalMesh(Nx, 0, L)
+x = mesh.coordinates()[:, 0]
 
 
-def update(val):
-    ti = (abs(times-val)).argmin()
-    cplot.set_ydata(c_array[ti])
-    ceplot.set_ydata(c_exact[ti])
+# error array
+error = np.zeros(len(nt_array))
 
 
-    err_plot.set_ydata(error[ti])
-    plt.draw()
+for i in range(0, len(nt_array)):
 
 
-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)
+    # exact solution
+    c_exact = np.zeros((nt_array[i]+1, Nx+1))
+    times = np.linspace(0, tmax, nt_array[i]+1)
+    for j in range(nt_array[i]+1):
+        c_exact[j] = 1 + 0.1 * np.cos(2*np.pi*x/L) * np.exp(-4*np.pi**2*D*times[j]/L**2)
+    c = advection_diffusion(Nx, L, nt_array[i], tmax, D)
+    error[i] = np.max(np.abs(c - c_exact))
 
 
+fig, ax = plt.subplots(1)
+ax.plot(dt_array, error, 'bo')
+ax.set_xlabel('$dt$')
+ax.set_ylabel('error')
 plt.show()
 plt.show()
+