<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Drift-diffusion simulation of electrochemical doping
This notebook reproduces published transient simulations of electrochemical doping.
Note that the simulation may not work for longer times than investigated in the paper because of the formation of an EDL which is poorly handled without further assumptions.
Also, all the profiles of electrical potential are inverted in the paper.
Step1:  Model and parameters
The model is multilayer, with two layers
Step2:  Doping of PEDOT
Step3:  BCs for Poisson's equation
Step4:  BCs for holes and ions
Step5:  Coupled model
Step6:  Parameters
Step7:  Results
Figure numbers given below correspond to Figures in the reference.
Fig. 2
Step8:  Fig. 3
Step9:  Fig. 4
Step10:  Fig. 5
Step11:  Fig. 6
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pylab as plt
from oedes import *
init_notebook()
Explanation: Drift-diffusion simulation of electrochemical doping
This notebook reproduces published transient simulations of electrochemical doping.
Note that the simulation may not work for longer times than investigated in the paper because of the formation of an EDL which is poorly handled without further assumptions.
Also, all the profiles of electrical potential are inverted in the paper.
End of explanation
electrolyte_mesh = fvm.mesh1d(
    25e-9,
    boundary_names=[
        'electrode0',
        'internal_boundary'])
pedot_mesh = fvm.mesh1d(
    175e-9,
    boundary_names=[
        'internal_boundary',
        'electrode1'])
mesh = fvm.multilayer1d(
    [('electrolyte', electrolyte_mesh), ('pedot', pedot_mesh)])
electrolyte_model = models.BaseModel()
electrolyte_model.name = 'electrolyte'
pedot_model = models.BaseModel()
pedot_model.name = 'pedot'
ions_electrolyte=models.add_ions(electrolyte_model,electrolyte_mesh)
ions_pedot=models.add_ions(pedot_model,pedot_mesh)
hole_eqn=models.add_transport(pedot_model,pedot_mesh,z=1,name='hole')
Explanation: Model and parameters
The model is multilayer, with two layers: electrolyte and pedot. In electrolyte, cations and anions are allowed. In pedot, cations, anions and holes are simulated.
End of explanation
pedot_model.species.append(models.FixedCharge(poisson=pedot_model.poisson))
Explanation: Doping of PEDOT:
End of explanation
electrolyte_model.poisson.bc = [
    models.AppliedVoltage('electrode0'),
    models.Equal(
        pedot_model.poisson,
        'internal_boundary')]
pedot_model.poisson.bc = [models.AppliedVoltage('electrode1')]
Explanation: BCs for Poisson's equation:
End of explanation
hole_eqn.bc = [models.DirichletFromParams('electrode1')]
for ion_electrolyte in ions_electrolyte:
    ion_electrolyte.bc = [models.DirichletFromParams('electrode0')]
for ion_pedot,ion_electrolyte in zip(ions_pedot,ions_electrolyte):
    ion_pedot.bc = [models.Equal(ion_electrolyte,'internal_boundary')]
