1D growing domain diffusion

growing_domain/given_v_get_rho_c/1D_growth/exp_growth.py

+import dolfin as df
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+import progressbar
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+def advection_diffusion(Nx, L, Nt, tmax, D, alpha):
+    # mesh, function space, function, test function
+    mesh = df.IntervalMesh(Nx, 0, L)
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    c = df.Function(SFS)
+    tc = df.TestFunction(SFS)
+
+    # x and t arrays
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+
+    # initial condition
+    c0 = df.Function(SFS)
+    c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=L, degree=1))
+
+    # arrays
+    c_array = np.zeros((Nt+1, Nx+1))
+    x_array = np.zeros((Nt+1, Nx+1))
+    x_array[0] = mesh.coordinates()[:, 0]
+    c_array[0] = c0.compute_vertex_values(mesh)
+
+    # velocity
+    v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
+    u = df.interpolate(v, SFS)
+
+    # form
+    cform = (df.inner((c - c0)/dt, tc) 
+            + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+            + df.inner((u*c).dx(0), tc) )* df.dx
+
+    # solve
+    for i in progressbar.progressbar(range(1, Nt+1)):
+        df.solve(cform == 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))
+        x_array[i] = mesh.coordinates()[:,0]
+
+    return c_array, x_array
+
+# plot c(x,t) numerical and analytical for given dt
+dx, L, dt, tmax, D, alpha = 0.01, 1, 0.01, 20, 0.1, 0.01
+Nx = int(L/dx)
+Nt = int(tmax/dt)
+x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
+
+# exact solution
+c_exact = np.zeros((Nt+1, Nx+1))
+times = np.linspace(0, tmax, Nt+1)
+
+for j in range(Nt+1):
+    l = L * np.exp(alpha * times[j])
+    if alpha == 0:
+        beta = -D * np.pi**2 * times[j]/L**2
+    else:
+        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))
+
+# numerical solution
+c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
+
+fig, (ax, ax1) = plt.subplots(2,1,figsize=(8,6))
+ax.set_xlabel(r'$x$')
+ax.set_ylabel(r'$c(x,t)$')
+ax.set_xlim(np.min(x), np.max(x)+2)
+ax.set_ylim(min(np.min(c), np.min(c_exact))-1, max(np.max(c), np.max(c_exact))+1)
+
+cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
+cexactplot, = ax.plot(x[0], c_exact[0], label = 'Exact solution') 
+
+error = np.abs(c - c_exact)
+print('dt = %6.3f, max error = %6.5f' % (tmax/Nt, np.max(error)))
+errorplot, = ax1.plot(x[0], error[0], 'bo', mfc = 'none')
+ax1.set_ylabel(r'$|c(x,t) - c_{exact}(x,t)|$')
+ax1.set_xlabel('$x$')
+ax1.set_xlim(np.min(x), np.max(x)+2)
+ax1.set_ylim([np.min(error)-1, np.max(error)+1])
+
+def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_xdata(x[ti])
+        cplot.set_ydata(c[ti])
+        cexactplot.set_xdata(x[ti])
+        cexactplot.set_ydata(c_exact[ti])
+        errorplot.set_xdata(x[ti])
+
+        errorplot.set_ydata(error[ti])
+        plt.draw()
+        
+sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+        valinit=min(times), valfmt='%3.1f',
+        fc='#999999')
+slider.drawon = False
+slider.on_changed(update)
+ax.legend(loc=0)
+
+fig, ax2 = plt.subplots(1,1,figsize=(8,6))
+ax2.semilogy(times, c[:,0], label='$c_{num}(0,t)$')
+if alpha == 0:
+    ax2.semilogy(times, 1 + 0.2 * np.exp(-D *np.pi**2 * times/L**2), label = '$c_{an}(0,t)$')
+else:
+    ax2.semilogy(times, np.exp(-alpha * times) * (1 + 0.2 * np.exp(-D * (np.pi**2/(2*alpha*L**2)) * (1 - np.exp(-2 * alpha * times)))), label = '$c_{an}(0,t)$' )
+ax2.legend()
+plt.show()
+

growing_domain/given_v_get_rho_c/1D_growth/linear_growth.py

