Nernst-Planck example scripts are broken (both 2d and 3d)
How did this happen? Don't we test the `scripts/` folder?
diff --git a/scripts/example_nernst_planck2d.py b/scripts/example_nernst_planck2d.py
index 91aac6fa6..fbdbbcdb7 100644
--- a/scripts/example_nernst_planck2d.py
+++ b/scripts/example_nernst_planck2d.py
@@ -12,13 +12,11 @@ from openpnm.phases import mixtures
 import numpy as np
 
 
+np.random.seed(0)
 ws = op.Workspace()
 proj = ws.new_project()
-export = False
-
-# network, geometry, phase
-np.random.seed(0)
 
+# Create network, geometry, phase
 net = op.network.Cubic(shape=[23, 15, 1], spacing=1e-6, project=proj)
 prs = (net['pore.back'] * net['pore.right'] + net['pore.back']
        * net['pore.left'] + net['pore.front'] * net['pore.right']
@@ -34,8 +32,7 @@ np.random.seed(0)
 op.topotools.reduce_coordination(net, 3)
 
 np.random.seed(0)
-geo = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
-
+geo = op.geometry.CirclesAndRectangles(network=net, pores=net.Ps, throats=net.Ts)
 
 sw = mixtures.SalineWater(network=net)
 # Retrieve handles to each species for use below
@@ -43,26 +40,26 @@ Na = sw.components['Na_' + sw.name]
 Cl = sw.components['Cl_' + sw.name]
 H2O = sw.components['H2O_' + sw.name]
 
-# physics
+# Create physics
 phys = op.physics.GenericPhysics(network=net, phase=sw, geometry=geo)
 
-flow = op.models.physics.hydraulic_conductance.hagen_poiseuille_2d
+flow = op.models.physics.hydraulic_conductance.hagen_poiseuille
 phys.add_model(propname='throat.hydraulic_conductance',
                pore_viscosity='pore.viscosity',
                throat_viscosity='throat.viscosity',
                model=flow, regen_mode='normal')
 
-current = op.models.physics.ionic_conductance.electroneutrality_2d
+current = op.models.physics.ionic_conductance.electroneutrality
 phys.add_model(propname='throat.ionic_conductance', ions=[Na.name, Cl.name],
                model=current, regen_mode='normal')
 
-eA_dif = op.models.physics.diffusive_conductance.ordinary_diffusion_2d
+eA_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
 phys.add_model(propname='throat.diffusive_conductance.' + Na.name,
                pore_diffusivity='pore.diffusivity.' + Na.name,
                throat_diffusivity='throat.diffusivity.' + Na.name,
                model=eA_dif, regen_mode='normal')
 
-eB_dif = op.models.physics.diffusive_conductance.ordinary_diffusion_2d
+eB_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
 phys.add_model(propname='throat.diffusive_conductance.' + Cl.name,
                pore_diffusivity='pore.diffusivity.' + Cl.name,
                throat_diffusivity='throat.diffusivity.' + Cl.name,
@@ -79,12 +76,12 @@ phys.add_model(propname='throat.ad_dif_mig_conductance.' + Cl.name,
                pore_pressure='pore.pressure', model=ad_dif_mig_Cl,
                ion=Cl.name, s_scheme=scheme)
 
-# settings for algorithms
+# Settings for algorithms
 setts1 = {'solver_max_iter': 5, 'solver_tol': 1e-08, 'solver_rtol': 1e-08,
           'nlin_max_iter': 10, 'cache_A': False, 'cache_b': False}
-setts2 = {'g_tol': 1e-4, 'g_max_iter': 100}
+setts2 = {'g_tol': 1e-4, 'g_max_iter': 50}
 
-# algorithms
+# Algorithms
 sf = op.algorithms.StokesFlow(network=net, phase=sw, settings=setts1)
 sf.set_value_BC(pores=net.pores('back'), values=11)
 sf.set_value_BC(pores=net.pores('front'), values=10)
@@ -94,7 +91,7 @@ sw.update(sf.results())
 p = op.algorithms.IonicConduction(network=net, phase=sw, settings=setts1)
 p.set_value_BC(pores=net.pores('left'), values=0.02)
 p.set_value_BC(pores=net.pores('right'), values=0.01)
-p.settings['charge_conservation'] = 'electroneutrality_2D'
+p.settings['charge_conservation'] = 'poisson'
 
 eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion=Na.name,
                                 settings=setts1)
@@ -108,7 +105,8 @@ eB.set_value_BC(pores=net.pores('front'), values=10)
 
 pnp = op.algorithms.NernstPlanckMultiphysicsSolver(network=net, phase=sw,
                                                    settings=setts2)
-pnp.setup(potential_field=p.name, ions=[eA.name, eB.name])
+pnp.settings['potential_field'] = p.name
+pnp.settings['ions'] = [eA.name, eB.name]
 pnp.run()
 
 sw.update(sf.results())
@@ -116,6 +114,5 @@ sw.update(p.results())
 sw.update(eA.results())
 sw.update(eB.results())
 
-# output data to Paraview
-if export:
-    proj.export_data(phases=[sw], filename='out', filetype='xdmf')
+# Output data to Paraview
+proj.export_data(phases=[sw], filename='out', filetype='xdmf')
diff --git a/scripts/example_nernst_planck3d.py b/scripts/example_nernst_planck3d.py
index d7a3cd6b7..9108e7447 100644
--- a/scripts/example_nernst_planck3d.py
+++ b/scripts/example_nernst_planck3d.py
@@ -14,9 +14,8 @@ import numpy as np
 
 ws = op.Workspace()
 proj = ws.new_project()
-export = False
 
-# network, geometry, phase
+# Create network, geometry, phase
 np.random.seed(0)
 
 net = op.network.Cubic(shape=[23, 15, 1], spacing=1e-6, project=proj)
@@ -36,7 +35,6 @@ op.topotools.reduce_coordination(net, 3)
 np.random.seed(0)
 geo = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
 
-
 sw = mixtures.SalineWater(network=net)
 # Retrieve handles to each species for use below
 Na = sw.components['Na_' + sw.name]
@@ -79,12 +77,12 @@ phys.add_model(propname='throat.ad_dif_mig_conductance.' + Cl.name,
                pore_pressure='pore.pressure', model=ad_dif_mig_Cl,
                ion=Cl.name, s_scheme=scheme)
 
-# settings for algorithms
+# Settings for algorithms
 setts1 = {'solver_max_iter': 5, 'solver_tol': 1e-08, 'solver_rtol': 1e-08,
           'nlin_max_iter': 10, 'cache_A': False, 'cache_b': False}
-setts2 = {'g_tol': 1e-4, 'g_max_iter': 100}
+setts2 = {'g_tol': 1e-4, 'g_max_iter': 50}
 
-# algorithms
+# Algorithms
 sf = op.algorithms.StokesFlow(network=net, phase=sw, settings=setts1)
 sf.set_value_BC(pores=net.pores('back'), values=11)
 sf.set_value_BC(pores=net.pores('front'), values=10)
@@ -94,7 +92,7 @@ sw.update(sf.results())
 p = op.algorithms.IonicConduction(network=net, phase=sw, settings=setts1)
 p.set_value_BC(pores=net.pores('left'), values=0.02)
 p.set_value_BC(pores=net.pores('right'), values=0.01)
-p.settings['charge_conservation'] = 'electroneutrality'
+p.settings['charge_conservation'] = 'poisson'
 
 eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion=Na.name,
                                 settings=setts1)
@@ -108,7 +106,8 @@ eB.set_value_BC(pores=net.pores('front'), values=10)
 
 pnp = op.algorithms.NernstPlanckMultiphysicsSolver(network=net, phase=sw,
                                                    settings=setts2)
-pnp.setup(potential_field=p.name, ions=[eA.name, eB.name])
+pnp.settings['potential_field'] = p.name
+pnp.settings['ions'] = [eA.name, eB.name]
 pnp.run()
 
 sw.update(sf.results())
@@ -117,5 +116,4 @@ sw.update(eA.results())
 sw.update(eB.results())
 
 # output data to Paraview
-if export:
-    proj.export_data(phases=[sw], filename='out', filetype='xdmf')
+proj.export_data(phases=[sw], filename='out', filetype='xdmf')
['scripts/example_nernst_planck3d.py', 'scripts/example_nernst_planck2d.py']