matscipy.electrochemistry.poisson_nernst_planck_solver_fenics

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

Copyright 2019 IMTEK Simulation University of Freiburg

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.

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

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

  3. 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)

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()

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