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.