Hello,
can you try by defining before entering the loop
u_func = Function(VD)
rho_func = Function(VD)
dMu = Function(VD)
c = Constant(1)
and then inside the loop do
u_func.vector().set_local(u)
rho_func.vector().set_local(rho)
c.assign(Constant(assemble(rho_func*u_func*u_func*dx)))
dmu = (dot(grad(u_func), grad(u_func)) - mu*u_func*u_func)/c
dMu.assign(project(dmu, VD)
also if mu is changing, define it as a Constant and assign it its new value inside the loop