mth mode cosine initial condition

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, tstop):
+def advection_diffusion(Nx, L, Nt, tmax, D, tstop, 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)
     SFS = df.FunctionSpace(mesh, 'P', 1)
     SFS = df.FunctionSpace(mesh, 'P', 1)
@@ -20,7 +20,7 @@ def advection_diffusion(Nx, L, Nt, tmax, D, tstop):
 
 
     # initial condition
     # initial condition
     c0 = df.Function(SFS)
     c0 = df.Function(SFS)
-    c0.interpolate(df.Expression('1 + 0.2 * cos(pi*x[0]/L)', pi=np.pi, L=L, degree=1))
+    c0.interpolate(df.Expression('1 + 0.2 * cos(m * pi*x[0]/L)', m = m, pi=np.pi, L=L, degree=1))
 
 
     # arrays
     # arrays
     c_array = np.zeros((Nt+1, Nx+1))
     c_array = np.zeros((Nt+1, Nx+1))
@@ -50,7 +50,7 @@ def advection_diffusion(Nx, L, Nt, tmax, D, tstop):
     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.001, 1, 0.001, 50, 0.01
+dx, L, dt, tmax, D, m = 0.01, 1, 0.01, 100, 0.01, 2
 alpha0, tstop = 0.01, 3
 alpha0, tstop = 0.01, 3
 Nx = int(L/dx)
 Nx = int(L/dx)
 Nt = int(tmax/dt)
 Nt = int(tmax/dt)
@@ -58,7 +58,7 @@ Nt = int(tmax/dt)
 times = np.linspace(0, tmax, Nt+1)
 times = np.linspace(0, tmax, Nt+1)
 
 
 # numerical solution
 # numerical solution
-c, x = advection_diffusion(Nx, L, Nt, tmax, D, tstop)
+c, x = advection_diffusion(Nx, L, Nt, tmax, D, tstop, m)
 
 
 # analytical solution
 # analytical solution
 c_exact = np.zeros((Nt+1, Nx+1))
 c_exact = np.zeros((Nt+1, Nx+1))
@@ -67,13 +67,13 @@ for j in range(0, len(times)):
         # diffusion-advection on moving domain with velocity alpha0*x
         # 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]))
-        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))
     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 * tstop)
         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))
+        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))
 fig, ax = plt.subplots(1,1,figsize=(8,6))