### Integrate derivative of known function

To calculate the velocity in an advection problem I first calculate pressure and then

`vel = (-w**2/(12*mu))*grad(p)`

where w and mu are constants, If I solve the pure advection equation in the case of a constant velocity vector

$\frac{\partial c}{\partial t}+u_i\frac{\partial c}{\partial x_i}=0$∂

`c`∂

`t`+

`u`

_{i}∂

`c`∂

`x`

_{i}=0

where c is the concentration to solve and u is the velocity vector. After integration by parts to implement Neumann BC, the following code works

```
q = TestFunction(V)
c = TrialFunction(V)
F = (1/dt)*w*c*q*dx - (1/dt)*w*c_1*q*dx \
-inner(c,inner(vel,grad(q)))*dx \
+inner(c,dot(vel,n))*q*ds(1)
a, L = lhs(F),rhs(F)
c = Function(V)
t = 0.0
for n in range(num_steps):
solve(a==L,c,bc)
c_1.assign(c)
print('t = %.4f' % (t))
t += dt
```

But when I try to solve the case where the velocity changes with x

$\frac{\partial c}{\partial t}+u_i\frac{\partial c}{\partial x_i}+c\frac{\partial u_i}{\partial x_i}=0$∂`c`∂`t` +`u`_{i}∂`c`∂`x`_{i} +`c`∂`u`_{i}∂`x`_{i} =0

the code below pops up the error "Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ()." How could I improve the code? thank you very much for your comments

```
F = (1/dt)*w*c*q*dx - (1/dt)*w*c_1*q*dx \
+c*Dx(vel,0)*q*dx+c*Dx(vel,1)*q*dx\
-inner(c,inner(vel,grad(q)))*dx \
+inner(c,dot(vel,n))*q*ds(1)
#I also try
F = (1/dt)*w*c*q*dx - (1/dt)*w*c_1*q*dx \
+inner(c,grad(vel))*q*dx\
-inner(c,inner(vel,grad(q)))*dx \
+inner(c,dot(vel,n))*q*ds(1)
```

### 1 Answer

`+c*Dx(vel,0)*q*dx+c*Dx(vel,1)*q*dx\`

translates to a vectorial expression. Both terms contain derivatives of the whole vector w.r.t. one spatial variable yielding a vectorial expression, which does not represent the divergence you gave in your equation.

Similarly the second line of your second version contains the term `grad(vel)`

which gives a second order tensor.

The divergence $\frac{\partial u_i}{\partial x_i}$∂`u`_{i}∂`x`_{i} can be written as `div(u)`

or as `u[i].dx(i)`

.