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:
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:
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.