Interactive simulation of PN junction

import oedes
from oedes import models
from oedes.ad import where
import numpy as np

Specification of device model and physical parameters

def doping_profile(mesh, ctx, eq):
    return where(mesh.x<mesh.length*0.5,ctx.param(eq, 'Nd'),- ctx.param(eq,'Na'))
poisson = models.PoissonEquation()
temperature = models.ConstTemperature()
electron = models.BandTransport(poisson=poisson, name='electron', z=-1, thermal=temperature)
hole = models.BandTransport(poisson=poisson, name='hole', z=1, thermal=temperature)
doping = models.FixedCharge(poisson, density=doping_profile)
semiconductor = models.Electroneutrality([electron, hole, doping],name='semiconductor')
recombination = models.DirectRecombination(semiconductor)
anode = models.OhmicContact(poisson, semiconductor, 'electrode0')
cathode = models.OhmicContact(poisson, semiconductor, 'electrode1')
current = models.RamoShockleyCurrentCalculation([poisson])
equations=[ poisson, temperature, electron, hole,
            doping, current, semiconductor,
            anode, cathode, recombination ]
Na=1e24
Nd=1e24

params={
    'T':300,
    'epsilon_r':12,
    'Na':1e24,
    'Nd':1e24,
    'hole.mu':1,
    'electron.mu':1,
    'hole.energy':-1.1,
    'electron.energy':0,
    'electrode0.voltage':0,
    'electrode1.voltage':0,
    'hole.N0':1e27,
    'electron.N0':1e27,
    'beta':1e-9
}

Interactive simulation

# Create discrete model
mesh = oedes.fvm.mesh1d(100e-9)
model = oedes.fvm.discretize(equations, mesh)

# Precalculate equilibrium solution to save time in interative simulation
c_eq=oedes.context(model)
c_eq.solve(params)
import ipywidgets
import pylab as plt

fig=plt.figure()
ax=plt.subplot()

voltage_selection = ipywidgets.FloatSlider(min=-2.,max=2.,step=0.025,description='Voltage')
plot_selection = ipywidgets.Select(options=['energy','potential','current density','concentrations'],description='Plot')

def simulation_at_voltage(v=voltage_selection,
                plot_what=plot_selection):
    params['electrode1.voltage']=v
    c=oedes.context(model,x=np.asarray(c_eq.x,dtype=np.longdouble))
    c.solve(params)
    ax.clear()
    p=c.mpl(fig,ax)
    if plot_what == 'energy':
        p.plot(['electron.Eband'],label='$E_c$')
        p.plot(['hole.Eband'],label='$E_v$')
        p.plot(['electron.Ef'],linestyle='--',label='$E_{Fn}$')
        p.plot(['hole.Ef'],linestyle='-.',label='$E_{Fp}$')
    elif plot_what == 'potential':
        p.plot(['potential'],label='$\psi$')
        p.plot(['electron.phi_f'],label='$\phi_n$')
        p.plot(['hole.phi_f'],label='$\phi_p$')
    elif plot_what == 'current density':
        p.plot(['electron.J'],label='$J_n$')
        p.plot(['hole.J'],label='$J_p$')
    elif plot_what == 'concentrations':
        p.plot(['electron.c'])
        p.plot(['hole.c'])
        p.apply_settings({'yscale':'log'})
    p.apply_settings({'xunit':'n','xlabel':'nm'})
    ax.legend(loc=0,frameon=False)
    display(fig)
ipywidgets.interact(simulation_at_voltage);