Commit 9f634beb by Jigyasa Watwani

first attempt

parent c3ff475c
Showing with 85 additions and 0 deletions
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
# parameters
tmax, dt, L, dx, b, m, D = 10, 0.01, 1, 0.01, 0.01, 2, 0.01
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/(L + b *t)*x[0]', 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 progressbar.progressbar(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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment