Hi all, So given this Darcy problem: a = (dot(v, w) - p*div(w) - div(v)*q)*dx L = f*q*dx where f = "12*pi*pi*sin(pi*x[0]*2)*sin(pi*x[1]*2)*sin(2*pi*x[2])" and with these solver options: 'ksp_type': 'gmres', #'ksp_monitor_true_residual': True, 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'pc_fieldsplit_schur_fact_type': 'upper', 'pc_fieldsplit_schur_precondition': 'selfp', 'fieldsplit_0_ksp_type': 'preonly', 'fieldsplit_0_pc_type': 'bjacobi', 'fieldsplit_0_sub_pc_type': 'ilu', 'fieldsplit_1_ksp_type': 'preonly', 'fieldsplit_1_pc_type': 'hypre', 'fieldsplit_1_pc_hypre_type': 'boomerang' And without specifying any boundaries (my forcing function f is chosen so that the pressure is homogeneous on the boundary), RT0 and BDM works beautifully. However, with Taylor-Hood elements and VMS (which is equal order CG1 plus a least squares stabilization), the velocity works but the pressure solution is screwed up. Why is that? Thanks, Justin