working numerical solution, non-jumpy analytical solution

moving_domain/moving_time_dependent_alpha.py

@@ -7,7 +7,7 @@ import progressbar
 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, D):
+def advection_diffusion(Nx, L, Nt, tmax, D, tstop):
     # mesh, function space, function, test function
     # mesh, function space, function, test function
     mesh = df.IntervalMesh(Nx, 0, L)
     mesh = df.IntervalMesh(Nx, 0, L)
     SFS = df.FunctionSpace(mesh, 'P', 1)
     SFS = df.FunctionSpace(mesh, 'P', 1)
@@ -29,70 +29,69 @@ def advection_diffusion(Nx, L, Nt, tmax, D):
     c_array[0] = c0.compute_vertex_values(mesh)
     c_array[0] = c0.compute_vertex_values(mesh)
 
 
     # velocity
     # velocity
-    alpha = df.Expression('times < 1 ? alpha0 : 0', times = 0, alpha0 = alpha0, degree=0)
-
-    v = df.Expression('alpha*x[0]', alpha=alpha, degree=1)
-    u = df.interpolate(v, SFS)
+    u = df.Expression('(t < tstop ? alpha0 : 0)*x[0]', alpha0 = alpha0, tstop=tstop, t=0, degree=0)
+    uh = df.project(u, SFS)
 
 
     # form
     # form
     cform = (df.inner((c - c0)/dt, tc) 
     cform = (df.inner((c - c0)/dt, tc) 
             + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
             + D * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
-            + df.inner((u*c).dx(0), tc) )* df.dx
+            + df.inner((uh*c).dx(0), tc) )* df.dx
 
 
     # solve
     # solve
     for i in progressbar.progressbar(range(1, Nt+1)):
     for i in progressbar.progressbar(range(1, Nt+1)):
+        u.t = times[i]
+        uh.assign(df.project(u, SFS))
         df.solve(cform == 0, c)
         df.solve(cform == 0, c)
         c_array[i] = c.compute_vertex_values(mesh)
         c_array[i] = c.compute_vertex_values(mesh)
         c0.assign(c)       
         c0.assign(c)       
-        df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
+        df.ALE.move(mesh, df.project(uh*dt, SFS))
         x_array[i] = mesh.coordinates()[:,0]
         x_array[i] = mesh.coordinates()[:,0]
-        alpha.times = times[i]
 
 
     return c_array, x_array
     return c_array, x_array
 
 
 # plot c(x,t) numerical and analytical for given dt
 # plot c(x,t) numerical and analytical for given dt
-dx, L, dt, tmax, D = 0.01, 1, 0.01, 10, 0.1
-alpha0, tc = 1, 1
+dx, L, dt, tmax, D = 0.001, 1, 0.001, 50, 0.01
+alpha0, tstop = 0.01, 3
 Nx = int(L/dx)
 Nx = int(L/dx)
 Nt = int(tmax/dt)
 Nt = int(tmax/dt)
 
 
 times = np.linspace(0, tmax, Nt+1)
 times = np.linspace(0, tmax, Nt+1)
-x = advection_diffusion(Nx, L, Nt, tmax, D)[1]
 
 
 # numerical solution
 # numerical solution
-c = advection_diffusion(Nx, L, Nt, tmax, D)[0]
+c, x = advection_diffusion(Nx, L, Nt, tmax, D, tstop)
 
 
 # analytical solution
 # analytical solution
 c_exact = np.zeros((Nt+1, Nx+1))
 c_exact = np.zeros((Nt+1, Nx+1))
 for j in range(0, len(times)):
 for j in range(0, len(times)):
-    if times[j] <= tc:
-        # diffusion on moving domain with velocity alpha0*x
+    if times[j] <= tstop:
+        # diffusion-advection on moving domain with velocity alpha0*x
         alpha = alpha0
         alpha = alpha0
         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]))
         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))
         c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
     else: 
     else: 
         # diffusion on fixed domain of length L = L0*exp(alpha0 * tc) with initial condition to be the profile at tc
         # 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 * tc)
-        beta = -D * np.pi**2*times[j]/l**2 - np.pi**2 * D/(2*l**2*alpha0) * (1 - np.exp(-2 *alpha0 *tc))
-        c_exact[j] = np.exp(-alpha0 * tc) * (1 + 0.2 * np.exp(beta) * np.cos(np.pi * x[j]/l))
+        l = L * np.exp(alpha0 * tstop)
+        beta = -D * np.pi**2*(times[j] - tstop)/l**2 - np.pi**2 * D/(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(np.pi * x[j]/l))
 
 
 
 
 fig, ax = plt.subplots(1,1,figsize=(8,6))
 fig, ax = plt.subplots(1,1,figsize=(8,6))
 ax.set_xlabel(r'$x$')
 ax.set_xlabel(r'$x$')
 ax.set_ylabel(r'$c(x,t)$')
 ax.set_ylabel(r'$c(x,t)$')
-ax.set_xlim(np.min(x), np.max(x)+2)
+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.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')
 cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
-# c_exactplot, = ax.plot(x[0], c_exact[0],label = 'Exact solution')
+c_exactplot, = ax.plot(x[0], c_exact[0],label = 'Exact solution')
 
 
 def update(value):
 def update(value):
         ti = np.abs(times-value).argmin()
         ti = np.abs(times-value).argmin()
         cplot.set_xdata(x[ti])
         cplot.set_xdata(x[ti])
         cplot.set_ydata(c[ti])
         cplot.set_ydata(c[ti])
-        # c_exactplot.set_xdata(x[ti])
-        # c_exactplot.set_ydata(c_exact[ti])
+        c_exactplot.set_xdata(x[ti])
+        c_exactplot.set_ydata(c_exact[ti])
         plt.draw()
         plt.draw()
         
         
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])