matscipy.electrochemistry.poisson_nernst_planck_solver_fenics

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

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

Classes

`Boundary`([x0, tol])	Mark a point to be within the domain boundary for fenics.
`PoissonNernstPlanckSystemFEniCS`(args, *kwargs)	Describes and solves a 1D Poisson-Nernst-Planck system, using log concentrations internally

class matscipy.electrochemistry.poisson_nernst_planck_solver_fenics.Boundary(x0=0, tol=1e-14)

Bases: SubDomain

Mark a point to be within the domain boundary for fenics.

Methods

`get_property`(self, arg0)
`inside`(x, on_boundary)	Mark a point to be within the domain boundary for fenics.
`map`(self, arg0, arg1)
`mark`(args, *kwargs)	Overloaded function.
`mark_cells`(self, mesh, sub_domain[, ...])
`mark_facets`(self, mesh, sub_domain[, ...])
`set_property`(self, arg0, arg1)
`snap`(self, arg0)

__init__(self: dolfin.cpp.mesh.SubDomain, map_tol: float = 1e-10) → None

inside(x, on_boundary): Mark a point to be within the domain boundary for fenics.

get_property(self: dolfin.cpp.mesh.SubDomain, arg0: str) → float

map(self: dolfin.cpp.mesh.SubDomain, arg0: numpy.ndarray[numpy.float64[m, 1]], arg1: numpy.ndarray[numpy.float64[m, 1], flags.writeable]) → None

mark(*args, **kwargs)

Overloaded function.

mark(self: dolfin.cpp.mesh.SubDomain, meshfunction: dolfin.cpp.mesh.MeshFunctionSizet, marker: int, check_midpoint: bool = True) -> None
mark(self: dolfin.cpp.mesh.SubDomain, meshfunction: dolfin.cpp.mesh.MeshFunctionDouble, marker: float, check_midpoint: bool = True) -> None
mark(self: dolfin.cpp.mesh.SubDomain, meshfunction: dolfin.cpp.mesh.MeshFunctionBool, marker: bool, check_midpoint: bool = True) -> None

mark_cells(self: dolfin.cpp.mesh.SubDomain, mesh: dolfin.cpp.mesh.Mesh, sub_domain: int, check_midpoint: bool = True) → None

mark_facets(self: dolfin.cpp.mesh.SubDomain, mesh: dolfin.cpp.mesh.Mesh, sub_domain: int, check_midpoint: bool = True) → None

set_property(self: dolfin.cpp.mesh.SubDomain, arg0: str, arg1: float) → None

snap(self: dolfin.cpp.mesh.SubDomain, arg0: numpy.ndarray[numpy.float64[m, 1], flags.writeable]) → None

class matscipy.electrochemistry.poisson_nernst_planck_solver_fenics.PoissonNernstPlanckSystemFEniCS(*args, **kwargs)

Bases: PoissonNernstPlanckSystem

Describes and solves a 1D Poisson-Nernst-Planck system, using log concentrations internally

Attributes:

X
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
`apply_central_reference_concentration_constraint`(k, c0)	FEniCS Dirichlet BC c0 for k'th ion species at right boundary.
`apply_left_concentration_dirichlet_bc`(k, c0)	FEniCS Dirichlet BC c0 for k'th ion species at left boundary.
`apply_left_potential_dirichlet_bc`(u0)	FEniCS Dirichlet BC u0 for potential at left boundary.
`apply_left_potential_robin_bc`(u0, lam0)
`apply_number_conservation_constraint`(k, c0)	Enforce number conservation constraint via Lagrange multiplier.
`apply_potential_dirichlet_bc`(u0, u1)	Potential Dirichlet BC u0 and u1 on left and right boundary.
`apply_potential_robin_bc`(u0, u1, lam0, lam1)	Potential Robin BC on left and right boundary.
`apply_right_concentration_dirichlet_bc`(k, c0)	FEniCS Dirichlet BC c0 for k'th ion species at right boundary.
`apply_right_potential_dirichlet_bc`(u0)	FEniCS Dirichlet BC u0 for potential at right boundary.
`boundary_C`(x, on_boundary)	Mark domain center.
`boundary_L`(x, on_boundary)	Mark left boundary.
`boundary_R`(x, on_boundary)	Mark right boundary.
`discretize`()	Builds function space, call again after introducing constraints
`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)
`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`()	Evoke FEniCS FEM 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_central_reference_concentration_based_cell_bc`()	Interfaces at left hand side and right hand side, species-wise concentration fixed at cell center.
`use_standard_cell_bc`()	Interfaces at left hand side and right hand side, species-wise number conservation within interval.
`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`()	Interfaces at left hand side and right hand side, species-wise number conservation within interval.

apply_right_potential_robin_bc
left_potential_dirichlet_bc
left_potential_robin_bc
right_dirichlet_bc
right_potential_dirichlet_bc
right_potential_robin_bc

property X

solve()

Evoke FEniCS FEM 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

boundary_L(x, on_boundary): Mark left boundary. Returns True if x on left boundary.

boundary_R(x, on_boundary): Mark right boundary. Returns True if x on right boundary.

boundary_C(x, on_boundary): Mark domain center.

apply_left_potential_dirichlet_bc(u0): FEniCS Dirichlet BC u0 for potential at left boundary.

apply_right_potential_dirichlet_bc(u0): FEniCS Dirichlet BC u0 for potential at right boundary.

apply_left_concentration_dirichlet_bc(k, c0): FEniCS Dirichlet BC c0 for k’th ion species at left boundary.

apply_right_concentration_dirichlet_bc(k, c0): FEniCS Dirichlet BC c0 for k’th ion species at right boundary.

apply_central_reference_concentration_constraint(k, c0): FEniCS Dirichlet BC c0 for k’th ion species at right boundary.

apply_left_potential_robin_bc(u0, lam0)

apply_right_potential_robin_bc(u0, lam0)

apply_number_conservation_constraint(k, c0): Enforce number conservation constraint via Lagrange multiplier. See https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/neumann-poisson/python/documentation.html

apply_potential_dirichlet_bc(u0, u1): Potential Dirichlet BC u0 and u1 on left and right boundary.

apply_potential_robin_bc(u0, u1, lam0, lam1): Potential Robin BC on left and right boundary.

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, species-wise number conservation within interval.

use_central_reference_concentration_based_cell_bc(): Interfaces at left hand side and right hand side, species-wise concentration fixed at cell center.

use_stern_layer_cell_bc(): Interfaces at left hand side and right hand side, species-wise number conservation within interval.

discretize(): Builds function space, call again after introducing constraints

__init__(*args, **kwargs): Same parameters as PoissonNernstPlanckSystem.init

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 charge_density

property concentration

property grid

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

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

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

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

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

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

left_potential_dirichlet_bc(x, u0=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

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

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

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

property potential

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

right_dirichlet_bc(x, k, x0=0)

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

right_potential_dirichlet_bc(x, x0=0)

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.

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_dirichlet_bc(): Dirichlet BC for all variables at all boundaries

property x1_scaled