changed range of dt + log log plot

euler_error/euler_error_moving_domain_diffusion_advection.py

 # moving domain diffusion advection equation
+# advection velocity v = alpha*x. domain also moves by thsi velocity
 # boundary condition: diffusive flux at the boundaries [0,1] is zero
-# initial condition c(x,0) = 1 + 0.2*cos(pi x) ; satisfies the boundary condition
+# initial condition c(x,0) = 1 + 0.2*cos(pi x/L) ; satisfies the boundary condition
+# exact solution: c(x,t) = exp(-alpha * t)(1 + 0.2 cos(pi * x/L0 *exp(-alpha*t)) * exp(-pi**2*D*(1-exp(-2*alpha*t)/(2*alpha*L0**2))))
 
 import numpy as np
 import dolfin as df
@@ -12,14 +14,22 @@ df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
 
 def advection_diffusion(Nx, L, Nt, tmax, v, D):
+    # mesh, function space, function, test function
     mesh = df.IntervalMesh(Nx, 0, L)
     function_space = df.FunctionSpace(mesh, 'P', 1)
     c, tc = df.Function(function_space), df.TestFunction(function_space)
-    times = np.linspace(0, tmax, Nt+1)
-    dt = times[1] - times[0]
 
+    # initial condition
     c0 = df.interpolate(df.Expression('1 + 0.2*cos(pi*x[0]/L)', pi = np.pi, L = L, degree=1), function_space)
 
+    # arrays
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+    x_array = np.zeros((Nt+1, Nx+1))
+    x_array[0] = mesh.coordinates()[:,0]
+    c_array = np.zeros((Nt+1, Nx+1))
+    c_array[0] = c0.compute_vertex_values(mesh)
+
     # form
     u = df.interpolate(v, function_space)
 
@@ -27,32 +37,42 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
             + 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))):
+    for i in progressbar.progressbar(range(1, len(times))):
         df.solve(form == 0, c)
-        df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
-        x = mesh.coordinates()
+        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))
+        x_array[i] = mesh.coordinates()[:,0]
+        
 
-    return [c.compute_vertex_values(mesh), x]
+    return [c_array,x_array]
 
-Nx, L, tmax, D, k = 2000, 1, 1, 1, 1
-v = df.Expression('k*x[0]', k=k, degree=1)
-nt_array = np.array([50, 100, 200, 400, 800, 1600, 3200, 6400, 12800])
-dt_array = tmax/(nt_array - 1)
+Nx, L, tmax, D, alpha = 100, 1, 1, 1, 1
+v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
+nt_array = np.array([5, 10, 15, 20])
+# , 400, 800, 1600, 3200, 6400, 12800])
+dt_array = tmax/(nt_array)
 error_in_c = np.zeros(len(nt_array))
 
 for i in range(len(nt_array)):
     x = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[-1]
-    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))
+    c = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[0]
+    times = np.linspace(0, tmax, nt_array[i]+1)
+    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)
+        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)))
 
-    error_in_c[i] = np.max(np.abs(ss_c - ss_c_exact))
+    error_in_c[i] = np.max(np.abs(c - c_exact))
 
-plt.scatter(dt_array, error_in_c)
+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.xlabel('log(dt)')
+plt.ylabel('log(error)')
 plt.show()
-
-# bug: tmax = 1 -- trend ok, not for any other tmax