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.
make it const. also, I think the name is poorly chosen. Maybe local_dof_index?