DG steady advection-diffusion-reaction with a skew velocity

import numpy as np
from ufl import SpatialCoordinate
from lucifex.mesh import rectangle_mesh
from lucifex.fem import Function, Constant
from lucifex.solver import bvp, BoundaryConditions, BoundaryValueProblem
from lucifex.plt import plot_colormap, save_figure
from lucifex.utils.fenicsx_utils import is_continuous_lagrange
from lucifex.pde.advection_diffusion import (
    steady_advection_diffusion, 
    dg_steady_advection_diffusion,
)

def create_solver(
    Lx: float,
    Ly: float,
    Nx: int,
    Ny: int,
    a,
    u_inflow,
    y_inflow,
    theta,
    d: float, 
    k: float,
    element: tuple[str, int],
    bc_type_upper_right: str,
    alpha: float = 10.0,
) -> BoundaryValueProblem:
    mesh = rectangle_mesh(Lx, Ly, Nx, Ny)
    x = SpatialCoordinate(mesh)
    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')
    r = 1 + k * x[0] * x[1]
    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_upper_right, lambda x: x[0] - Lx, 0.0),
        (bc_type_upper_right, lambda x: x[1] - Ly, 0.0),
    )
    u = Function((mesh, *element), name='u')
    if is_continuous_lagrange(u.function_space):
        return bvp(steady_advection_diffusion, bcs)(u, a, d, r=r)
    else:
        return bvp(dg_steady_advection_diffusion)(u, alpha, a, d, r=r, bcs=bcs)
    

Lx = 1.0
Ly = 1.0
Nx = 50
Ny = 50

a = 1.0
u_inflow = 1.0
y_inflow = 0.5
theta = 30.0
d = 5e-4
k = 2.0

elem_opts = (
    ('DP', 1),
    ('DP', 2),
    ('P', 1),
    ('P', 2),
)
bc_upper_right = 'neumann'

solutions: dict[tuple[str, int]: Function] = {}
for elem in elem_opts:
    solver = create_solver(Lx, Ly, Nx, Ny, a, u_inflow, y_inflow, theta, d, k, elem, bc_upper_right)
    solver.solve()
    solutions[elem] = solver.solution

for elem, u in solutions.items():
    fam, deg = elem
    fig, ax = plot_colormap(u, title=f'{fam}$_{deg}$\n$u$')
    save_figure(f'u(x,y)_{fam}{deg}', thumbnail=(elem == ('DP', 1)))(fig)