Diffusion along vector field (Nernst eq.)

You need to provide the definition of all your functions, or we will have
to guess what you are doing.

Martin
Den 11. apr. 2012 19:15 skrev "Charl" <email address hidden>
følgende:

> New question #193347 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/193347
>
> Hello all,
>
> I am attempting to solve the Nernst-Planck equation, i.e. diffusion of
> charged particles along an electrochemical gradient. Purely chemical
> diffusion works, but incorporating the electric field into particle motion
> is causing trouble.
>
> The basic Nernst-Planck equation is: dc/dt = div(a*E + b*grad(c)), where c
> is the (charged) particle concentration, E is the electric field, and a and
> b are scalar constants.
>
> Assembly of the problem is coded as follows:
>
> --------
> trialfunc = TrialFunction(mesh_funcspace_scalar)
> testfunc = TestFunction(mesh_funcspace_scalar)
> LHS = trialfunc*testfunc*dx + \
> dt*diffusivity*inner(nabla_grad(trialfunc),
> nabla_grad(testfunc))*dx + \
> dt*nernst_constant*nabla_div(efield)*testfunc*dx
> RHS = (concentration_prev + dt*concentration_injection)*testfunc*dx
> solve(LHS == RHS, concentration, bc_diffusion)
> --------
>
> The second term of the LHS is chemical diffusion (working), which is of
> rank 2.
> The third term of the LHS is electrical diffusion, which is of rank 1.
>
> Obviously, the Nernst term should be of the same rank. However, I cannot
> seem to figure out where it is going wrong.
>
> Much obliged,
> Charl
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

It should be:

   dc/dt = div(c*a*E + b*grad(c))

However, discretizing this using standard FEM is infamous for generating
unstable solutions. Upwind stabilization can help, but if E comes from a
static potential you can discretize and solve another problem instead.
Have a look at:

   http://dx.doi.org/10.1529/biophysj.106.102533

For a thorough description of how to set up the equation and solve it
using FEM.

Johan

On 04/11/2012 07:15 PM, Charl wrote:
> New question #193347 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/193347
>
> Hello all,
>
> I am attempting to solve the Nernst-Planck equation, i.e. diffusion of charged particles along an electrochemical gradient. Purely chemical diffusion works, but incorporating the electric field into particle motion is causing trouble.
>
> The basic Nernst-Planck equation is: dc/dt = div(a*E + b*grad(c)), where c is the (charged) particle concentration, E is the electric field, and a and b are scalar constants.
>
> Assembly of the problem is coded as follows:
>
> --------
> trialfunc = TrialFunction(mesh_funcspace_scalar)
> testfunc = TestFunction(mesh_funcspace_scalar)
> LHS = trialfunc*testfunc*dx + \
> dt*diffusivity*inner(nabla_grad(trialfunc), nabla_grad(testfunc))*dx + \
> dt*nernst_constant*nabla_div(efield)*testfunc*dx
> RHS = (concentration_prev + dt*concentration_injection)*testfunc*dx
> solve(LHS == RHS, concentration, bc_diffusion)
> --------
>
> The second term of the LHS is chemical diffusion (working), which is of rank 2.
> The third term of the LHS is electrical diffusion, which is of rank 1.
>
> Obviously, the Nernst term should be of the same rank. However, I cannot seem to figure out where it is going wrong.
>
> Much obliged,
> Charl
>

Gentlemen,

A belated thank you for your support. Johan, I have gone over the series of papers by Cheng, Song, Holst, McCammon et al., but the only reference I can find to numerical instabilities is their use of either backward Euler or Crank-Nicholson for time stepping.

Based on the tutorial at http://fenicsproject.org/documentation/tutorial/time-dependent.html I believe I now have the equation correct (using backward Euler):

-------------
LHS = trialfunc_diffusion*testfunc_diffusion*dx \
                 + dt*diffusivity*nernst_constant*trialfunc_diffusion*inner(efield, nabla_grad(testfunc_diffusion))*dx
RHS = concentration_prev*testfunc_diffusion*dx
solve(LHS == RHS, concentration, bc_diffusion)
-------------

diffusivity is a Function(), nernst_constant = Expression('80') and dt = 0.1. For the full script and a demo mesh, see http://www.turingbirds.com/upload/files/demo_nernst.tar

The LHS contains only the Nernst term; I left out the diffusion term to facilitate troubleshooting. A further simplification is a constant E-field rather than feedback via charged ion species. Starting off the simulation with a delta peak concentration, the result quickly becomes 'ragged', and the concentration begins to take on negative values, as illustrated:

http://turingbirds.com/upload/files/nernst_planck.png

Could this be due to a poor or lacking interpolation somewhere? Or is it a result of numerical errors, which might be resolved by using a considerably lower dt?

Thanks in advance,
Charl

On 04/26/2012 05:15 PM, Charl wrote:
> Question #193347 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/193347
>
> Status: Answered => Open
>
> Charl is still having a problem:
> Gentlemen,
>
> A belated thank you for your support. Johan, I have gone over the series
> of papers by Cheng, Song, Holst, McCammon et al., but the only reference
> I can find to numerical instabilities is their use of either backward
> Euler or Crank-Nicholson for time stepping.

The point is that they avoid the numerical instabilities by solving a
slightly different PDE, but see below.

> Based on the tutorial at http://fenicsproject.org/documentation/tutorial
> /time-dependent.html I believe I now have the equation correct (using
> backward Euler):
>
> -------------
> LHS = trialfunc_diffusion*testfunc_diffusion*dx \
> + dt*diffusivity*nernst_constant*trialfunc_diffusion*inner(efield, nabla_grad(testfunc_diffusion))*dx
> RHS = concentration_prev*testfunc_diffusion*dx
> solve(LHS == RHS, concentration, bc_diffusion)
> -------------
>
> diffusivity is a Function(), nernst_constant = Expression('80') and dt =
> 0.1. For the full script and a demo mesh, see
> http://www.turingbirds.com/upload/files/demo_nernst.tar
>
> The LHS contains only the Nernst term; I left out the diffusion term to
> facilitate troubleshooting. A further simplification is a constant
> E-field rather than feedback via charged ion species.

Having a non-constant field is making the problem even more troublesome
problem. Then you cannot use the method as suggested by Cheng et al.
Then you need to solve a poisson bolzmann equation coupled with the
electro diffusion problem.

> Starting off the
> simulation with a delta peak concentration, the result quickly becomes
> 'ragged', and the concentration begins to take on negative values, as
> illustrated:
>
> http://turingbirds.com/upload/files/nernst_planck.png
>
> Could this be due to a poor or lacking interpolation somewhere? Or is it
> a result of numerical errors, which might be resolved by using a
> considerably lower dt?

That or using a more refined mesh.

Johan

> Thanks in advance,
> Charl
>