New to FEniCS. Trying to port my MATLAB codes. Porting the basic laminar incompressible was straightforward and I used the driven 2D cavity for verification at Re = 10000 using P2-P1 elements. The problem quickly approaches reference data when a graded mesh with corner refining is done. So for the test I typically use an unstructured mesh. In my experience this smudges the solution for the centerline profiles a bit. With SUPG activated the solution fits the expected profile exactly. It seems I am not allowed to post the images as a new user. So here is the matlab snippet the Re is varied only through mu here. This works as expected for me..
```
% interpolate velocities in x and y directions with shape function N
a1 = Ni.' * U1el; a2 = Ni.' * U2el;
% Fill up matrix A: symmetric shear terms + non-linear velocity
% terms + SUPG terms
supg_test = a1*dNx + a2*dNy;
lap_u_gp = (d2Nx2 + d2Ny2); lap_v_gp = (d2Nx2 + d2Ny2); shap_i = Ni;
visc_xx = mu*((dNx*dNx.') + dNy*dNy.') ...
- taus1*mu*(supg_test*lap_u_gp.');
visc_yy = mu*(dNx*dNx.' + (dNy*dNy.')) ...
- taus1*mu*(supg_test*lap_v_gp.');
conv_x = rho*Re*(shap_i*(a1*dNx + a2*dNy).') + taus1*(supg_test*(a1*dNx + a2*dNy).');
conv_y = rho*Re*(shap_i*(a1*dNx + a2*dNy).') + taus1*(supg_test*(a1*dNx + a2*dNy).');
Kx1e = Kx1e + (visc_xx + conv_x)*wdet;
Kx2e = Kx2e + (visc_xy )*wdet;
Ky1e = Ky1e + (visc_yx )*wdet;
Ky2e = Ky2e + (visc_yy + conv_y)*wdet;
% Fill up matrix B^T: pressure terms + SUPG terms
Kx3e = Kx3e - (dNx*Np.' - taus1*(supg_test*(dNxp).')) * wdet;
Ky3e = Ky3e - (dNy*Np.' - taus1*(supg_test*(dNyp).')) * wdet;
For FEniCS I attempted to exactly replicate it as below in the provided snippet
```
# ----------------------- Stabilization (SUPG/PSPG/LSIC) ------------------
def stabilization_terms(u, p, v, q):
if not (enable_supg or enable_pspg or enable_lsic):
return 0
u_mag = sqrt(dot(u,u) + Constant(u_mag_eps))
A_expr = CellVolume(mesh)
h_lin = 2.0*sqrt(A_expr/np.pi)
h_loc = h_lin / (vel_degree + 1.0)
term_conv = (2.0*u_mag / h_loc)**2
term_diff = (C_I*nu / h_loc**2)**2
tau_m = tau_m_scale / sqrt(term_conv + term_diff)
tau_p = tau_p_scale * tau_m if enable_pspg else Constant(0.0)
tau_c = tau_c_scale * h_loc / u_mag if enable_lsic else Constant(0.0)
R_m = rho*grad(u)*u + grad(p) - nu*div(grad(u))
supg_term = tau_m * inner(grad(v)*u, R_m) if enable_supg else 0
pspg_term = tau_p * inner(grad(q), R_m) if enable_pspg else 0
lsic_term = tau_c * div(v)*div(u) if enable_lsic else 0
return (supg_term + pspg_term + lsic_term)*dx
F += stabilization_terms(u, p, v, q)
For constructing the convective term I had used inner(grad u*u,v) and that worked for me. But for the SUPG term it seems to be making no difference and I get almost the same result as without SUPG (no benefit as I got from my MATLAB code). The difference is slight in reality but it makes a big difference in my actual reacting flows simulations. You can ignore the PSPG or LSIC part since I am not using it at the moment.
In both cases the meshes are nearly the same and they behave equally when it comes to the non stabilized P2-P1 formulation but in the stabilized formulation they are mismatching. Sorry if I made any stupid mistakes, and I did use some ChatGPT while coding FEniCS because I found it too high level
