fix direct solver with nonlinear schemes

We need to compute and return the correct initial_residual in solve_stokes() when using a direct solver for the nonlinear methods like iterated IMPES to work.
modify comment at the top to indicate what you really want to test here? E.g.,
  # modification of compressibility_iterated_stokes.prm, but using a direct solver. etc bla bla.