Implementing a very simple 1D advection-diffusion demo

Question

I am just starting to learn FEM and FEniCS by constructing some simple problems (although in the past I have written my own finite volume code in python so I'm not a total beginner to this area). I have taken the Poisson example code and successfully modified it to solved a steady-state diffusion problem. I am now trying to extend this by adding an advection term. I'd like to solve.
$$
0 = \nabla \cdot \left(d\nabla u\right) - \nabla \cdot\left(\mathbf{v} u\right)
$$
subject to the inflow/outflow boundary conditions,
$$
u(0) = 1 \quad u(1) = 0
$$

So in weak form,
$$
0 = \int_{\Omega} d\nabla u \nabla\nu dx + \int_{\Omega}\nabla\cdot\left( \mathbf{v} u \right) \nu dx
$$

I think you can write that as the following in FEniCS,

d = Constant(1.0)    # diffusion coefficient
vel = Constant(0.1)  # advection velocity
a = inner(d*grad(u), grad(v))*dx + div(vel*grad(u))*v*dx

In particular I am stuck on how to implement the advection term in FEniCS. How can an specify a fixed vector field $\mathbf{v}$ for the advection velocity?

The script I have written is obviously wrong (a python dolfin script is on gist as 1d_advection_diffusion.py if anybody is interested in taking a look) because the solution doesn't seem to be dependent on the coefficients. This problem has the analytical solution,
$$
u(x) = \frac{e^{\lambda} - e^{x\lambda} }{e^{\lambda} - 1}
$$
where $\lambda=a/d$. So the ratio of the coefficient should effect the width of the boundary layer at x=1.

Can anybody guide me with the implementation of this example problem?

Comments

You are missing a $-$ sign in the weak formulation.

---

Corrected above, thanks. I can't run my code here (only installed on work computer). I'm not too sure that the sign error was the problem because nothing seemed to effect the solution in the script that I wrote, but I will check as soon as I can.

Answer 1

The second term in your equation differs from the second term in your program code. Is it supposed to be div(vel * u)*v or div(vel * grad(u))*v?

To make vel a velocity field rather than a constant, you have two options. Either write it as an expression

V = VectorFunctionSpace(your_mesh, 'CG', 1)
# Example for vel(x) = x
vel = Expression(('x[0]', 'x[1]', 'x[2]'), element=V.ufl_element())

or as a function defined by its degrees of freedom:

V = VectorFunctionSpace(your_mesh, 'CG', 1)
vel = Function(V)

# now assign degrees of freedom to vel.vector()

# Example for vel(x) = x:
all_coords = V.dofmap().tabulate_all_coordinates(V.mesh())
(xs, ys, zs) = all_coords.reshape((-1,3)).T
vel_v = vel.vector()

indices = V.sub(0).dofmap().dofs()
vel_v[indices] = xs

indices = V.sub(1).dofmap().dofs()
vel_v[indices] = ys

indices = V.sub(2).dofmap().dofs()
vel_v[indices] = zs