# Drift-Diffusion: Neumann BC with du/dx instead of du/dn

**URL:** <https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478>\
**Category:** variational formulation\
**Created:** [October 17, 2023, 2:20pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478 "2023-10-17T14:20:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![ausler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ausler/32/8877_2.png) [@ausler](https://fenicsproject.discourse.group/u/ausler)\
**Post date:** [October 17, 2023, 2:20pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/1 "2023-10-17T14:20:16Z")

</div>

Hello all,

i want to impose a Neumann boundary condition to a drift–diffusion-problem that simply lets a flux pass through the boundaries. I do this by adding a boundary integral (ds) over the flux term to the variational form. To test this approach, I have implemented a simple example where a concentration peak moves through the cell with a constant drift velocity \mathbf{b} and, eventually, leaves the cell.

So the drift–diffusion equation looks like this:

\frac{\partial u}{\partial t} = -\nabla\mathbf{j}\\ \mathbf{j} = -D\left(\nabla u - \mathbf{b}\cdot u\right)

And, initially, I thought that the Neumann boundary condition can be set in a simple way, like this:  
\int\limits\_\Omega\left(u(t\_{n+1})-u(t\_n)\right)\cdot v\,\mathrm{d}x -\Delta t\cdot\int\limits\_\Omega\mathbf{j}\cdot\nabla v\,\mathrm{d}x+\Delta t\cdot\int\limits\_{\partial\Omega}\mathbf{j}\,v\,\mathrm{d}\mathbf{s} = 0

Minimal code example:

```auto
import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

mesh = IntervalMesh(1000, -5, 5)
V = FunctionSpace(mesh, "CG", 1)
u, u_copy = Function(V), Function(V)
v = TestFunction(V)

# Drift-Diffusion equation
D, dt, drift = Constant(1.), Constant(.01), Constant(+10.)
j = -D*(u.dx(0)-u*drift)
F = (u-u_copy)*v*dx - dt*j*v.dx(0)*dx + dt*j*v*ds

# initial concentration profile
u0 = Expression('exp(-(x[0]-m)*(x[0]-m)/(2*l*l))', m=0., l=.1, degree=1)
u.interpolate(u0)

# solve time propagation
for i in range(100):
    u_copy.vector()[:] = u.vector()[:]
    solve(F==Constant(0.), u, [])
    plot(u)
plt.show()

```

For b=+10, this appears to work, but, if I switch the sign to b=-10, the outcome of the simulation is completely different:

 ![driftdiff_bm10_bp10](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d4c3695add91d1e126cca45658dfe57ae29a0fb7.png)

Basically, the Neumann boundary condition does what I want it to do on the right boundary, but not on the left one. I think this is because I actually want to prescribe gradients of the form \mathrm{d}/\mathrm{d}x, rather than \mathrm{d}/\mathrm{d}\mathbf{n} (with \mathbf{n} being the normal vector of the outer boundary).

For the 1-dimensional case, I could fix this by having separate terms (`-dt*j*ds(1)` for the right boundary, `+dt*j*ds(2)` for the left boundary) for the two boundaries. I am wondering, though, if there is a simpler and more general solution to this that would also work for more complex geometries in higher dimensions.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [October 17, 2023, 2:35pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/2 "2023-10-17T14:35:16Z")

</div>

Why doesn’t `n = ufl.FacetNormal(mesh)` accommodate what you need? In 1D it will be n = -1 on the left and n = +1 on the right as you require.

---

<div class="post-metadata">

**Author:** ![ausler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ausler/32/8877_2.png) [@ausler](https://fenicsproject.discourse.group/u/ausler)\
**Post date:** [October 17, 2023, 2:45pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/3 "2023-10-17T14:45:51Z")

</div>

There seems to be a problem with the dimensions in the variational expression then. I tried formulating the boundary term as `+ dt*n*j*v*ds`, with `n = FacetNormal(mesh)`

But then, I get this error:  
`ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (1,) and free indices with labels ().`

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [October 17, 2023, 2:47pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/4 "2023-10-17T14:47:28Z")

</div>

Looks like a generalisation of the tensor formulated nature of the facet normal. You could just try `n[0]`. In 1D it will always have shape `(1,)`.

---

<div class="post-metadata">

**Author:** ![ausler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ausler/32/8877_2.png) [@ausler](https://fenicsproject.discourse.group/u/ausler)\
**Post date:** [October 17, 2023, 2:48pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/5 "2023-10-17T14:48:22Z")

</div>

Nice, thanks a lot! In case of higher dimensions, would I have to do the same?

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [October 17, 2023, 2:50pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/6 "2023-10-17T14:50:40Z")

</div>

I’ve only glanced at your formulation. The flux is typically written:

\frac{\partial u}{\partial \vec{n}} = \nabla u \cdot \vec{n}.

You can formulate these terms as necessary in UFL, it all depends on the weak formulation and boundary conditions you’re attempting to discretise.

---

<div class="post-metadata">

**Author:** ![ausler](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ausler/32/8877_2.png) [@ausler](https://fenicsproject.discourse.group/u/ausler)\
**Post date:** [October 17, 2023, 2:54pm UTC](https://fenicsproject.discourse.group/t/drift-diffusion-neumann-bc-with-du-dx-instead-of-du-dn/12478/7 "2023-10-17T14:54:56Z")

</div>

Alright, thanks. Guess I’ll try it out in 2D to see if something like `inner(j, n)` works.