Resistor (3D) - Quasi Fermi Method

   ###################################################################
   #                          Constants                              #
   ###################################################################
   
   set T 300.0
   set k 1.38066e-23
   set q 1.619e-19
   set Vt [expr {$k*$T/$q}]
   set ni 1.1e10
   set esi [expr 11.8 * 8.85418e-14]
   set eps [expr $esi / $q]
   set Emob 350.0
   set Hmob 150.0
   set small 1.0e-10
   
   ###################################################################
   #                     Define the Grid                             #
   ###################################################################
   set MeshSpac 0.10
   set Side 2.0
   line x loc=0.0   spac=$MeshSpac   tag=Top
   line x loc=$Side spac=$MeshSpac   tag=Bottom
   
   line y loc=0.0   spac=$MeshSpac   tag=Lside
   line y loc=$Side spac=$MeshSpac   tag=Rside
   
   line z loc=0.0   spac=$MeshSpac   tag=Zlo
   line z loc=1.0   spac=$MeshSpac*2 tag=Zhi
   
   mater add name=Silicon
   region Silicon xlo=Top xhi=Bottom ylo=Lside yhi=Rside zlo=Zlo zhi=Zhi
   init
   
   ###################################################################
   #                  Contacts and Doping                            #
   ###################################################################
   # place the contacts on the corners of the square
   set Contact 0.2
   set buf [ expr $MeshSpac/2]
   contact name=GND Silicon xlo=$MeshSpac-$buf xhi=$Contact+$buf ylo=$MeshSpac-$buf yhi=$Contact+$buf zlo=-$buf zhi=$buf add
   contact name=VSS Silicon xlo=$Side-$Contact-$buf xhi=$Side+$buf-$MeshSpac ylo=$Side-$Contact-$buf yhi=$Side+$buf-$MeshSpac zlo=-$buf zhi=$buf add
  
   sel z=1.0e20 name=Nd
   sel z=1.0e17 name=Na
   sel z=(Nd-Na) name=Doping
   #Plot the doping concentration
   sel z=log10(abs(Doping+1.0))
   window xwin col=2
   plot2d zvalue=0.2 levels=05 graph= Doping\ Concentration\ at\ z=0.2
    
   ###################################################################
   #                       Equations                                 #
   ###################################################################
   solution add name=DevPsi pde solve negative damp
   solution add name=Elec   pde solve !negative
   solution add name=Hole   pde solve !negative
   solution add name=Circuit solve negative damp circuit const val=0.0
   set eqnP "$eps * grad( DevPsi ) + Doping - Elec + Hole"
   set eqnE "ddt(Elec) - ($Emob) * $Vt * sgrad(Elec, DevPsi/$Vt)"
   set eqnH "ddt(Hole) - ($Hmob) * $Vt * sgrad(Hole, -DevPsi/$Vt)"
   pdbSetString Silicon DevPsi Equation $eqnP
   pdbSetString Silicon Elec   Equation $eqnE
   pdbSetString Silicon Hole   Equation $eqnH
   
   pdbSetDouble Silicon DevPsi DampValue $Vt
   pdbSetDouble Silicon DevPsi Abs.Error 1.0e-9
   pdbSetDouble Silicon Elec   Abs.Error 1.0e-5
   pdbSetDouble Silicon Hole   Abs.Error 1.0e-5
   
   #Contact Equations
   proc OhmicContact {Contact} {
       global Vt ni
       pdbSetBoolean $Contact Elec Flux 1
       pdbSetBoolean $Contact Hole Flux 1
       pdbSetBoolean $Contact DevPsi Flux 1
       pdbSetBoolean $Contact Elec Fixed 1
       pdbSetBoolean $Contact Hole Fixed 1
       pdbSetBoolean $Contact DevPsi Fixed 1
       pdbSetDouble $Contact Elec Flux.Scale 1.619e-19
       pdbSetDouble $Contact Hole Flux.Scale 1.619e-19
       pdbSetString $Contact DevPsi Equation "Nd - Na - Elec + Hole"
       pdbSetString $Contact Elec Equation "DevPsi - $Vt*log((Elec)/$ni) -$Contact"
       pdbSetString $Contact Hole Equation "DevPsi + $Vt*log((Hole)/$ni) -$Contact"
   }
   OhmicContact VSS
   OhmicContact GND
   
   #Bias contacts
   contact name=VSS voltage supply=0.0
   contact name=GND voltage supply=0.0
   
   ###################################################################
   #                         Solve                                   #
   ###################################################################
   #Initial Guess
   sel z= {(Doping>0.0)
       ?  (0.025*log( (Doping + 1e-10) / 1.1e10))
       :  (-0.025*log(-(Doping + 1e-10) / 1.1e10))} name = DevPsi
   sel z=$ni*exp(DevPsi/$Vt)  name=Elec
   sel z=$ni*exp(-DevPsi/$Vt) name=Hole
  
   #First DC solve at Equilibrium (0V)
   device
   
   puts "Electron Flux [contact name=VSS sol=Elec flux]"
   puts "Hole Flux [contact name=VSS sol=Hole flux]"
   
  
   #initialize an array, bias (Vb) and counter (i), and create and new graph window
   array set curr {}
   set i 1
   #Ramp from 0.0-1.2V and Plot
   for {set bias 0.0} {$bias < 1.2} {set bias [expr $bias+0.05]} {
       contact name=VSS supply=$bias
       device init
       set curr($i) [expr ([contact name=VSS sol=Elec flux] - [contact name=VSS sol=Hole flux])]
       chart graph=IV curve=IV xval=$bias yval=$curr($i)
       incr i
   }