rewriting analytical solution

euler_error/euler_error_moving_domain_diffusion_advection.py

@@ -7,20 +7,20 @@
 import numpy as np
 import numpy as np
 import dolfin as df
 import dolfin as df
 import matplotlib.pyplot as plt
 import matplotlib.pyplot as plt
-import progressbar
+# import progressbar
 from scipy.optimize import curve_fit
 from scipy.optimize import curve_fit
 
 
 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
 
 
-def advection_diffusion(Nx, L, Nt, tmax, v, D):
+def advection_diffusion(Nx, L, Nt, tmax, v, D, m):
     # mesh, function space, function, test function
     # mesh, function space, function, test function
     mesh = df.IntervalMesh(Nx, 0, L)
     mesh = df.IntervalMesh(Nx, 0, L)
     function_space = df.FunctionSpace(mesh, 'P', 1)
     function_space = df.FunctionSpace(mesh, 'P', 1)
     c, tc = df.Function(function_space), df.TestFunction(function_space)
     c, tc = df.Function(function_space), df.TestFunction(function_space)
 
 
     # initial condition
     # initial condition
-    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(m * pi*x[0]/L)', pi = np.pi, L = L, m = m, degree=1), function_space)
 
 
     # arrays
     # arrays
     times = np.linspace(0, tmax, Nt+1)
     times = np.linspace(0, tmax, Nt+1)
@@ -37,7 +37,7 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
             + df.inner((u * c).dx(0), tc)
             + df.inner((u * c).dx(0), tc)
             + D * df.inner(tc.dx(0), c.dx(0))) * df.dx(mesh)
             + D * df.inner(tc.dx(0), c.dx(0))) * df.dx(mesh)
 
 
-    for i in progressbar.progressbar(range(1, len(times))):
+    for i in range(1, len(times)):
         df.solve(form == 0, c)
         df.solve(form == 0, c)
         c_array[i] = c.compute_vertex_values(mesh)
         c_array[i] = c.compute_vertex_values(mesh)
         c0.assign(c)
         c0.assign(c)
@@ -47,26 +47,25 @@ def advection_diffusion(Nx, L, Nt, tmax, v, D):
 
 
     return [c_array,x_array]
     return [c_array,x_array]
 
 
-Nx, L, tmax, D, alpha = 3200, 1, 1, 1, 0.1
+Nx, L, tmax, D, alpha, m = 3200, 1, 1, 0.1, 0.01, 2
 v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
 v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
-nt_array = np.array([5,6,7,8,9,10, 20, 30])
+nt_array = np.array([5,6,7,8,9,10, 20, 30, 40, 50])
 # , 400, 800, 1600, 3200, 6400, 12800])
 # , 400, 800, 1600, 3200, 6400, 12800])
 dt_array = tmax/(nt_array)
 dt_array = tmax/(nt_array)
 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]
-    c = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[0]
+    c, x = advection_diffusion(Nx, L, nt_array[i], tmax, v, D, m)
     times = np.linspace(0, tmax, nt_array[i]+1)
     times = np.linspace(0, tmax, nt_array[i]+1)
     c_exact = np.zeros((nt_array[i]+1, Nx+1))
     c_exact = np.zeros((nt_array[i]+1, Nx+1))
     for j in range(0, nt_array[i]+1):
     for j in range(0, nt_array[i]+1):
         if alpha == 0:
         if alpha == 0:
-            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)
+            beta = -D * m**2 * np.pi**2 * times[j]/L**2
+            c_exact[j] = 1 + 0.2 * np.cos(m * np.pi * x[j]/L) * np.exp(beta)
         else:
         else:
             l = L * np.exp(alpha * times[j])
             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] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
+            beta = (-D * m**2 * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
+            c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(m * np.pi*x[j]/l) * np.exp(beta))
             
             
     error_in_c[i] = np.max(np.abs(c[-1] - c_exact[-1]))
     error_in_c[i] = np.max(np.abs(c[-1] - c_exact[-1]))