Merge branch 'master' of gitlab.icts.res.in:jigyasa.watwani/growth-pattern-control

euler_error/linear_growth_euler_errors.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/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
+import matplotlib.pyplot as plt
+import progressbar
+from scipy.optimize import curve_fit
+
+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)
+
+    # initial condition
+    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
+    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.project(v, function_space)
+
+    form = (df.inner((c - c0)/dt, tc)
+            + df.inner((u * c0).dx(0), tc)
+            + D * df.inner(tc.dx(0), c.dx(0))) * df.dx(mesh)
+
+    for i in progressbar.progressbar(range(1, len(times))):
+        v.t = times[i]
+        u.assign(df.project(v, function_space))
+        df.solve(form == 0, c)
+        c_array[i] = c.compute_vertex_values(mesh)
+        c0.assign(c)
+        df.ALE.move(mesh, df.project(u*dt, function_space))
+        x_array[i] = mesh.coordinates()[:,0]
+        
+
+    return [c_array,x_array]
+
+Nx, L, tmax, D, b, m = 100, 1, 1, 0.1, 0.01, 2
+v = df.Expression('b*x[0]/(L + b*t)', b=b, L=L, t=0, degree=1)
+Nt_array = np.array([5, 10, 15, 20, 25, 30])
+# , 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)):
+    c, x = advection_diffusion(Nx, L, Nt_array[i], tmax, v, D)
+    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):
+        l = L + b * times[j]
+        c_exact[j] = (L/l) * (1 + 0.2 * np.cos(m * np.pi * x[j]/l) * np.exp(-m**2 * np.pi**2 * D * times[j]/(L * l)))
+    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')
+plt.xlabel('log(dt)')
+plt.ylabel('log(error)')
+plt.show()