with the correct exact solution

euler_error/euler_error_moving_domain_diffusion_advection.py

@@ -18,15 +18,16 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
     times = np.linspace(0, tmax, Nt+1)
     times = np.linspace(0, tmax, Nt+1)
     dt = times[1] - times[0]
     dt = times[1] - times[0]
 
 
-    c0 = df.interpolate(df.Expression('1 + 0.2*cos(pi*x[0]/L)', pi = np.pi, L=L, degree=1), function_space)
+    c0 = df.interpolate(df.Expression('1 + 0.2*cos(pi*x[0]/L)', pi = np.pi, L = L, degree=1), function_space)
 
 
     # form
     # form
     u = df.interpolate(v, function_space)
     u = df.interpolate(v, function_space)
-    form = (df.inner((c-c0)/dt, tc)
-            + df.inner((u*c).dx(0), tc)
-            + D * df.inner(tc.dx(0), c.dx(0)))*df.dx(mesh)
 
 
-    for _ in progressbar.progressbar(range(1, len(times)+1)):
+    form = (df.inner((c - c0)/dt, tc)
+            + df.inner((u * c).dx(0), tc)
+            + D * df.inner(tc.dx(0), c.dx(0))) * df.dx(mesh)
+
+    for _ in progressbar.progressbar(range(1, len(times))):
         df.solve(form == 0, c)
         df.solve(form == 0, 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 = mesh.coordinates()
         x = mesh.coordinates()
@@ -36,16 +37,22 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
 
 
 Nx, L, tmax, D, k = 2000, 1, 1, 1, 1
 Nx, L, tmax, D, k = 2000, 1, 1, 1, 1
 v = df.Expression('k*x[0]', k=k, degree=1)
 v = df.Expression('k*x[0]', k=k, degree=1)
-nt_array = np.array([50, 100, 200, 400, 800])
-# , 1600, 3200, 6400, 12800])
+nt_array = np.array([50, 100, 200, 400, 800, 1600, 3200, 6400, 12800])
 dt_array = tmax/(nt_array - 1)
 dt_array = tmax/(nt_array - 1)
 error_in_c = np.zeros(len(nt_array))
 error_in_c = np.zeros(len(nt_array))
 
 
 for i in range(len(nt_array)):
 for i in range(len(nt_array)):
     x = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[-1]
     x = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[-1]
-    ss_c_exact = 2 * np.exp(-k*tmax) * (1 + 0.2*np.cos(np.pi*x*np.exp(-k*tmax)) * np.exp(-np.pi**2 * (D/(2*k)) *(1-np.exp(-2*k*tmax))))
     ss_c = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[0]
     ss_c = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[0]
+
+    xprime = x / (L * np.exp(k * tmax))
+    tprime = ( D / (2 * k * L**2)) * ( 1 - np.exp(-2 * k * tmax))
+    int = np.exp(-k * tmax)
+    ss_c_exact = int * (1 + 0.2 * np.cos(np.pi * xprime) * np.exp(-np.pi**2 * tprime))
+
     error_in_c[i] = np.max(np.abs(ss_c - ss_c_exact))
     error_in_c[i] = np.max(np.abs(ss_c - ss_c_exact))
 
 
 plt.scatter(dt_array, error_in_c)
 plt.scatter(dt_array, error_in_c)
 plt.show()
 plt.show()
+
+# bug: tmax = 1 -- trend ok, not for any other tmax