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 ${c}_{i}^{n+1}={c}_{i}^{n}+\Delta t\phantom{\rule{0.167em}{0ex}}\nabla \cdot \left({K}_{i}{c}_{i}^{n}\nabla {\varphi }^{n}\right)$ where ${c}_{i}$ is the volume density of the $i$-specie, ${K}_{i}$ its volume electric conductivity and $\varphi$ the electric potential.

``````extern scalar φ;

struct Species {
scalar * c; // A list of the species concentration and their corresponding
int * z;    // valences
double dt;
// optional
vector * K; // electric mobility (default the valence)
};

void ohmic_flux (struct Species sp)
{``````

If the volume conductivity is not provided it is set to the value of the valence.

``````  if (!sp.K) { // fixme: this does not work yet
int i = 0;
for (scalar s in sp.c) {
const face vector kc[] = {sp.z[i], sp.z[i]}; i++;
sp.K = vectors_append (sp.K, kc);
}
}

scalar c;
(const) face vector K;
for (c, K in sp.c, sp.K) {``````

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

``````    face vector f[];
foreach_face()
f.x[] = K.x[]*(c[] + c[-1,0])*(φ[] - φ[-1,0])/(2.*Δ);
boundary_flux ({f});``````

The specie concentration is updated using the net amount of that specie leaving/entering each cell through the face in the interval $dt$

``````    foreach()
foreach_dimension()
c[] += sp.dt*(f.x[1,0] - f.x[])/Δ;
boundary ({c});
}
}``````