Define time dependent vector valued function used in bilinear form

General
1

I want to define a vector valued function which will be used in bilinear form.
The problem I consider is u_{t}+dot(b,grad(u))+epsilon*div(grad(u))+cu=f
where b=(b_1(x,t),b_2(x,t)) in 2D and c=c(x,t).
I try to define them as

def u_exact(mod, t):
            return lambda x: mod.exp(mod.sin(2*mod.pi*t))*2*mod.sin(x[0])*(1-mod.exp(-2*(1-x[0])/epsilon.value))*(1-mod.exp(-(1-x[1])/epsilon.value))*x[1]**2

 #define b
 def bx(t):
      return lambda x:  2+(2*x[0]-x[1])*ufl.cos(2*ufl.pi*t)
def by(t):
      return  lambda x: 3+(-x[0]+2*x[1])*ufl.cos(2*ufl.pi*t)
def c(t):
       return lambda x: 3+t**2

 x=ufl.SpatialCoordinate(domain)
 t_var=ufl.variable(t)

exact_soln_ufl=u_exact(ufl,t_var)  #used for linear form
exact_soln_np=u_exact(np,t.value)    #used for interpolation

 bx=bx(t_var)
 by=by(t_var)
 c=c(t_var)

and I try to use them in

a=fem.form(u*v*dx+dt*epsilon*inner(grad(u),grad(v))*dx+dt*(bx*u.dx(0)*v*dx+by*u.dx(1))*v*dx+ dt*c*u*v*dx +
           dt*delta*(-epsilon*div(grad(u))+(bx*u.dx(0)+by*u.dx(1))+c*u)*(bx*v.dx(0)+by*v.dx(1))*dx+ 
           (u*delta*(bx*v.dx(0)+by*v.dx(1))*dx))

However , I receive error as

File "/root/shared/Fenics Codes/fenics/FSI_Solver/SUPG_TEST.py", line 82, in <module>
    a=fem.form(u*v*dx+dt*epsilon*inner(grad(u),grad(v))*dx+dt*(bx*u.dx(0)*v*dx+by*u.dx(1))*v*dx+ dt*c*u*v*dx +\
TypeError: unsupported operand type(s) for *: 'function' and 'Indexed'

So how to fix it ? And is there any other way to define a vector valued function used in bilinear form ?
Please help! Thanks in advance.

The complete code is

from dolfinx import mesh,fem,plot
import numpy as np
from ufl import div,grad,diff,dot,inner,dx,ds,CellDiameter,TrialFunction,TestFunction,sin,pi,exp,Measure
import ufl
from mpi4py import MPI
from petsc4py import PETSc

T=1
ny=nx
num_steps=nx*nx   
domain=mesh.create_unit_square(MPI.COMM_WORLD,nx,ny,mesh.CellType.triangle)
V=fem.FunctionSpace(domain,("CG",1))
k=T/num_steps
t=fem.Constant(domain,PETSc.ScalarType(0))
dt=fem.Constant(domain,PETSc.ScalarType(k))
epsilon=fem.Constant(domain,PETSc.ScalarType(1e-8))

def u_exact(mod, t):
    return lambda x: mod.exp(mod.sin(2*mod.pi*t))*2*mod.sin(x[0])*(1-mod.exp(-2*(1-x[0])/epsilon.value))*(1-mod.exp(-(1-x[1])/epsilon.value))*x[1]**2

#define b
def bx(t):
    return lambda x:  2+(2*x[0]-x[1])*ufl.cos(2*ufl.pi*t)
def by(t):
    return  lambda x: 3+(-x[0]+2*x[1])*ufl.cos(2*ufl.pi*t)
def c(t):
    return lambda x: 3+t**2


x=ufl.SpatialCoordinate(domain)
t_var=ufl.variable(t)

exact_soln_ufl=u_exact(ufl,t_var)  #used for linear form
exact_soln_np=u_exact(np,t.value)    #used for interpolation

bx=bx(t_var)
by=by(t_var)
c=c(t_var)

#Apply homogeneous Dirichlet boundary condition
boundary_fun=fem.Function(V)
boundary_fun.interpolate(exact_soln_np)

tdim=domain.topology.dim
fdim=tdim-1
domain.topology.create_connectivity(fdim,tdim)
boundary_facets=mesh.exterior_facet_indices(domain.topology)
bc=fem.dirichletbc(boundary_fun,fem.locate_dofs_topological(V,fdim,boundary_facets))

#Initial condition
u_old=fem.Function(V)
u_old.interpolate(exact_soln_np)

#Define right hand side f
f=diff(exact_soln_ufl(x),t_var)-epsilon*div(grad(exact_soln_ufl(x)))+bx(x)*exact_soln_ufl(x).dx(0)+by(x)*exact_soln_ufl(x).dx(1)+c(x)*exact_soln_ufl(x)

u,v=TrialFunction(V),TestFunction(V)
h=CellDiameter(domain)
delta=h/40

a=fem.form(u*v*dx+dt*epsilon*inner(grad(u),grad(v))*dx+dt*(bx*u.dx(0)*v*dx+by*u.dx(1))*v*dx+ dt*c*u*v*dx +\
            dt*delta*(-epsilon*div(grad(u))+(bx*u.dx(0)+by*u.dx(1))+c*u)*(bx*v.dx(0)+by*v.dx(1))*dx+ \
            (u*delta*(bx*v.dx(0)+by*v.dx(1))*dx))
L=fem.form(u_old*v*dx+dt*f*(v+delta*(bx*v.dx(0)+by*v.dx(1)))*dx+u_old*delta*(bx*v.dx(0)+by*v.dx(1))*dx)
2

You really need to put some work into simplifying your code to pin down errors.
For instance, you should figure out what term of your form that causes the issue (as your form is rather long). By reducing it term by term, you will realize that the issue is with the definition of bx and by.

You can for instance use the following definition:

from dolfinx import mesh,fem,plot
import numpy as np
from ufl import div,grad,diff,dot,inner,dx,ds,CellDiameter,TrialFunction,TestFunction,sin,pi,exp,Measure
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from functools import partial
nx = 1
T=1
ny=nx
num_steps=nx*nx   
domain=mesh.create_unit_square(MPI.COMM_WORLD,nx,ny,mesh.CellType.triangle)
V=fem.FunctionSpace(domain,("CG",1))
k=T/num_steps
t=fem.Constant(domain,PETSc.ScalarType(0))
dt=fem.Constant(domain,PETSc.ScalarType(k))

def u_exact(mod, t):
    return lambda x: mod.exp(mod.sin(2*mod.pi*t))*2*mod.sin(x[0])*(1-mod.exp(-2*(1-x[0])/epsilon.value))*(1-mod.exp(-(1-x[1])/epsilon.value))*x[1]**2

#define b
def bx(x, t):
    return 2+(2*x[0]-x[1])*ufl.cos(2*ufl.pi*t)
def by(x, t):
    return 3+(-x[0]+2*x[1])*ufl.cos(2*ufl.pi*t)
def c(t):
    return 3+t**2


x=ufl.SpatialCoordinate(domain)
t_var=ufl.variable(t)

bx=partial(bx, t=t_var)
by=partial(by, t=t_var)
c=c(t_var)


u,v=TrialFunction(V),TestFunction(V)

a = dt*(bx(x)*u.dx(0)+by(x)*u.dx(1))*v*dx

which makes the first part of the code run.
You can then extend it term by term.