Flow of a Darcy fluid through a lens of low permeability

\[\begin{split} \mathbb{S}_{\textbf{u},p} \begin{cases} \Omega = [-L_x/2, L_x/2] \times [-L_y/2, L_y/2] \\ \mathsf{K}=\begin{cases} K^-(x,y)\mathsf{I} & (x/a)^2 + (y/b)^2 < 1 \\ K^+(x,y)\mathsf{I} & \text{otherwise} \end{cases} & \text{lens of low permeability} \\ \mu = 1 \\ u_{\text{E}}(x,y=\pm L_y/2) = 0 & \text{no-penetration on upper and lower boundaries} \\ p_{\text{N}}(x=-L_x/2,y) = p_{\text{in}} & \text{high-pressure left boundary} \\ p_{\text{N}}(x=L_x/2,y) = 0 & \text{low-pressure left boundary} \\ \end{cases} \end{split}\]

import numpy as np
from lucifex.mesh import rectangle_mesh, mesh_boundary
from lucifex.fem import Function, Constant
from lucifex.solver import bvp, BoundaryConditions, OptionsPETSc
from lucifex.utils import get_component_fem_functions
from lucifex.viz import plot_colormap, plot_streamlines, plot_quiver, save_figure
from lucifex.pde.darcy import darcy_incompressible

Lx = 2.0
Ly = 1.0
Nx = 200
Ny = 200
mesh = rectangle_mesh((-Lx/2 , Lx/2), (-Ly/2, Ly/2), Nx, Ny, cell='right')
boundary = mesh_boundary(
    mesh, 
    {
        "left": lambda x: x[0] + Lx/2,
        "right": lambda x: x[0] - Lx/2,
        "lower": lambda x: x[1] + Ly/2,
        "upper": lambda x: x[1] - Ly/2,
    },
)

mu = Constant(mesh, 1.0, 'mu')

u_deg = 1
u_elem = ('BDM', u_deg)
p_elem = ('DP', u_deg - 1)
u_p_elem = [u_elem, p_elem]

u_bcs = BoundaryConditions(
    ('essential', boundary['upper', 'lower'], (0.0, 0.0), 0),
)
p_in = 1.0
p_bcs = BoundaryConditions(
    ('natural', boundary['left'], p_in, 1),
    ('natural', boundary['right'], 0.0, 1),
)

petsc = OptionsPETSc(
    ksp_type='gmres', 
    pc_type='lu', 
    pc_factor_mat_solver_type='mumps',
)

k_minus, k_plus = 0.05, 1.0
a, b = (0.2 * Lx, 0.2 * Ly)
r2 = lambda x: (x[0] / a)**2 + (x[1] / b)**2
k_dcts = Function(
    (mesh, 'DP', 0), 
    lambda x: k_plus + (k_minus - k_plus) * (r2(x) < 1), 
    name='kDP',
)
eps = 0.2
k_cts = Function(
    (mesh, 'P', 1), 
    lambda x: 0.5 * (k_plus - k_minus) * np.tanh((r2(x) - 1) / eps) + 0.5 * (k_plus + k_minus), 
    name='kP',
)

solutions: dict[str, Function] = {}

for k in (k_dcts, k_cts):
    up = Function((mesh, u_p_elem), name="up")
    up_solver = bvp(darcy_incompressible, u_bcs, petsc)(up, k, mu, bcs=p_bcs)
    up_solver.solve()
    solutions[k.name] = up

for k in (k_dcts, k_cts):
    up = solutions[k.name]
    u, p = up.split(('u', 'p'), collapse=True)
    ux, uy = get_component_fem_functions(('P', 1), u, names=('ux', 'uy'))
    fig, ax = plot_colormap(p, title='$p$', x_label='$x$', y_label='$y$')
    plot_streamlines(fig, ax, (ux, uy), density=0.75, color='cyan')
    save_figure(f'{k.name}_p(x,y)', thumbnail=(k is k_dcts))(fig)
    fig, ax = plot_quiver((ux, uy), arrow_slc=7, title="$\mathbf{u}=(u_x, u_y)$")
    save_figure(f'{k.name}_u(x,y)')(fig)
    fig, ax = plot_colormap(k, title='$K$')

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.