Add compressible Navier-Stokes operator with BR1 viscous fluxes

[Copilot]

Adds a compressible Navier-Stokes (CNS) operator to grudge/models/, analogous to the existing Euler operator, using Bassi-Rebay 1 (BR1) fluxes for the viscous terms. Includes inflow, outflow, and no-slip boundary conditions, and a Poiseuille flow spatial convergence test.

New module: grudge/models/navierstokes.py

BR1 gradient (compute_viscous_gradient):

• Implements σ = −weak_local_grad(q) + M⁻¹ face_mass({q}⊗n)
• Interior faces use averages {q}; boundary faces use prescribed BC values

Viscous flux (compute_ns_viscous_flux):

• Newtonian stress tensor τ_{ij} = μ(∂ⱼuᵢ + ∂ᵢuⱼ − ²⁄₃δᵢⱼ∇·u)
• Fourier heat flux q = −κ∇T, with κ = μcₚ/Pr
• Returns a ConservedEulerField with momentum=τ, energy=τ·u − q

CNSOperator discretizes ∂Q/∂t + ∇·(F_inv − σ_vis) = 0:

cns = CNSOperator(
    dcoll,
    bdry_conditions={
        "left":   InflowBC(exact_state),
        "right":  OutflowBC(),
        "walls":  NoSlipBC(T_wall=T0),
    },
    mu=1e-3, Pr=0.72, gamma=1.4,
)
rhs = cns.operator(actx, t, q)

Boundary conditions: NoSlipBC (isothermal or adiabatic), InflowBC (prescribed callable), OutflowBC (zero-order extrapolation). All implement both the inviscid trace pair and the BR1 gradient BC data interface.

poiseuille_flow(x_vec, *, mu, ...): exact steady Poiseuille solution for use as IC/BC.

New test: test/test_navierstokes_model.py

test_poiseuille_convergence evaluates the CNS operator residual at the exact Poiseuille solution on a channel mesh with all three BC types. Checks that the mass and x-momentum residuals (both analytically zero for this flow) converge at ≥ O(h^{p+0.5}) for polynomial orders 2 and 3, directly measuring spatial discretization accuracy of the viscous operator.