Unable to define linear variational problem a(u, v) == L(v) for all v

Hey guys,
below you find my attempt to implement the navier stokes eq. step by step. now i am just trying to solve two of t hose and work my way up:

from __future__ import print_function
import fenics
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mshr import *
from fenics import XDMFFile
from fenics import derivative


# Mesh und Umgebung erstellen
l1 = 3
l2 = 1
domain = Rectangle(Point(0, 0), Point(l1, l2))
mesh = generate_mesh(domain, 128)


# Elemente, Funktionen und Funktionsraum erstellen
P2 = VectorElement("P", triangle, 2)
P1 = FiniteElement("P", triangle, 1)
element = MixedElement([P2, P1])


V = FunctionSpace(mesh, element)

solution = Function(V)
u, rho = solution.split()
solution_n = Function(V)
u_n, rho_n = solution_n.split()
w_1, w_2 = TestFunctions(V)



# Zeitdiskretisierung
T = 5.0                        # final time
f = 1/T
num_steps = 100                # number of time steps
dt = T / num_steps             # time step size


# erste vereinfachte Gleichung
F1 = inner((u-u_n), w_1) * dx \
    + dt * inner(nabla_grad(u), nabla_grad(w_1)) * dx

F2 = inner((rho-rho_n), w_2) * dx \
    + dt * inner(div(u), w_2) * dx

F = F1 + F2
a = lhs(F)
L = rhs(F)


# Anfangsbedingungen
u_init = Constant((0.0, 0.0))
u_1_init = Constant(0.0)
u_2_init = Constant(0.0)
rho_init = Constant(998.0)
solution_init = Expression(('u_1_init', 'u_2_init', 'rho_init'), u_1_init = u_1_init, u_2_init = u_2_init, rho_init = rho_init, degree=2)


# Randbedingungen + Schallwelle
t = 0
amplitudeu = Constant((1.5, 0.0))
Schallwelle = Expression(('u_init[0] + amplitudeu[0] * sin(2 * pi * f * t)', 'u_init[1] + amplitudeu[1] * sin(2 * pi * f * t)'), u_init=u_init, amplitudeu=amplitudeu,  f=f, t=t, degree=2, domain=mesh)
Wand = Constant((0.0, 0.0))


tol = 1E-14

def boundary_Dl(x, on_boundary):
    if on_boundary:
        if near(x[0], 0, tol):
            return True
        else:
            return False
    else:
        return False

def boundary_Dr(x, on_boundary):
    if on_boundary:
        if near(x[0], l1, tol):
            return True
        else:
            return False
    else:
        return False


bcl = DirichletBC(V.sub(0), Schallwelle, boundary_Dl)
bcr = DirichletBC(V.sub(0), Wand, boundary_Dr)
bc = [bcl, bcr]


# pvd files erstellen´
xdmf_file_u = XDMFFile('RAWHEATNAVIERtest/u_1.xdmf')
xdmf_file_rho = XDMFFile('RAWHEATNAVIERtest/rho_1.xdmf')


solutionX = Function(V)
solution = Function(V)                         **LINE 100**
t = 0


for n in range(num_steps):


    # Solve for u

    solve(a == L, solution, bc)


    # Assign values to the storage function

    solutionX.assign(solution)


    # Write to XDMF file

    xdmf_file_u.write(u_solution, t)


    # Update for the next time step


    t += float(dt)
    Schallwelle.t = t
    bc = [bcl, bcr]
    u_n.assign(u)

This leads to the following error:

Form has no parts with arity 2.
Traceback (most recent call last):
File “RAWHEATNAVIER.py”, line 115, in
solve(a == L, solution, bc)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 242, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/problem.py”, line 59, in init
cpp.fem.LinearVariationalProblem.init(self, a, L, u._cpp_object, bcs)
RuntimeError:
*** Error: Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason: Expecting the left-hand side to be a bilinear form (not rank 0).
*** Where: This error was encountered inside LinearVariationalProblem.cpp.
*** Process: 0

I read that the problem might be the mixed element space but I dont really know what to do differently. Thank you very much for reading and helping!

Since you are solving a linear problem, this should be u,rho =TrialFunctions(V)

Thanks man!
That helped