Diffusion of a Gaussian in an interval#
\[\begin{split}
\mathbb{S}
\begin{cases}
\Omega = [0, L_x] \\
u_{\text{D}}(x=0) = 0 \\
u_{\text{D}}(x=L_x) = 0 \\
\mathsf{D}=D\mathsf{I} \\
\end{cases}
\end{split}\]
import numpy as np
from lucifex.mesh import interval_mesh, mesh_boundary
from lucifex.fdm import BE, FunctionSeries, diffusive_timestep
from lucifex.fem import Constant
from lucifex.solver import ibvp, BoundaryConditions
from lucifex.viz import plot_line, create_animation, save_figure, display_animation
from lucifex.pde.diffusion import diffusion
Nx = 100
Lx = 2.0
mesh = interval_mesh(Lx, 100)
boundary = mesh_boundary(
mesh,
{
"left": lambda x: x[0],
"right": lambda x: x[0] - Lx,
},
)
dt = Constant(mesh, 0.01, name='dt')
d = Constant(mesh, 1.0, name='d')
D_diff = BE
h = Lx/ Nx
dt_diffusive = diffusive_timestep(d, h)
courant = dt.value / dt_diffusive
ics = lambda x: np.exp(-(x[0] - Lx/2)**2 / (0.01 * Lx))
bcs = BoundaryConditions(
("dirichlet", boundary['left', 'right'], 0.0),
)
u = FunctionSeries((mesh, 'P', 1), name='u', store=1)
w = FunctionSeries((mesh, 'P', 1), name='w', store=1)
u_solver = ibvp(diffusion, ics, bcs)(u, dt, d, D_diff)
t = 0.0
t_stop = 10.0
n = 0
n_stop = 10
while n < n_stop and t < t_stop:
u_solver.solve()
u.forward(t)
t += dt.value
n += 1
time_slice = slice(0, None)
legend_title = f"$\mathcal{{D}}_{{\mathsf{{D}}, u}}=\\mathrm{{{str(D_diff)}}}$\n$C_{{\mathsf{{D}}}}={courant}$\n\n$t$"
legend_labels = [f'{t:.2f}' for t in u.time_series[time_slice]]
fig, ax = plot_line(
u.series[time_slice], legend_labels, legend_title,
cyc='jet', x_label='$x$', y_label='$u$',
)
save_figure('u(x,t)', thumbnail=True)(fig)
title_series = [f'$t={t:.2f}$' for t in u.time_series]
anim = create_animation(
plot_line,
x_label='$x$',
y_label='$u(x,t)$',
)(u.series, title=title_series)
anim_path = save_figure('u(x,t)', get_path=True)(anim)
display_animation(anim_path)
The Kernel crashed while executing code in the current cell or a previous cell.
Please review the code in the cell(s) to identify a possible cause of the failure.
Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info.
View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.