Need help solving Convective and Conductive heat equation.

Dear All.
I am trying to solve a heat transfer equation with a conductive and convective term with the form of

k*grad^2*T-Cp*u*grad(T)=0

over a cube space with 2 boundary conditions, a Dirichlet boundary on one surface and a numan boundary conditions on an other where -dT/dx=Q(T-Text)
Here k is conductive heat transfer coefficient, Cp is heat capacity and u is velocity profile defined by an function for laminar flow so I do not solve Navia stokes here.

The problem above is linear, but my final goal is to solve a non-linear problem (as Cp and k are functions of T) so the example code will use a nonlinera solver in it, and the function is written as such as well.
I did integration by parts on the first term of the function above and got the following variational problem as written in Fenics UFL format
Here T is the trial function, and v is the test function

F=inner(k*grad(T),grad(v))*dx-inner(Cp*dot(grad(T),u),v)*dx+Q*T*v*ds(1)-Q*T_ext*v*ds(1)

Here T_ext, Cp, and k are constants, and the dot product is needed as velocity function is a vector as well as grad(T) should be a vector, and needs to be multiplied out, and ds(1) is the surface perpendicular to the flow direction.

I am able to solve this now with the code bellow, and would appreciate if some one could double check my implementation for any errors. Thank you.

from dolfin import *
# Create mesh and define function space
class lf_function(Expression):
    def setup(self,FH=0.25,FL=0.25,Fxoff=0,Fyoff=0):
        self.l=FH #y length
        self.L=FL #x length
        self.A=self.l/self.L
        self.n=100
        self.xOffset=Fxoff
        self.yOffset=Fyoff

    def m_0(self,n):

        a=lambda n:(128.0*(-1.0)**n)/(self.A*pi**5.0*(2.0*n+1.0)**5.0)
        b=lambda n:np.tanh((2.0*n+1.0)*pi*self.A/2.0)
        c=lambda n:np.sin((2.0*n+1.0)*pi/2.0)
        sumabc=lambda n:a(n)*b(n)*c(n)

        return -2.0/3.0+np.sum(sumabc(n))
    def U(self,x):
        n=np.array(range(self.n))
        X=(x[0]-self.xOffset)/self.L
        Y=(x[1]-self.yOffset)/self.l
        sumTop =lambda n:32.0*((-1.0)**n)*np.cosh((2.0*n+1.0)*pi*self.A*Y/2.0)
        sumBottom=lambda n:((2.0*n+1.0)**3.0)*(pi**3.0)*np.cosh((2.0*n+1)*pi*self.A/2.0)
        sumMultiplier=lambda n:np.cos((2.0*n+1)*pi*X/2.0)
        sumTotal=lambda n:(sumTop(n)/sumBottom(n))*sumMultiplier(n)
        m_0=self.m_0(n)
        return (1.0/m_0)*((X**2.0-1.0)+np.sum(sumTotal(n)))

    def eval(self, value, x):
        if x[0]<=self.L+self.xOffset and x[0]>=self.xOffset-self.L \
                and x[1]<=self.l+self.yOffset and \
                x[1]>=self.yOffset-self.l:
            value[2]=self.U(x)
            value[0]=0
            value[1]=0
        else:
            value[2]=0.0
            value[1]=0.0
            value[0]=0.0
    def value_shape(self):
        return (3,)
cellMesh = UnitCubeMesh(5, 5,5)
V = FunctionSpace(cellMesh, "Lagrange", 2)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[2] < DOLFIN_EPS
class LowerRobinBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14   # tolerance for coordinate comparisons
        #if abs(x[1]) <= tol:
            #print x[0],x[1],x[2]
        return on_boundary and abs(x[1]) <= tol

# Define DiricletBC condition
u0 = Constant(300)
bc = DirichletBC(FunctionSpace(cellMesh, "Lagrange", 2), u0, boundary)
bcs=[bc]
# Defining Numan Boundary

boundary_parts = \
MeshFunction("size_t", cellMesh, cellMesh.topology().dim()-1)
boundary_parts.set_all(0)
Gamma_R = LowerRobinBoundary()
Gamma_R.mark(boundary_parts, 1)

uflow=lf_function(degree=2)
uflow.setup(FH=0.5,FL=0.5,Fyoff=0.5,Fxoff=0.5)
# Define variational problem
ds=Measure('ds')[boundary_parts]
T = TrialFunction(V)
v = TestFunction(V)
dT = TrialFunction(V)
k=Constant(0.6)
Cp=Constant(200)
u=uflow
Q=Constant(100)
T_ext=Constant(270)
#uvec=interpolate(u,VectorFunctionSpace(cellMesh, "Lagrange", 2))
F=inner(k*nabla_grad(T),nabla_grad(v))*dx\
    -inner(Cp*dot(grad(T),u),v)*dx(2)\
    +Q*T*v*ds(1,subdomain_data = boundary_parts)\
    -Q*T_ext*v*ds(1,subdomain_data = boundary_parts)

T_=Function(V)

F=action(F,T_)

J=derivative(F,T_,dT)
problem = NonlinearVariationalProblem(F,T_,J=J,bcs=bcs) #,bcs=BC3d
solver = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-8
prm['newton_solver']['maximum_iterations'] = 10000
prm['newton_solver']['relaxation_parameter'] = 1.0

set_log_level(PROGRESS)
solver.solve()

# Save solution in VTK format
file = File("convectionTest.pvd")
file << T_