Explanation: BCs for holes and ions:
End of explanation
model = models.CompositeModel()
model.sub = [electrolyte_model, pedot_model]
model.other = [models.RamoShockleyCurrentCalculation([electrolyte_model.poisson,pedot_model.poisson])]
model.setUp()
Explanation: Coupled model:
End of explanation
params = {
    'electrolyte.epsilon_r': 3.,
    'pedot.epsilon_r': 3,
    'electrolyte.T': 300,
    'pedot.T': 300,
    'electrolyte.electrode0.voltage': 0,
    'electrolyte.electrode0.workfunction': 0,
    'pedot.electrode1.voltage': 0,
    'pedot.electrode1.workfunction': 0,
    'pedot.hole.N0': 1e27,
    'pedot.hole.electrode1': 1e25,
    'pedot.hole.mu': 1e-7,
    'pedot.hole.energy': 0,
    'pedot.Na': 1e25,
    'pedot.Nd': 0,
    'electrolyte.cation.electrode0': 1e25,
    'electrolyte.cation.mu': 1e-7,
    'electrolyte.cation.N0': 1e27,
    'electrolyte.cation.energy': 0.,
    'electrolyte.anion.electrode0': 1e25,
    'electrolyte.anion.mu': 1e-7,
    'electrolyte.anion.N0': 1e27,
    'electrolyte.anion.energy': 0.,
    'pedot.cation.mu': 1e-9,
    'pedot.cation.N0': 1e27,
    'pedot.cation.energy': 0.,
    'pedot.anion.mu': 1e-9,
    'pedot.anion.N0': 1e27,
    'pedot.anion.energy': 0.
}
def store(out):
    testing.store(out['electrolyte.potential'])
    testing.store(out['electrolyte.anion.c'], atol=1e12)
    testing.store(out['electrolyte.cation.c'], atol=1e12)
    testing.store(out['pedot.potential'], atol=1e12)
    testing.store(out['pedot.hole.c'], atol=1e12)
    testing.store(out['pedot.anion.c'], atol=1e12)
    testing.store(out['pedot.cation.c'], atol=1e12)
Explanation: Parameters:
End of explanation
x = model.X.copy()
x[hole_eqn.idx] = params['pedot.Na']
c = context(model, x=x)
params['electrolyte.electrode0.voltage'] = 0.
params['pedot.electrode1.voltage'] = -1
c.transient(params, 1e-5, 1e-9)
from oedes import mpl
fig, (ion_plot, hole_plot, potential_plot) = plt.subplots(
    nrows=3, figsize=(6, 8))
for time in [2e-6, 4e-6, 6e-6, 8e-6, 10e-6]:
    ts = c.attime(time)
    store(ts.output())
    ts.mpl(fig, ion_plot).species('cation', settings=mpl.Settings.c_x_nm)
    ts.mpl(fig, ion_plot).species('anion', settings=mpl.Settings.c_x_nm)
    o = ts.output()
    ts.mpl(fig, hole_plot).species('hole')
    ts.mpl(fig, potential_plot).potential()
fig.tight_layout()
Explanation: Results
Figure numbers given below correspond to Figures in the reference.
Fig. 2
End of explanation
def length(out):
    store(out)
    c = out['pedot.hole.c']
    i = np.argwhere(c > 0.5e25)[0]
    return pedot_mesh.cells['center'][i] - pedot_mesh.cells['center'][0]
times = np.asarray([2e-6, 4e-6, 6e-6, 8e-6, 10e-6])
plt.plot(times * 1e6, [length(c.attime(t).output())
                       ** 2 * 1e18 for t in times], '-o')
plt.xlabel('t [$\mu s$]')
plt.ylabel('$\ell^2$ [$nm%2$]')
t, j = c.teval('time', 'electrode1.J')
t, j = t[t > 0], j[t > 0]
plt.plot(1. / np.sqrt(t * 1e6), j * 1e-4, 'o-')
plt.xlim([0.25, 1.1])
plt.ylim([0.75, 4])
plt.xlabel('$t^{-1/2}$ [$\mu s^{-1/2}$]')
plt.ylabel('j [$A/cm^2$]')
Explanation: Fig. 3
End of explanation
x0 = model.X.copy()
x0[hole_eqn.idx] = params['pedot.Na']
by_v = {}
voltages = [0.01, 0.05, 0.1, 0.5, 1]
for v in voltages:
    params['pedot.electrode1.voltage'] = -v
    c = context(model, x=x0)
    c.transient(params, 1e-5, 1e-9)
    by_v[v] = c
