Hello, I'm quite new to using Firedrake (although I have used FEniCS in the past). Based on some of the Firedrake demos I wrote a Stokes code in which I'm using a mixed BDM/DG or RT/DG function space (see attached). In this code I compare my numerical solution to an exact solution on different grids and computed the rates of convergence. This runs fine in serial but in parallel the rates of convergence are not as I expect. For example, in serial, rates of convergence using BDM(order 2)/DG (order 1) are: Err_u r Err_p r div(u) r 8.21e-05 2.97 9.98e-03 1.95 5.17e-09 -0.80 1.04e-05 2.98 2.54e-03 1.97 7.75e-09 -0.59 1.31e-06 2.99 6.41e-04 1.99 5.43e-09 0.51 but in parallel (mpiexec -np 4 python stokes_hdiv.py) Err_u r Err_p r div(u) r 8.18e-05 3.48 9.99e-03 2.07 2.99e-09 0.22 5.07e-05 0.69 2.91e-03 1.78 5.69e-09 -0.93 3.48e-05 0.54 1.34e-03 1.12 7.88e-09 -0.47 Am I missing something in my code to deal with runs in parallel? Using a standard Taylor-Hood P2-P1 finite element (code is also attached), I am getting the correct rates of convergence for this Stokes problem. Thanks, Sander