Help Needed with Transient Simulation for BJT Example

models
1

Hello DEVSIM Community,

I am working on performing a transient simulation for the BJT example provided in the DEVSIM BJT Example Repository. My goal is to simulate transient behavior where the collector voltage Vc​ is incrementally increased over time.

Below is the code I am using:

import devsim
import bjt_common

def transient_simulation():
    # Uncomment the following lines for extended mode
    # devsim.set_parameter(name="extended_solver", value=True) 
    # devsim.set_parameter(name="extended_model", value=True) 
    # devsim.set_parameter(name="extended_equation", value=True)

    bjt_common.run()
    
    # Initial DC solution
    devsim.solve(type="transient_dc", absolute_error=1e6, relative_error=1e-1, maximum_iterations=40)
    devsim.circuit_alter(name="Vb", value=0.6)
    devsim.solve(type="dc", absolute_error=1e6, relative_error=1e-1, maximum_iterations=40)
    devsim.solve(type="transient_dc", absolute_error=1e6, relative_error=1e-1, maximum_iterations=40)

    # Transient simulation settings
    time_step = 1e-9  # Time step
    total_time = 1e-7  # Total simulation time
    current_time = 0.0
    target_vc = 1.0  # Final Vc value
    vc_step = target_vc / (total_time / time_step)

    while current_time < total_time:
        new_vc = current_time / time_step * vc_step  # Increment Vc
        devsim.circuit_alter(name="Vc", value=new_vc)
        
        devsim.solve(type="transient_bdf1",
                     absolute_error=1e6,
                     relative_error=1e-1,
                     maximum_iterations=100,
                     tdelta=time_step,
                     charge_error=1)
        
        current_time += time_step

if __name__ == "__main__":
    transient_simulation()

Issue: The script executes, but I am encountering different behaviors depending on the mode I use:

  • In normal mode, the simulation fails due to overflow or convergence failure.
  • In extended mode (with extended_solver, extended_model, and extended_equation parameters enabled), the simulation runs but stops progressing at a certain point without explicitly reporting a convergence failure.

Observations:

  • I am using the mesh and setup from the DEVSIM BJT Example Repository.
  • The transient simulation starts but fails to complete as expected.
  • Adding a collector resistor did not help improve convergence or stability.
  • The way I increment Vc​ over time might be contributing to instability.

Questions:

  1. Are there specific solver settings or parameters to improve the stability and convergence for this transient simulation?
  2. Is my method of incrementing Vc​ in the loop correct, or does it need adjustment?
  3. Are there additional strategies or configurations I should consider to enhance solver robustness?

I have attached the log file with further details for reference. Thank you for your time and assistance!

Best regards,
AC

Log File Output

