matscipy.electrochemistry.poisson_nernst_planck_solver

Compute ion concentrations with general Poisson-Nernst-Planck (PNP) equations.

Authors:: Johannes Hoermann <johannes.hoermann@imtek-uni-freiburg.de>

Functions

`B`(x)	Bernoulli function.
`jacobian`(f, x0[, dx])	Naive way to construct N x N Jacobin Fij from N-valued function f of N-valued vector x0.

Classes

PoissonNernstPlanckSystem(*args, **kwargs)

Describes and solves a 1D Poisson-Nernst-Planck system

matscipy.electrochemistry.poisson_nernst_planck_solver.B(x): Bernoulli function.

matscipy.electrochemistry.poisson_nernst_planck_solver.jacobian(f, x0, dx=nan)

Naive way to construct N x N Jacobin Fij from N-valued function f of N-valued vector x0.

Parameters:

f (callable) – N-valued function of N-valued variable vector
x0 ((N,) ndarray) – N-valued variable vector
dx (float (default: np.nan)) – Jacobian built with finite difference scheme of spacing dx. If np.nan, then use machine precision.

Returns:

F – NxN-valued 2nd order finite difference scheme approximate of Jacobian convention: F_ij = dfidxj, where i are array rows, j are array columns

Return type:

(N,N) ndarray

class matscipy.electrochemistry.poisson_nernst_planck_solver.PoissonNernstPlanckSystem(*args, **kwargs)

Bases: object

Describes and solves a 1D Poisson-Nernst-Planck system

Attributes:

charge_density
concentration
grid
ionic_strength: Compute the system’s ionic strength from charges and concentrations.
lambda_D: Compute the system’s Debye length.
potential
x1_scaled

Methods

`G`(x)	Non-linear system
`discretize`()	Sets up discretization scheme and initial value
`init`([c, z, L, lambda_S, x0, T, delta_u, ...])	Initializes a 1D Poisson-Nernst-Planck system description.
`initial_values`()	Solves decoupled linear system to get initial potential distribution.
`left_controlled_volume_scheme_flux_bc`(x, k)	Compute left hand side flux boundary condition residual in accord with controlled volume scheme.
`left_dirichlet_bc`(x, k[, x0])	Construct Dirichlet BC at left boundary
`left_finite_difference_scheme_flux_bc`(x, k)	param x: N-valued variable vector
`left_robin_bc`(x, k, lam[, x0])	Compute left hand side Robin (u + lam*dudx = u0 ) BC at in accord with 2nd order finite difference scheme.
`nernst_planck_pde`(x)	Returns Nernst-Planck equation residual by applying controlled volume scheme
`newton`(f, xij, **kwargs)	Newton solver expects system f and initial value xij
`number_conservation_constraint`(x, k, N0)	N0: total amount of species, k: ion species
`poisson_pde`(x)	Returns Poisson equation residual by applying 2nd order FD scheme
`right_controlled_volume_scheme_flux_bc`(x, k)	Compute right hand side flux boundary condition residual in accord with controlled volume scheme.
`right_finite_difference_scheme_flux_bc`(x, k)	See `left_finite_difference_scheme_flux_bc`
`right_robin_bc`(x, k, lam[, x0])	Construct Robin (u + lam*dudx = u0 ) BC at right boundary.
`solve`()	Evokes solver
`solver_callback`(xij, *_)	Callback function that can be used by optimizers of scipy.optimize. The second argument "*_" makes sure that it still works when the optimizer calls the callback function with more than one argument. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html.
`use_standard_cell_bc`()	Interfaces at left hand side and right hand side
`use_standard_dirichlet_bc`()	Dirichlet BC for all variables at all boundaries
`use_standard_interface_bc`()	Interface at left hand side and open bulk at right hand side
`use_stern_layer_cell_bc`([implicit])	Interfaces at left hand side and right hand side, Stern layer either by prescribing linear potential regime between cell boundary and outer Helmholtz plane (OHP), or by applying Robin BC; zero flux BC on all ion species.

left_potential_dirichlet_bc
left_potential_robin_bc
right_dirichlet_bc
right_potential_dirichlet_bc
right_potential_robin_bc

property grid

property potential

property concentration

property charge_density

property x1_scaled

newton(f, xij, **kwargs)

Newton solver expects system f and initial value xij

Parameters:

f (callable) – N-valued function of N-valued vector
xij ((N,) ndarray) – N-valued initial value vector

Returns:

xij – N-valued solution vector

Return type:

(N,) ndarray

solver_callback(xij, *_): Callback function that can be used by optimizers of scipy.optimize. The second argument “*_” makes sure that it still works when the optimizer calls the callback function with more than one argument. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

discretize(): Sets up discretization scheme and initial value

initial_values(): Solves decoupled linear system to get initial potential distribution.

solve()

Evokes solver

Returns:

uij ((Ni,) ndarray) – potential at Ni grid points
nij ((M,Nij) ndarray) – concentrations of M species at Ni grid points
lamj ((L,) ndarray) – value of L Lagrange multipliers (not implemented, empty)

use_standard_interface_bc(): Interface at left hand side and open bulk at right hand side

use_standard_cell_bc(): Interfaces at left hand side and right hand side

use_stern_layer_cell_bc(implicit=False)

Interfaces at left hand side and right hand side, Stern layer either by prescribing linear potential regime between cell boundary and outer Helmholtz plane (OHP), or by applying Robin BC; zero flux BC on all ion species.

Parameters:: implicit (bool, optional) – If true, then true Robin BC are applied. Attention: if desired, domain must be cropped by manually by twice the Stern layer thickness lambda_S. Otherwise, enforces constant potential gradient across Stern layer region of thickness lambda_S. (default:False)

