SUPG steady advection-diffusion in an interval

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

\[\begin{split} \mathbb{S} \begin{cases} \Omega = [0, 1] & \text{unit interval} \\ u_{\text{D}}(x=0,1)=0 \\ \textbf{a}=a\,\textbf{e}_x \\ \mathsf{D}=D\mathsf{I} \\ R=0 \\ J=1 \\ u_{\text{e}}(x)=\frac{1}{a}\frac{(x-e^{ax/D})}{1-e^{a/D}} & \text{exact solution} \\ Pe= \frac{a}{2DN_x} & \text{local Peclet number} \\ \end{cases} \end{split}\]

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

def create_solver(
    supg: str | None,
    Nx: int,
    a: float,
    d: float, 
    j: float | None = 1.0,
    Lx: float = 1.0,
    dirichlet: tuple[float, float] = (0.0, 0.0),
) -> BoundaryValueProblem:
    mesh = interval_mesh(Lx, Nx)
    a = Constant(mesh, (a, ), name='a')
    d = Constant(mesh, d, name='d')
    if j is not None:
        j = Constant(mesh, j, name='s')
    u = Function((mesh, 'P', 1), name='u')
    bcs = BoundaryConditions(
        ('dirichlet', lambda x: x[0], dirichlet[0]),
        ('dirichlet', lambda x: x[0] - Lx, dirichlet[1])
    )
    u_solver = bvp(steady_advection_diffusion, bcs=bcs)(
        u, a, d, j=j, supg=supg,
    )
    return u_solver

Lx = 1.0
Nx = 10
h = Lx / Nx
a = 1.0
j = 1.0

supg_opts = (None, 'coth', 'upwind')
d_opts = (0.2, 1/18, 1/100)
solutions = nested_dict(dict[float, dict[str | None, Function]], depth=2)
for d in d_opts:
    for supg in supg_opts:
        h = Lx / Nx
        u_solver = create_solver(supg, Nx, a, d, j, Lx)
        u_solver.solve()
        solutions[d][supg] = u_solver.solution

def exact_solution(
    x: np.ndarray, 
    a: float, 
    d: float,
) -> np.ndarray:    
    return (1 / a) * (x - (1 - np.exp(a * x / d)) / (1.0 - np.exp(a / d)))

x = np.linspace(0, Lx, num=500)

for d in d_opts:
    lines = [(x, exact_solution(x, a, d))]
    legend_labels = ['exact']
    for supg in supg_opts:
        lines.append(solutions[d][supg])
        if supg is None:
            supg = 'unstabilized'
        legend_labels.append(supg)
    Pe = peclet(h, a, d)
    fig, ax = plot_line(lines, legend_labels, 'SUPG', x_lims=x, x_label='$x$', y_label='$u$', title=f'$Pe={Pe:.2f}$')
    save_figure(f'Pe={Pe:.2f}', thumbnail=(d == d_opts[-1]))(fig)