Adding 14773 nodes
Adding 43698 edges with 0 duplicates removed
Adding 28926 triangles with 0 duplicate removed
Contact base in region bjt with 63 nodes
Contact collector in region bjt with 345 nodes
Contact emitter in region bjt with 63 nodes
Replacing Node Model Acceptors in region bjt of material Silicon
Replacing Node Model DEG in region bjt of material Silicon
Replacing Edge Model DField in region bjt of material Silicon
Replacing Edge Model DField:Potential@n0 in region bjt of material Silicon
Replacing Edge Model DField:Potential@n1 in region bjt of material Silicon
Replacing Node Model Donors in region bjt of material Silicon
Replacing Node Model EC in region bjt of material Silicon
Replacing Node Model EC:Potential in region bjt of material Silicon
Replacing Node Model EFN in region bjt of material Silicon
Replacing Node Model EFN:Electrons in region bjt of material Silicon
Replacing Node Model EFN:Potential in region bjt of material Silicon
Replacing Node Model EFP in region bjt of material Silicon
Replacing Node Model EFP:Holes in region bjt of material Silicon
Replacing Node Model EFP:Potential in region bjt of material Silicon
Replacing Edge Model EField in region bjt of material Silicon
Replacing Edge Model EField:Potential@n0 in region bjt of material Silicon
Replacing Edge Model EField:Potential@n1 in region bjt of material Silicon
Replacing Triangle Edge Model EField_x in region bjt of material Silicon
Replacing Triangle Edge Model EField_y in region bjt of material Silicon
Replacing Node Model EG in region bjt of material Silicon
Replacing Node Model EI in region bjt of material Silicon
Replacing Node Model EI:Potential in region bjt of material Silicon
Replacing Node Model EV in region bjt of material Silicon
Replacing Node Model EV:Potential in region bjt of material Silicon
Replacing Node Model ElectronGeneration in region bjt of material Silicon
Replacing Node Model ElectronGeneration:Electrons in region bjt of material Silicon
Replacing Node Model ElectronGeneration:Holes in region bjt of material Silicon
Replacing Edge Model Electrons@n0 in region bjt of material Silicon
Replacing Edge Model Electrons@n1 in region bjt of material Silicon
Replacing Triangle Edge Model Emag in region bjt of material Silicon
Replacing Edge Model Epar_n in region bjt of material Silicon
Replacing Edge Model Epar_n:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Epar_n:Potential@n1 in region bjt of material Silicon
Replacing Edge Model Epar_p in region bjt of material Silicon
Replacing Edge Model Epar_p:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Epar_p:Potential@n1 in region bjt of material Silicon
Replacing Node Model HoleGeneration in region bjt of material Silicon
Replacing Node Model HoleGeneration:Electrons in region bjt of material Silicon
Replacing Node Model HoleGeneration:Holes in region bjt of material Silicon
Replacing Edge Model Holes@n0 in region bjt of material Silicon
Replacing Edge Model Holes@n1 in region bjt of material Silicon
Replacing Node Model IntrinsicCharge in region bjt of material Silicon
Replacing Node Model IntrinsicCharge:Potential in region bjt of material Silicon
Replacing Node Model IntrinsicElectrons in region bjt of material Silicon
Replacing Node Model IntrinsicElectrons:Potential in region bjt of material Silicon
Replacing Node Model IntrinsicHoles in region bjt of material Silicon
Replacing Node Model IntrinsicHoles:Potential in region bjt of material Silicon
Replacing Edge Model Jn in region bjt of material Silicon
Replacing Edge Model Jn:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model Jn:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model Jn:Holes@n0 in region bjt of material Silicon
Replacing Edge Model Jn:Holes@n1 in region bjt of material Silicon
Replacing Edge Model Jn:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Jn:Potential@n1 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Holes@n0 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Holes@n1 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Jn_arora_lf:Potential@n1 in region bjt of material Silicon
Replacing Triangle Edge Model Jn_x in region bjt of material Silicon
Replacing Triangle Edge Model Jn_y in region bjt of material Silicon
Replacing Triangle Edge Model Jnmag in region bjt of material Silicon
Replacing Edge Model Jp in region bjt of material Silicon
Replacing Edge Model Jp:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model Jp:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model Jp:Holes@n0 in region bjt of material Silicon
Replacing Edge Model Jp:Holes@n1 in region bjt of material Silicon
Replacing Edge Model Jp:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Jp:Potential@n1 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Holes@n0 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Holes@n1 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Potential@n0 in region bjt of material Silicon
Replacing Edge Model Jp_arora_lf:Potential@n1 in region bjt of material Silicon
Replacing Triangle Edge Model Jp_x in region bjt of material Silicon
Replacing Triangle Edge Model Jp_y in region bjt of material Silicon
Replacing Triangle Edge Model Jpmag in region bjt of material Silicon
Replacing Node Model LogNetDoping in region bjt of material Silicon
Replacing Node Model NC in region bjt of material Silicon
Replacing Node Model NCharge in region bjt of material Silicon
Replacing Node Model NCharge:Electrons in region bjt of material Silicon
Replacing Node Model NIE in region bjt of material Silicon
Replacing Node Model NTOT in region bjt of material Silicon
Replacing Node Model NV in region bjt of material Silicon
Replacing Node Model NetDoping in region bjt of material Silicon
Replacing Node Model PCharge in region bjt of material Silicon
Replacing Node Model PCharge:Holes in region bjt of material Silicon
Replacing Edge Model Potential@n0 in region bjt of material Silicon
Replacing Edge Model Potential@n1 in region bjt of material Silicon
Replacing Node Model PotentialIntrinsicCharge in region bjt of material Silicon
Replacing Node Model PotentialIntrinsicCharge:Potential in region bjt of material Silicon
Replacing Node Model PotentialNodeCharge in region bjt of material Silicon
Replacing Node Model PotentialNodeCharge:Electrons in region bjt of material Silicon
Replacing Node Model PotentialNodeCharge:Holes in region bjt of material Silicon
Replacing Node Model Tn in region bjt of material Silicon
Replacing Node Model USRH in region bjt of material Silicon
Replacing Node Model USRH:Electrons in region bjt of material Silicon
Replacing Node Model USRH:Holes in region bjt of material Silicon
Replacing Node Model V_t in region bjt of material Silicon
Replacing Edge Model V_t_edge in region bjt of material Silicon
Replacing Node Model basenodeelectrons in region bjt of material Silicon
Replacing Node Model basenodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model basenodeholes in region bjt of material Silicon
Replacing Node Model basenodeholes:Holes in region bjt of material Silicon
Replacing Node Model basenodemodel in region bjt of material Silicon
Replacing Node Model basenodemodel:Potential in region bjt of material Silicon
Replacing Edge Model beta_n in region bjt of material Silicon
Replacing Edge Model beta_p in region bjt of material Silicon
Replacing Node Model collectornodeelectrons in region bjt of material Silicon
Replacing Node Model collectornodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model collectornodeholes in region bjt of material Silicon
Replacing Node Model collectornodeholes:Holes in region bjt of material Silicon
Replacing Node Model collectornodemodel in region bjt of material Silicon
Replacing Node Model collectornodemodel:Potential in region bjt of material Silicon
Replacing Node Model contactcharge_node in region bjt of material Silicon
Replacing Node Model emitternodeelectrons in region bjt of material Silicon
Replacing Node Model emitternodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model emitternodeholes in region bjt of material Silicon
Replacing Node Model emitternodeholes:Holes in region bjt of material Silicon
Replacing Node Model emitternodemodel in region bjt of material Silicon
Replacing Node Model emitternodemodel:Potential in region bjt of material Silicon
Replacing Edge Model mu_arora_n_lf in region bjt of material Silicon
Replacing Node Model mu_arora_n_node in region bjt of material Silicon
Replacing Edge Model mu_arora_p_lf in region bjt of material Silicon
Replacing Node Model mu_arora_p_node in region bjt of material Silicon
Replacing Edge Model mu_n in region bjt of material Silicon
Replacing Edge Model mu_n:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model mu_n:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model mu_n:Holes@n0 in region bjt of material Silicon
Replacing Edge Model mu_n:Holes@n1 in region bjt of material Silicon
Replacing Edge Model mu_n:Potential@n0 in region bjt of material Silicon
Replacing Edge Model mu_n:Potential@n1 in region bjt of material Silicon
Replacing Edge Model mu_p in region bjt of material Silicon
Replacing Edge Model mu_p:Electrons@n0 in region bjt of material Silicon
Replacing Edge Model mu_p:Electrons@n1 in region bjt of material Silicon
Replacing Edge Model mu_p:Holes@n0 in region bjt of material Silicon
Replacing Edge Model mu_p:Holes@n1 in region bjt of material Silicon
Replacing Edge Model mu_p:Potential@n0 in region bjt of material Silicon
Replacing Edge Model mu_p:Potential@n1 in region bjt of material Silicon
Replacing Edge Model vsat_n in region bjt of material Silicon
Replacing Edge Model vsat_p in region bjt of material Silicon
Replacing Node Model basenodemodel in region bjt of material Silicon
Replacing Node Model basenodemodel:Potential in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: base, Equation: PotentialEquation
Replacing Node Model basenodeelectrons in region bjt of material Silicon
Replacing Node Model basenodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model basenodeholes in region bjt of material Silicon
Replacing Node Model basenodeholes:Holes in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: base, Equation: ElectronContinuityEquation
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: base, Equation: HoleContinuityEquation
Replacing Node Model emitternodemodel in region bjt of material Silicon
Replacing Node Model emitternodemodel:Potential in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: emitter, Equation: PotentialEquation
Replacing Node Model emitternodeelectrons in region bjt of material Silicon
Replacing Node Model emitternodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model emitternodeholes in region bjt of material Silicon
Replacing Node Model emitternodeholes:Holes in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: emitter, Equation: ElectronContinuityEquation
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: emitter, Equation: HoleContinuityEquation
Replacing Node Model collectornodemodel in region bjt of material Silicon
Replacing Node Model collectornodemodel:Potential in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: collector, Equation: PotentialEquation
Replacing Node Model collectornodeelectrons in region bjt of material Silicon
Replacing Node Model collectornodeelectrons:Electrons in region bjt of material Silicon
Replacing Node Model collectornodeholes in region bjt of material Silicon
Replacing Node Model collectornodeholes:Holes in region bjt of material Silicon
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: collector, Equation: ElectronContinuityEquation
Warning: Replacing Contact Equation with Contact Equation of the same name.
Contact: collector, Equation: HoleContinuityEquation
creating solution dcop with 6 nodes
number of equations 44325
creating solution dcop_prev with 6 nodes
Iteration: 0
Device: “bjt” RelError: 2.79806e-01 AbsError: 1.13014e+07
Region: “bjt” RelError: 2.79806e-01 AbsError: 1.13014e+07
Equation: “ElectronContinuityEquation” RelError: 1.24746e-02 AbsError: 8.73575e+06
Equation: “HoleContinuityEquation” RelError: 2.67332e-01 AbsError: 2.56566e+06
Equation: “PotentialEquation” RelError: 2.43219e-08 AbsError: 5.39830e-12
Circuit: RelError: 2.83219e-06 AbsError: 2.84249e-16
Iteration: 1
Device: “bjt” RelError: 6.41556e-02 AbsError: 8.28313e+03
Region: “bjt” RelError: 6.41556e-02 AbsError: 8.28313e+03
Equation: “ElectronContinuityEquation” RelError: 5.46127e-02 AbsError: 8.17775e+03
Equation: “HoleContinuityEquation” RelError: 9.54290e-03 AbsError: 1.05385e+02
Equation: “PotentialEquation” RelError: 1.84650e-12 AbsError: 2.28823e-16
Circuit: RelError: 2.84794e-06 AbsError: 2.88152e-16
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 4.93139e-02 AbsError: 8.27283e+03
Region: “bjt” RelError: 4.93139e-02 AbsError: 8.27283e+03
Equation: “ElectronContinuityEquation” RelError: 4.92892e-02 AbsError: 8.17775e+03
Equation: “HoleContinuityEquation” RelError: 2.47177e-05 AbsError: 9.50865e+01
Equation: “PotentialEquation” RelError: 4.34556e-12 AbsError: 1.38444e-16
Circuit: RelError: 2.33512e-08 AbsError: 3.05173e-18
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 1.04642e+01 AbsError: 7.61305e+17
Region: “bjt” RelError: 1.04642e+01 AbsError: 7.61305e+17
Equation: “ElectronContinuityEquation” RelError: 9.58639e-01 AbsError: 4.38268e+17
Equation: “HoleContinuityEquation” RelError: 9.64430e-01 AbsError: 3.23037e+17
Equation: “PotentialEquation” RelError: 8.54112e+00 AbsError: 6.00000e-01
Circuit: RelError: 1.00000e+00 AbsError: 6.00000e-01
Iteration: 1
Device: “bjt” RelError: 2.73680e+00 AbsError: 6.37560e+17
Region: “bjt” RelError: 2.73680e+00 AbsError: 6.37560e+17
Equation: “ElectronContinuityEquation” RelError: 1.00000e+00 AbsError: 3.19413e+17
Equation: “HoleContinuityEquation” RelError: 1.00000e+00 AbsError: 3.18147e+17
Equation: “PotentialEquation” RelError: 7.36795e-01 AbsError: 1.32042e-01
Circuit: RelError: 1.00000e+00 AbsError: 3.66509e-04
Iteration: 2
Device: “bjt” RelError: 2.13024e+00 AbsError: 1.93793e+17
Region: “bjt” RelError: 2.13024e+00 AbsError: 1.93793e+17
Equation: “ElectronContinuityEquation” RelError: 9.94521e-01 AbsError: 9.28883e+16
Equation: “HoleContinuityEquation” RelError: 9.95892e-01 AbsError: 1.00904e+17
Equation: “PotentialEquation” RelError: 1.39829e-01 AbsError: 5.00367e-02
Circuit: RelError: 4.59849e+00 AbsError: 4.95548e-04
Iteration: 3
Device: “bjt” RelError: 1.38090e+00 AbsError: 7.28651e+15
Region: “bjt” RelError: 1.38090e+00 AbsError: 7.28651e+15
Equation: “ElectronContinuityEquation” RelError: 4.95951e-01 AbsError: 3.14751e+15
Equation: “HoleContinuityEquation” RelError: 8.48025e-01 AbsError: 4.13900e+15
Equation: “PotentialEquation” RelError: 3.69213e-02 AbsError: 1.26512e-02
Circuit: RelError: 1.23915e+01 AbsError: 1.25466e-04
Iteration: 4
Device: “bjt” RelError: 1.41460e-01 AbsError: 2.07027e+14
Region: “bjt” RelError: 1.41460e-01 AbsError: 2.07027e+14
Equation: “ElectronContinuityEquation” RelError: 7.35700e-02 AbsError: 1.44573e+14
Equation: “HoleContinuityEquation” RelError: 6.48800e-02 AbsError: 6.24539e+13
Equation: “PotentialEquation” RelError: 3.01010e-03 AbsError: 1.04471e-03
Circuit: RelError: 7.93485e+02 AbsError: 1.95994e-05
Iteration: 5
Device: “bjt” RelError: 1.55914e-03 AbsError: 2.12314e+12
Region: “bjt” RelError: 1.55914e-03 AbsError: 2.12314e+12
Equation: “ElectronContinuityEquation” RelError: 7.27669e-04 AbsError: 1.50960e+12
Equation: “HoleContinuityEquation” RelError: 8.08492e-04 AbsError: 6.13541e+11
Equation: “PotentialEquation” RelError: 2.29807e-05 AbsError: 7.93527e-06
Circuit: RelError: 1.15915e+00 AbsError: 1.70627e-07
Iteration: 6
Device: “bjt” RelError: 1.90023e-04 AbsError: 2.07918e+10
Region: “bjt” RelError: 1.90023e-04 AbsError: 2.07918e+10
Equation: “ElectronContinuityEquation” RelError: 1.87659e-04 AbsError: 1.10948e+10
Equation: “HoleContinuityEquation” RelError: 2.28603e-06 AbsError: 9.69698e+09
Equation: “PotentialEquation” RelError: 7.73106e-08 AbsError: 2.15189e-08
Circuit: RelError: 9.39505e-05 AbsError: 2.23722e-09
Iteration: 7
Device: “bjt” RelError: 1.95439e-08 AbsError: 3.93944e+07
Region: “bjt” RelError: 1.95439e-08 AbsError: 3.93944e+07
Equation: “ElectronContinuityEquation” RelError: 1.73213e-08 AbsError: 3.01255e+07
Equation: “HoleContinuityEquation” RelError: 2.14447e-09 AbsError: 9.26888e+06
Equation: “PotentialEquation” RelError: 7.81739e-11 AbsError: 2.43572e-11
Circuit: RelError: 1.36799e-07 AbsError: 7.93968e-14
Iteration: 8
Device: “bjt” RelError: 7.62995e-11 AbsError: 2.32146e+04
Region: “bjt” RelError: 7.62995e-11 AbsError: 2.32146e+04
Equation: “ElectronContinuityEquation” RelError: 5.74589e-11 AbsError: 2.02093e+04
Equation: “HoleContinuityEquation” RelError: 1.87928e-11 AbsError: 3.00526e+03
Equation: “PotentialEquation” RelError: 4.78202e-14 AbsError: 1.35556e-14
Circuit: RelError: 9.42600e-11 AbsError: 1.25978e-17
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 1.41734e-11 AbsError: 8.25129e+03
Region: “bjt” RelError: 1.41734e-11 AbsError: 8.25129e+03
Equation: “ElectronContinuityEquation” RelError: 7.61290e-12 AbsError: 8.17103e+03
Equation: “HoleContinuityEquation” RelError: 6.56013e-12 AbsError: 8.02604e+01
Equation: “PotentialEquation” RelError: 4.15251e-16 AbsError: 7.65172e-17
Circuit: RelError: 3.48785e-11 AbsError: 5.14189e-18
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 3.05222e-11 AbsError: 8.31617e+03
Region: “bjt” RelError: 3.05222e-11 AbsError: 8.31617e+03
Equation: “ElectronContinuityEquation” RelError: 8.84944e-12 AbsError: 8.17103e+03
Equation: “HoleContinuityEquation” RelError: 2.16717e-11 AbsError: 1.45136e+02
Equation: “PotentialEquation” RelError: 1.12779e-15 AbsError: 1.53100e-16
Circuit: RelError: 1.73927e-08 AbsError: 1.76881e-15
Charge Relative Error 1.63468e-02
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 1.72092e-11 AbsError: 9.56686e+03
Region: “bjt” RelError: 1.72092e-11 AbsError: 9.56686e+03
Equation: “ElectronContinuityEquation” RelError: 8.10094e-12 AbsError: 8.17103e+03
Equation: “HoleContinuityEquation” RelError: 9.09306e-12 AbsError: 1.39583e+03
Equation: “PotentialEquation” RelError: 1.52031e-14 AbsError: 2.03374e-15
Circuit: RelError: 3.64276e-09 AbsError: 5.79067e-16
Charge Relative Error 1.63350e-02
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 2.68240e-10 AbsError: 1.18980e+05
Region: “bjt” RelError: 2.68240e-10 AbsError: 1.18980e+05
Equation: “ElectronContinuityEquation” RelError: 5.26740e-12 AbsError: 6.44312e+04
Equation: “HoleContinuityEquation” RelError: 2.62367e-10 AbsError: 5.45484e+04
Equation: “PotentialEquation” RelError: 6.05485e-13 AbsError: 8.10606e-14
Circuit: RelError: 1.28537e-07 AbsError: 1.81937e-14
Charge Relative Error 1.63336e-02
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 1.19562e-08 AbsError: 5.32026e+06
Region: “bjt” RelError: 1.19562e-08 AbsError: 5.32026e+06
Equation: “ElectronContinuityEquation” RelError: 1.78224e-10 AbsError: 2.87982e+06
Equation: “HoleContinuityEquation” RelError: 1.17509e-08 AbsError: 2.44044e+06
Equation: “PotentialEquation” RelError: 2.70383e-11 AbsError: 3.61992e-12
Circuit: RelError: 5.72876e-06 AbsError: 8.25063e-13
Iteration: 1
Device: “bjt” RelError: 4.77850e-11 AbsError: 8.38070e+03
Region: “bjt” RelError: 4.77850e-11 AbsError: 8.38070e+03
Equation: “ElectronContinuityEquation” RelError: 2.65621e-11 AbsError: 8.18802e+03
Equation: “HoleContinuityEquation” RelError: 2.12214e-11 AbsError: 1.92677e+02
Equation: “PotentialEquation” RelError: 1.53896e-15 AbsError: 2.05789e-16
Circuit: RelError: 2.62286e-08 AbsError: 2.41003e-15
Charge Relative Error 1.62717e-02
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 5.26411e-07 AbsError: 2.34613e+08
Region: “bjt” RelError: 5.26411e-07 AbsError: 2.34613e+08
Equation: “ElectronContinuityEquation” RelError: 7.87554e-09 AbsError: 1.27022e+08
Equation: “HoleContinuityEquation” RelError: 5.17343e-07 AbsError: 1.07591e+08
Equation: “PotentialEquation” RelError: 1.19261e-09 AbsError: 1.59666e-10
Circuit: RelError: 2.52674e-04 AbsError: 3.63479e-11
Iteration: 1
Device: “bjt” RelError: 6.48790e-11 AbsError: 8.30848e+03
Region: “bjt” RelError: 6.48790e-11 AbsError: 8.30848e+03
Equation: “ElectronContinuityEquation” RelError: 3.27627e-11 AbsError: 8.18600e+03
Equation: “HoleContinuityEquation” RelError: 3.21154e-11 AbsError: 1.22486e+02
Equation: “PotentialEquation” RelError: 8.51962e-16 AbsError: 1.19510e-16
Circuit: RelError: 1.14065e-07 AbsError: 1.35255e-14
Charge Relative Error 1.03878e-02
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 2.33788e-05 AbsError: 1.04140e+10
Region: “bjt” RelError: 2.33788e-05 AbsError: 1.04140e+10
Equation: “ElectronContinuityEquation” RelError: 3.49575e-07 AbsError: 5.63772e+09
Equation: “HoleContinuityEquation” RelError: 2.29763e-05 AbsError: 4.77629e+09
Equation: “PotentialEquation” RelError: 5.29320e-08 AbsError: 7.08653e-09
Circuit: RelError: 1.10922e-02 AbsError: 1.61360e-09
Iteration: 1
Device: “bjt” RelError: 4.30693e-09 AbsError: 2.85163e+05
Region: “bjt” RelError: 4.30693e-09 AbsError: 2.85163e+05
Equation: “ElectronContinuityEquation” RelError: 1.20274e-11 AbsError: 1.55496e+05
Equation: “HoleContinuityEquation” RelError: 4.29344e-09 AbsError: 1.29667e+05
Equation: “PotentialEquation” RelError: 1.46048e-12 AbsError: 1.95510e-13
Circuit: RelError: 2.12032e-07 AbsError: 6.32542e-14
Charge Relative Error 4.12999e-03
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 3.07237e+03 AbsError: 4.67169e+16
Region: “bjt” RelError: 3.07237e+03 AbsError: 4.67169e+16
Equation: “ElectronContinuityEquation” RelError: 2.05051e+03 AbsError: 2.52709e+16
Equation: “HoleContinuityEquation” RelError: 1.02155e+03 AbsError: 2.14459e+16
Equation: “PotentialEquation” RelError: 3.09367e-01 AbsError: 3.17885e-02
Circuit: RelError: 1.05240e+00 AbsError: 7.23519e-03
Iteration: 1
Device: “bjt” RelError: 4.88095e+03 AbsError: 2.34092e+16
Region: “bjt” RelError: 4.88095e+03 AbsError: 2.34092e+16
Equation: “ElectronContinuityEquation” RelError: 1.66681e+02 AbsError: 1.27705e+16
Equation: “HoleContinuityEquation” RelError: 4.71410e+03 AbsError: 1.06387e+16
Equation: “PotentialEquation” RelError: 1.71631e-01 AbsError: 3.73252e-02
Circuit: RelError: 1.11064e+00 AbsError: 6.10344e-02
Iteration: 2
Device: “bjt” RelError: 9.43829e+02 AbsError: 1.69650e+16
Region: “bjt” RelError: 9.43829e+02 AbsError: 1.69650e+16
Equation: “ElectronContinuityEquation” RelError: 4.03964e+00 AbsError: 7.55082e+15
Equation: “HoleContinuityEquation” RelError: 9.39643e+02 AbsError: 9.41418e+15
Equation: “PotentialEquation” RelError: 1.46287e-01 AbsError: 4.37047e-02
Circuit: RelError: 6.57858e-01 AbsError: 2.28269e-02
Iteration: 3
Device: “bjt” RelError: 2.76754e+02 AbsError: 1.65791e+16
Region: “bjt” RelError: 2.76754e+02 AbsError: 1.65791e+16
Equation: “ElectronContinuityEquation” RelError: 4.96005e+00 AbsError: 7.52971e+15
Equation: “HoleContinuityEquation” RelError: 2.71674e+02 AbsError: 9.04938e+15
Equation: “PotentialEquation” RelError: 1.19791e-01 AbsError: 3.77931e-02
Circuit: RelError: 2.70849e-01 AbsError: 8.40170e-03
Iteration: 4
Device: “bjt” RelError: 9.31479e+02 AbsError: 1.10955e+16
Region: “bjt” RelError: 9.31479e+02 AbsError: 1.10955e+16
Equation: “ElectronContinuityEquation” RelError: 1.90627e+00 AbsError: 5.20046e+15
Equation: “HoleContinuityEquation” RelError: 9.29505e+02 AbsError: 5.89504e+15
Equation: “PotentialEquation” RelError: 6.71850e-02 AbsError: 2.47877e-02
Circuit: RelError: 8.13139e-03 AbsError: 9.73656e-05
Iteration: 5
Device: “bjt” RelError: 8.48824e+03 AbsError: 2.64936e+16
Region: “bjt” RelError: 8.48824e+03 AbsError: 2.64936e+16
Equation: “ElectronContinuityEquation” RelError: 1.10308e+00 AbsError: 1.22364e+16
Equation: “HoleContinuityEquation” RelError: 8.48697e+03 AbsError: 1.42572e+16
Equation: “PotentialEquation” RelError: 1.67272e-01 AbsError: 6.69905e-02
Circuit: RelError: 3.30267e-01 AbsError: 1.23252e-02
Iteration: 6
Device: “bjt” RelError: 8.21682e+03 AbsError: 6.93979e+17
Region: “bjt” RelError: 8.21682e+03 AbsError: 6.93979e+17
Equation: “ElectronContinuityEquation” RelError: 7.78065e+00 AbsError: 2.76293e+17
Equation: “HoleContinuityEquation” RelError: 8.20464e+03 AbsError: 4.17686e+17
Equation: “PotentialEquation” RelError: 4.40276e+00 AbsError: 1.82916e+00
Circuit: RelError: 2.74308e+00 AbsError: 1.24827e-02
Iteration: 7
Device: “bjt” RelError: 9.86477e+04 AbsError: 1.40023e+18
Region: “bjt” RelError: 9.86477e+04 AbsError: 1.40023e+18
Equation: “ElectronContinuityEquation” RelError: 6.00597e+04 AbsError: 5.85128e+17
Equation: “HoleContinuityEquation” RelError: 3.08589e+03 AbsError: 8.15103e+17
Equation: “PotentialEquation” RelError: 3.55022e+04 AbsError: 3.29294e+00
Circuit: RelError: 1.40239e+00 AbsError: 9.90147e-02
Iteration: 8
Device: “bjt” RelError: 8.17687e+04 AbsError: 2.17717e+18
Region: “bjt” RelError: 8.17687e+04 AbsError: 2.17717e+18
Equation: “ElectronContinuityEquation” RelError: 5.57840e+03 AbsError: 9.52540e+17
Equation: “HoleContinuityEquation” RelError: 2.96513e+03 AbsError: 1.22463e+18
Equation: “PotentialEquation” RelError: 7.32252e+04 AbsError: 5.63977e+00
Circuit: RelError: 1.69114e+00 AbsError: 3.60203e-01
Iteration: 9
Device: “bjt” RelError: 7.63745e+03 AbsError: 6.19079e+20
Region: “bjt” RelError: 7.63745e+03 AbsError: 6.19079e+20
Equation: “ElectronContinuityEquation” RelError: 7.26125e+03 AbsError: 2.82408e+20
Equation: “HoleContinuityEquation” RelError: 8.34334e+01 AbsError: 3.36671e+20
Equation: “PotentialEquation” RelError: 2.92767e+02 AbsError: 5.30890e+02
Circuit: RelError: 1.00130e+00 AbsError: 2.10281e+02
while evaluating node model IntrinsicElectrons on Device: bjt on Region: bjt
There was a Overflow floating point exception while evaluating exp((Potential * pow(V_t,(-1))))
while evaluating node model IntrinsicElectrons on Device: bjt on Region: bjt
There was a Overflow floating point exception while evaluating exp((Potential * pow(V_t,(-1))))
while evaluating model IntrinsicElectrons: expression (NIE * exp((Potential * pow(V_t,(-1))))) evaluates to invalid

