Finalize the Nernst-Planck algs

example_PNP.py

@@ -1,9 +1,8 @@
-"""
 import openpnm as op
 import numpy as np
 ws = op.Workspace()
 proj = ws.new_project()
-# ws.settings['loglevel'] = 20
+
 
 # network, geometry, phase
 np.random.seed(0)
@@ -14,10 +13,6 @@
  net['pore.front'] * net['pore.left'])
 thrts = net['throat.surface']
 op.topotools.trim(network=net, pores=net.Ps[prs], throats=net.Ts[thrts])
-net['pore.coords'][:, 0] -= 0.00045
-net['pore.coords'][:, 1] -= 0.00045
-net['pore.coords'][:, 2] = 0
-# np.savetxt('coords.txt', net['pore.coords'])
 
 
 geo = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
@@ -28,8 +23,6 @@
 geo.regenerate_models()
 
 sw = op.phases.Mixtures.SalineWater(network=net)
-sw['pore.diffusivity.Na'] = sw['pore.diffusivity.Cl']
-sw['throat.diffusivity.Na'] = sw['throat.diffusivity.Cl']
 
 # physics
 phys = op.physics.GenericPhysics(network=net, phase=sw, geometry=geo)
@@ -67,17 +60,17 @@
 p = op.algorithms.OhmicConduction(network=net, phase=sw)
 p.settings['conductance'] = 'throat.ionic_conductance'
 p.settings['quantity'] = 'pore.potential'
-p.set_value_BC(pores=net.pores('left'), values=0.05)
+p.set_value_BC(pores=net.pores('left'), values=0.01)
 p.set_value_BC(pores=net.pores('right'), values=0.00)
 p.settings['rxn_tolerance'] = 1e-12
-sw.update(p.results())
 
-eA = op.algorithms.NernstPlanck(network=net, phase=sw, electrolyte='Na')
+
+eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion='Na')
 eA.set_value_BC(pores=net.pores('back'), values=100)
 eA.set_value_BC(pores=net.pores('front'), values=90)
 eA.settings['rxn_tolerance'] = 1e-12
 
-eB = op.algorithms.NernstPlanck(network=net, phase=sw, electrolyte='Cl')
+eB = op.algorithms.NernstPlanck(network=net, phase=sw, ion='Cl')
 eB.set_value_BC(pores=net.pores('back'), values=100)
 eB.set_value_BC(pores=net.pores('front'), values=90)
 eB.settings['rxn_tolerance'] = 1e-12
@@ -93,27 +86,14 @@
  model=ad_dif_mig_Cl, ion='Cl',
  s_scheme='exponential')
 
-pnp = op.algorithms.PoissonNernstPlanck(network=net, phase=sw)
-pnp.setup(potential_field=p, electrolytes=[eA, eB])
-pnp.settings['max_iter'] = 400
+pnp = op.algorithms.ChargeConservationNernstPlanck(network=net, phase=sw)
+pnp.setup(potential_field=p, ions=[eA, eB])
+pnp.settings['max_iter'] = 10
 pnp.settings['tolerance'] = 1e-04
-# pnp.settings['charge_conservation'] = 'laplace'
-pnp.settings['charge_conservation'] = 'electroneutrality'
+pnp.settings['charge_conservation'] = 'laplace'
 pnp.run()
 
-# Comsol results
-# cNa_cmsl = np.genfromtxt('cNa_at_coords.txt', skip_header=9)
-# sw['pore.concentration.cmsl'] = cNa_cmsl[:, 3]
-
-# phi_cmsl = np.genfromtxt('phi_at_coords.txt', skip_header=9)
-# sw['pore.potential_cmsl'] = phi_cmsl[:, 3]
-
-# p_cmsl = np.genfromtxt('p_at_coords.txt', skip_header=9)
-# sw['pore.pressure_cmsl'] = p_cmsl[:, 3]
-
 sw.update(sf.results())
 sw.update(p.results())
 sw.update(eA.results())
 sw.update(eB.results())
-# op.io.VTK.save(network=net, phases=[sw], filename='OUTPUT_PNP')
-"""

