This is a three-dimensional problem with billions of degrees of freedom
Then MUMPS and direct solvers are clearly not feasible. See, e.g., How to choose the optimal solver for a PDE problem? - #2 by nate.
See also the Stokes demo which shows some basic iterative solution methods. I also have a very basic example here.
The specific construction of an iterative solver and appropriate preconditioner will derive from your problem’s formulation. This is a vast research problem and would likely require consultation of research literature.