2

Hi @AC_hm

I think there is a bug, but I am not sure where it is. It may be somewhere in the time-dependent models for the bjt example. Changing “transient_bdf1” to “dc” in the loop causes every thing to run, but unfortunately it is not time dependent. This example is very old, and may take a while to debug.

3

This script works by using the models packaged in devsim, instead of the bjt example. It also loads the mesh from gmsh and applies the doping profile. Loading from the format made it confusing because it was loading the previous physics along with the structure data.

Please let me know if this is useful to you:

from devsim import (
    circuit_alter,
    circuit_element,
    equation,
    get_contact_list,
    get_equation_command,
    set_node_values,
    set_parameter,
    solve,
    )
from devsim.python_packages.model_create import CreateSolution
from devsim.python_packages.simple_physics import (
    CreateSiliconDriftDiffusion,
    CreateSiliconDriftDiffusionAtContact,
    CreateSiliconPotentialOnly,
    CreateSiliconPotentialOnlyContact,
    GetContactBiasName,
    SetSiliconParameters,
    )
import read_gmsh
import netdoping

device="bjt"
region="bjt"

# read in gmsh mesh
# we are avoiding the devsim mesh format, since it loads new physics models
# which may hae a bug
read_gmsh.run("bjt.msh", device, region, "Silicon", ("base", "collector", "emitter"))

# set net doping
netdoping.run(device, region)


####
#### Initial DC solution
####
# here we are avoiding new physical models and using the ones from devsim
CreateSolution(device, region, "Potential")
# start with temperature as a model and not a parameter
set_parameter(device=device, name="T", value=300)
SetSiliconParameters(device, region, 300)
#set_parameter(device=device, name="taun", value=1e-10)
#set_parameter(device=device, name="taup", value=1e-10)



CreateSiliconPotentialOnly(device, region)

fixed_contacts = []
circuit_contacts = ['base', 'collector', 'emitter']

for c in fixed_contacts:
    CreateSiliconPotentialOnlyContact(device, region, c)
    set_parameter(device=device, region=region, name=GetContactBiasName(c), value=0.0)
for c in circuit_contacts:
    circuit_element(name="V%s" % c[0], n1=GetContactBiasName(c), n2="0", value=0.0)
    CreateSiliconPotentialOnlyContact(device, region, c, is_circuit=True)

solve(type="dc", absolute_error=1.0, relative_error=1e-9, maximum_iterations=40)

###
### Drift diffusion and circuit
###
CreateSolution(device, region, "Electrons")
CreateSolution(device, region, "Holes")
set_node_values(device=device, region=region, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=region, name="Holes",     init_from="IntrinsicHoles")

CreateSiliconDriftDiffusion(device, region)

#eqcmd = get_equation_command(device="bjt", region="bjt", name="PotentialEquation")
#eqcmd['variable_update'] = 'default'
#equation(**eqcmd)


for c in fixed_contacts:
    CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=False)
for c in circuit_contacts:
    #set_parameter(device=device, region=region, name=GetContactBiasName(c), value=0.0)
    CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=True)
    # use first initial of each contact name

circuit_alter(name="Vb", value=0.6)

