ValueError in my code

Hello, i have task:
\frac{ d \rho f}{dt} \, + \, \nabla \cdot (\rho f \; \vec{v}) \, = \, \nabla \cdot (\rho D \nabla f)
My geometry:

Boundary condition:
f|_{inlet_{top}} = \rho_1
f|_{inlet_{bottom}} = \rho_2
\frac{df}{dn}|_{wall} = 0
\rho - scalar function, f - scalar function, D - diffusion constant, \vec{v} - velocity vector

F = FunctionSpace(mesh, 'P', 1)
Rho = FunctionSpace(mesh, 'P', 1)
Mu = FunctionSpace(mesh, 'P', 1)
...
f = TrialFunction(F)
fi = TestFunction(F)

rho = Function(Rho)
mu = Function(Mu)
f_n = Function(F)
f_ = Function(F)

Variational formulation:
\int \rho \frac{f -f_n}{dt} fi \; dx + \int \nabla \cdot \, (\rho \, f \, \vec{v}) \, fi \, dx \;= \; \int (D \, \rho \nabla f) \cdot \nabla fi \, dx

Is there a mistake here?
Code:

D = Constant(0.05)
F4 = rho*dot((f - f_n) / k, fi)*dx \
    + div(rho * f * u_)*fi*dx \
    - dot(D*rho * grad(f), grad(fi))*dx

a4 = lhs(F4)
L4 = rhs(F4)

Rho and mu - constant functions (=1)

Tthe calculation is proceeding normally,
but when i try to output field f_n an error occurs:

plot(f_n)

Error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-120-5564c15e22f8> in <module>
----> 1 plot(f_n)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py in plot(object, *args, **kwargs)
    436     # Plot
    437     if backend == "matplotlib":
--> 438         return _plot_matplotlib(object, mesh, kwargs)
    439     elif backend == "x3dom":
    440         return _plot_x3dom(object, kwargs)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py in _plot_matplotlib(obj, mesh, kwargs)
    302 
    303     if isinstance(obj, cpp.function.Function):
--> 304         return mplot_function(ax, obj, **kwargs)
    305     elif isinstance(obj, cpp.function.Expression):
    306         return mplot_expression(ax, obj, mesh, **kwargs)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/common/plotting.py in mplot_function(ax, f, **kwargs)
    150             if mode == "contourf":
    151                 levels = kwargs.pop("levels", 40)
--> 152                 return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)
    153             elif mode == "color":
    154                 shading = kwargs.pop("shading", "gouraud")

/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py in tricontourf(ax, *args, **kwargs)
    273 def tricontourf(ax, *args, **kwargs):
    274     kwargs['filled'] = True
--> 275     return TriContourSet(ax, *args, **kwargs)
    276 
    277 

/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py in __init__(self, ax, *args, **kwargs)
     33         are described in the docstring of `tricontour`.
     34         """
---> 35         ContourSet.__init__(self, ax, *args, **kwargs)
     36 
     37     def _process_args(self, *args, **kwargs):

/usr/lib/python3/dist-packages/matplotlib/contour.py in __init__(self, ax, levels, filled, linewidths, linestyles, alpha, origin, extent, cmap, colors, norm, vmin, vmax, extend, antialiased, *args, **kwargs)
    853         self._transform = kwargs.pop('transform', None)
    854 
--> 855         kwargs = self._process_args(*args, **kwargs)
    856         self._process_levels()
    857 

/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py in _process_args(self, *args, **kwargs)
     45         else:
     46             from matplotlib import _tri
---> 47             tri, z = self._contour_args(args, kwargs)
     48             C = _tri.TriContourGenerator(tri.get_cpp_triangulation(), z)
     49             self._mins = [tri.x.min(), tri.y.min()]

/usr/lib/python3/dist-packages/matplotlib/tri/tricontour.py in _contour_args(self, args, kwargs)
     92                              'triangulation')
     93         if not np.isfinite(z_check).all():
---> 94             raise ValueError('z array must not contain non-finite values '
     95                              'within the triangulation')
     96 

ValueError: z array must not contain non-finite values within the triangulation

I don’t understand what(

If i ask for a task.
Rho and mu:
\mu = \mu_1 \, f + (1-f)\, \mu_2
\rho = \frac{\rho_1 \rho_2}{f\rho_2 + (1-f)\rho_1}
Code:

rho_1 = Constant(1)
rho_2 = Constant(2)
mu_1 = Constant(1)
mu_2 = Constant(2)

mu = project(mu_1 * f_n + (1 - f_n)*mu_2, Mu)
rho = project(rho_1*rho_2/(f_n*rho_2 + (1 - f_n)*rho_1), Rho)

My solution is nan from 1 iteration.

I can post the solution full.

I found error in my program, my solution is nan.