openpnm/algorithms/PoissonNernstPlanck.py → openpnm/algorithms/ChargeConservationNernstPlanck.py

@@ -5,14 +5,15 @@
 logger = logging.getLogger(__name__)
 
 
-class PoissonNernstPlanck(ReactiveTransport):
+class ChargeConservationNernstPlanck(ReactiveTransport):
  r"""
- A subclass of GenericTransport to solve the Poisson Nernst-Planck equations
+ A subclass of GenericTransport to solve the charge conservation and
+ Nernst-Planck equations.
  """
  def __init__(self, settings={}, phase=None, **kwargs):
  def_set = {'phase': None,
  'potential_field': None,
- 'electrolytes': [],
+ 'ions': [],
  'charge_conservation': 'electroneutrality',
  'tolerance': 1e-4,
  'max_iter': 10}
@@ -22,7 +23,7 @@ def __init__(self, settings={}, phase=None, **kwargs):
  if phase is not None:
  self.setup(phase=phase)
 
- def setup(self, phase=None, potential_field=None, electrolytes=[],
+ def setup(self, phase=None, potential_field=None, ions=[],
  charge_conservation=None, tolerance=None,
  max_iter=None, **kwargs):
  r"""
@@ -31,8 +32,8 @@ def setup(self, phase=None, potential_field=None, electrolytes=[],
  self.settings['phase'] = phase.name
  if potential_field:
  self.settings['potential_field'] = potential_field
- if electrolytes:
- self.settings['electrolytes'] = electrolytes
+ if ions:
+ self.settings['ions'] = ions
  if charge_conservation:
  self.settings['charge_conservation'] = charge_conservation
  if tolerance:
@@ -42,25 +43,32 @@ def setup(self, phase=None, potential_field=None, electrolytes=[],
  super().setup(**kwargs)
 
  def run(self):
- # Phase, potential and electrolytes algorithms
+ # Phase, potential and ions algorithms
  phase = self.project.phases()[self.settings['phase']]
  p_alg = self.settings['potential_field']
- e_alg = self.settings['electrolytes']
+ e_alg = self.settings['ions']
 
  # Define initial conditions (if not defined by the user)
  try:
  p_alg[p_alg.settings['quantity']]
  except KeyError:
- p_alg[p_alg.settings['quantity']] = np.zeros(shape=[p_alg.Np, ],
- dtype=float)
+ try:
+ p_alg[p_alg.settings['quantity']] = (
+ phase[p_alg.settings['quantity']])
+ except KeyError:
+ p_alg[p_alg.settings['quantity']] = np.zeros(
+ shape=[p_alg.Np, ], dtype=float)
  for e in e_alg:
  try:
  e[e.settings['quantity']]
  except KeyError:
- e[e.settings['quantity']] = np.zeros(shape=[e.Np, ],
- dtype=float)
+ try:
+ e[e.settings['quantity']] = phase[e.settings['quantity']]
+ except KeyError:
+ e[e.settings['quantity']] = np.zeros(shape=[e.Np, ],
+ dtype=float)
 
- # Define source term for Poisson equation
+ # Source term for Poisson or charge conservation (electroneutrality) eq
  Ps = (p_alg['pore.all'] * np.isnan(p_alg['pore.bc_value']) *
  np.isnan(p_alg['pore.bc_rate']))
  mod = gst.charge_conservation
@@ -80,8 +88,7 @@ def run(self):
  # Iterate until solutions converge
  for itr in range(int(self.settings['max_iter'])):
  r = str([float(format(i, '.3g')) for i in res.values()])[1:-1]
- logger.info('Iter: ' + str(itr) + ', Residuals: ' + r)
- print('Iter: ' + str(itr) + ', Residuals: ' + r)
+ logger.info('Iter: ' + str(itr+1) + ', Residuals: ' + r)
  convergence = max(i for i in res.values()) < tol
  if not convergence:
  # Poisson eq
@@ -95,7 +102,7 @@ def run(self):
  phase.update(p_alg.results())
  phys[0].regenerate_models()
 
- # Electrolytes
+ # Ions
  for e in e_alg:
  c_old = e[e.settings['quantity']].copy()
  e._run_reactive(x=c_old)

openpnm/algorithms/NernstPlanck.py

@@ -5,13 +5,16 @@
 
 class NernstPlanck(ReactiveTransport):
  r"""
+ A class to simulate transport of charged species (such as ions) in dilute
+ solutions.
+
  """
- def __init__(self, settings={}, phase=None, electrolyte='', **kwargs):
+ def __init__(self, settings={}, phase=None, ion='', **kwargs):
  def_set = {'phase': None,
- 'quantity': 'pore.concentration.'+electrolyte,
- 'conductance': 'throat.ad_dif_mig_conductance.'+electrolyte}
+ 'quantity': 'pore.concentration.'+ion,
+ 'conductance': 'throat.ad_dif_mig_conductance.'+ion}
  super().__init__(**kwargs)
- self.name = electrolyte
+ self.name = ion
  self.settings.update(def_set)
  self.settings.update(settings)
  if phase is not None:

openpnm/algorithms/__init__.py

@@ -27,4 +27,4 @@
 from .MixedInvasionPercolation import MixedInvasionPercolation
 from .Porosimetry import Porosimetry
 from .NernstPlanck import NernstPlanck
-from .PoissonNernstPlanck import PoissonNernstPlanck
+from .ChargeConservationNernstPlanck import ChargeConservationNernstPlanck

openpnm/models/physics/ad_dif_mig_conductance.py

@@ -233,9 +233,15 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
  T = phase[throat_temperature][throats]
  except KeyError:
  T = phase.interpolate_data(propname=pore_temperature)[throats]
-
- P = phase[pore_pressure]
- V = phase[pore_potential]
+ # Check if pressure and potential values exist, otherwise, assign zeros
+ try:
+ P = phase[pore_pressure]
+ except KeyError:
+ P = _sp.zeros(shape=[phase.Np, ], dtype=float)
+ try:
+ V = phase[pore_potential]
+ except KeyError:
+ V = _sp.zeros(shape=[phase.Np, ], dtype=float)
  gh = phase[throat_hydraulic_conductance]
  gd = phase[throat_diffusive_conductance]
  gd = _sp.tile(gd, 2)
@@ -244,16 +250,19 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
  F = 96485.3329
  R = 8.3145
 
- S = (A1*L1+A2*L2+At*Lt)/(L1+L2+Lt)
- L = L1 + Lt + L2
-
  # Advection
  Qij = -gh*_sp.diff(P[cn], axis=1).squeeze()
  Qij = _sp.append(Qij, -Qij)
 
  # Migration
- grad_V = _sp.diff(V[cn], axis=1).squeeze() / L
- mig = ((z*F*D*S)/(R*T)) * grad_V
+ gm1, gm2, gmt = _sp.zeros((3, len(Lt)))
+ gm1[~m1] = gm2[~m2] = gmt[~mt] = _sp.inf
+ gm1[m1] = ((z*F*D*A1)/(R*T))[m1] / L1[m1]
+ gm2[m2] = ((z*F*D*A2)/(R*T))[m2] / L2[m2]
+ gmt[mt] = ((z*F*D*At)/(R*T))[mt] / Lt[mt]
+ gm = (1/gm1 + 1/gm2 + 1/gmt)**(-1)
+ delta_V = _sp.diff(V[cn], axis=1).squeeze()
+ mig = gm * delta_V
  mig = _sp.append(mig, -mig)
 
  # Advection-migration

openpnm/models/physics/diffusive_conductance.py

@@ -337,7 +337,7 @@ def classic_ordinary_diffusion(target,
  # Get shape factor
  try:
  sf = network[shape_factor]
- except:
+ except KeyError:
  sf = _sp.ones(network.num_throats())
  sf[_sp.isnan(sf)] = 1.0
  gt = (1 / sf) * ct * DABt * tarea / tlen