pnp.h

Go to the documentation of this file.

/** @file pnp.h
 */
/**
# Ohmic conduction flux of charged species
 
This function computes the fluxes due to ohmic conduction appearing in
the [Nernst--Planck
equation](http://en.wikipedia.org/wiki/Nernst%E2%80%93Planck_equation). The
species charge concentrations are then updated using the explicit
scheme
\f[
c^{n+1}_i = c^n_i +\Delta t \, \nabla \cdot( K_i c^n_i \nabla \phi^n)
\f] 
where \f$c_i\f$ is the volume density of the \f$i\f$-specie, \f$K_i\f$ its volume
electric conductivity and \f$\phi\f$ the electric potential. */
 
extern scalar phi;
 
void ohmic_flux (scalar * c,         // A list of the species concentration...
     int * z,            // ... and their corresponding valences
     double dt,
     vector * K = NULL)  // electric mobility (default the valence)
{
  /**
  If the volume conductivity is not provided it is set to the value of
  the valence. */
  
  if (!K) { // fixme: this does not work yet
    int i = 0;
    for (int _s = 0; _s < 1; _s++) /* scalar in c */ {
      const vector kc[] = {z[i], z[i]}; i++;
      K = vectors_append (K, kc); // fixme: K should be freed eventually
    }
  }
 
  scalar s;
  const vector k;
  for (s, k in c, K) {
 
    /**
    The fluxes of each specie through each face due to ohmic transport
    are */
 
    vector f[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      f.x[] = k.x[]*(s[] + s[-1])*(phi[] - phi[-1])/(2.*Delta);
 
    /**
    The specie concentration is updated using the net amount of that
    specie leaving/entering each cell through the face in the interval
    \f$dt\f$ */
 
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++)
        s[] += dt*(f.x[1] - f.x[])/Delta;
  }
}