solve(type="dc", absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)
solve(type="transient_dc", absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

# Transient simulation settings
time_step = 1e-11  # Time step
total_time = 1e-9  # Total simulation time
current_time = 0.0
target_vc = 1.0  # Final Vc value
vc_step = target_vc / (total_time / time_step)
new_vc = 0.0
while current_time < total_time:
    current_time += time_step
    new_vc += vc_step
    circuit_alter(name="Vc", value=new_vc)
    
    try:
        solve(type="transient_bdf1",
                     absolute_error=1e6,
                     relative_error=1e-5,
                     maximum_iterations=40,
                     tdelta=time_step,
                     charge_error=1e3)
    except:
        print((new_vc,time_step))
        raise
4

Hi @Juan ,

Thank you for providing the script that utilizes DEVSIM’s built-in models and loads the mesh from Gmsh. This approach is indeed helpful, as it ensures a clean setup by avoiding the inadvertent loading of previous physics along with the structural data.

Regarding the transient simulation, I have observed discrepancies in the current calculations compared to the DC analysis. Specifically, the base current in the transient simulation is significantly larger and negative, and Kirchhoff’s Current Law (KCL) is not satisfied within the specified tolerance. This contrasts with the DC simulation, where KCL is satisfied.

Transient Simulation Results:

  • Base Voltage: 0.6 V
  • Collector Voltage: 2.0 V
  • Emitter Voltage: ~0 V
  • Base Current: -4.177e-03 A
  • Collector Current: 6.076e-04 A
  • Emitter Current: 3.569e-03 A
  • KCL Status: Not satisfied

DC Simulation Results:

  • Base Voltage: 0.6 V
  • Collector Voltage: 2.0 V
  • Emitter Voltage: ~0 V
  • Base Current: 1.594e-07 A
  • Collector Current: 1.580e-04 A
  • Emitter Current: -1.582e-04 A
  • KCL Status: Satisfied

I have ensured that the mesh is correctly loaded from Gmsh and that the doping profile is appropriately applied. The script utilizes DEVSIM’s built-in models, and I have verified that the initial conditions and biasing are consistent between simulations.

  1. Current Calculation Methodology: Is there a difference in how DEVSIM calculates currents in transient versus DC simulations that might explain these discrepancies?

  2. Current densities: In a 2D DEVSIM simulation, is it necessary to adjust the calculated current densities by multiplying them with the respective contact widths to obtain the total current, or does DEVSIM inherently account for the contact width when computing the total current?

  3. is there a solution for transient simulation like the “dcop” or a better way to do it ?

def print_circuit_solution():
    """Prints the circuit node values for DC operating point."""
    for node in devsim.get_circuit_node_list():
        r = devsim.get_circuit_node_value(solution="dcop", node=node)
        print(f"{node}\t{r:1.15e}")

The widths of these contacts are as follows:

Contact Width (cm)
Emitter 5 × 10⁻⁴
Base 5 × 10⁻⁴
Collector 3 × 10⁻³

I have attached the log files and scripts for both simulations for your reference. Your insights into resolving these issues would be greatly appreciated.

Best regards,
AC

Transient_simulation

from devsim import (
circuit_alter,
circuit_element,
get_contact_current,
get_circuit_node_value,
set_node_values,
set_parameter,
solve,
)
from devsim.python_packages.model_create import CreateSolution
from devsim.python_packages.simple_physics import (
CreateSiliconDriftDiffusion,
CreateSiliconDriftDiffusionAtContact,
CreateSiliconPotentialOnly,
CreateSiliconPotentialOnlyContact,
GetContactBiasName,
SetSiliconParameters,
)
import read_gmsh
import netdoping

device = “bjt”
region = “bjt”

Read in Gmsh mesh

read_gmsh.run(“bjt.msh”, device, region, “Silicon”, (“base”, “collector”, “emitter”))

Set net doping

netdoping.run(device, region)

Initial DC solution setup

CreateSolution(device, region, “Potential”)
set_parameter(device=device, name=“T”, value=300)
SetSiliconParameters(device, region, 300)

CreateSiliconPotentialOnly(device, region)

fixed_contacts =
circuit_contacts = [‘base’, ‘collector’, ‘emitter’]

Setup for contacts

for c in fixed_contacts:
CreateSiliconPotentialOnlyContact(device, region, c)
set_parameter(device=device, region=region, name=GetContactBiasName(c), value=0.0)

for c in circuit_contacts:
circuit_element(name=f"V{c[0]}", n1=GetContactBiasName(c), n2=“0”, value=0.0)
CreateSiliconPotentialOnlyContact(device, region, c, is_circuit=True)

solve(type=“dc”, absolute_error=1.0, relative_error=1e-9, maximum_iterations=40)

Drift-diffusion and circuit setup

CreateSolution(device, region, “Electrons”)
CreateSolution(device, region, “Holes”)
set_node_values(device=device, region=region, name=“Electrons”, init_from=“IntrinsicElectrons”)
set_node_values(device=device, region=region, name=“Holes”, init_from=“IntrinsicHoles”)

CreateSiliconDriftDiffusion(device, region)

for c in fixed_contacts:
CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=False)

for c in circuit_contacts:
CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=True)

Function to print contact currents

def print_contact_currents(device, contacts):
for contact in contacts:
electron_current = get_contact_current(device=device, contact=contact, equation=“ElectronContinuityEquation”)
hole_current = get_contact_current(device=device, contact=contact, equation=“HoleContinuityEquation”)
total_current = electron_current + hole_current
print(f"Total current at {contact} contact: {total_current} A")

def print_contact_voltages_dc(device, contacts):
for contact in contacts:
circuit_node = GetContactBiasName(contact)
if circuit_node:
voltage = get_circuit_node_value(solution=‘dcop’, node=circuit_node)
print(f"Voltage at {contact} contact (DC): {voltage} V")
else:
print(f"Error: No circuit node found for contact {contact}")

Function for Printing Contact Voltages during Transient Simulations:

def print_contact_voltages_transient(device, contacts):
for contact in contacts:
circuit_node = GetContactBiasName(contact)
if circuit_node:
voltage = get_circuit_node_value(solution=‘dcop’, node=circuit_node)
print(f"Voltage at {contact} contact (Transient): {voltage} V")
else:
print(f"Error: No circuit node found for contact {contact}")

Function to verify Kirchhoff’s Current Law (KCL) using relative error

def verify_kcl(device, contacts, relative_error_threshold=1e-5):
currents =
for contact in contacts:
electron_current = get_contact_current(device=device, contact=contact, equation=“ElectronContinuityEquation”)
hole_current = get_contact_current(device=device, contact=contact, equation=“HoleContinuityEquation”)
total_current = electron_current + hole_current
currents.append(total_current)

ic, ie, ib = currents
numerator = abs(ic + ie + ib)
denominator = abs(ic) + abs(ie) + abs(ib)
relative_error = numerator / denominator

print(f"Sum of currents (numerator): {numerator} A")
print(f"Total absolute currents (denominator): {denominator} A")
print(f"Relative error: {relative_error}")

if relative_error < relative_error_threshold:
    print("KCL is satisfied within the specified tolerance.")
    return True
else:
    print("KCL is NOT satisfied within the specified tolerance.")
    return False

Apply biases and solve

circuit_alter(name=“Vb”, value=0.6)
solve(type=“dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

circuit_alter(name=“Vc”, value=1.0)
solve(type=“dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)
solve(type=“transient_dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

Print contact currents

print_contact_voltages_dc(device, circuit_contacts)
print_contact_currents(device, circuit_contacts)

Verify Kirchhoff’s Current Law

kcl_satisfied = verify_kcl(device, circuit_contacts)
if not kcl_satisfied:
print(“Investigate discrepancies in current calculations.”)

Transient simulation settings

time_step = 1e-11 # Time step
total_time = 1e-9 # Total simulation time
current_time = 0.0
target_vc = 2.0 # Final Vc value
vc_step = target_vc / (total_time / time_step)
new_vc = 1.0

while current_time < total_time:
current_time += time_step
new_vc += vc_step
if new_vc > target_vc: # Prevent exceeding the target voltage
new_vc = target_vc

circuit_alter(name="Vc", value=new_vc)

try:
    solve(type="transient_bdf1",
          absolute_error=1e6,
          relative_error=1e-5,
          maximum_iterations=40,
          tdelta=time_step,
          charge_error=1e3)
    print_contact_voltages_transient(device, circuit_contacts)
    print_contact_currents(device, circuit_contacts)
    kcl_satisfied = verify_kcl(device, circuit_contacts)
except Exception as e:
    print(f"Error at Vc={new_vc}, time_step={time_step}: {e}")
    raise
DC_simulation

from devsim import (
circuit_alter,
circuit_element,
get_contact_current,
set_node_values,
set_parameter,
solve,
get_circuit_node_value,

)
from devsim.python_packages.model_create import CreateSolution
from devsim.python_packages.simple_physics import (
CreateSiliconDriftDiffusion,
CreateSiliconDriftDiffusionAtContact,
CreateSiliconPotentialOnly,
CreateSiliconPotentialOnlyContact,
GetContactBiasName,
SetSiliconParameters,
)
import read_gmsh
import netdoping

device = “bjt”
region = “bjt”

Read in Gmsh mesh

read_gmsh.run(“bjt.msh”, device, region, “Silicon”, (“base”, “collector”, “emitter”))

Set net doping

netdoping.run(device, region)

Initial DC solution setup

CreateSolution(device, region, “Potential”)
set_parameter(device=device, name=“T”, value=300)
SetSiliconParameters(device, region, 300)

CreateSiliconPotentialOnly(device, region)

fixed_contacts =
circuit_contacts = [‘base’, ‘collector’, ‘emitter’]

Setup for contacts

for c in fixed_contacts:
CreateSiliconPotentialOnlyContact(device, region, c)
set_parameter(device=device, region=region, name=GetContactBiasName(c), value=0.0)

for c in circuit_contacts:
circuit_element(name=f"V{c[0]}", n1=GetContactBiasName(c), n2=“0”, value=0.0)
CreateSiliconPotentialOnlyContact(device, region, c, is_circuit=True)

solve(type=“dc”, absolute_error=1.0, relative_error=1e-9, maximum_iterations=40)

Drift-diffusion and circuit setup

CreateSolution(device, region, “Electrons”)
CreateSolution(device, region, “Holes”)
set_node_values(device=device, region=region, name=“Electrons”, init_from=“IntrinsicElectrons”)
set_node_values(device=device, region=region, name=“Holes”, init_from=“IntrinsicHoles”)

CreateSiliconDriftDiffusion(device, region)

for c in fixed_contacts:
CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=False)

for c in circuit_contacts:
CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=True)

Function to print contact currents

def print_contact_currents(device, contacts):
for contact in contacts:
electron_current = get_contact_current(device=device, contact=contact, equation=“ElectronContinuityEquation”)
hole_current = get_contact_current(device=device, contact=contact, equation=“HoleContinuityEquation”)
total_current = electron_current + hole_current
print(f"Total current at {contact} contact: {total_current} A")

def print_contact_voltages_dc(device, contacts):
for contact in contacts:
circuit_node = GetContactBiasName(contact)
if circuit_node:
voltage = get_circuit_node_value(solution=‘dcop’, node=circuit_node)
print(f"Voltage at {contact} contact (DC): {voltage} V")
else:
print(f"Error: No circuit node found for contact {contact}")

Function to verify Kirchhoff’s Current Law (KCL) using relative error

def verify_kcl(device, contacts, relative_error_threshold=1e-10):
currents =
for contact in contacts:
electron_current = get_contact_current(device=device, contact=contact, equation=“ElectronContinuityEquation”)
hole_current = get_contact_current(device=device, contact=contact, equation=“HoleContinuityEquation”)
total_current = electron_current + hole_current
currents.append(total_current)

ic, ie, ib = currents
numerator = abs(ic + ie + ib)
denominator = abs(ic) + abs(ie) + abs(ib)
relative_error = numerator / denominator

print(f"Sum of currents (numerator): {numerator} A")
print(f"Total absolute currents (denominator): {denominator} A")
print(f"Relative error: {relative_error}")

if relative_error < relative_error_threshold:
    print("KCL is satisfied within the specified tolerance.")
    return True
else:
    print("KCL is NOT satisfied within the specified tolerance.")
    return False

Apply biases and solve

