when computing the RHS with do_mom_diff==1, multiply U_new by rho_new (rather than rho_old); i.e. use
rho^(n+1) Unew = rho^(n+1) [(rho^n U^n - dt U dot grad U + ...)/rho^(n+1)]
pass a better initial guess for the solution to the tensor solver. This can change results up to the tolerance of the viscous solve (default relative tol = 1e-10)
This makes two changes: