SUPG steady advection-diffusion with a skew velocity

Donea, J. & Huerta, A. (2003). Finite Element Methods for Flow Problems. \(\S 2.6.4\)

\[\begin{split} \mathbb{S}_u \begin{cases} \Omega = [0, 1]\times[0,1] \\ u_{\text{D}}(x=0,y<y_{\text{in}}) = 0 \\ u_{\text{D}}(x=0,y \geq y_{\text{in}}) = u_{\text{in}} \\ u_{\text{D}}(x,y=0) = 0 \\ u_{\text{N}}(x,y=1)=u_{\text{D}}(x=1,y) = 0 ~\text{or}~ u_{\text{N}}(x,y=1)=u_{\text{N}}(x=1,y) = 0\\ \textbf{a}=a\cos\theta\textbf{e}_x + a\sin\theta\textbf{e}_y \\ \mathsf{D}=D\mathsf{I} \\ R=0 \\ J=0 \\ \end{cases} \end{split}\]

import numpy as np
from lucifex.mesh import rectangle_mesh
from lucifex.fem import Function, Constant
from lucifex.fdm import peclet
from lucifex.solver import BoundaryValueProblem, bvp, BoundaryConditions
from lucifex.plt import plot_colormap, save_figure
from lucifex.utils.py_utils import nested_dict
from lucifex.pde.advection_diffusion import steady_advection_diffusion

def create_solver(
    tau: str | None,
    Lx: float,
    Ly: float,
    Nx: int,
    Ny: int,
    u_inflow: float,
    y_inflow: float,
    bc_type: str,
    d: float, 
    a: float,
    theta: float,
) -> BoundaryValueProblem:
    mesh = rectangle_mesh(Lx, Ly, Nx, Ny)
    ax = a * np.cos(np.radians(theta))
    ay = a * np.sin(np.radians(theta))
    a = Constant(mesh, (ax, ay), name='a')
    d = Constant(mesh, d, name='d')
    u = Function((mesh, 'P', 1), name='u')
    bcs = BoundaryConditions(
        ('dirichlet', lambda x: x[1], 0.0),
        ('dirichlet', lambda x: np.isclose(x[0], 0) & (x[1] < y_inflow), 0.0),
        ('dirichlet', lambda x: np.isclose(x[0], 0) & (x[1] >= y_inflow), u_inflow),
        (bc_type, lambda x: x[0] - Lx, 0.0),
        (bc_type, lambda x: x[1] - Ly, 0.0),
    )
    u_solver = bvp(steady_advection_diffusion, bcs=bcs)(
        u, a, d, j=Constant(mesh, 0.0), tau=tau,
    )
    return u_solver


Lx = 1.0
Ly = 1.0
Nx = 50
Ny = 50
u_inflow = 1.0
y_inflow = 0.2
d = 1e-3
a = 1.0
theta = 30.0

hx = Lx / Nx
hy = Ly / Ny
Pe_x = peclet(hx, a, d)
Pe_y = peclet(hy, a, d)

tau_opts = (None, 'coth', 'coth_infty')
bc_opts = ('neumann', 'dirichlet')
solutions: dict[str | None, dict[str, Function]] = nested_dict(depth=2)
for tau in tau_opts:
    for bc in bc_opts:
        u_solver = create_solver(tau, Lx, Ly, Nx, Ny, u_inflow, y_inflow, bc, d, a, theta)
        u_solver.solve()
        solutions[tau][bc] = u_solver.solution

for bc in bc_opts:
    for tau in tau_opts:
        u = solutions[tau][bc]
        if tau is None:
            tau = 'unstabilized'
        if tau == 'coth_infty':
            tau = 'coth$_\infty$'
        title = f'{bc.capitalize()}-{tau}\n$u$'
        fig, ax = plot_colormap(u, title=title)
        thumbnail = (bc == 'neumann' and tau == 'coth')
        save_figure(f'u(x,y)_{bc}_{tau}', thumbnail=thumbnail)(fig)