Ohmic conduction flux of charged species

    This function computes the fluxes due to ohmic conduction appearing in the Nernst–Planck equation. The species charge concentrations are then updated using the explicit scheme \displaystyle c^{n+1}_i = c^n_i +\Delta t \, \nabla \cdot( K_i c^n_i \nabla \phi^n) where c_i is the volume density of the i-specie, K_i its volume electric conductivity and \phi 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 (scalar s in c) {
          const face vector kc[] = {z[i], z[i]}; i++;
          K = vectors_append (K, kc); // fixme: K should be freed eventually
      scalar s;
      (const) face vector k;
      for (s, k in c, K) {

    The fluxes of each specie through each face due to ohmic transport are

        face vector f[];
          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 dt

            s[] += dt*(f.x[1] - f.x[])/Delta;