Error in determining a linear variational problem

variational formulation
1

Hi, everyone.

I am trying to solve the problem of two-phase filtering, but when I solve it, I get an error:

*** 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.

Can you help me solve this problem?
from dolfin import *
mesh = UnitSquareMesh(50, 50, “crossed”)

order = 1
BDM = FiniteElement("BDM", triangle, order)
DG = FiniteElement("DG", triangle, order - 1)
mixed = MixedElement(BDM, DG, DG)

W = FunctionSpace(mesh, mixed)

dw = TrialFunction(W)

wt = TestFunction(W)
v, q, r = split(wt)

w = Function(W)
w0 = Function(W)

u, p, s = split(w)
u0, p0, s0 = split(w0)

sm = 0.5*(s0+s)

class PressureBC(UserExpression):
    def eval_cell(self, values, x, cell):
        values[0] = 1.0-x[0]

class SaturationBC(UserExpression):
    def eval(self, values, x):
        if x[0] < DOLFIN_EPS:
            values[0] = 1.0

pbar = PressureBC()
sbar = SaturationBC()

dt = Constant(0.01)

kinv = Expression("1 / max (exp(-pow(10 * (x[1] - 0.5 - 0.1 * sin(10 * x[0])), 2)), 0.01)",     degree=2)
zero = Constant(0.0)
Kinv = as_matrix(((kinv, zero), (zero, kinv)))

mu_rel = 0.2

def linv(s):
    return 1.0 / ((1.0 / mu_rel)*s**2 + (1.0 - s)**2)

def F(s):
    return s**2 / (s**2 + mu_rel * (1.0 - s)**2)

n = FacetNormal(mesh)

un   = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))
stab = inner(jump(r), un('+') * F(s0)('+') - un('-') * F(s0)('-')) * dS + r * un_h * sbar * ds

F1 = inner(v, linv(s0) * Kinv * u) * dx - div(v) * p * dx + inner(v, pbar * n) * ds
F2 = q * div(u0) * dx
F3 = r * (s - s0) * dx - dt * inner(grad(r), F(s0) * u) * dx + dt * r * F(s0) * un * ds + dt * stab

F = F1 + F2 + F3

a = lhs(F)
L = rhs(F)

w = Function(W)

folder = "res/examp-liner/"

u_file = File(folder + "velocity.pvd")
p_file = File(folder + "pressure.pvd")
s_file = File(folder + "saturation.pvd")

max_time_iters = 100
for iter in range(max_time_iters + 1):
    print("Time iteration = %d" % iter) 
    solve(a == L, w)
    u, p, s = w.split() 

    u_file << u
    p_file << p
    s_file << s

    w0.assign(w)
2

You are not using your trial functions in the variational form, thus it can not be solved as a==L.
You need to replace