  A parallel implementation of a compatible discretization scheme for
steady-state Stokes problems is presented in this work. The scheme uses
generalized moving least squares to generate differential operators and apply
boundary conditions. This meshless scheme allows a high-order convergence for
both the velocity and pressure, while also incorporates finite-difference-like
sparse discretization. Additionally, the method is inherently scalable: the
stencil generation process requires local inversion of matrices amenable to GPU
acceleration, and the divergence-free treatment of velocity replaces the
traditional saddle point structure of the global system with elliptic diagonal
blocks amenable to algebraic multigrid. The implementation in this work uses a
variety of Trilinos packages to exploit this local and global parallelism, and
benchmarks demonstrating high-order convergence and weak scalability are
provided.
