Flow of a Darcy fluid across an anticline of heterogeneous permeability

\[\begin{split} \mathbb{S}_{\psi} \begin{cases} \Omega = \{(x,y)~:~~R_{\text{inner}}^2 < x^2 + y^2 < R_{\text{outer}}^2~,~y>0\} \\ \partial\Omega_{\text{outer}} = \{(x,y)~:~x^2 + y^2 < R_{\text{outer}}^2\} \\ \partial\Omega_{\text{inner}} = \{(x,y)~:~x^2 + y^2 < R_{\text{inner}}^2\} \\ \partial\Omega_{\text{left}} = \{(x,y)~:~y=0~,~-R_{\text{outer}}<x<-R_{\text{inner}}\} \\ \partial\Omega_{\text{right}} = \{(x,y)~:~y=0~,~R_{\text{inner}}<x<R_{\text{outer}}\} \\ \psi_{\text{D}}\vert_{\partial\Omega_{\text{outer}}}=0 & \text{no-penetration on outer boundary}\\ \psi_{\text{D}}\vert_{\partial\Omega_{\text{left}}}=-u_{\text{left}}(x+R_{\text{outer}}) & \text{inflow on left boundary} \\ \psi_{\text{D}}\vert_{\partial\Omega_{\text{right}}}=u_{\text{left}}(x-R_{\text{outer}}) & \text{outflow on right boundary} \\ \psi_{\text{D}}\vert_{\partial\Omega_{\text{inner}}}=-u_{\text{left}}(R_{\text{outer}}-R_{\text{inner}}) & \text{no-penetration on inner boundary} \\ \mathsf{K}(\textbf{x})=\mathcal{N}(\textbf{x})\mathsf{I} & \text{random isotropic permeability} \\ \mu=1 & \text{constant viscosity} \\ \end{cases} \end{split}\]

import numpy as np
from lucifex.fem import Function, Constant, DofsPerturbation
from lucifex.mesh import annulus_sector_mesh, mesh_boundary
from lucifex.solver import BoundaryConditions, bvp
from lucifex.plt import plot_contours, plot_mesh, plot_colormap, save_figure
from lucifex.pde.darcy import darcy_streamfunction


Rinner = 1.0
Router = 2.0
Nradial = 10
dr = (Router - Rinner) / Nradial
mesh = annulus_sector_mesh(dr, 'triangle', 'anticline')(Rinner, Router, 180)

r2 = lambda x: x[0]**2 + x[1]**2
boundary = mesh_boundary(
    mesh, 
    {
        "inner": lambda x: r2(x) - Rinner**2,
        "outer": lambda x: r2(x) - Router**2,
        "left": lambda x: np.logical_and(np.isclose(x[1], 0.0), x[0] < 0),
        "right": lambda x: np.logical_and(np.isclose(x[1], 0.0), x[0] > 0),
    },
)

mu = Constant(mesh, 1.0)
k_lims = (0.05, 1.0)
k_noise = DofsPerturbation(
    0.0,
    1234,
    k_lims,
    (16, 16)
)
k = k_noise.combine_base_noise((mesh, 'P', 1), name='k')

psi = Function((mesh, 'P', 2), name='psi')
u_in = 1.0
psi_bcs = BoundaryConditions(
    ('dirichlet', boundary['outer'], 0.0),
    ('dirichlet', boundary['left'], lambda x: -u_in * (x[0] + Router)),
    ('dirichlet', boundary['inner'], -u_in * (Router - Rinner)),
    ('dirichlet', boundary['right'], lambda x: u_in * (x[0] - Router)),
)
psi_solver = bvp(darcy_streamfunction, psi_bcs)(psi, k, mu)
psi_solver.solve()

fig, ax = plot_colormap(k, title='$K$', x_label='$x$', y_label='$y$')
plot_mesh(fig, ax, mesh, color='cyan', linewidth=0.5)
save_figure('K(x,y)_mesh')(fig)

fig, ax = plot_colormap(k, title='$K$', x_label='$x$', y_label='$y$')
plot_contours(fig, ax, psi, colors='cyan', levels=10)
save_figure('K(x,y)_streamlines', thumbnail=True)(fig)