Steady diffusion-reaction with mesh refinement

\[\begin{split} \mathbb{S}_u \begin{cases} \Omega = [0, L_x] \times [0, L_y] \\ u_{\text{D}}(x,y=L_y)=1 \\ u_{\text{D}}(x,y=0)=0 \\ u_{\text{D}}(x=0,y) \\ u_{\text{D}}(x=L_x,y)=0 \\ \mathsf{D}=D\mathsf{I} \\ R = -1 \\ J = 0 \end{cases} \end{split}\]

from typing import Callable
import numpy as np
from lucifex.fem import Function, Constant, cross_section_grid
from lucifex.mesh import rectangle_mesh, refine_mesh, deform_mesh
from lucifex.utils.fenicsx_utils import cell_sizes, number_of_entities
from lucifex.solver import bvp, BoundaryConditions, BoundaryValueProblem
from lucifex.pde.diffusion import steady_diffusion_reaction
from lucifex.plt import plot_mesh, plot_colormap, plot_line, create_multifigure, save_figure

from dolfinx.mesh import create_mesh


def stretch_transform(
    s: float,
    Lx: float | tuple[float, float],
) -> Callable[[float], float]:
    """
    `(x - x₋) / (x₊ - x₋) –> (x - x₋)ˢ / (x₊ - x₋)ˢ`
    """
    if isinstance(Lx, float):
        Lx = (0.0, Lx)
    x_min, x_max = Lx
    Lx = x_max - x_min
    return lambda x: x_min + Lx * ((x - x_min) / Lx) ** s


def shishkin_transform(
    xi: float,
    x0: float,
    Lx: float | tuple[float, float],
) -> Callable[[float], float]:
    """
    `xi ≷ x0` to refine `x ≷ x0`
    """
    if isinstance(Lx, float):
        Lx = (0.0, Lx)
    x_min, x_max = Lx
    x_step = x_min + x0
    m_lt = xi / x0
    m_gt = (x_max - xi) / (x_max - x0)
    return lambda x: (
        (x_min + m_lt * x) * (x <= x_step) 
        +  (x_min + m_lt * x0  + m_gt * (x - x_step)) * (x > x_step)
    )


def create_solver(
    Lx: float,
    Ly: float,
    Nx: int,
    Ny: int,
    cell: str,
    d: float, 
    r: float,
    j: float,
    refinement: str | int | None,
) -> BoundaryValueProblem:
    mesh = rectangle_mesh(Lx, Ly, Nx, Ny, cell)
    bcs = BoundaryConditions(
        ("dirichlet", lambda x: x[1] - Ly, 1.0),
        ("dirichlet", lambda x: x[1], 0.0),  
        ("dirichlet", lambda x: x[0] - Lx, 0.0),  
        ("dirichlet", lambda x: x[0], 0.0),  
    )
    degree = 1
    if refinement is None:
        pass
    elif isinstance(refinement, int):
        degree = refinement
    elif refinement == 'refine':
        marker = lambda x: x[1] > 0.75 * Ly
        mesh = refine_mesh(mesh, marker, n_stop=2)
    elif refinement == 'stretch':
        mesh = deform_mesh(mesh, stretch_transform(0.5, Ly), axis=1)
    elif refinement == 'shishkin':
        mesh = deform_mesh(mesh, shishkin_transform(0.75 * Ly, 0.5 * Ly, Ly), axis=1)
    else:
        raise ValueError(refinement)

    d = Constant(mesh, d)
    r = Constant(mesh, r)
    j = Constant(mesh, j)
    u = Function((mesh, 'P', degree), name='u')
    return bvp(steady_diffusion_reaction, bcs)(u, d, r, j)


Lx = 1.0
Ly = 1.0
Nx = 10
Ny = 10
cell = 'right'
elem = ('P', 1)
d = 1e-3
r = -1.0
j = 0.0

refinement_opts = (None, 'shishkin', 'refine', 'stretch', 2)
solutions: dict[str, Function] = {}

for rfn in refinement_opts:
    u_solver = create_solver(Lx, Ly, Nx, Ny, cell, d, r, j, rfn)
    u_solver.solve()
    solutions[rfn] = u_solver.solution

dy_bl = np.sqrt(d)
y_lims = (-0.3, 1.0)

for rfn, u in solutions.items():
    mesh = u.function_space.mesh
    n_cells = number_of_entities(mesh, 2)
    hmin = np.min(cell_sizes(mesh, 'hmin'))
    hmax = np.max(cell_sizes(mesh, 'hmax'))
    print(f"Refinement option '{rfn}':\n hmin, hmax, ncells = {hmin:.2f}, {hmax:.2f}, {n_cells}")
    u_cs, x_value = cross_section_grid(u, 'x', 0.5)
    deg = u.function_space.ufl_element().degree()
    mfig, axs, _ = create_multifigure(n_cols=2)
    plot_colormap(mfig, axs[0], u, title='$u$')
    plot_mesh(mfig, axs[0], mesh, color='cyan', linewidth=0.75)
    plot_line(mfig, axs[1], u_cs, x_label='$y$', y_label=f'$u(x={x_value})$', title=f'P$_{deg}$')
    axs[1].vlines([Ly - dy_bl], *y_lims, color='gray', ls='dashed', lw=0.75),
    axs[1].set_ylim(*y_lims)
    save_figure(f'u(x,y)_{rfn}', thumbnail=(rfn == 'shishkin'))(mfig)

Refinement option 'None':
 hmin, hmax, ncells = 0.10, 0.14, 200
Refinement option 'shishkin':
 hmin, hmax, ncells = 0.05, 0.18, 200
Refinement option 'refine':
 hmin, hmax, ncells = 0.02, 0.14, 910
Refinement option 'stretch':
 hmin, hmax, ncells = 0.05, 0.33, 200
Refinement option '2':
 hmin, hmax, ncells = 0.10, 0.14, 200

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.