increasing Nx not affecting error values

euler_error/euler_error_moving_domain_diffusion_advection.py

@@ -41,15 +41,15 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
         df.solve(form == 0, c)
         c_array[i] = c.compute_vertex_values(mesh)
         c0.assign(c)
-        # df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
+        df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
         x_array[i] = mesh.coordinates()[:,0]
         
 
     return [c_array,x_array]
 
-Nx, L, tmax, D, alpha = 100, 1, 1, 1, 1
+Nx, L, tmax, D, alpha = 512, 1, 0.25, 1, 1
 v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
-nt_array = np.array([5, 10, 15, 20])
+nt_array = np.array([5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 30, 35, 40, 45, 50])
 # , 400, 800, 1600, 3200, 6400, 12800])
 dt_array = tmax/(nt_array)
 error_in_c = np.zeros(len(nt_array))
@@ -61,18 +61,21 @@ for i in range(len(nt_array)):
     c_exact = np.zeros((nt_array[i]+1, Nx+1))
     for j in range(0, nt_array[i]+1):
         if alpha == 0:
-            c_exact[j] = 1 + 0.2*np.cos(np.pi*x[j]/L)*np.exp(-np.pi**2*D/L**2)
+            beta = -D * np.pi**2 * times[j]/L**2
+            c_exact[j] = 1 + 0.2 * np.cos(np.pi * x[j]/L) * np.exp(beta)
         else:
-            c_exact[j] = np.exp(-alpha*times[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)))
+            l = L * np.exp(alpha * times[j])
+            beta = (-D * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
+            c_exact[j] = 1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(-alpha * times[j]) * np.exp(beta)
             
-    error_in_c[i] = np.max(np.abs(c - c_exact))
+    error_in_c[i] = np.max(np.abs(c[-1] - c_exact[-1]))
 
 def linear_func(x, a, b):
     return a * x + b
     
 popt, pcov = curve_fit(linear_func, np.log(dt_array), np.log(error_in_c))
 print('The slope of log(error) vs log(dt) is', popt[0])
-plt.loglog(dt_array, error_in_c, 'bo', ms=3)
+plt.loglog(dt_array, error_in_c, 'bo')
 plt.xlabel('log(dt)')
 plt.ylabel('log(error)')
 plt.show()