The problem was that WV_map and WQ_map reference (cpu-)local dof numbers, whereas PETSc’s index-set requires global numbers. Secondly, the ghost dofs need to be thrown out. The following tweak fixes the parallel runnability:
W_indexmap = W.dofmap.index_map
globaldofs_u = np.setdiff1d( W_indexmap.local_to_global(np.array(WV_map)), W_indexmap.ghosts ) # getting the global velocity dofs that live on this core
globaldofs_p = np.setdiff1d( W_indexmap.local_to_global(np.array(WQ_map)), W_indexmap.ghosts ) # getting the global pressure dofs that live on this core
IS_u = PETSc.IS().createGeneral(np.array(globaldofs_u,dtype=np.int32), comm=domain.comm).sort()
IS_p = PETSc.IS().createGeneral(np.array(globaldofs_p,dtype=np.int32), comm=domain.comm).sort()
I’ll post a small scaling study when I find the time