correct analytical solution

moving_domain/moving_heat_equation_analytical.py

@@ -7,17 +7,15 @@ df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
 df.parameters['form_compiler']['optimize'] = True
 
 
 # parameters
 # parameters
-k = 1.0
+alpha = 1.0
 T = 1
 T = 1
-
 dt = 0.001
 dt = 0.001
-L = 1
+L0 = 1
 D = 1.0
 D = 1.0
 Nx = 2000
 Nx = 2000
 Nt = 1000
 Nt = 1000
 t = np.linspace(0,T,Nt)
 t = np.linspace(0,T,Nt)
 
 
-
 # diffusion and advection
 # diffusion and advection
 def diffusion(c, tc):
 def diffusion(c, tc):
         return (D * df.inner(c.dx(0), tc.dx(0)))
         return (D * df.inner(c.dx(0), tc.dx(0)))
@@ -27,17 +25,17 @@ def advection(c, tc, v):
         return (df.inner((u*c).dx(0),tc))
         return (df.inner((u*c).dx(0),tc))
 
 
 # create mesh
 # create mesh
-mesh = df.IntervalMesh(Nx, 0, L)
+mesh = df.IntervalMesh(Nx, 0, L0)
 x = mesh.coordinates()
 x = mesh.coordinates()
 # v = df.Constant(1.0)
 # v = df.Constant(1.0)
-v = df.Expression('k*x[0]',k=k,degree=1)
+v = df.Expression('alpha*x[0]',alpha = alpha,degree=1)
 
 
 # create function space
 # create function space
 conc_element = df.FiniteElement('P', mesh.ufl_cell(), 1)
 conc_element = df.FiniteElement('P', mesh.ufl_cell(), 1)
 function_space = df.FunctionSpace(mesh,conc_element)
 function_space = df.FunctionSpace(mesh,conc_element)
 
 
 # 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(pi*x[0]/L0)', pi = np.pi, L0 = L0,degree=1),function_space)
 c0_array = c0.compute_vertex_values(mesh)
 c0_array = c0.compute_vertex_values(mesh)
 
 
 # define variational problem
 # define variational problem
@@ -71,10 +69,10 @@ for n in progressbar.progressbar(range(1, len(t))):
 c_exact= np.zeros((len(t), len(x)))
 c_exact= np.zeros((len(t), len(x)))
 
 
 for i in range(len(t)):
 for i in range(len(t)):
-        xprime = x_array[0]*np.exp(-k*t[i])/L
-        tprime = (D/(2*k*L**2))*(1-np.exp(-2*k*t[i]))
-        # int = ((2*k**2*L**2)/(D**2))*(np.exp(-4*k*t[i])-1)
-        c_exact[i] = np.exp(-k*t[i]) * (1 + 0.2*np.cos(np.pi*xprime)*np.exp(-np.pi**2*tprime))
+        xprime = x_array[0] / (L0 * np.exp(alpha * t[i]))
+        tprime = ( D / (2 * alpha * L0**2)) * ( 1 - np.exp(-2 * alpha * t[i]))
+        int = np.exp(-alpha * t[i])
+        c_exact[i] = int * (1 + 0.2 * np.cos(np.pi * xprime) * np.exp(-np.pi**2 * tprime))
 
 
 # plot c(x,t) computed numerically
 # plot c(x,t) computed numerically
 fig, ax_comp = plt.subplots(1,1,figsize=(8,6))
 fig, ax_comp = plt.subplots(1,1,figsize=(8,6))
@@ -83,7 +81,7 @@ ax_comp.set_ylabel(r'$c(x,t)$')
 ax_comp.set_xlim(np.min(x_array)-1, np.max(x_array)+1)
 ax_comp.set_xlim(np.min(x_array)-1, np.max(x_array)+1)
 ax_comp.set_ylim(np.min(c_array)-1, np.max(c_array)+1)
 ax_comp.set_ylim(np.min(c_array)-1, np.max(c_array)+1)
 cplot, = ax_comp.plot(x_array[0], c0_array)
 cplot, = ax_comp.plot(x_array[0], c0_array)
-c_exactplot, = ax_comp.plot(x_array[0], c_exact[0], 'ro', markevery=50)
+c_exactplot, = ax_comp.plot(x_array[0], c_exact[0], 'ro', markersize=3, markevery=50)
 
 
 def update(value):
 def update(value):
         ti = np.abs(t-value).argmin()
         ti = np.abs(t-value).argmin()