Diffusion of a Gaussian in an interval

Diffusion of a Gaussian in an interval#

\[\begin{split} \mathbb{S}_u \begin{cases} \Omega = [0, L_x] \\ u_0(x) = \exp\left(-\frac{(x-x_0)^2}{L_x/100}\right) \\ 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.plt 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)
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)
../../_images/10cb0364c42e8f253d78b1faae8ba1c9a4e7f30e795df7d127d757c35d4b5588.png
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.
y_label_series = [f'$u(t={t:.2f})$' for t in u.time_series]

anim = create_animation(
    plot_line,
    x_label='$x$',
)(u.series, y_label=y_label_series)
anim_path = save_figure('u(x,t)', return_path=True)(anim)

display_animation(anim_path)