works for some h, region of stability?

euler_error/euler_error_moving_domain_diffusion_advection.py

+# moving domain diffusion advection equation
+# boundary condition: diffusive flux at the boundaries [0,1] is zero
+# initial condition c(x,0) = 1 + 0.2*cos(pi x) ; satisfies the boundary condition
+# exact solution at steady state: c(x) = c(0) e^(vx/D)
+
+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 = df.IntervalMesh(Nx, 0, L)
+    function_space = df.FunctionSpace(mesh, 'P', 1)
+    c, tc = df.Function(function_space), df.TestFunction(function_space)
+    times = np.linspace(0, tmax, Nt+1)
+    dt = times[1] - times[0]
+
+    c0 = df.interpolate(df.Expression('1 + 0.2*cos(pi*x[0]/L)', pi = np.pi, L=L, degree=1), function_space)
+
+    # form
+    u = df.interpolate(v, function_space)
+    form = (df.inner((c-c0)/dt, tc)
+            + df.inner((u*c).dx(0), tc)
+            + D * df.inner(tc.dx(0), c.dx(0)))*df.dx(mesh)
+
+    for _ in progressbar.progressbar(range(1, len(times)+1)):
+        df.solve(form == 0, c)
+        df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
+        x = mesh.coordinates()
+        c0.assign(c)
+
+    return [c.compute_vertex_values(mesh), x]
+
+Nx, L, tmax, D, k = 2000, 1, 5, 1, 1
+v = df.Expression('k*x[0]', k=k, degree=1)
+nt_array = np.array([50, 100, 200, 400, 800, 1600, 3200, 6400, 12800])
+dt_array = tmax/(nt_array - 1)
+error_in_c = np.zeros(len(nt_array))
+
+for i in range(len(nt_array)):
+    x = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[-1]
+    ss_c_exact = 2 * np.exp(-k*tmax) * (1 + 0.2*np.cos(np.pi*x*np.exp(-k*tmax)) * np.exp(-np.pi**2 * (D/(2*k)) *(1-np.exp(-2*k*tmax))))
+    ss_c = advection_diffusion(Nx, L, nt_array[i], tmax, v, D)[0]
+    error_in_c[i] = np.max(np.abs(ss_c - ss_c_exact))
+
+plt.scatter(dt_array, error_in_c)
+plt.show()

moving_domain/moving_domain_heat_equation_2.py

-import dolfin as df
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider
-import progressbar
-
-# bother me only if there is an error
-df.set_log_level(df.LogLevel.ERROR)                             
-df.parameters['form_compiler']['optimize'] = True
-
-# parameters
-Nx = 2000
-L = 2*np.pi
-dt = 0.005
-T = 5
-D = 1.0
-k = 0.5
-times = np.arange(0, T, dt)
-
-# diffusion and advection
-def diffusion(func, testfunc, D):
-    return D * df.inner(func.dx(0), testfunc.dx(0))
-
-def advection(func, testfunc, vel):
-    return df.inner((vel * func).dx(0), testfunc)
-
-# mesh
-mesh = df.IntervalMesh(Nx, 0, L)
-x = mesh.coordinates()
-
-# function space
-conc_element = df.FiniteElement('P', mesh.ufl_cell(), 1)
-function_space = df.FunctionSpace(mesh, conc_element)
-
-# define velocity
-v = df.Constant(0.5)
-# v = df.interpolate(df.Expression('1.0', degree=1), function_space)
-
-# initial condition
-c0 = df.interpolate(df.Expression('1 + 0.2*cos(x[0])',degree=1), function_space)
-
-# define function, test function
-c = df.Function(function_space)
-tc = df.TestFunction(function_space)
-
-# weak form
-form = (df.inner((c - c0)/dt, tc) 
-        + diffusion(c, tc, D)
-        + advection(c, tc, v)
-        )*df.dx
-
-# define the arrays
-c_array = np.zeros((len(times), len(x)))
-c_array[0] = c0.compute_vertex_values(mesh)
-
-x_array = np.zeros((len(times), mesh.num_vertices()))
-x_array[0] = mesh.coordinates()[:,0]
-
-c_tot = np.zeros_like(times)
-c_tot[0] = df.assemble(c0 * df.dx(mesh))
-
-# time stepping
-for i in progressbar.progressbar(range(1,len(times))):
-    df.solve(form == 0, c)
-    c_tot[i] = df.assemble(c * df.dx(mesh))
-    c_array[i] = c.compute_vertex_values(mesh)
-    df.ALE.move(mesh, df.Expression('v*dt', v=v, dt=dt, degree=1))
-    x_array[i] = mesh.coordinates()[:,0]
-    c0.assign(c)
-
-
-# plotting
-
-# plot c(x,t) vs x for all t
-fig, axc = plt.subplots(1, 1, figsize=(8,6))
-axc.set_xlabel(r'$x$')
-axc.set_ylabel(r'$c(x,t)$')
-
-cplot, = axc.plot(x_array[0], c_array[0], 'r')
-axc.set_xlim(np.min(x_array)-2, np.max(x_array)+2)
-axc.set_ylim(np.min(c_array)-2, np.max(c_array)+2)
-
-def update(value):
-    ti = np.abs(times-value).argmin()
-    cplot.set_xdata(x_array[ti])
-
-    cplot.set_ydata(c_array[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)
-
-# plot ctot vs time
-fig1, ax1 = plt.subplots(1, figsize=(8,6))
-fig1.subplots_adjust(left=0.1, bottom=0.1, top=0.9, right=0.9,
-        wspace=0.0, hspace=0.0)
-ax1.plot(times, c_tot)
-ax1.set_ylim([0,np.max(c_tot)+1])
-ax1.set_xlabel(r'$t$')
-ax1.set_title('Total concentration vs time')
-ax1.set_ylabel(r'$\int_{\Gamma} c$', color='#ff5b00')
-plt.show()
-