fig, (ion_plot, potential_plot) = plt.subplots(ncols=2, figsize=(12, 4))
for v in voltages:
    o = by_v[v].attime(1e-5).output()
    store(o)
    ion_plot.plot(
        electrolyte_mesh.cells['center'] *
        1e9,
        o['electrolyte.cation.c'] -
        o['electrolyte.anion.c'])
    ion_plot.plot(
        pedot_mesh.cells['center'] *
        1e9,
        o['pedot.cation.c'] -
        o['pedot.anion.c'])
    ion_plot.set_xlabel('x [nm]')
    ion_plot.set_ylabel('P-N [1/m^3]')
    potential_plot.plot(
        electrolyte_mesh.cells['center'] * 1e9, o['electrolyte.potential'] / (-v))
    potential_plot.plot(
        pedot_mesh.cells['center'] * 1e9, o['pedot.potential'] / (-v))
    potential_plot.set_xlabel('x [nm]')
    potential_plot.set_ylabel('$V/V_{appl}$')
Explanation: Fig. 4
End of explanation
x0 = model.X.copy()
x0[hole_eqn.idx] = params['pedot.Na']
params['pedot.electrode1.voltage'] = -1.
concentrations = [1e24, 1e25, 1e26]
by_c = {}
for conc in concentrations:
    params['electrolyte.cation.electrode0'] = conc
    params['electrolyte.anion.electrode0'] = conc
    c = context(model, x=x0)
    c.transient(params, 1e-5, 1e-9)
    by_c[conc] = c
fig, (ion_plot, potential_plot) = plt.subplots(ncols=2, figsize=(12, 4))
for conc in concentrations:
    o = by_c[conc].attime(1e-5).output()
    store(o)
    ion_plot.plot(
        electrolyte_mesh.cells['center'] *
        1e9,
        o['electrolyte.cation.c'] -
        o['electrolyte.anion.c'])
    ion_plot.plot(
        pedot_mesh.cells['center'] *
        1e9,
        o['pedot.cation.c'] -
        o['pedot.anion.c'])
    ion_plot.set_xlabel('x [nm]')
    ion_plot.set_ylabel('P-N [1/m^3]')
    potential_plot.plot(
        electrolyte_mesh.cells['center'] * 1e9, -o['electrolyte.potential'])
    potential_plot.plot(
        pedot_mesh.cells['center'] * 1e9, -o['pedot.potential'])
    potential_plot.set_xlabel('x [nm]')
    potential_plot.set_ylabel('- electrostatic potential [V]')
Explanation: Fig. 5
End of explanation
x0 = model.X.copy()
x0[hole_eqn.idx] = params['pedot.Na']
params['pedot.electrode1.voltage'] = -1.
params['electrolyte.cation.electrode'] = 1e25
params['electrolyte.anion.electrode'] = 1e25
mu_values = [1e-10, 1e-9, 1e-8]
by_mu = {}
for mu in mu_values:
    params['pedot.cation.mu'] = mu
    params['pedot.anion.mu'] = mu
    c = context(model, x=x0)
    c.transient(params, 1e-5, 1e-9)
    by_mu[mu] = c
fig, (ion_plot, potential_plot) = plt.subplots(ncols=2, figsize=(12, 4))
time = {1e-10: 6e-6, 1e-9: 6e-6, 1e-8: 1e-6}
for mu in mu_values:
    o = by_mu[mu].attime(time[mu]).output()
    store(o)
    ion_plot.plot(
        electrolyte_mesh.cells['center'] *
        1e9,
        o['electrolyte.cation.c'] -
        o['electrolyte.anion.c'])
    ion_plot.plot(
        pedot_mesh.cells['center'] *
        1e9,
        o['pedot.cation.c'] -
        o['pedot.anion.c'])
    ion_plot.set_xlabel('x [nm]')
    ion_plot.set_ylabel('P-N [1/m^3]')
    potential_plot.plot(
        electrolyte_mesh.cells['center'] * 1e9, -o['electrolyte.potential'] / v)
    potential_plot.plot(
        pedot_mesh.cells['center'] * 1e9, -o['pedot.potential'] / v)
    potential_plot.set_xlabel('x [nm]')
    potential_plot.set_ylabel('- electrostatic potential [V]')
Explanation: Fig. 6
End of explanation
<END_TASK>