+import dolfin as df
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+# parameters
+tmax, dt, L, dx, b, m, D = 20, 0.01, 1, 0.01, 0.1, 2, 0.1
+Nt = int(tmax/dt)
+Nx = int(L/dx)
+
+# mesh, function space, functions
+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)
+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)
+
+# velocity 
+v = df.Expression('b*x[0]/(L + b *t)', b=b, L=L, t=0, degree=0)
+vh = df.project(v, function_space)
+
+# form
+cform = (df.inner((c - c0)/dt, tc) 
+        + D * df.inner((c).dx(0), (tc).dx(0))
+        + df.inner((vh*c0).dx(0), tc) 
+        )* df.dx
+
+# time stepping
+for i in range(1, Nt+1):
+    v.t = times[i]
+    vh.assign(df.project(v, function_space))
+    df.solve(cform == 0, c)
+    c_array[i] = c.compute_vertex_values(mesh)
+    c0.assign(c)
+    df.ALE.move(mesh, df.project(vh*dt, function_space))
+    x_array[i] = mesh.coordinates()[:, 0]
+
+
+# exact solution
+c_exact = np.zeros((Nt + 1, Nx + 1))
+for j in range(0, Nt+1):
+    l = L + b * times[j]
+    c_exact[j] = (L/l) * (1 + 0.2 * np.cos(m * np.pi * x_array[j]/l) * np.exp(-m**2 * np.pi**2 * D * times[j]/(L * l)))
+
+
+# plotting
+fig, ax = plt.subplots(1,1,figsize=(8,6))
+ax.set_xlabel(r'$x$')
+ax.set_ylabel(r'$c(x,t)$')
+ax.set_xlim(np.min(x_array), np.max(x_array))
+ax.set_ylim(min(np.min(c_exact), np.min(c_array)), max(np.max(c_array), np.max(c_exact)))
+ax.grid(True)
+
+cplot, = ax.plot(x_array[0], c_array[0], '--',label = 'Numerical solution')
+c_exactplot, = ax.plot(x_array[0], c_exact[0],label = 'Exact solution')
+
+def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_xdata(x_array[ti])
+        cplot.set_ydata(c_array[ti])
+        c_exactplot.set_xdata(x_array[ti])
+        c_exactplot.set_ydata(c_exact[ti])
+        plt.draw()
+        
+sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+        valinit=min(times), valfmt='%3.1f',
+        fc='#999999')
+slider.drawon = False
+slider.on_changed(update)
+ax.legend(loc=0)
+plt.show()
\ No newline at end of file

growing_domain/given_v_get_rho_c/1D_growth/time_dependent_alpha_exp_growth.py.py

+import dolfin as df
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+import progressbar
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+def advection_diffusion(Nx, L, Nt, tmax, D, tstop, m):
+    # mesh, function space, function, test function
+    mesh = df.IntervalMesh(Nx, 0, L)
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    c = df.Function(SFS)
+    tc = df.TestFunction(SFS)
+
+    # x and t arrays
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+
+    # initial condition
+    c0 = df.Function(SFS)
+    c0.interpolate(df.Expression('1 + 0.2 * cos(m * pi*x[0]/L)', m = m, pi=np.pi, L=L, degree=1))
+
+    # arrays
+    c_array = np.zeros((Nt+1, Nx+1))
+    x_array = np.zeros((Nt+1, Nx+1))
+    x_array[0] = mesh.coordinates()[:, 0]
+    c_array[0] = c0.compute_vertex_values(mesh)
+
+    # velocity
+    u = df.Expression('(t < tstop ? alpha0 : 0)*x[0]', alpha0 = alpha0, tstop=tstop, t=0, degree=0)
+    uh = df.project(u, SFS)
+
+    # form
+    cform = (df.inner((c - c0)/dt, tc) 
+            + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+            + df.inner((uh*c).dx(0), tc) )* df.dx
+
+    # solve
+    for i in progressbar.progressbar(range(1, Nt+1)):
+        u.t = times[i]
+        uh.assign(df.project(u, SFS))
+        df.solve(cform == 0, c)
+        c_array[i] = c.compute_vertex_values(mesh)
+        c0.assign(c)       
+        df.ALE.move(mesh, df.project(uh*dt, SFS))
+        x_array[i] = mesh.coordinates()[:,0]
+
+    return c_array, x_array
+
+# plot c(x,t) numerical and analytical for given dt
+dx, L, dt, tmax, D, m = 0.01, 1, 0.01, 100, 0.01, 2
+alpha0, tstop = 0.01, 3
+Nx = int(L/dx)
+Nt = int(tmax/dt)
+
+times = np.linspace(0, tmax, Nt+1)
+
+# numerical solution
+c, x = advection_diffusion(Nx, L, Nt, tmax, D, tstop, m)
+
+# analytical solution
+c_exact = np.zeros((Nt+1, Nx+1))
+for j in range(0, len(times)):
+    if times[j] <= tstop:
+        # diffusion-advection on moving domain with velocity alpha0*x
+        alpha = alpha0
+        l = L * np.exp(alpha * times[j])
+        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))
+    else: 
+        # diffusion on fixed domain of length L = L0*exp(alpha0 * tc) with initial condition to be the profile at tc
+        l = L * np.exp(alpha0 * tstop)
+        beta = -D * m**2 * np.pi**2*(times[j] - tstop)/l**2 - np.pi**2 * D * m**2/(2*L**2*alpha0) * (1 - np.exp(-2 *alpha0 *tstop))
+        c_exact[j] = np.exp(-alpha0 * tstop) * (1 + 0.2 * np.exp(beta) * np.cos(m * np.pi * x[j]/l))
+
+
+fig, ax = plt.subplots(1,1,figsize=(8,6))
+ax.set_xlabel(r'$x$')
+ax.set_ylabel(r'$c(x,t)$')
+ax.set_xlim(np.min(x), np.max(x))
+ax.set_ylim(min(np.min(c), np.min(c_exact)), max(np.max(c), np.max(c_exact)))
+ax.grid(True)
+
+cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
+c_exactplot, = ax.plot(x[0], c_exact[0],label = 'Exact solution')
+
+def update(value):
+        ti = np.abs(times-value).argmin()
+        cplot.set_xdata(x[ti])
+        cplot.set_ydata(c[ti])
+        c_exactplot.set_xdata(x[ti])
+        c_exactplot.set_ydata(c_exact[ti])
+        plt.draw()
+        
+sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+        valinit=min(times), valfmt='%3.1f',
+        fc='#999999')
+slider.drawon = False
+slider.on_changed(update)
+ax.legend(loc=0)
+plt.show()
+