use_standard_dirichlet_bc(): Dirichlet BC for all variables at all boundaries

left_finite_difference_scheme_flux_bc(x, k, j0=0)

Parameters:

x ((Ni,) ndarray) – N-valued variable vector
k (int) – ion species (-1 for potential)
j0 (float) – flux of ion species k at left hand boundary

Returns:

float

Return type:

boundary condition residual

right_finite_difference_scheme_flux_bc(x, k, j0=0): See `left_finite_difference_scheme_flux_bc`

left_controlled_volume_scheme_flux_bc(x, k, j0=0)

Compute left hand side flux boundary condition residual in accord with controlled volume scheme.

Parameters:

x ((Ni,) ndarray) – N-valued variable vector
k (int) – ion species (-1 for potential)
j0 (float) – flux of ion species k at left hand boundary

Returns:

float

Return type:

boundary condition residual

right_controlled_volume_scheme_flux_bc(x, k, j0=0): Compute right hand side flux boundary condition residual in accord with controlled volume scheme. See left_controlled_volume_scheme_flux_bc

left_potential_dirichlet_bc(x, u0=0)

left_dirichlet_bc(x, k, x0=0): Construct Dirichlet BC at left boundary

right_potential_dirichlet_bc(x, x0=0)

right_dirichlet_bc(x, k, x0=0)

left_potential_robin_bc(x, lam, u0=0)

left_robin_bc(x, k, lam, x0=0)

Compute left hand side Robin (u + lam*dudx = u0 ) BC at in accord with 2nd order finite difference scheme.

Parameters:

x ((Ni,) ndarray) – N-valued variable vector
k (int) – ion species (-1 for potential)
lam (float) – BC coefficient, corresponds to Stern layer thickness if applied to potential variable in PNP problem. Here, this steric layer is assumed to constitute a region of uniform charge density and thus linear potential drop across the interface.
x0 (float) – right hand side value of BC, corresponds to potential beyond Stern layer if applied to potential variable in PNP system.

Returns:

float

Return type:

boundary condition residual

right_potential_robin_bc(x, lam, u0=0)

right_robin_bc(x, k, lam, x0=0): Construct Robin (u + lam*dudx = u0 ) BC at right boundary.

number_conservation_constraint(x, k, N0): N0: total amount of species, k: ion species

poisson_pde(x): Returns Poisson equation residual by applying 2nd order FD scheme

nernst_planck_pde(x): Returns Nernst-Planck equation residual by applying controlled volume scheme

G(x)

Non-linear system

Discretization of Poisson-Nernst-Planck system with M ion species. Implements “controlled volume” method as found in

Selbherr, Analysis and Simulation of Semiconductor Devices, Springer 1984

Parameters:: x (((M+1)*Ni,) ndarray) – system variables. 1D array of (M+1)*Ni values, where M is number of ion species, Ni number of spatial discretization points. First Ni entries are expected to contain potential, following M*Ni points contain ion concentrations.
Returns:: residual
Return type:: ((M+1)*Ni,) ndarray

property ionic_strength

Compute the system’s ionic strength from charges and concentrations.

Returns:: ionic_strength – ionic strength ( 1/2 * sum(z_i^2*c_i) ) [concentration unit, i.e. mol m^-3]
Return type:: float

property lambda_D

Compute the system’s Debye length.

Returns:: lambda_D – Debye length, sqrt( epsR*eps*R*T/(2*F^2*I) ) [length unit, i.e. m]
Return type:: float

init(c=array([0.1, 0.1]), z=array([1, -1]), L=1e-07, lambda_S=0, x0=0, T=298.15, delta_u=0.05, relative_permittivity=79, vacuum_permittivity=8.8541878128e-12, R=8.314462618, F=96485.33212, N=200, e=1e-10, maxit=20, solver=None, options=None, potential0=None, concentration0=None, **kwarsg)

Initializes a 1D Poisson-Nernst-Planck system description.

Expects quantities in SI units per default.

Parameters:

c ((M,) ndarray, optional) – bulk concentrations of each ionic species [mol/m^3] (default: [ 0.1, 0.1 ])
z ((M,) ndarray, optional) – charge of each ionic species [1] (default: [ +1, -1 ])
x0 (float, optional) – left hand side reference position (default: 0)
L (float, optional) – 1D domain size [m] (default: 100e-9)
lambda_S (float, optional) – Stern layer thickness in case of Robin BC [m] (default: 0)
T (float, optional) – temperature of the solution [K] (default: 298.15)
delta_u (float, optional) – potential drop across 1D cell [V] (default: 0.05)
relative_permittivity (float, optional) – relative permittivity of the ionic solution [1] (default: 79)
vacuum_permittivity (float, optional) – vacuum permittivity [F m^-1] (default: 8.854187817620389e-12 )
R (float, optional) – molar gas constant [J mol^-1 K^-1] (default: 8.3144598)
F (float, optional) – Faraday constant [C mol^-1] (default: 96485.33289)
N (int, optional) – number of discretization grid segments (default: 200)
e (float, optional) – absolute tolerance for Newton solver convergence (default: 1e-10)
maxit (int, optional) – maximum number of Newton iterations (default: 20)
solver (func( func(x), x0), optional) – solver to use (default: None, will use own simple Newton solver)
potential0 ((N+1,) ndarray, optional (default: None)) – potential initial values
concentration0 ((M,N+1) ndarray, optional (default: None)) – concentration initial values

__init__(*args, **kwargs): Constructor, see init doc string for arguments.