circuit_alter(name=“Vb”, value=0.6)
solve(type=“dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

circuit_alter(name=“Vc”, value=1.0)
solve(type=“dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

Print contact currents

print_contact_currents(device, circuit_contacts)

Verify Kirchhoff’s Current Law

kcl_satisfied = verify_kcl(device, circuit_contacts)
if not kcl_satisfied:
print(“Investigate discrepancies in current calculations.”)

Additional solves

circuit_alter(name=“Vc”, value=2.0)
solve(type=“dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)
solve(type=“transient_dc”, absolute_error=1e6, relative_error=1e-10, maximum_iterations=40)

Print final contact currents

print_contact_voltages_dc(device, circuit_contacts)
print_contact_currents(device, circuit_contacts)
kcl_satisfied = verify_kcl(device, circuit_contacts)

LogFile_Transient

Iteration: 0
Device: “bjt” RelError: 5.69178e-02 AbsError: 5.16664e+14
Region: “bjt” RelError: 5.69178e-02 AbsError: 5.16664e+14
Equation: “ElectronContinuityEquation” RelError: 2.83395e-02 AbsError: 2.77942e+14
Equation: “HoleContinuityEquation” RelError: 2.60882e-02 AbsError: 2.38722e+14
Equation: “PotentialEquation” RelError: 2.49014e-03 AbsError: 3.43039e-04
Circuit: RelError: 4.21915e-02 AbsError: 3.66993e-04
Iteration: 1
Device: “bjt” RelError: 2.69822e-04 AbsError: 3.36630e+11
Region: “bjt” RelError: 2.69822e-04 AbsError: 3.36630e+11
Equation: “ElectronContinuityEquation” RelError: 7.96143e-05 AbsError: 1.66795e+11
Equation: “HoleContinuityEquation” RelError: 1.89682e-04 AbsError: 1.69835e+11
Equation: “PotentialEquation” RelError: 5.25627e-07 AbsError: 7.25311e-08
Circuit: RelError: 7.96281e-05 AbsError: 1.27676e-07
Iteration: 2
Device: “bjt” RelError: 1.84350e-09 AbsError: 1.11674e+05
Region: “bjt” RelError: 1.84350e-09 AbsError: 1.11674e+05
Equation: “ElectronContinuityEquation” RelError: 1.07660e-10 AbsError: 4.10965e+04
Equation: “HoleContinuityEquation” RelError: 1.73576e-09 AbsError: 7.05771e+04
Equation: “PotentialEquation” RelError: 8.43567e-14 AbsError: 2.08496e-14
Circuit: RelError: 3.07594e-09 AbsError: 2.00836e-12
Charge Relative Error 9.60282e-03
Voltage at base contact (Transient): 0.6 V
Voltage at collector contact (Transient): 2.0 V
Voltage at emitter contact (Transient): -2.8597602890957234e-33 V
Total current at base contact: -0.004523244987649031 A
Total current at collector contact: 0.0006450019873405258 A
Total current at emitter contact: 0.0038775338339720232 A
Sum of currents (numerator): 7.091663364821275e-07 A
Total absolute currents (denominator): 0.009045780808961582 A
Relative error: 7.83974707611268e-05
KCL is NOT satisfied within the specified tolerance.
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 5.57415e-02 AbsError: 4.96560e+14
Region: “bjt” RelError: 5.57415e-02 AbsError: 4.96560e+14
Equation: “ElectronContinuityEquation” RelError: 2.75304e-02 AbsError: 2.67149e+14
Equation: “HoleContinuityEquation” RelError: 2.58114e-02 AbsError: 2.29411e+14
Equation: “PotentialEquation” RelError: 2.39969e-03 AbsError: 3.29788e-04
Circuit: RelError: 4.22790e-02 AbsError: 3.52810e-04
Iteration: 1
Device: “bjt” RelError: 2.48687e-04 AbsError: 3.10835e+11
Region: “bjt” RelError: 2.48687e-04 AbsError: 3.10835e+11
Equation: “ElectronContinuityEquation” RelError: 7.35661e-05 AbsError: 1.54016e+11
Equation: “HoleContinuityEquation” RelError: 1.74635e-04 AbsError: 1.56818e+11
Equation: “PotentialEquation” RelError: 4.85983e-07 AbsError: 6.69001e-08
Circuit: RelError: 7.56371e-05 AbsError: 1.17727e-07
Iteration: 2
Device: “bjt” RelError: 1.56979e-09 AbsError: 9.50880e+04
Region: “bjt” RelError: 1.56979e-09 AbsError: 9.50880e+04
Equation: “ElectronContinuityEquation” RelError: 9.17493e-11 AbsError: 3.49181e+04
Equation: “HoleContinuityEquation” RelError: 1.47797e-09 AbsError: 6.01700e+04
Equation: “PotentialEquation” RelError: 7.13341e-14 AbsError: 1.78783e-14
Circuit: RelError: 1.37934e-09 AbsError: 1.99484e-12
Charge Relative Error 8.32832e-03
Voltage at base contact (Transient): 0.6 V
Voltage at collector contact (Transient): 2.0 V
Voltage at emitter contact (Transient): -3.0591541061011293e-32 V
Total current at base contact: -0.004346801087934734 A
Total current at collector contact: 0.0006259327487453262 A
Total current at emitter contact: 0.0037201866089124004 A
Sum of currents (numerator): 6.817302770074744e-07 A
Total absolute currents (denominator): 0.00869292044559246 A
Relative error: 7.842361853813233e-05
KCL is NOT satisfied within the specified tolerance.
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 5.45489e-02 AbsError: 4.77228e+14
Region: “bjt” RelError: 5.45489e-02 AbsError: 4.77228e+14
Equation: “ElectronContinuityEquation” RelError: 2.67221e-02 AbsError: 2.56770e+14
Equation: “HoleContinuityEquation” RelError: 2.55146e-02 AbsError: 2.20458e+14
Equation: “PotentialEquation” RelError: 2.31226e-03 AbsError: 3.17038e-04
Circuit: RelError: 4.23690e-02 AbsError: 3.39165e-04
Iteration: 1
Device: “bjt” RelError: 2.28996e-04 AbsError: 2.87006e+11
Region: “bjt” RelError: 2.28996e-04 AbsError: 2.87006e+11
Equation: “ElectronContinuityEquation” RelError: 6.79739e-05 AbsError: 1.42212e+11
Equation: “HoleContinuityEquation” RelError: 1.60573e-04 AbsError: 1.44794e+11
Equation: “PotentialEquation” RelError: 4.49301e-07 AbsError: 6.17077e-08
Circuit: RelError: 7.18248e-05 AbsError: 1.08556e-07
Iteration: 2
Device: “bjt” RelError: 3.01938e-09 AbsError: 8.10935e+04
Region: “bjt” RelError: 3.01938e-09 AbsError: 8.10935e+04
Equation: “ElectronContinuityEquation” RelError: 7.84776e-11 AbsError: 2.97338e+04
Equation: “HoleContinuityEquation” RelError: 2.94084e-09 AbsError: 5.13597e+04
Equation: “PotentialEquation” RelError: 6.04886e-14 AbsError: 1.51782e-14
Circuit: RelError: 8.49244e-10 AbsError: 5.40534e-13
Charge Relative Error 9.14359e-03
Voltage at base contact (Transient): 0.6 V
Voltage at collector contact (Transient): 2.0 V
Voltage at emitter contact (Transient): 8.20389997066352e-33 V
Total current at base contact: -0.004177181640944649 A
Total current at collector contact: 0.0006076078739598118 A
Total current at emitter contact: 0.003568918398199972 A
Sum of currents (numerator): 6.553687848655058e-07 A
Total absolute currents (denominator): 0.008353707913104434 A
Relative error: 7.845244191952546e-05
KCL is NOT satisfied within the specified tolerance.

LogFile_DC

Iteration: 13
Device: “bjt” RelError: 2.80335e+04 AbsError: 1.87212e+14
Region: “bjt” RelError: 2.80335e+04 AbsError: 1.87212e+14
Equation: “ElectronContinuityEquation” RelError: 2.05169e+00 AbsError: 1.35653e+14
Equation: “HoleContinuityEquation” RelError: 2.80314e+04 AbsError: 5.15582e+13
Equation: “PotentialEquation” RelError: 1.21145e-02 AbsError: 2.82814e-02
Circuit: RelError: 1.28352e-02 AbsError: 1.22428e-08
Iteration: 14
Device: “bjt” RelError: 1.31168e+03 AbsError: 9.44289e+13
Region: “bjt” RelError: 1.31168e+03 AbsError: 9.44289e+13
Equation: “ElectronContinuityEquation” RelError: 1.00000e+00 AbsError: 7.24969e+13
Equation: “HoleContinuityEquation” RelError: 1.31067e+03 AbsError: 2.19321e+13
Equation: “PotentialEquation” RelError: 9.77010e-03 AbsError: 2.30334e-02
Circuit: RelError: 1.09794e-03 AbsError: 2.85018e-09
Iteration: 15
Device: “bjt” RelError: 1.67054e+03 AbsError: 8.69035e+13
Region: “bjt” RelError: 1.67054e+03 AbsError: 8.69035e+13
Equation: “ElectronContinuityEquation” RelError: 2.40664e-01 AbsError: 6.52564e+13
Equation: “HoleContinuityEquation” RelError: 1.67030e+03 AbsError: 2.16471e+13
Equation: “PotentialEquation” RelError: 3.62118e-04 AbsError: 8.25373e-04
Circuit: RelError: 1.47921e-03 AbsError: 2.64726e-09
Iteration: 16
Device: “bjt” RelError: 2.95501e-01 AbsError: 4.92981e+11
Region: “bjt” RelError: 2.95501e-01 AbsError: 4.92981e+11
Equation: “ElectronContinuityEquation” RelError: 1.17988e-03 AbsError: 4.18893e+11
Equation: “HoleContinuityEquation” RelError: 2.94319e-01 AbsError: 7.40879e+10
Equation: “PotentialEquation” RelError: 1.89134e-06 AbsError: 4.31092e-06
Circuit: RelError: 2.40703e-08 AbsError: 2.15010e-12
Iteration: 17
Device: “bjt” RelError: 6.81948e-07 AbsError: 1.64152e+07
Region: “bjt” RelError: 6.81948e-07 AbsError: 1.64152e+07
Equation: “ElectronContinuityEquation” RelError: 3.13281e-08 AbsError: 1.46390e+07
Equation: “HoleContinuityEquation” RelError: 6.50561e-07 AbsError: 1.77625e+06
Equation: “PotentialEquation” RelError: 5.92172e-11 AbsError: 1.34974e-10
Circuit: RelError: 2.75873e-12 AbsError: 2.04696e-17
Iteration: 18
Device: “bjt” RelError: 1.54549e-15 AbsError: 1.88462e-02
Region: “bjt” RelError: 1.54549e-15 AbsError: 1.88462e-02
Equation: “ElectronContinuityEquation” RelError: 1.26283e-15 AbsError: 1.72272e-02
Equation: “HoleContinuityEquation” RelError: 2.82596e-16 AbsError: 1.61901e-03
Equation: “PotentialEquation” RelError: 6.43335e-20 AbsError: 1.46635e-19
Circuit: RelError: 4.89100e-18 AbsError: 1.11173e-21
number of equations 44325
Iteration: 0
Device: “bjt” RelError: 2.73929e-25 AbsError: 7.24246e-15
Region: “bjt” RelError: 2.73929e-25 AbsError: 7.24246e-15
Equation: “ElectronContinuityEquation” RelError: 2.72571e-25 AbsError: 7.10234e-15
Equation: “HoleContinuityEquation” RelError: 1.35821e-27 AbsError: 1.40123e-16
Equation: “PotentialEquation” RelError: 7.52529e-34 AbsError: 1.53433e-33
Circuit: RelError: 4.89100e-18 AbsError: 1.11173e-21
Voltage at base contact (DC): 0.6 V
Voltage at collector contact (DC): 2.0 V
Voltage at emitter contact (DC): 9.9743548256647e-42 V
Total current at base contact: 1.5943467399288255e-07 A
Total current at collector contact: 0.00015805397985680026 A
Total current at emitter contact: -0.00015821341453079314 A
Sum of currents (numerator): 0.0 A
Total absolute currents (denominator): 0.0003164268290615863 A
Relative error: 0.0
KCL is satisfied within the specified tolerance.

5

I appreciate your patience. I will be able to take a look in the next few days

6

Hello @AC_hm

I have been on vacation, but I am now looking at your scripts. It will take me some more time to get you an answer concerning your simulation. I think there may be some charge models that are in the bjt example equations, which are not correct and I will update you after I see if removing them makes a difference.

At the end of each step, the current and charge are stored on the contact.
Using get_contact_current does not account for any transient charge. For an ohmic contact to a semiconductor, it should be sufficient. For a MOS capacitor, you can use get_contact_charge to get the charge on the insulator.

To get the net current, including the time derivative of charge, you should use a circuit node. If your voltage source is Vb, then the current through it will be Vb.I, and will account for the net current resulting from electron current and charge from the contact. At the end of each transient step, it will be stored on dcop, which is the same storage as the dc solution.

The flux into the circuit is from the current and charge and goes into this option of the contact_equation command.

    circuit_node : str, optional
       Name of the circuit we integrate the flux into

The models for charge and current integrated into the circuit node are:

    edge_charge_model : str, optional
       Name of the edge model used to determine the charge at this contact
    edge_current_model : str, optional
       Name of the edge model used to determine the current flowing out of this contact
    element_charge_model : str, optional
       Name of the element edge model used to determine the charge at this contact
    element_current_model : str, optional
       Name of the element edge model used to determine the current flowing out of this contact
    node_charge_model : str, optional
       Name of the node model used to determine the charge at this contact
    node_current_model : str, optional
       Name of the node model used to determine the current flowing out of this

For an ohmic contact, I would expect that only edge_current_model would need to be defined. You can use get_contact_equation_command to verify this. I suspect for historical reasons, the BJT example may be defining edge_charge_model and node_charge_model which I now believe is incorrect. I will check this and see about removing it using:

eqcmd = get_contact_equation_command(.....)
eqcmd['edge_charge_model'] = ''
eqcmd['node_charge_model'] = ''
contact_equation(**eqcmd)

DEVSIM accounts for the contact width by integrating the edge current going into the contact by its EdgeCouple model.

I = \sum_i \sum_j J_{i,j} \text{EdgeCouple}_{i,j}

where J_{i,j} is the edge current between a node on the contact and another node inside the bulk of the device.

image
image1270×1127 28.2 KB

dcop is the correct solution to use for both transient and dc simulation.

7

@AC_hm

Please let me know if you are still having a simulation issue.

8

@Juan
Thank you, Juan. The simulation now works as expected, and I’m currently working on other simulations. If any issues arise, I’ll be sure to reach out. I appreciate your support!

1 Like
9

Subject: Issues with Base Gradient and Potential Hole Injection in NPN Transistor Simulation

Hello @Juan,
I’m currently simulating an NPN transistor in forward-active mode using devsim, and I’ve encountered some unexpected behaviors regarding the minority carrier distributions in the base region. I’ve attached a photo that shows the theoretical profile I expect in the base. Below, I provide my analytical results and simulation setup details, along with several questions.

1. Background and Analytical Results (Forward-Active Case)

  • Built-in Voltages:

    • V_bi1 = 0.9172 V
    • V_bi2 = 0.7980 V

  • Depletion Widths:

    • W1 = 0.120 μm
    • W2 = 0.716 μm

  • Junction Edge Positions (abrupt):

    • x_e = 0.00 μm
    • x1 = 2.00 μm
    • x2 = 3.00 μm
    • x_c = 9.00 μm
    • x_n1 = 2.00 μm
    • x_p1 = 2.12 μm
    • x_n2 = 3.07 μm
    • x_p2 = 2.35 μm

  • Minority Carrier Concentrations:

    • p(x_e) = 4.10 cm⁻³
    • p(x_n1) = 1.08×10¹⁴ cm⁻³
    • n(x_p1) = 1.08×10¹⁷ cm⁻³
    • n(x_p2) = 1.92×10⁻⁴⁷ cm⁻³
    • p(x_n2) = 1.92×10⁻⁴⁸ cm⁻³
    • p(x_c) = 4.10×10² cm⁻³

  • Diffusion Lengths:

    • L_pe = 111.47 μm
    • L_nb = 186.94 μm
    • L_pc = 111.47 μm
      Current Values:

  • I_en = 2.642×10⁻¹ A

  • I_ep = 1.076×10⁻⁵ A

  • I_e = 2.642×10⁻¹ A

  • I_cn = 2.642×10⁻¹ A

  • I_cp = 1.374×10⁻¹⁷ A

  • I_c = 2.642×10⁻¹ A

  • I_b = 1.095×10⁻⁵ A

  • V_ce = 3.800 V

  • V1 = 0.800 V

  • V2 = -3.000 V

  • γ = 1.000, B = 1.000, α = 1.000, β = 2.412×10⁴

2. Simulation Results and Setup Details

The simulation yields the following contact currents and bias values:

  • Currents:
    • Vb.I = -0.02497 A
    • Vc.I = -2.95778 A
    • Ve.I = 2.98275 A
  • Bias Voltages:
    • Base bias = 0.0 V
    • Collector bias = 3.0 V
    • Emitter bias = -0.8 V

Code Highlights:

  • Carrier Initialization:
set_node_values(device=device, region=region, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=region, name="Holes", init_from="IntrinsicHoles")

Is this approach sufficient for initializing the carrier values from the potential solution?

  • Parameter Definitions:
    For mobility and carrier lifetimes, I set the following:
set_parameter(device=device, region=region, name="mu_n", value=1350)
set_parameter(device=device, region=region, name="mu_p", value=480)
set_parameter(device=device, region=region, name="taun_emitter", value=1e-5)
set_parameter(device=device, region=region, name="taun_base", value=1e-5)
set_parameter(device=device, region=region, name="taun_collector", value=1e-5)

Waht is the optimal way to define these parameters in devsim for each region diffrently? also the mobilities (“mu_n”, “mu_p”) be defined differently per region?

  • Mesh and Doping Profile:
    The mesh is generated from a .geo file refined using bjt_refine.py. My simulation uses abrupt doping profiles, and I observe similar issues even with a smooth erfc profile in the bjt_example.

3. Specific Questions

  1. Base Gradient:
    Why does the electron concentration gradient in the base not approach the expected nb₀ (as depicted in the attached photo here case D)? Could this be due to an unexpected hole injection from the base into the emitter that does not decay as expected?
  2. Carrier Initialization:
    Is initializing the carriers via set_node_values(..., init_from="IntrinsicElectrons"/"IntrinsicHoles") from the potential solution sufficient, or are additional initialization steps required?
  3. Parameter Setup in devsim:
    What is the recommended approach for setting parameters like mobilities (mu_n, mu_p) and carrier lifetimes in different regions? How should one properly define these (and possibly related parameters like mun/mup) in devsim?
  4. Simulation of Abrupt vs. Smooth Doping Profiles:
    Are there known improvements or best practices for simulating structures with abrupt doping profiles (and even with smooth profiles) to achieve results that better match theoretical expectations?

Attachments:

  • A photo showing the theoretical hole and electron concentration profile in the base.
  • The complete simulation code, including the calculate_and_print_parameters(device, region, case_name, vb, vc, ve) function and related parts of the setup.

I appreciate any insights or suggestions on how to address these issues. Thank you in advance for your help!

bjt_2d.py 
# -------------------------------------------------------------------------
# Import necessary functions from DEVSIM and associated Python packages
# -------------------------------------------------------------------------
from print_module import*

import math
from devsim import (
    create_gmsh_mesh,
    add_gmsh_region,
    add_gmsh_contact,
    finalize_mesh,
    create_device,
    write_devices,
    create_2d_mesh,
    reset_devsim,
    circuit_alter,
    circuit_element,
    get_contact_current,
    set_node_values,
    get_circuit_node_list,
    set_parameter,
    solve,
    get_circuit_node_value,
    node_model,
    get_node_model_values,
    print_node_values,
    get_contact_list,
    edge_from_node_model,
    edge_model,
    get_parameter,
    get_edge_model_values
)
from devsim.python_packages.model_create import CreateSolution, CreateNodeModel
from devsim.python_packages.simple_physics import (
    CreateSiliconDriftDiffusion,
    CreateSiliconDriftDiffusionAtContact,
    CreateSiliconPotentialOnly,
    CreateSiliconPotentialOnlyContact,
    GetContactBiasName,
    SetSiliconParameters,
    PrintCurrents
)

# -------------------------------------------------------------------------

# Define device and region names
device = "bjt"
region = "bjt"
filename = "bjt.msh"
region_material = "Silicon"

# Set physical parameters assuming constants 
n_i = 6.405e9 #1.0e10  # #/cm^3 
q = 1.6e-19  # coul
k = 1.3806503e-23  # J/K
eps_0 = 8.85e-14  # F/cm^2
eps_si = 11.1
T = 300

reset_devsim()

def initialize_simulation():
    
    # -------------------------------------------------------------------------
    # Load mesh from Gmsh file
    # -------------------------------------------------------------------------
    create_gmsh_mesh(mesh=device, file=filename)

    # -------------------------------------------------------------------------
    # Define region:
    # The Physical Surface defined in the Gmsh file ("bjt") is assigned to region
    # "bjt" and material "Silicon" is set.
    # -------------------------------------------------------------------------
    add_gmsh_region(mesh=device, gmsh_name=region, region=region, material=region_material)

    # -------------------------------------------------------------------------
    # Define contacts
    # -------------------------------------------------------------------------
    add_gmsh_contact(mesh=device, gmsh_name="base", region=region, material="metal", name="base")
    add_gmsh_contact(mesh=device, gmsh_name="collector", region=region, material="metal", name="collector")
    add_gmsh_contact(mesh=device, gmsh_name="emitter", region=region, material="metal", name="emitter")

    # -------------------------------------------------------------------------
    # Finalize mesh, create device and save
    # -------------------------------------------------------------------------
    finalize_mesh(mesh=device)
    create_device(mesh=device, device=device)

    set_parameter(device=device, region=region,name='emitter_doping' ,value= 1e19)
    set_parameter(device=device, region=region,name='base_doping' ,value= 1e16)
    set_parameter(device=device, region=region,name='collector_doping' ,value= 1e17)
    set_parameter(device=device, region=region,name='x1' ,value= 2e-4)
    set_parameter(device=device, region=region,name='x2' ,value= 3e-4)
    set_parameter(device=device, region=region,name='xcl' ,value= 9e-4)

    # 2.5 Define doping profiles using step functions
    node_model(device=device, region=region, name="Acceptors",equation="base_doping * step(x - x1) * step(x2 - x)")

    node_model(device=device, region=region, name="Donors",equation="emitter_doping * step(x1 - x) + collector_doping * step(x - x2)")

    # NetDoping: Difference between Donors and Acceptors
    node_model(device=device, region=region, name="NetDoping" ,equation= "Donors - Acceptors")
    node_model(device=device, region=region, name="LogNetDoping", equation="asinh(Donors-Acceptors/2)/log(10)")

    write_devices(file="netdoping.msh", type="devsim")
    # 2.6 Create solution for electric potential

    CreateSolution(device, region, "Potential")
    CreateNodeModel(device, region, "IntrinsicElectrons", "ni")
    CreateNodeModel(device, region, "IntrinsicHoles", "ni")

    set_parameter(device=device, region=region, name="Permittivity", value=eps_si * eps_0)
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)
    set_parameter(device=device, region=region, name="n_i", value=n_i)
    set_parameter(device=device, region=region, name="T", value=T)
    set_parameter(device=device, region=region, name="kT", value=k * T)
    set_parameter(device=device, region=region, name="V_t", value=k * T / q)
    set_parameter(device=device, region=region, name="mu_n", value=1350)
    set_parameter(device=device, region=region, name="mu_p", value=480)

    # default SRH parameters
    set_parameter(device=device, region=region, name="n1", value=n_i)
    set_parameter(device=device, region=region, name="p1", value=n_i)

    set_parameter(device=device, region=region, name="taun_emitter", value=1e-5)
    set_parameter(device=device, region=region, name="taun_base", value=1e-5)
    set_parameter(device=device, region=region, name="taun_collector", value=1e-5)

    # Create Node Model for tau_n as function of x
    node_model(device=device, region=region, name="taun",
           equation="taun_emitter * step(x1 - x) + taun_base * step(x - x1) * step(x2 - x) + taun_collector * step(x - x2)")

    set_parameter(device=device, region=region, name="taup", value=1e-5)

    # 2.7 Create physical model considering only potential
    CreateSiliconPotentialOnly(device, region)

    # -------------------------------------------------------------------------
    # 3. Connect contacts to circuit (Potential-only contacts)
    # -------------------------------------------------------------------------
    circuit_contacts = get_contact_list(device=device)
    for c in circuit_contacts:
        # Create circuit element. Name is based on first letter of contact.
        circuit_element(name=f"V{c[0]}", n1=GetContactBiasName(c), n2="0", value=0.0)
        # Create "Potential Only" contact model.
        CreateSiliconPotentialOnlyContact(device, region, c, is_circuit=True)

    # -------------------------------------------------------------------------
    # 4. First DC solution (Potential-only solution)
    # -------------------------------------------------------------------------
    solve(type="dc", absolute_error=1e-3, relative_error=1e-12, maximum_iterations=100)

    # -------------------------------------------------------------------------
    # 5. Drift-Diffusion simulation: Charge carrier densities and electric fields
    # -------------------------------------------------------------------------
    CreateSolution(device, region, "Electrons")
    CreateSolution(device, region, "Holes")

    set_node_values(device=device, region=region, name="Electrons", init_from="IntrinsicElectrons")
    set_node_values(device=device, region=region, name="Holes", init_from="IntrinsicHoles")

    edge_from_node_model(device=device, region=region, node_model="Potential")

    edge_model(device=device, region=region, name="EField",
                  equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")

    edge_model(device=device, region=region, name="EField:Potential@n0",
                 equation="EdgeInverseLength")  

    edge_model(device=device, region=region, name="EField:Potential@n1",
                  equation="-EdgeInverseLength")  

    CreateSiliconDriftDiffusion(device, region)
    for c in circuit_contacts:
        CreateSiliconDriftDiffusionAtContact(device, region, c, is_circuit=True)


def ramp_voltage(contact, initial_voltage, final_voltage, steps=10, abs_err=1e6, rel_err=1e-3, max_iter=40):
    """Ramps the voltage for a given contact from initial_voltage to final_voltage in specified steps."""
    dv = (final_voltage - initial_voltage) / steps
    voltage = initial_voltage
    for step in range(steps):
        voltage += dv
        circuit_alter(name=contact, value=voltage)
        solve(type="dc", absolute_error=abs_err, relative_error=rel_err, maximum_iterations=max_iter)
        print(f"Ramping {contact}: Step {step+1}/{steps}, voltage = {voltage:.3f} V")


def run_simulation_case(case_name, vb, vc, ve=0.0):
    """Runs simulation for specific operating condition"""
    ramp_voltage("Ve", initial_voltage=0.0, final_voltage=ve, steps=4)
    ramp_voltage("Vb", initial_voltage=0.0, final_voltage=vb, steps=4)
    ramp_voltage("Vc", initial_voltage=0.0, final_voltage=vc, steps=4)
    
    # File output
    write_devices(
        file=f"bjt_{case_name}.tec", 
        type="tecplot"
    )
    
    # Result logging
    print_circuit_solution()
    calculate_and_print_parameters(device, region, case_name, vb, vc, ve)


def print_circuit_solution():
    """Prints values of circuit nodes for DC operating point."""
    for node in get_circuit_node_list():
        r = get_circuit_node_value(solution="dcop", node=node)
        print(f"{node}\t{r:1.15e}")

if __name__ == "__main__":
    # Initialization
    initialize_simulation()
    solve(type="dc", absolute_error=1e6, relative_error=1e-3, maximum_iterations=100)
    
    # Case 1: Equilibrium state
    run_simulation_case(
        case_name="equilibrium",
        vb=0,
        vc=0,
        ve=0
    )
    
    # Case 2: Active forward operation
    run_simulation_case(
        case_name="forward_active",
        vb=0,
        vc=3,
        ve=-0.8
    )
def calculate 
from print_module import *
import math
from devsim import (
    create_gmsh_mesh,
    add_gmsh_region,
    add_gmsh_contact,
    finalize_mesh,
    create_device,
    write_devices,
    create_2d_mesh,
    reset_devsim,
    circuit_alter,
    circuit_element,
    get_contact_current,
    set_node_values,
    get_circuit_node_list,
    set_parameter,
    solve,
    get_circuit_node_value,
    node_model,
    get_node_model_values,
    print_node_values,
    get_contact_list,
    edge_from_node_model,
    edge_model,
    get_parameter,
    get_edge_model_values
)
from devsim.python_packages.model_create import CreateSolution, CreateNodeModel
from devsim.python_packages.simple_physics import (
    CreateSiliconDriftDiffusion,
    CreateSiliconDriftDiffusionAtContact,
    CreateSiliconPotentialOnly,
    CreateSiliconPotentialOnlyContact,
    GetContactBiasName,
    SetSiliconParameters,
    PrintCurrents
)

def calculate_and_print_parameters(device, region, case_name, vb, vc, ve):

    # --- Kontaktfläche (in DEVSIM in cm²) ---
    # Hier wird z. B. angenommen, dass Kontaktbreite und Tiefe bereits in cm vorliegen.
    contact_width_cm = 0.1e-4   # kontaktbreite 0.1 um
    depth_cm = 1.0           # Tiefe in cm
    A_cb_cm2 = contact_width_cm * depth_cm  # in cm²
    # Für Stromberechnung: Umrechnung in m² (1 cm² = 1e-4 m²)
    A_cb_SI = A_cb_cm2 * 1e-4
    set_parameter(device=device, region=region, name="A_cb", value=A_cb_SI)
    
    # --- Abfrage der Parameter (alle Längen in cm, Doping in cm⁻³ etc.) ---
    q  = get_parameter(device=device, region=region, name="ElectronCharge")  # [C]
    k  = 1.38064852e-23  # [J/K]
    T  = get_parameter(device=device, region=region, name="T")            # [K]
    V_t = get_parameter(device=device, region=region, name="V_t")          # [V]
    n_i = get_parameter(device=device, region=region, name="n_i")          # in cm⁻³

    # Dopingwerte (in cm⁻³)
    N_de = get_parameter(device=device, region=region, name="emitter_doping")
    N_ab = get_parameter(device=device, region=region, name="base_doping")
    N_dc = get_parameter(device=device, region=region, name="collector_doping")
    
    # Positionsparameter (in cm)
    x1_cm  = get_parameter(device=device, region=region, name="x1")
    x2_cm  = get_parameter(device=device, region=region, name="x2")
    xcl_cm = get_parameter(device=device, region=region, name="xcl")
    # Für die Simulation: xe = 0 cm, xc = xcl_cm
    x_e_cm = 0.0
    x_c_cm = xcl_cm

    # Beweglichkeiten (in cm²/(V·s))
    mu_n = get_parameter(device=device, region=region, name="mu_n")
    mu_p = get_parameter(device=device, region=region, name="mu_p")
    
    # Permittivität (in F/cm)
    eps = get_parameter(device=device, region=region, name="Permittivity")
    
    # Zeitkonstanten (in s)
    tau_pe = get_parameter(device=device, region=region, name="taup")
    tau_nb = get_parameter(device=device, region=region, name="taun_base")
    tau_pc = get_parameter(device=device, region=region, name="taup")
    
    # --- Diffusionskoeffizienten in cm²/s (D = mu * V_t) ---
    D_nb = mu_n * V_t
    D_pe = mu_p * V_t
    D_pc = D_pe

    # Spannungdifferenzen
    V1_val = vb - ve
    V2_val = vb - vc

    # --- (1) Built-in-Spannungen ---
    V_bi1 = (k * T / q) * math.log((N_de * N_ab) / (n_i ** 2))
    V_bi2 = (k * T / q) * math.log((N_dc * N_ab) / (n_i ** 2))

    # --- (2) Depletionsbreiten (ergibt cm) ---
    W1_cm = math.sqrt((2 * eps * (N_de + N_ab) * (V_bi1 - V1_val)) / (q * N_de * N_ab))
    W2_cm = math.sqrt((2 * eps * (N_dc + N_ab) * (V_bi2 - V2_val)) / (q * N_dc * N_ab))

    # --- (3) Sperrkantenpositionen (abrupt) in cm ---
    xn1_cm = x1_cm - ((N_ab * W1_cm) / (N_ab + N_de))
    xp1_cm = x1_cm + ((N_de * W1_cm) / (N_ab + N_de))
    xn2_cm = x2_cm + ((N_ab * W2_cm) / (N_ab + N_dc))
    xp2_cm = x2_cm - ((N_dc * W2_cm) / (N_ab + N_dc))
    
    # --- (4) Minoritätskonzentrationen (in cm⁻³) ---
    pxe   = (n_i ** 2) / N_de
    pxn1  = N_ab * math.exp(-q * (V_bi1 - V1_val) / (k * T))
    nxp1  = N_de * math.exp(-q * (V_bi1 - V1_val) / (k * T))
    nxp2  = N_dc * math.exp(-q * (V_bi2 - V2_val) / (k * T))
    pxn2  = N_ab * math.exp(-q * (V_bi2 - V2_val) / (k * T))
    pxc   = (n_i ** 2) / N_dc

    # --- (5) Diffusionslängen (ergibt cm) ---
    L_pe_cm = math.sqrt(D_pe * tau_pe)
    L_nb_cm = math.sqrt(D_nb * tau_nb)
    L_pc_cm = math.sqrt(D_pc * tau_pc)

    # --- (6) Koeffizienten der Minoritätsverteilungen ---
    A_e = (pxe * math.exp(x_e_cm / L_pe_cm) - pxn1 * math.exp(xn1_cm / L_pe_cm)) / (math.exp(2 * x_e_cm / L_pe_cm) - math.exp(2 * xn1_cm / L_pe_cm))*1e6
    B_e = (pxe * math.exp(-x_e_cm / L_pe_cm) - pxn1 * math.exp(-xn1_cm / L_pe_cm)) / (math.exp(-2 * x_e_cm / L_pe_cm) - math.exp(-2 * xn1_cm / L_pe_cm))*1e6
    A_b = (nxp1 * math.exp(xp1_cm / L_nb_cm) - nxp2 * math.exp(xp2_cm / L_nb_cm)) / (math.exp(2 * xp1_cm / L_nb_cm) - math.exp(2 * xp2_cm / L_nb_cm))*1e6
    B_b = (nxp1 * math.exp(-xp1_cm / L_nb_cm) - nxp2 * math.exp(-xp2_cm / L_nb_cm)) / (math.exp(-2 * xp1_cm / L_nb_cm) - math.exp(-2 * xp2_cm / L_nb_cm))*1e6
    A_c = (pxn2 * math.exp(xn2_cm / L_pc_cm) - pxc * math.exp(x_c_cm / L_pc_cm)) / (math.exp(2 * xn2_cm / L_pc_cm) - math.exp(2 * x_c_cm / L_pc_cm))*1e6
    B_c = (pxn2 * math.exp(-xn2_cm / L_pc_cm) - pxc * math.exp(-x_c_cm / L_pc_cm)) / (math.exp(-2 * xn2_cm / L_pc_cm) - math.exp(-2 * x_c_cm / L_pc_cm))*1e6

    # --- (7) Ableitungen der Minoritätskonzentrationen (ergibt 1/cm) ---
    dp_e_dx_at_xn1 = (A_e / L_pe_cm) * math.exp(xn1_cm / L_pe_cm) - (B_e / L_pe_cm) * math.exp(-xn1_cm / L_pe_cm)
    dn_b_dx_at_xp1 = (A_b / L_nb_cm) * math.exp(xp1_cm / L_nb_cm) - (B_b / L_nb_cm) * math.exp(-xp1_cm / L_nb_cm)
    dn_b_dx_at_xp2 = (A_b / L_nb_cm) * math.exp(xp2_cm / L_nb_cm) - (B_b / L_nb_cm) * math.exp(-xp2_cm / L_nb_cm)
    dp_c_dx_at_xn2 = (A_c / L_pc_cm) * math.exp(xn2_cm / L_pc_cm) - (B_c / L_pc_cm) * math.exp(-xn2_cm / L_pc_cm)

    # --- (8) Ströme ---
    # Für die Stromberechnung: Diffusionskoeffizienten in m²/s, Ableitungen in 1/m
    D_nb_SI = D_nb * 1e-4
    D_pe_SI = D_pe * 1e-4
    D_pc_SI = D_pc * 1e-4
    dn_b_dx_at_xp1_SI = dn_b_dx_at_xp1 * 100
    dn_b_dx_at_xp2_SI = dn_b_dx_at_xp2 * 100
    dp_e_dx_at_xn1_SI = dp_e_dx_at_xn1 * 100
    dp_c_dx_at_xn2_SI = dp_c_dx_at_xn2 * 100

    I_en = -A_cb_SI * q * D_nb_SI * dn_b_dx_at_xp1_SI
    I_ep =  A_cb_SI * q * D_pe_SI * dp_e_dx_at_xn1_SI
    I_e  = I_en + I_ep
    I_cn = -A_cb_SI * q * D_nb_SI * dn_b_dx_at_xp2_SI
    I_cp =  A_cb_SI * q * D_pc_SI * dp_c_dx_at_xn2_SI
    I_c  = I_cn + I_cp
    I_b  = I_e - I_c
    V_ce = V1_val - V2_val

    gamma = I_en / (I_en + I_ep) if (I_en + I_ep) != 0 else float('nan')
    B_val = I_c / I_en if I_en != 0 else float('nan')
    alpha = I_c / I_e if I_e != 0 else float('nan')
    beta  = alpha / (1 - alpha) if (1 - alpha) != 0 else float('nan')

    # --- Ausgabe der Ergebnisse ---
    # Für geometrische Größen: Umrechnung in μm (1 cm = 1e4 μm)
    print(f"\n——— {case_name} ———")
    print("(1) Built-in-Spannungen:")
    print(f"    V_bi1 = (k*T/q)*ln((N_de*N_ab)/(n_i²))")
    print(f"          = ({k:.3e}*{T:.1f}/{q:.3e})*ln(({N_de:.3e}*{N_ab:.3e})/({n_i:.3e}²)) = {V_bi1:.4f} V")
    print(f"    V_bi2 = (k*T/q)*ln((N_dc*N_ab)/(n_i²))")
    print(f"          = ({k:.3e}*{T:.1f}/{q:.3e})*ln(({N_dc:.3e}*{N_ab:.3e})/({n_i:.3e}²)) = {V_bi2:.4f} V\n")

    print("(2) Depletionsbreiten:")
    print(f"    W1 = {W1_cm*1e4:.3f} μm")
    print(f"    W2 = {W2_cm*1e4:.3f} μm\n")

    print("(3) Sperrkantenpositionen (abrupt):")
    print(f"    x_e  = {x_e_cm*1e4:.2f} μm")
    print(f"    x1   = {x1_cm*1e4:.2f} μm")
    print(f"    x2   = {x2_cm*1e4:.2f} μm")
    print(f"    x_c  = {x_c_cm*1e4:.2f} μm")
    print(f"    x_n1 = x1 - ((N_ab*W1)/(N_ab+N_de))")
    print(f"         = {x1_cm*1e4:.2f} μm - (({N_ab:.3e}*{W1_cm*1e4:.2f} μm)/({(N_ab+N_de):.3e})) = {xn1_cm*1e4:.2f} μm")
    print(f"    x_p1 = x1 + ((N_de*W1)/(N_ab+N_de))")
    print(f"         = {x1_cm*1e4:.2f} μm + (({N_de:.3e}*{W1_cm*1e4:.2f} μm)/({(N_ab+N_de):.3e})) = {xp1_cm*1e4:.2f} μm")
    print(f"    x_n2 = x2 + ((N_ab*W2)/(N_ab+N_dc))")
    print(f"         = {x2_cm*1e4:.2f} μm + (({N_ab:.3e}*{W2_cm*1e4:.2f} μm)/({(N_ab+N_dc):.3e})) = {xn2_cm*1e4:.2f} μm")
    print(f"    x_p2 = x2 - ((N_dc*W2)/(N_ab+N_dc))")
    print(f"         = {x2_cm*1e4:.2f} μm - (({N_dc:.3e}*{W2_cm*1e4:.2f} μm)/({(N_ab+N_dc):.3e})) = {xp2_cm*1e4:.2f} μm\n")

    print("(4) Minoritätskonzentrationen (an den Kontaktknoten):")
    print(f"    p(x_e)   = n_i²/N_de = ({n_i:.3e}²)/{N_de:.3e} = {pxe:.2e} cm⁻³")
    print(f"    p(x_n1)  = {N_ab:.3e}*exp(-q*(V_bi1 - (vb-ve))/(k*T)) = {pxn1:.2e} cm⁻³")
    print(f"    n(x_p1)  = {N_de:.3e}*exp(-q*(V_bi1 - (vb-ve))/(k*T)) = {nxp1:.2e} cm⁻³")
    print(f"    n(x_p2)  = {N_dc:.3e}*exp(-q*(V_bi2 - (vb-vc))/(k*T)) = {nxp2:.2e} cm⁻³")
    print(f"    p(x_n2)  = {N_ab:.3e}*exp(-q*(V_bi2 - (vb-vc))/(k*T)) = {pxn2:.2e} cm⁻³")
    print(f"    p(x_c)   = n_i²/N_dc = ({n_i:.3e}²)/{N_dc:.3e} = {pxc:.2e} cm⁻³\n")

    print("(5) Diffusionslängen:")
    print(f"    L_pe = {L_pe_cm*1e4:.2f} μm")
    print(f"    L_nb = {L_nb_cm*1e4:.2f} μm")
    print(f"    L_pc = {L_pc_cm*1e4:.2f} μm\n")

    print("(6) Koeffizienten der Minoritätsverteilungen:")
    print("    Im Emitter:")
    print(f"      A_e = {A_e:.3e}")
    print(f"      B_e = {B_e:.3e}\n")
    print("    In der Base:")
    print(f"      A_b = {A_b:.3e}")
    print(f"      B_b = {B_b:.3e}\n")
    print("    Im Collector:")
    print(f"      A_c = {A_c:.3e}")
    print(f"      B_c = {B_c:.3e}\n")

    print("(7) Ableitungen der Minoritätskonzentrationen (in 1/cm):")
    print(f"    dp_e/dx (bei x=x_n1) = {dp_e_dx_at_xn1:.3e} 1/cm")
    print(f"    dn_b/dx (bei x=x_p1) = {dn_b_dx_at_xp1:.3e} 1/cm")
    print(f"    dn_b/dx (bei x=x_p2) = {dn_b_dx_at_xp2:.3e} 1/cm")
    print(f"    dp_c/dx (bei x=x_n2) = {dp_c_dx_at_xn2:.3e} 1/cm\n")

    print("(8) Ströme:")
    print(f"    I_en = {I_en:.3e} A")
    print(f"    I_ep = {I_ep:.3e} A")
    print(f"    I_e  = {I_e:.3e} A")
    print(f"    I_cn = {I_cn:.3e} A")
    print(f"    I_cp = {I_cp:.3e} A")
    print(f"    I_c  = {I_c:.3e} A")
    print(f"    I_b  = {I_b:.3e} A")
    print(f"    V_ce = {V_ce:.3f} V")
    print(f"    V_1 = {V1_val:.3f} V")
    print(f"    V_2 = {V2_val:.3f} V")
    print(f"    γ = {gamma:.3f}")
    print(f"    B = {B_val:.3f}")
    print(f"    α = {alpha:.3f}")
    print(f"    β = {beta:.3e}")
abrupt_bjt_2d.geo 
/***********************************************
 * 2D-BJT Geometrie – Parametrisierte .geo-Datei
 * =============================================
 * Kontaktbreite modifiziert auf 0.1 µm
 ***********************************************/

/* 1. Mesh-Parameter und allgemeine Einstellungen */
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.Algorithm = 5; // Frontal-Delaunay Algorithmus
Mesh.CharacteristicLengthMax = 1e-5; // Maximale Gittergröße
Mesh.MshFileVersion = 2.2; // Gmsh-Dateiversion

/* 2. Skalierung und charakteristische Längen */
sf = 1.0e-4; // Skalierungsfaktor (1 µm = 1e-4 cm)
cl_emitter   = 2.5e-5; // Gittergröße für Emitter
cl_base      = 2.5e-5; // Gittergröße für Basis
cl_collector = 2.5e-5; // Gittergröße für Collector
cl_noncontact = cl_base; // Gittergröße für nicht-kontaktierte Bereiche

/* 3. Geräteabmessungen */
device_width = 12 * sf; // Gesamtweite des Geräts
device_depth = 5 * sf; // Tiefe des Geräts

/* 4. Kontaktsegmente (0.1 µm breit) */
contact_width = 0.1 * sf; // Neue Kontaktbreite (0.1 µm)
x_e_mid = 1 * sf; // Mittelpunkt des Emitters
x_e_contact_start = x_e_mid - contact_width / 2; // Startpunkt des Emitters
x_e_contact_end   = x_e_mid + contact_width / 2; // Endpunkt des Emitters

x_base_mid = 2.5 * sf; // Mittelpunkt der Basis
x_base_contact_start = x_base_mid - contact_width / 2; // Startpunkt der Basis
x_base_contact_end   = x_base_mid + contact_width / 2; // Endpunkt der Basis

x_c_mid = 7 * sf; // Mittelpunkt des Collectors
x_c_contact_start = x_c_mid - contact_width / 2; // Startpunkt des Collectors
x_c_contact_end   = x_c_mid + contact_width / 2; // Endpunkt des Collectors

x0 = 0; // Startpunkt des Geräts
x_end = device_width; // Endpunkt des Geräts

/* 5. Punkte */
Point(1) = {x0, 0, 0, cl_noncontact}; // linker unterer Eckpunkt
Point(2) = {x_e_contact_start, 0, 0, cl_emitter}; // Emitter-Kontaktstart
Point(3) = {x_e_contact_end, 0, 0, cl_emitter}; // Emitter-Kontaktende
Point(4) = {x_base_contact_start, 0, 0, cl_base}; // Basis-Kontaktstart
Point(5) = {x_base_contact_end, 0, 0, cl_base}; // Basis-Kontaktende
Point(6) = {x_c_contact_start, 0, 0, cl_collector}; // Collector-Kontaktstart
Point(7) = {x_c_contact_end, 0, 0, cl_collector}; // Collector-Kontaktende
Point(8) = {x_end, 0, 0, cl_noncontact}; // rechter unterer Eckpunkt
Point(9) = {x_end, device_depth, 0, cl_noncontact}; // rechter oberer Eckpunkt
Point(10) = {x0, device_depth, 0, cl_noncontact}; // linker oberer Eckpunkt

/* 6. Linien */
Line(1) = {1, 2}; // von linker Ecke zum Emitter-Kontaktstart
Line(2) = {2, 3}; // Emitter-Kontakt
Line(3) = {3, 4}; // von Emitter-Kontaktende zum Basis-Kontaktstart
Line(4) = {4, 5}; // Basis-Kontakt
Line(5) = {5, 6}; // von Basis-Kontaktende zum Collector-Kontaktstart
Line(6) = {6, 7}; // Collector-Kontakt
Line(7) = {7, 8}; // von Collector-Kontaktende zum rechten Rand
Line(8) = {8, 9}; // rechte Seite des Geräts
Line(9) = {9, 10}; // obere Seite des Geräts
Line(10) = {10, 1}; // linke Seite des Geräts

/* 7. Fläche */
Line Loop(20) = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; // geschlossener Linienloop
Plane Surface(21) = {20}; // definierte Fläche

/* 8. Physikalische Gruppen */
Physical Line("emitter") = {2}; // Emitter-Kontakt
Physical Line("base") = {4}; // Basis-Kontakt
Physical Line("collector") = {6}; // Collector-Kontakt
Physical Surface("bjt") = {21}; // gesamte BJT-Fläche

Capture d'écran 2025-03-02 051811
Capture d'écran 2025-03-02 0518112929×1155 300 KB

10

Below is my response to your specific questions. I apologize if I am not able to run your simulations at this time. Please consider my answers as being from the perspective of a qualitative discussion and not specific quantitative reasoning. Online, it appears that there are multiple resources available with a thoughtful analytical theory of BJT operation.

What I am seeing in your figure is a comparison of a logarithmic plot being compared to a theoretical result, which I am interpreting as being linear. Perhaps a linear y-axis plot may be more revealing. I also don’t believe the carrier distribution in the base to be a straight line, but more of a reasonable approximation.

The simulation result is always governed by a discretized solution, whereas an analytic result may be based on the depletion approximation. There may also be other factors, concerning the quality of the mesh or other simulation inputs.

This is a reasonable initial guess. Without it, it is likely that the drift diffusion simulation will not converge. After used for the first bias point, this guess is no longer needed for the rest of the simulation.

The set_parameter command can be used to set a parameter globally, per device, or per region. If you set mu_n, mu_p on a region, then the global or per-device value is ignored.

If your structure has only 1 region, then you could use an edge_model to set the mobility on a per-edge basis. You can then either use set_edge_values, or even have a position dependent equation.

This introduces some complexity in ensuring that the mu_n, mu_p values are treated as being models, and not parameters.

If you wish for abrupt doping profiles, it is best to either:

  1. Create a structure with a different region for each dopant density. You would then need to connect the regions with interface equations.
  2. Create a single-region device which has adequate mesh refinement to recreate an abrupt junction.