Improve diffuse tensor velocity (AMReX-Fluids#160)

This makes two changes:
* 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)

Source/Diffusion.cpp

@@ -816,7 +816,7 @@ Diffusion::diffuse_tensor_velocity (Real                   dt,
        const Box& bx  = mfi.tilebox();
        auto const& rhs      = Rhs.array(mfi);
        auto const& unew     = U_new.array(mfi,Xvel);
-       auto const& rho      = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_old_data(State_Type).array(mfi,Density);
+       auto const& rho      = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_new_data(State_Type).array(mfi,Density);
        auto const& deltarhs = (has_delta_rhs) ? delta_rhs->array(mfi,rhsComp) : U_new.array(mfi);
        amrex::ParallelFor(bx, [rhs, unew, rho, deltarhs, has_delta_rhs, dt]
        AMREX_GPU_DEVICE(int i, int j, int k) noexcept
@@ -827,6 +827,9 @@ Diffusion::diffuse_tensor_velocity (Real                   dt,
              if ( has_delta_rhs ) {
                 rhs(i,j,k,n) += deltarhs(i,j,k,n) * dt;
              }
+             // Put unew back since it's a reference to the MF that
+             // ultimately gets used as an initial guess for the solve
+             unew(i,j,k,n) /= rho(i,j,k);
           }
        });
     }