The heat generation calculation in th MOSFET

Hello everyone, I am trying to calculate the Joule heat only and its influence to temperature distribution in the MOSFET. For now, there are two problems bothering me. The equation is :

First, the heat conduction equation has same form with Potential equation. So I modified the Poisson equation to get the heat conduction equation. For simplicity, I use the uniform heat source to test the code. But the results shows that the heat source didn’t work. My modified code is shown below：

The gmsh_mos2d.py :

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create

device = "mos2d"
silicon_regions=("gate", "bulk")
oxide_regions=("oxide",)
regions = ("gate", "bulk", "oxide")
interfaces = ("bulk_oxide", "gate_oxide")

for i in regions:
CreateSolution(device, i, "Potential")

for i in silicon_regions:
SetSiliconParameters(device, i, 300)
CreateSiliconPotentialOnly(device, i)

for i in oxide_regions:
SetOxideParameters(device, i, 300)
CreateOxidePotentialOnly(device, i, "log_damp")

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
tmp = get_region_list(device=device, contact=i)
r = tmp[0]
print("%s %s" % (r, i))
CreateSiliconPotentialOnlyContact(device, r, i)
# set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
set_parameter(device=device, name="drain_bias", value=0.0)
set_parameter(device=device, name="source_bias", value=0.0)
set_parameter(device=device, name="gate_bias", value=0.0)
set_parameter(device=device, name="body_bias", value=0.0)
for i in interfaces:
CreateSiliconOxideInterface(device, i)
# for d in get_device_list():
#  for gn in get_parameter_list():
#    print("{0} {1}").format(gn, get_parameter(device=d, name=gn))
#  for gn in get_parameter_list(device=d):
#    print("{0} {1} {2}").format(d, gn, get_parameter(device=d, nasme=gn))
#  for r in get_region_list(device=d):
#    for gn in get_parameter_list(device=d, region=r):
#      print("{0} {1} {2} {3}").format(d, r, gn, get_parameter(device=d, region=r, name=gn))
# write_devices(file="foo.msh", type="devsim")
#

CreateSolution(device, "bulk", "Temperature")
CreateTemperaturefiled(device,"bulk")
set_node_values(device=device, region="bulk", name="Temperature", init_from="init_temperature")
CreateThermalcontact(device,"bulk","body")
CreateThermalcontact(device,"bulk","drain")
CreateThermalcontact(device,"bulk","source")
set_parameter(device=device, name="drain_biasT", value=300)
set_parameter(device=device, name="source_biasT", value=300)
set_parameter(device=device, name="body_biasT", value=300)

solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
# solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
# write_devices -file gmsh_mos2d_potentialonly.flps -type floops
# write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

for i in silicon_regions:
CreateSolution(device, i, "Electrons")
CreateSolution(device, i, "Holes")
set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
CreateSiliconDriftDiffusion(device, i, "mu_n", "mu_p")

for c in contacts:
tmp = get_region_list(device=device, contact=c)
r = tmp[0]
CreateSiliconDriftDiffusionAtContact(device, r, c)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)

for r in silicon_regions:
node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data

for r in silicon_regions:
element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)

#
# printAllCurrents(device)
# rampbias(device, "drain",  1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)    ##current unit is A/cm
# rampbias(device, "gate", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

``````

The simple_physics.py:

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
#T      = 300 # K
eps_si = 11.1
eps_ox = 3.9
# TODO: make this temperature dependent
n_i    = 1.0e10 # #/cm^3
# constant in our approximation
mu_n   = 400
mu_p   = 200
Sithermalconducitivity = 19100  ##W/K*cm
heatsource = 10000     ##w/cm^2
init_T = 300
boundarytemperature = 250
def GetContactBiasName(contact):
return "{0}_bias".format(contact)

def GetContactTemperatureName(contact):
return "{0}_biasT".format(contact)

def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)

def GetContactTemperatureNodeModelName(contact):
return "{0}Tnodemodel".format(contact)

def PrintCurrents(device, contact):
'''
print out contact currents
'''
contact_bias_name = GetContactBiasName(contact)
electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
total_current   = electron_current + hole_current
voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
conteq="Permittivity*ElectricField"
contact_bias_name  = GetContactBiasName(contact)
contact_model_name = GetContactNodeModelName(contact)
eq = "Potential - {0}".format(contact_bias_name)
CreateContactNodeModel(device, contact, contact_model_name, eq)
CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

#TODO: make everyone use dfield
if not InEdgeModelList(device, region, contactcharge_edge):
CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_charge_model= contactcharge_edge)

#####
##### Constants
#####
def SetOxideParameters(device, region, T):
'''
Sets physical parameters
'''
set_parameter(device=device, region=region, name="Permittivity",   value=eps_ox * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)

def SetSiliconParameters(device, region, T):
'''
Sets physical parameters assuming constants
'''
#### TODO: make T a free parameter and T dependent parameters as models
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=mu_n)
set_parameter(device=device, region=region, name="mu_p",           value=mu_p)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Sithermalconducitivity)
set_parameter(device=device, region=region, name="heatsource",     value=heatsource)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)
set_parameter(device=device, region=region, name="boundaryTemperature",     value=boundarytemperature)
#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", value=1e-5)
set_parameter(device=device, region=region, name="taup", value=1e-5)

def CreateSiliconPotentialOnly(device, region):
'''
Creates the physical models for a Silicon region
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
elec_i = "n_i*exp(Potential/V_t)"
hole_i = "n_i^2/IntrinsicElectrons"
charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
pcharge_i = "-ElectronCharge * IntrinsicCharge"

# require NetDoping
for i in (
("IntrinsicElectrons", elec_i),
("IntrinsicHoles", hole_i),
("IntrinsicCharge", charge_i),
("PotentialIntrinsicCharge", pcharge_i)
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Potential")

### TODO: Edge Average Model
for i in (
("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
("PotentialEdgeFlux", "Permittivity * ElectricField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Potential")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
'''
Creates the potential equation at the contact
if is_circuit is true, than use node given by GetContactBiasName
'''
# Means of determining contact charge
# Same for all contacts
if not InNodeModelList(device, region, "contactcharge_node"):
CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contactcharge_edge"):
CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
-V_t*log({1}/n_i), \
V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)

contact_model_name = GetContactNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")

def CreateSRH(device, region):
USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
Gn = "-ElectronCharge * USRH"
Gp = "+ElectronCharge * USRH"
CreateNodeModel(device, region, "USRH", USRH)
CreateNodeModel(device, region, "ElectronGeneration", Gn)
CreateNodeModel(device, region, "HoleGeneration", Gp)
for i in ("Electrons", "Holes"):
CreateNodeModelDerivative(device, region, "USRH", USRH, i)
CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)

def CreateECE(device, region, mu_n):
CreateElectronCurrent(device, region, mu_n)

NCharge = "-ElectronCharge * Electrons"
CreateNodeModel(device, region, "NCharge", NCharge)
CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model = "NCharge",
edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region, mu_p):
CreateHoleCurrent(device, region, mu_p)
PCharge = "ElectronCharge * Holes"
CreateNodeModel(device, region, "PCharge", PCharge)
CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
time_node_model = "PCharge",
edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
CreateNodeModel(device, region, "PotentialNodeCharge", pne)
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
time_node_model="", variable_update="log_damp")

def CreateSiliconDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p"):
CreatePE(device, region)
CreateBernoulli(device, region)
CreateSRH(device, region)
CreateECE(device, region, mu_n)
CreateHCE(device, region, mu_p)

def CreateSiliconDriftDiffusionAtContact(device, region, contact, is_circuit=False):
'''
Restrict electrons and holes to their equilibrium values
Integrates current into circuit
'''
contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
contact_electrons_name = "{0}nodeelectrons".format(contact)
contact_holes_name = "{0}nodeholes".format(contact)

CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
#TODO: The simplification of the ifelse statement is time consuming
#   CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

#TODO: keyword args
if is_circuit:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

else:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent")

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent")

def CreateOxidePotentialOnly(device, region, update_type="default"):
'''
Create electric field model in oxide
Creates Potential solution variable if not available
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")

efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
# this needs to remove derivatives w.r.t. independents
CreateEdgeModel(device, region, "ElectricField", efield)
CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
dfield="Permittivity*ElectricField"
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="PotentialEdgeFlux", variable_update=update_type)

def CreateSiliconOxideInterface(device, interface):
'''
continuous potential at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiliconSiliconInterface(device, interface):
'''
Enforces potential, electron, and hole continuity across the interface
'''
CreateSiliconOxideInterface(device, interface)
ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
hname = CreateContinuousInterfaceModel(device, interface, "Holes")
interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")

def CreateTemperaturefiled(device, region):
'''
Creates a thermal Equation in a region
'''
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")

for i in (
("Source", "heatsource"),
("init_temperature", "initialTemperature"),
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Temperature")

### TODO: Edge Average Model
for i in (
("Temperaturefield",     "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("ThermalEdgeflux", "thermalconducitivity * Temperaturefield")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")

equation(device=device, region=region, name="ThermalEquation", variable_name="Temperature",
node_model="Source", edge_model="ThermalEdgeflux", variable_update="log_damp")

def CreateThermalcontact(device, region, contact):   ##constant temperature boundary
# Boundary_name_temperature="%sbiasT" % contact
# format_dict = { "contact" : contact}
# set_parameter(device=device, region=region, name=Boundary_name_temperature, value=300)
# for name, equation in (
#     ("%(contact)snodemodel",  "Temperature-%(contact)sbiasT"),
#     ("%(contact)snodemodel:Potential", "1"),
# ):
#     name_sub = name % format_dict
#     equation_sub = equation % format_dict
#     contact_node_model(device=device, contact=contact, name=name_sub, equation=equation_sub)

# contact_equation(device=device, contact=contact, name="ThermalEquation",
#                         node_model="%snodemodel" % contact)

if not InNodeModelList(device, region, "contacttemperature_node"):
CreateNodeModel(device, region, "contacttemperature_node", "-heatsource")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contacttemperature_edge"):
CreateEdgeModel(device, region, "contacttemperature_edge", "thermalconducitivity*Temperaturefield")
CreateEdgeModelDerivatives(device, region, "contacttemperature_edge", "thermalconducitivity*Temperaturefield", "Temperature")
contact_model = "Temperature -{0}".format(GetContactTemperatureName(contact))

contact_model_name = GetContactTemperatureNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Temperature"), "1")
contact_equation(device=device, contact=contact, name="ThermalEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contacttemperature_node", edge_charge_model="contacttemperature_edge",
node_current_model="", edge_current_model="")

``````

I didn’t change the simple_dd.py and gmsh_mos2d_create.py. And the temperature boundary I set just for test, which didn’t have actual physical background. The code works fine, but it turns out that the heat source isn’t doing the trick. Could anyone tell me where I am wrong?

The second problem is that the current and electric field both belongs to the edge model. But to calculate the heat generation, it seem to I need transfer them to a node model. Has anyone done something similar ？Please give me some advice.

It should be possible to integrate the Joule heating term on the `EdgeNodevolume` shown in:
https://devsim.net/models.html#id4

this is by specifying the `edge_volume_model` option in the `equation` command:
https://devsim.net/CommandReference.html#devsim.equation

The model specified in this option is integrated into the subvolume of the edge into both of the nodes. You would only need to put `(Jn + Jp) * E` into your equation. You would also need to provide the derivatives of these models with respect to the solution variables. For example, this would be:

``````H:Potential@n0 = (Jn:Potential@n0 + Jp:Potential@n0) * E + (Jn + Jp) * E:Potential@n0
H:Potential@n1
H:Electrons@n0
H:Electrons@n1
H:Holes@n0
H:Holes@n1
H:Temperature@n0
H:Temperature@n1
``````

where I show how the chain rule can be applied for to get the equation for the first term.

This approach is used for the density gradient formulation showed in this document:
and is referred to as the ‘EV’ component in chapter 4.

Hi @Juan
Thank you for your suggestion. I will try it later. But now, I can’t even calculate the result for a simple uniform heat source. Could you please take a moment to help me check where is the wrong with equation assembly ? The results shows “Convergence failure!”. Actully, I don’t think the equation is solved at all. I’m following the Poisson equation, but I haven’t checked out what exactly is the problem.I changed the code as follow:
gmsh_mos2d.py

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create

device = "mos2d"
silicon_regions=("gate", "bulk")
oxide_regions=("oxide",)
regions = ("gate", "bulk", "oxide")
interfaces = ("bulk_oxide", "gate_oxide")

for i in regions:
CreateSolution(device, i, "Potential")
#T field
CreateSolution(device,i,"Temperature")

for i in oxide_regions:
SetOxideParameters(device, i, 300)

for i in silicon_regions:
SetSiliconParameters(device, i, 300)
CreateSiliconPotentialOnly(device, i)
#T field
CreateTemperaturefield(device,i)

for i in oxide_regions:
# SetOxideParameters(device, i, 300)
CreateOxidePotentialOnly(device, i, "log_damp")

#T field
CreateTemperaturefield(device,i)

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
tmp = get_region_list(device=device, contact=i)
r = tmp[0]
print("%s %s" % (r, i))
CreateSiliconPotentialOnlyContact(device, r, i)

#T field
CreateTemperatureboundary(device,r,i)

# set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
set_parameter(device=device, name="drain_bias", value=1.0)
set_parameter(device=device, name="source_bias", value=0.0)
set_parameter(device=device, name="gate_bias", value=1.0)
set_parameter(device=device, name="body_bias", value=0.0)

#T field
set_parameter(device=device, name="drain_biasT", value=300)
set_parameter(device=device, name="source_biasT", value=300)
set_parameter(device=device, name="gate_biasT", value=300)
set_parameter(device=device, name="body_biasT", value=300)

for i in interfaces:
CreateSiliconOxideInterface(device, i)

#T field
CreateTemperatureinterface(device,i)

# for d in get_device_list():
#  for gn in get_parameter_list():
#    print("{0} {1}").format(gn, get_parameter(device=d, name=gn))
#  for gn in get_parameter_list(device=d):
#    print("{0} {1} {2}").format(d, gn, get_parameter(device=d, nasme=gn))
#  for r in get_region_list(device=d):
#    for gn in get_parameter_list(device=d, region=r):
#      print("{0} {1} {2} {3}").format(d, r, gn, get_parameter(device=d, region=r, name=gn))
# write_devices(file="foo.msh", type="devsim")
#

solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=100)
# solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
# write_devices -file gmsh_mos2d_potentialonly.flps -type floops
# write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

for i in silicon_regions:
CreateSolution(device, i, "Electrons")
CreateSolution(device, i, "Holes")
set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
CreateSiliconDriftDiffusion(device, i, "mu_n", "mu_p")

for c in contacts:
tmp = get_region_list(device=device, contact=c)
r = tmp[0]
CreateSiliconDriftDiffusionAtContact(device, r, c)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=100)

for r in silicon_regions:
node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data

for r in silicon_regions:
element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)

#
# printAllCurrents(device)
# rampbias(device, "drain",  1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)    ##current unit is A/cm
# rampbias(device, "gate", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

``````

simple_physics.py

``````#Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
#T      = 300 # K
eps_si = 11.1
eps_ox = 3.9
# TODO: make this temperature dependent
n_i    = 1.0e10 # #/cm^3
# constant in our approximation
mu_n   = 400
mu_p   = 200
Sithermalconducitivity = 19100  ##W/K*cm
Oxidethermalconducitivity = 14000
heatsource = 10000     ##w/cm^2
init_T = 200
boundarytemperature = 250
def GetContactBiasName(contact):
return "{0}_bias".format(contact)

def GetContactTemperatureName(contact):
return "{0}_biasT".format(contact)

def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)

def GetContactTemperatureNodeModelName(contact):
return "{0}Tnodemodel".format(contact)

def PrintCurrents(device, contact):
'''
print out contact currents
'''
contact_bias_name = GetContactBiasName(contact)
electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
total_current   = electron_current + hole_current
voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
conteq="Permittivity*ElectricField"
contact_bias_name  = GetContactBiasName(contact)
contact_model_name = GetContactNodeModelName(contact)
eq = "Potential - {0}".format(contact_bias_name)
CreateContactNodeModel(device, contact, contact_model_name, eq)
CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

#TODO: make everyone use dfield
if not InEdgeModelList(device, region, contactcharge_edge):
CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_charge_model= contactcharge_edge)

#####
##### Constants
#####
def SetOxideParameters(device, region, T):
'''
Sets physical parameters
'''
set_parameter(device=device, region=region, name="Permittivity",   value=eps_ox * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)
set_parameter(device=device, region=region, name="heatsource",     value=heatsource)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Oxidethermalconducitivity)

def SetSiliconParameters(device, region, T):
'''
Sets physical parameters assuming constants
'''
#### TODO: make T a free parameter and T dependent parameters as models
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=mu_n)
set_parameter(device=device, region=region, name="mu_p",           value=mu_p)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Sithermalconducitivity)
set_parameter(device=device, region=region, name="heatsource",     value=heatsource)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)
set_parameter(device=device, region=region, name="boundaryTemperature",     value=boundarytemperature)
#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", value=1e-5)
set_parameter(device=device, region=region, name="taup", value=1e-5)

def CreateSiliconPotentialOnly(device, region):
'''
Creates the physical models for a Silicon region
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
elec_i = "n_i*exp(Potential/V_t)"
hole_i = "n_i^2/IntrinsicElectrons"
charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
pcharge_i = "-ElectronCharge * IntrinsicCharge"

# require NetDoping
for i in (
("IntrinsicElectrons", elec_i),
("IntrinsicHoles", hole_i),
("IntrinsicCharge", charge_i),
("PotentialIntrinsicCharge", pcharge_i)
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Potential")

### TODO: Edge Average Model
for i in (
("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
("PotentialEdgeFlux", "Permittivity * ElectricField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Potential")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
'''
Creates the potential equation at the contact
if is_circuit is true, than use node given by GetContactBiasName
'''
# Means of determining contact charge
# Same for all contacts
if not InNodeModelList(device, region, "contactcharge_node"):
CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contactcharge_edge"):
CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
-V_t*log({1}/n_i), \
V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)

contact_model_name = GetContactNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")

def CreateSRH(device, region):
USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
Gn = "-ElectronCharge * USRH"
Gp = "+ElectronCharge * USRH"
CreateNodeModel(device, region, "USRH", USRH)
CreateNodeModel(device, region, "ElectronGeneration", Gn)
CreateNodeModel(device, region, "HoleGeneration", Gp)
for i in ("Electrons", "Holes"):
CreateNodeModelDerivative(device, region, "USRH", USRH, i)
CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)

def CreateECE(device, region, mu_n):
CreateElectronCurrent(device, region, mu_n)

NCharge = "-ElectronCharge * Electrons"
CreateNodeModel(device, region, "NCharge", NCharge)
CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model = "NCharge",
edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region, mu_p):
CreateHoleCurrent(device, region, mu_p)
PCharge = "ElectronCharge * Holes"
CreateNodeModel(device, region, "PCharge", PCharge)
CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
time_node_model = "PCharge",
edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
CreateNodeModel(device, region, "PotentialNodeCharge", pne)
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
time_node_model="", variable_update="log_damp")

def CreateSiliconDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p"):
CreatePE(device, region)
CreateBernoulli(device, region)
CreateSRH(device, region)
CreateECE(device, region, mu_n)
CreateHCE(device, region, mu_p)

def CreateSiliconDriftDiffusionAtContact(device, region, contact, is_circuit=False):
'''
Restrict electrons and holes to their equilibrium values
Integrates current into circuit
'''
contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
contact_electrons_name = "{0}nodeelectrons".format(contact)
contact_holes_name = "{0}nodeholes".format(contact)

CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
#TODO: The simplification of the ifelse statement is time consuming
#   CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

#TODO: keyword args
if is_circuit:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

else:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent")

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent")

def CreateOxidePotentialOnly(device, region, update_type="default"):
'''
Create electric field model in oxide
Creates Potential solution variable if not available
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")

efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
# this needs to remove derivatives w.r.t. independents
CreateEdgeModel(device, region, "ElectricField", efield)
CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
dfield="Permittivity*ElectricField"
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="PotentialEdgeFlux", variable_update=update_type)

def CreateSiliconOxideInterface(device, interface):
'''
continuous potential at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiliconSiliconInterface(device, interface):
'''
Enforces potential, electron, and hole continuity across the interface
'''
CreateSiliconOxideInterface(device, interface)
ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
hname = CreateContinuousInterfaceModel(device, interface, "Holes")
interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")

def CreateTemperaturefield(device,region):
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")
Source = "-heatsource"
# Init_temperature = "initialTemperature"

for i in (
("InterSource", Source),

):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Temperature")

### TODO: Edge Average Model
for i in (
("TemperatureField",     "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("TemperatureEdgeFlux", "thermalconducitivity*TemperatureField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")

equation(device=device, region=region, name="HeatconductionEquation", variable_name="Temperature",
node_model="InterSource", edge_model="TemperatureEdgeFlux", variable_update="log_damp")

def CreateTemperatureboundary(device,region,contact):
# if not InNodeModelList(device, region, "contactTemperature_node"):
#     CreateNodeModel(device, region, "contactTemperature_node", "heatsource")
# if not InEdgeModelList(device, region, "contactTemperature_edge"):
#     CreateEdgeModel(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField")
#     CreateEdgeModelDerivatives(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField", "Temperature")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Temperature -{0}".format(GetContactTemperatureName(contact))

contact_model_name = GetContactTemperatureNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Temperature"), "1")

contact_equation(device=device, contact=contact, name="HeatconductionEquation",
node_model=contact_model_name, edge_model="")
# node_charge_model="contactTemperature_node", edge_charge_model="contactTemperature_edge",
# node_current_model="", edge_current_model="")

def CreateTemperatureinterface(device,interface):

'''
continuous Temperature at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Temperature")
interface_equation(device=device, interface=interface, name="HeatconductionEquation", interface_model=model_name, type="continuous")

``````

The simple_dd.py and gmsh_mos2d_create.py don’t be changed. The error:

I recommend moving your temperature equation into where the drift diffusion equations have been set up. The first `solve` command is for the equilibrium solution, so I would recommend after this step.

I recommend creating a node solution for temperature and initializing its value to 300 using `set_node_values`. Then have the electron and hole continuity equations depend on this value. You would then need to make sure that the electron and hole continuity equations have to be accounted for terms like `mu_n`, `mu_p` and `V_t` being models instead of fixed parameters.

To prevent conflicts between parameters and models, make sure to avoid using the same name for a model name and a parameter name.

Once this is working, you can then then add the temperature equation so that temperature is treated as a solution variable.

Hi @Juan
I don’t know if I understand your suggestion. I rewrite the code as bellow:
the gmesh_mos2d.py

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create

device = "mos2d"
silicon_regions=("gate", "bulk")
oxide_regions=("oxide",)
regions = ("gate", "bulk", "oxide")
interfaces = ("bulk_oxide", "gate_oxide")

for i in regions:
CreateSolution(device, i, "Potential")

# for i in oxide_regions:
#     SetOxideParameters(device, i, 300)

for i in silicon_regions:
SetSiliconParameters(device, i, 300)
CreateSiliconPotentialOnly(device, i)

for i in oxide_regions:
SetOxideParameters(device, i, 300)
CreateOxidePotentialOnly(device, i, "log_damp")

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
tmp = get_region_list(device=device, contact=i)
r = tmp[0]
print("%s %s" % (r, i))
CreateSiliconPotentialOnlyContact(device, r, i)
# set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
set_parameter(device=device, name="drain_bias", value=1.0)
set_parameter(device=device, name="source_bias", value=0.0)
set_parameter(device=device, name="gate_bias", value=1.0)
set_parameter(device=device, name="body_bias", value=0.0)

for i in interfaces:
CreateSiliconOxideInterface(device, i)

# for d in get_device_list():
#  for gn in get_parameter_list():
#    print("{0} {1}").format(gn, get_parameter(device=d, name=gn))
#  for gn in get_parameter_list(device=d):
#    print("{0} {1} {2}").format(d, gn, get_parameter(device=d, nasme=gn))
#  for r in get_region_list(device=d):
#    for gn in get_parameter_list(device=d, region=r):
#      print("{0} {1} {2} {3}").format(d, r, gn, get_parameter(device=d, region=r, name=gn))
# write_devices(file="foo.msh", type="devsim")
#

solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=100)
# solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
# write_devices -file gmsh_mos2d_potentialonly.flps -type floops
# write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

#T field
for i in regions:
CreateSolution(device,i,"Temperature")
CreateTinital(device,i)
set_node_values(device=device, region=i, name="Temperature", init_from="Init_temperature")
CreateTemperaturefield(device,i)

for i in silicon_regions:
CreateSolution(device, i, "Electrons")
CreateSolution(device, i, "Holes")
set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
CreateSiliconDriftDiffusion(device, i)

for c in contacts:
tmp = get_region_list(device=device, contact=c)
r = tmp[0]
CreateSiliconDriftDiffusionAtContact(device, r, c)
CreateTemperatureboundary(device, r, c)
set_parameter(device=device, name="drain_biasT", value=300)
set_parameter(device=device, name="source_biasT", value=300)
set_parameter(device=device, name="gate_biasT", value=300)
set_parameter(device=device, name="body_biasT", value=300)

for i in interfaces:
CreateTemperatureinterface(device, i)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=100)

for r in silicon_regions:
node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data

for r in silicon_regions:
element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)

#
# printAllCurrents(device)
# rampbias(device, "drain",  1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)    ##current unit is A/cm
# rampbias(device, "gate", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

``````

the simple_physics.py:

``````#Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
#T      = 300 # K
eps_si = 11.1
eps_ox = 3.9
# TODO: make this temperature dependent
n_i    = 1.0e10 # #/cm^3
# constant in our approximation
mu_n   = 400
mu_p   = 200
Sithermalconducitivity = 19100  ##W/K*cm
Oxidethermalconducitivity = 14000
heatsource = 10000     ##w/cm^2
init_T = 300
boundarytemperature = 250
def GetContactBiasName(contact):
return "{0}_bias".format(contact)

def GetContactTemperatureName(contact):
return "{0}_biasT".format(contact)

def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)

def GetContactTemperatureNodeModelName(contact):
return "{0}Tnodemodel".format(contact)

def PrintCurrents(device, contact):
'''
print out contact currents
'''
contact_bias_name = GetContactBiasName(contact)
electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
total_current   = electron_current + hole_current
voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
conteq="Permittivity*ElectricField"
contact_bias_name  = GetContactBiasName(contact)
contact_model_name = GetContactNodeModelName(contact)
eq = "Potential - {0}".format(contact_bias_name)
CreateContactNodeModel(device, contact, contact_model_name, eq)
CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

#TODO: make everyone use dfield
if not InEdgeModelList(device, region, contactcharge_edge):
CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_charge_model= contactcharge_edge)

#####
##### Constants
#####
def SetOxideParameters(device, region, T):
'''
Sets physical parameters
'''
set_parameter(device=device, region=region, name="Permittivity",   value=eps_ox * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)
set_parameter(device=device, region=region, name="source",     value=heatsource)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Oxidethermalconducitivity)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)

def SetSiliconParameters(device, region, T):
'''
Sets physical parameters assuming constants
'''
#### TODO: make T a free parameter and T dependent parameters as models
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="boltzmannConstant",              value=k)
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=mu_n)
set_parameter(device=device, region=region, name="mu_p",           value=mu_p)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Sithermalconducitivity)
set_parameter(device=device, region=region, name="source",     value=heatsource)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)
set_parameter(device=device, region=region, name="boundaryTemperature",     value=boundarytemperature)
#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", value=1e-5)
set_parameter(device=device, region=region, name="taup", value=1e-5)

def CreateSiliconPotentialOnly(device, region):
'''
Creates the physical models for a Silicon region
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
elec_i = "n_i*exp(Potential/V_t)"
hole_i = "n_i^2/IntrinsicElectrons"
charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
pcharge_i = "-ElectronCharge * IntrinsicCharge"

# require NetDoping
for i in (
("IntrinsicElectrons", elec_i),
("IntrinsicHoles", hole_i),
("IntrinsicCharge", charge_i),
("PotentialIntrinsicCharge", pcharge_i)
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Potential")

### TODO: Edge Average Model
for i in (
("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
("PotentialEdgeFlux", "Permittivity * ElectricField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Potential")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
'''
Creates the potential equation at the contact
if is_circuit is true, than use node given by GetContactBiasName
'''
# Means of determining contact charge
# Same for all contacts
if not InNodeModelList(device, region, "contactcharge_node"):
CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contactcharge_edge"):
CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
-V_t*log({1}/n_i), \
V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)

contact_model_name = GetContactNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")

def CreateSRH(device, region):
USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
Gn = "-ElectronCharge * USRH"
Gp = "+ElectronCharge * USRH"
CreateNodeModel(device, region, "USRH", USRH)
CreateNodeModel(device, region, "ElectronGeneration", Gn)
CreateNodeModel(device, region, "HoleGeneration", Gp)
for i in ("Electrons", "Holes"):
CreateNodeModelDerivative(device, region, "USRH", USRH, i)
CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)

def CreateECE(device, region):
CreateElectronCurrent(device, region)

NCharge = "-ElectronCharge * Electrons"
CreateNodeModel(device, region, "NCharge", NCharge)
CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model = "NCharge",
edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region):
CreateHoleCurrent(device, region)
PCharge = "ElectronCharge * Holes"
CreateNodeModel(device, region, "PCharge", PCharge)
CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
time_node_model = "PCharge",
edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
CreateNodeModel(device, region, "PotentialNodeCharge", pne)
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
time_node_model="", variable_update="log_damp")

def CreateSiliconDriftDiffusion(device, region):

mu_n_T = "mu_n*pow((Temperature@n0+Temperature@n1)/(2*T),-2.2)"
mu_p_T = "mu_p*pow((Temperature@n0+Temperature@n1)/(2*T),-2.2)"
V_t_T = "boltzmannConstant*(Temperature@n0+Temperature@n1)/(2*ElectronCharge)"
for i in (
("Nmobility", mu_n_T),
("Pmobility", mu_p_T),
("VtT", V_t_T),
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
for i in("Temperature@n0","Temperature@n1"):
CreateEdgeModelDerivatives(device, region, n, e, i)
CreatePE(device, region)
CreateBernoulli(device, region)
CreateSRH(device, region)
CreateECE(device, region)
CreateHCE(device, region)

def CreateSiliconDriftDiffusionAtContact(device, region, contact, is_circuit=False):
'''
Restrict electrons and holes to their equilibrium values
Integrates current into circuit
'''
contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
contact_electrons_name = "{0}nodeelectrons".format(contact)
contact_holes_name = "{0}nodeholes".format(contact)

CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
#TODO: The simplification of the ifelse statement is time consuming
#   CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

#TODO: keyword args
if is_circuit:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

else:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent")

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent")

def CreateOxidePotentialOnly(device, region, update_type="default"):
'''
Create electric field model in oxide
Creates Potential solution variable if not available
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")

efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
# this needs to remove derivatives w.r.t. independents
CreateEdgeModel(device, region, "ElectricField", efield)
CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
dfield="Permittivity*ElectricField"
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="PotentialEdgeFlux", variable_update=update_type)

def CreateSiliconOxideInterface(device, interface):
'''
continuous potential at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiliconSiliconInterface(device, interface):
'''
Enforces potential, electron, and hole continuity across the interface
'''
CreateSiliconOxideInterface(device, interface)
ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
hname = CreateContinuousInterfaceModel(device, interface, "Holes")
interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")

def CreateTinital(device,region):
initial_Temperature = "initialTemperature"
for i in (
("Init_temperature", initial_Temperature),

):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Temperature")

def CreateSource(device,region):
Heatsource ="source"
CreateNodeModel(device,region,"InterSource",Heatsource)
CreateNodeModelDerivative(device, region, "InterSource", Heatsource, "Temperature")

def CreateTemperaturefield(device,region):
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")
CreateSource(device,region)
for i in (
("TemperatureField",     "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("TemperatureEdgeFlux", "thermalconducitivity*TemperatureField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")

equation(device=device, region=region, name="HeatconductionEquation", variable_name="Temperature",
node_model="InterSource", edge_model="TemperatureEdgeFlux", variable_update="log_damp")

def CreateTemperatureboundary(device,region,contact):
# if not InNodeModelList(device, region, "contactTemperature_node"):
#     CreateNodeModel(device, region, "contactTemperature_node", "heatsource")
# if not InEdgeModelList(device, region, "contactTemperature_edge"):
#     CreateEdgeModel(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField")
#     CreateEdgeModelDerivatives(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField", "Temperature")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Temperature -{0}".format(GetContactTemperatureName(contact))

contact_model_name = GetContactTemperatureNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Temperature"), "1")

contact_equation(device=device, contact=contact, name="HeatconductionEquation",
node_model=contact_model_name, edge_model="")
# node_charge_model="contactTemperature_node", edge_charge_model="contactTemperature_edge",
# node_current_model="", edge_current_model="")

def CreateTemperatureinterface(device,interface):

'''
continuous Temperature at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Temperature")
interface_equation(device=device, interface=interface, name="HeatconductionEquation", interface_model=model_name, type="continuous")

``````

the simple_dd.py:

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .model_create import *
def CreateBernoulli (device, region):
'''
Creates the Bernoulli function for Scharfetter Gummel
'''
#### test for requisite models here
EnsureEdgeFromNodeModelExists(device, region, "Potential")
vdiffstr="(Potential@n0 - Potential@n1)/VtT"
CreateEdgeModel(device, region, "vdiff", vdiffstr)
CreateEdgeModel(device, region, "vdiff:Potential@n0",  "VtT^(-1)")
CreateEdgeModel(device, region, "vdiff:Potential@n1",  "-vdiff:Potential@n0")
CreateEdgeModelDerivatives(device, region,"vdiff",vdiffstr,"Temperature@n0")
CreateEdgeModelDerivatives(device, region,"vdiff",vdiffstr,"Temperature@n1")
CreateEdgeModel(device, region, "Bern01",              "B(vdiff)")
CreateEdgeModel(device, region, "Bern01:Potential@n0", "dBdx(vdiff) * vdiff:Potential@n0")
CreateEdgeModel(device, region, "Bern01:Potential@n1", "-Bern01:Potential@n0")
CreateEdgeModelDerivatives(device, region,"Bern01","B(vdiff)","Temperature@n0")
CreateEdgeModelDerivatives(device, region,"Bern01","B(vdiff)","Temperature@n1")
#identity of Bernoulli functions
#  CreateEdgeModel(device, region, "Bern10",              "Bern01 + vdiff")
#  CreateEdgeModel(device, region, "Bern10:Potential@n0", "Bern01:Potential@n0 + vdiff:Potential@n0")
#  CreateEdgeModel(device, region, "Bern10:Potential@n1", "Bern01:Potential@n1 + vdiff:Potential@n1")

def CreateElectronCurrent(device, region):
'''
Electron current
'''
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Electrons")
EnsureEdgeFromNodeModelExists(device, region, "Holes")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
# Make sure the bernoulli functions exist
if not InEdgeModelList(device, region, "Bern01"):
CreateBernoulli(device, region)
#### test for requisite models here
#  Jn = "ElectronCharge*{0}*EdgeInverseLength*V_t*(Electrons@n1*Bern10 - Electrons@n0*Bern01)".format(mu_n)
Jn = "ElectronCharge*Nmobility*EdgeInverseLength*VtT*kahan3(Electrons@n1*Bern01,  Electrons@n1*vdiff,  -Electrons@n0*Bern01)"
#  Jn = "ElectronCharge*{0}*EdgeInverseLength*V_t*((Electrons@n1-Electrons@n0)*Bern01 +  Electrons@n1*vdiff)".format(mu_n)

CreateEdgeModel(device, region, "ElectronCurrent", Jn)
for i in ("Electrons", "Potential", "Holes","Temperature"):
CreateEdgeModelDerivatives(device, region, "ElectronCurrent", Jn, i)

def CreateHoleCurrent(device, region):
'''
Hole current
'''
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Holes")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
# Make sure the bernoulli functions exist
if not InEdgeModelList(device, region, "Bern01"):
CreateBernoulli(device, region)
##### test for requisite models here
#  Jp ="-ElectronCharge*{0}*EdgeInverseLength*V_t*(Holes@n1*Bern01 - Holes@n0*Bern10)".format(mu_p)
Jp ="-ElectronCharge*Pmobility*EdgeInverseLength*VtT*kahan3(Holes@n1*Bern01, -Holes@n0*Bern01, -Holes@n0*vdiff)"
#  Jp ="-ElectronCharge*{0}*EdgeInverseLength*V_t*((Holes@n1 - Holes@n0)*Bern01 - Holes@n0*vdiff)".format(mu_p)
CreateEdgeModel(device, region, "HoleCurrent", Jp)
for i in ("Holes", "Potential", "Electrons","Temperature"):
CreateEdgeModelDerivatives(device, region, "HoleCurrent", Jp, i)

``````

The code can run when I set the contact temperature equal to the initial temperature. But the truth is that I still found that the source term didn’t work and the heat conduction equation didn’t be calculated. Because when the iteration increases, the error of heat conduction equation are constant:

Due to the influence of the heat source, the temperature must change, but the actual output does not change from the beginning to the end. So I’m guessing that the temperature equation didn’t work at all. But because of too few errors, I couldn’t find my problem at all. Could you help me?

Hi @Juan
I found that reducing the temperature value to single digits, which is one hundredth of the original. The program is then able to perform convergent calculations.But my results show that my heat source term is not doing anything, which makes me even more confused.The new code is as follows：
the gmsh_mos2d.py:

``````# Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *

import gmsh_mos2d_create

device = "mos2d"
silicon_regions=("gate", "bulk")
oxide_regions=("oxide",)
regions = ("gate", "bulk", "oxide")
interfaces = ("bulk_oxide", "gate_oxide")

for i in regions:
CreateSolution(device, i, "Potential")

# for i in oxide_regions:
#     SetOxideParameters(device, i, 300)

for i in silicon_regions:
SetSiliconParameters(device, i, 300)
CreateSiliconPotentialOnly(device, i)

for i in oxide_regions:
SetOxideParameters(device, i, 300)
CreateOxidePotentialOnly(device, i, "log_damp")

### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
tmp = get_region_list(device=device, contact=i)
r = tmp[0]
print("%s %s" % (r, i))
CreateSiliconPotentialOnlyContact(device, r, i)
# set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
set_parameter(device=device, name="drain_bias", value=1.0)
set_parameter(device=device, name="source_bias", value=0.0)
set_parameter(device=device, name="gate_bias", value=1.0)
set_parameter(device=device, name="body_bias", value=0.0)

for i in interfaces:
CreateSiliconOxideInterface(device, i)

# for d in get_device_list():
#  for gn in get_parameter_list():
#    print("{0} {1}").format(gn, get_parameter(device=d, name=gn))
#  for gn in get_parameter_list(device=d):
#    print("{0} {1} {2}").format(d, gn, get_parameter(device=d, nasme=gn))
#  for r in get_region_list(device=d):
#    for gn in get_parameter_list(device=d, region=r):
#      print("{0} {1} {2} {3}").format(d, r, gn, get_parameter(device=d, region=r, name=gn))
# write_devices(file="foo.msh", type="devsim")
#

solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=100)
# solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
# write_devices -file gmsh_mos2d_potentialonly.flps -type floops
# write_devices(file="gmsh_mos2d_potentialonly", type="vtk")

#T field
for i in regions:
CreateSolution(device,i,"Temperature")
CreateTinitial(device,i)
set_node_values(device=device, region=i, name="Temperature", init_from="Init_temperature")
CreateTemperaturefield(device,i)

for i in silicon_regions:
CreateSolution(device, i, "Electrons")
CreateSolution(device, i, "Holes")
set_node_values(device=device, region=i, name="Electrons", init_from="IntrinsicElectrons")
set_node_values(device=device, region=i, name="Holes",     init_from="IntrinsicHoles")
CreateSiliconDriftDiffusion(device, i)

for c in contacts:
tmp = get_region_list(device=device, contact=c)
r = tmp[0]
CreateSiliconDriftDiffusionAtContact(device, r, c)
CreateTemperatureboundary(device, r, c)
set_parameter(device=device, name="drain_biasT", value=3)
set_parameter(device=device, name="source_biasT", value=3)
set_parameter(device=device, name="gate_biasT", value=3)
set_parameter(device=device, name="body_biasT", value=2)

for i in interfaces:
CreateTemperatureinterface(device, i)

solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=100)

for r in silicon_regions:
node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data

for r in silicon_regions:
element_from_edge_model(edge_model="ElectricField",   device=device, region=r)
element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
element_from_edge_model(edge_model="HoleCurrent",     device=device, region=r)

#
# printAllCurrents(device)
# rampbias(device, "drain",  1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)    ##current unit is A/cm
# rampbias(device, "gate", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")

``````

the simple_physics.py:

``````#Copyright 2013 Devsim LLC
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"

q      = 1.6e-19 # coul
k      = 1.3806503e-23 # J/K
eps_0  = 8.85e-14 # F/cm^2
#T      = 300 # K
eps_si = 11.1
eps_ox = 3.9
# TODO: make this temperature dependent
n_i    = 1.0e10 # #/cm^3
# constant in our approximation
mu_n   = 400
mu_p   = 200
Sithermalconducitivity = 1.48  ##W/K*cm
Oxidethermalconducitivity = 0.014
heatsource = 10000     ##w/cm^2
init_T = 3
boundarytemperature = 250
def GetContactBiasName(contact):
return "{0}_bias".format(contact)

def GetContactTemperatureName(contact):
return "{0}_biasT".format(contact)

def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)

def GetContactTemperatureNodeModelName(contact):
return "{0}Tnodemodel".format(contact)

def PrintCurrents(device, contact):
'''
print out contact currents
'''
contact_bias_name = GetContactBiasName(contact)
electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
hole_current    = get_contact_current(device=device, contact=contact, equation=hce_name)
total_current   = electron_current + hole_current
voltage         = get_parameter(device=device, name=GetContactBiasName(contact))
print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))

#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
conteq="Permittivity*ElectricField"
contact_bias_name  = GetContactBiasName(contact)
contact_model_name = GetContactNodeModelName(contact)
eq = "Potential - {0}".format(contact_bias_name)
CreateContactNodeModel(device, contact, contact_model_name, eq)
CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")

#TODO: make everyone use dfield
if not InEdgeModelList(device, region, contactcharge_edge):
CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")

contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_charge_model= contactcharge_edge)

#####
##### Constants
#####
def SetOxideParameters(device, region, T):
'''
Sets physical parameters
'''
set_parameter(device=device, region=region, name="Permittivity",   value=eps_ox * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)
set_parameter(device=device, region=region, name="source",     value=heatsource)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Oxidethermalconducitivity)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)

def SetSiliconParameters(device, region, T):
'''
Sets physical parameters assuming constants
'''
#### TODO: make T a free parameter and T dependent parameters as models
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="boltzmannConstant",              value=k)
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=mu_n)
set_parameter(device=device, region=region, name="mu_p",           value=mu_p)
set_parameter(device=device, region=region, name="thermalconducitivity",     value=Sithermalconducitivity)
set_parameter(device=device, region=region, name="source",     value=heatsource)
set_parameter(device=device, region=region, name="initialTemperature",     value=init_T)
set_parameter(device=device, region=region, name="boundaryTemperature",     value=boundarytemperature)
#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", value=1e-5)
set_parameter(device=device, region=region, name="taup", value=1e-5)

def CreateSiliconPotentialOnly(device, region):
'''
Creates the physical models for a Silicon region
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
elec_i = "n_i*exp(Potential/V_t)"
hole_i = "n_i^2/IntrinsicElectrons"
charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
pcharge_i = "-ElectronCharge * IntrinsicCharge"

# require NetDoping
for i in (
("IntrinsicElectrons", elec_i),
("IntrinsicHoles", hole_i),
("IntrinsicCharge", charge_i),
("PotentialIntrinsicCharge", pcharge_i)
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Potential")

### TODO: Edge Average Model
for i in (
("ElectricField",     "(Potential@n0-Potential@n1)*EdgeInverseLength"),
("PotentialEdgeFlux", "Permittivity * ElectricField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Potential")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")

def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
'''
Creates the potential equation at the contact
if is_circuit is true, than use node given by GetContactBiasName
'''
# Means of determining contact charge
# Same for all contacts
if not InNodeModelList(device, region, "contactcharge_node"):
CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contactcharge_edge"):
CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
-V_t*log({1}/n_i), \
V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)

contact_model_name = GetContactNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")

if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")

def CreateSRH(device, region):
USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
Gn = "-ElectronCharge * USRH"
Gp = "+ElectronCharge * USRH"
CreateNodeModel(device, region, "USRH", USRH)
CreateNodeModel(device, region, "ElectronGeneration", Gn)
CreateNodeModel(device, region, "HoleGeneration", Gp)
for i in ("Electrons", "Holes"):
CreateNodeModelDerivative(device, region, "USRH", USRH, i)
CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)

def CreateECE(device, region):
CreateElectronCurrent(device, region)

NCharge = "-ElectronCharge * Electrons"
CreateNodeModel(device, region, "NCharge", NCharge)
CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")

equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model = "NCharge",
edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")

def CreateHCE(device, region):
CreateHoleCurrent(device, region)
PCharge = "ElectronCharge * Holes"
CreateNodeModel(device, region, "PCharge", PCharge)
CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")

equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
time_node_model = "PCharge",
edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")

def CreatePE(device, region):
pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
CreateNodeModel(device, region, "PotentialNodeCharge", pne)
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")

equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
time_node_model="", variable_update="log_damp")

def CreateSiliconDriftDiffusion(device, region):

mu_n_T = "mu_n*pow((Temperature@n0+Temperature@n1)/(2*T/100),-2.2)"
mu_p_T = "mu_p*pow((Temperature@n0+Temperature@n1)/(2*T/100),-2.2)"
V_t_T = "boltzmannConstant*(Temperature@n0+Temperature@n1)/(2*ElectronCharge)"
for i in (
("Nmobility", mu_n_T),
("Pmobility", mu_p_T),
("VtT", V_t_T),
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
for i in("Temperature@n0","Temperature@n1"):
CreateEdgeModelDerivatives(device, region, n, e, i)
CreatePE(device, region)
CreateBernoulli(device, region)
CreateSRH(device, region)
CreateECE(device, region)
CreateHCE(device, region)

def CreateSiliconDriftDiffusionAtContact(device, region, contact, is_circuit=False):
'''
Restrict electrons and holes to their equilibrium values
Integrates current into circuit
'''
contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
contact_electrons_name = "{0}nodeelectrons".format(contact)
contact_holes_name = "{0}nodeholes".format(contact)

CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
#TODO: The simplification of the ifelse statement is time consuming
#   CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")

CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")

#TODO: keyword args
if is_circuit:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))

else:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent")

contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent")

def CreateOxidePotentialOnly(device, region, update_type="default"):
'''
Create electric field model in oxide
Creates Potential solution variable if not available
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")

efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
# this needs to remove derivatives w.r.t. independents
CreateEdgeModel(device, region, "ElectricField", efield)
CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
dfield="Permittivity*ElectricField"
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="PotentialEdgeFlux", variable_update=update_type)

def CreateSiliconOxideInterface(device, interface):
'''
continuous potential at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")

#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiliconSiliconInterface(device, interface):
'''
Enforces potential, electron, and hole continuity across the interface
'''
CreateSiliconOxideInterface(device, interface)
ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
hname = CreateContinuousInterfaceModel(device, interface, "Holes")
interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")

def CreateTinitial(device,region):
initial_Temperature = "initialTemperature"
for i in (
("Init_temperature", initial_Temperature),

):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Temperature")

def CreateSource(device,region):
Heatsource ="source"
CreateNodeModel(device,region,"InterSource",Heatsource)
CreateNodeModelDerivative(device, region, "InterSource", Heatsource, "Temperature")

def CreateTemperaturefield(device,region):
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")
CreateSource(device,region)
for i in (
("TemperatureField",     "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("TemperatureEdgeFlux", "thermalconducitivity*TemperatureField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")

equation(device=device, region=region, name="HeatconductionEquation", variable_name="Temperature",
node_model="InterSource", edge_model="TemperatureEdgeFlux",
time_node_model="", variable_update="log_damp")

def CreateTemperatureboundary(device,region,contact):
# if not InNodeModelList(device, region, "contactTemperature_node"):
#     CreateNodeModel(device, region, "contactTemperature_node", "heatsource")
# if not InEdgeModelList(device, region, "contactTemperature_edge"):
#     CreateEdgeModel(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField")
#     CreateEdgeModelDerivatives(device, region, "contactTemperature_edge", "thermalconducitivity * TemperatureField", "Temperature")

#  set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)

contact_model = "Temperature -{0}".format(GetContactTemperatureName(contact))

contact_model_name = GetContactTemperatureNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Temperature"), "1")

contact_equation(device=device, contact=contact, name="HeatconductionEquation",
node_model=contact_model_name, edge_model="")
# node_charge_model="contactTemperature_node", edge_charge_model="contactTemperature_edge",
# node_current_model="", edge_current_model="")

def CreateTemperatureinterface(device,interface):

'''
continuous Temperature at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Temperature")
interface_equation(device=device, interface=interface, name="HeatconductionEquation", interface_model=model_name, type="continuous")

``````

the simple_dd.py:

``````
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

from .model_create import *
def CreateBernoulli (device, region):
'''
Creates the Bernoulli function for Scharfetter Gummel
'''
#### test for requisite models here
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
vdiffstr="(Potential@n0 - Potential@n1)/VtT"
CreateEdgeModel(device, region, "vdiff", vdiffstr)
CreateEdgeModel(device, region, "vdiff:Potential@n0",  "VtT^(-1)")
CreateEdgeModel(device, region, "vdiff:Potential@n1",  "-vdiff:Potential@n0")
CreateEdgeModelDerivatives(device, region,"vdiff",vdiffstr,"Temperature@n0")
CreateEdgeModelDerivatives(device, region,"vdiff",vdiffstr,"Temperature@n1")
CreateEdgeModel(device, region, "Bern01",              "B(vdiff)")
CreateEdgeModel(device, region, "Bern01:Potential@n0", "dBdx(vdiff) * vdiff:Potential@n0")
CreateEdgeModel(device, region, "Bern01:Potential@n1", "-Bern01:Potential@n0")
CreateEdgeModelDerivatives(device, region,"Bern01","B(vdiff)","Temperature@n0")
CreateEdgeModelDerivatives(device, region,"Bern01","B(vdiff)","Temperature@n1")
#identity of Bernoulli functions
#  CreateEdgeModel(device, region, "Bern10",              "Bern01 + vdiff")
#  CreateEdgeModel(device, region, "Bern10:Potential@n0", "Bern01:Potential@n0 + vdiff:Potential@n0")
#  CreateEdgeModel(device, region, "Bern10:Potential@n1", "Bern01:Potential@n1 + vdiff:Potential@n1")

def CreateElectronCurrent(device, region):
'''
Electron current
'''
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Electrons")
EnsureEdgeFromNodeModelExists(device, region, "Holes")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
# Make sure the bernoulli functions exist
if not InEdgeModelList(device, region, "Bern01"):
CreateBernoulli(device, region)
#### test for requisite models here
#  Jn = "ElectronCharge*{0}*EdgeInverseLength*V_t*(Electrons@n1*Bern10 - Electrons@n0*Bern01)".format(mu_n)
Jn = "ElectronCharge*Nmobility*EdgeInverseLength*VtT*kahan3(Electrons@n1*Bern01,  Electrons@n1*vdiff,  -Electrons@n0*Bern01)"
#  Jn = "ElectronCharge*{0}*EdgeInverseLength*V_t*((Electrons@n1-Electrons@n0)*Bern01 +  Electrons@n1*vdiff)".format(mu_n)

CreateEdgeModel(device, region, "ElectronCurrent", Jn)
for i in ("Electrons", "Potential", "Holes","Temperature"):
CreateEdgeModelDerivatives(device, region, "ElectronCurrent", Jn, i)

def CreateHoleCurrent(device, region):
'''
Hole current
'''
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Holes")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
# Make sure the bernoulli functions exist
if not InEdgeModelList(device, region, "Bern01"):
CreateBernoulli(device, region)
##### test for requisite models here
#  Jp ="-ElectronCharge*{0}*EdgeInverseLength*V_t*(Holes@n1*Bern01 - Holes@n0*Bern10)".format(mu_p)
Jp ="-ElectronCharge*Pmobility*EdgeInverseLength*VtT*kahan3(Holes@n1*Bern01, -Holes@n0*Bern01, -Holes@n0*vdiff)"
#  Jp ="-ElectronCharge*{0}*EdgeInverseLength*V_t*((Holes@n1 - Holes@n0)*Bern01 - Holes@n0*vdiff)".format(mu_p)
CreateEdgeModel(device, region, "HoleCurrent", Jp)
for i in ("Holes", "Potential", "Electrons","Temperature"):
CreateEdgeModelDerivatives(device, region, "HoleCurrent", Jp, i)

``````

I give the temperature distribution of this results below:

Next, I will continue to explore why the absolute value of the temperature field affects the results of convergence, and why the heat source term I set has no effect？

Hi @ghost
Would you mind creating a github issue and attach your files as a zip file?
https://github.com/devsim/devsim/issues

It should be easier to manage a discussion. Please post the link to the issue here for the benefit of others.

Hi @Juan
Thank you for your suggestion. I will do that later.

Hi @Juan
I have created the issue in th github. I hope you can have a look at it If you have time. Thank you for your kindness.

Thanks @ghost:
https://github.com/devsim/devsim/issues/93
I will try to look at it this weekend

Hi @Juan
I try to use “edgevolumemodel” to calculate the Joule heat. I don’t know if it is right.
The code:

``````def CreateJouleHeat(device,region):
CreateEdgeModel(device,region,"JouleHeat","(ElectronCurrent+HoleCurrent)*ElectricField")
EnsureEdgeFromNodeModelExists(device, region, "Potential")
EnsureEdgeFromNodeModelExists(device, region, "Temperature")
EnsureEdgeFromNodeModelExists(device, region, "Holes")
EnsureEdgeFromNodeModelExists(device, region, "Electrons")
for i in ("Holes", "Potential", "Electrons","Temperature"):
CreateEdgeModelDerivatives(device, region, "JouleHeat", "(ElectronCurrent+HoleCurrent)*ElectricField", i)
# Heatsource ="source"
# CreateNodeModel(device,region,"InterSource",Heatsource)
# CreateNodeModelDerivative(device, region, "InterSource", Heatsource, "Temperature")

def CreateSiTemperatureField(device,region):
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")
CreateJouleHeat(device,region)
for i in (
("TemperatureField",     "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("TemperatureEdgeFlux", "thermalconducitivity*TemperatureField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")

edge_model(device=device, region=region, name="EdgeNodeVolume", equation="0.5*EdgeCouple*EdgeLength")

set_parameter(name="edge_node0_volume_model", value="EdgeNodeVolume")
set_parameter(name="edge_node1_volume_model", value="EdgeNodeVolume")
equation(device=device, region=region, name="HeatconductionEquation", variable_name="Temperature",
node_model="", edge_model="TemperatureEdgeFlux",edge_volume_model="JouleHeat",time_node_model="", variable_update="log_damp")

``````

But it can’t converge.

Hello @ghost

I was able to get it to converge with regular absolute temperatures after making numerous changes. Your `JouleHeating` term as an `edge_volume_model` term is correct, but there were several other changes required to be make this converge. This isn’t in production ready form, but I think it is a good start:
https://github.com/devsim/devsim/issues/93#issuecomment-1295018631

As of now, the Mobility is constant. The temperature dependence can be added by making the mobility models edge models in terms of `edgeTemp`, which is a model I created which is the average of the temperature on the edge.

Hi @Juan
Thanks. You help me a lot. Although the code can run, I still have some problems.
I found that if I set the contact temperature boundary as adiabatic Boundary Conditions, the USRH will appears “overflow floating point exception” if the drain_bias is bigger than 0.1V. I guess it is because the drain current is huge and the temperature rise too fast so it can’t converge. However, if I use the function “rampbias”, this phenomenon won’t happen. This is probably because the voltage changes slowly so that the data doesn’t change too fast without causing a floating point exception. This phenomenon makes me want to ask, in devism, if the result of the previous calculation will be used as the initial value of the next calculation? Or could you describe the iteration mechanism in devism in detail?
The second question is about the calculation of vector. The results shows that the Jouleheat in the calculation will get negative value in some area. But this is non-physical. The current and electricfield has same direction in theroy. Besides, the current also gets negative value. How to understand the vector value in the devsim? And the calculation between vectors？

Thanks. You help me a lot. Although the code can run, I still have some problems.
I found that if I set the contact temperature boundary as adiabatic Boundary Conditions, the USRH will appears “overflow floating point exception” if the drain_bias is bigger than 0.1V. I guess it is because the drain current is huge and the temperature rise too fast so it can’t converge. However, if I use the function “rampbias”, this phenomenon won’t happen. This is probably because the voltage changes slowly so that the data doesn’t change too fast without causing a floating point exception. This phenomenon makes me want to ask, in devism, if the result of the previous calculation will be used as the initial value of the next calculation?

The previous solution is always used as an initial guess for the next bias point in a dc solution. If it fails to converge, the simulator will restore the previous solution values. The rampbias can then try a smaller step.

Or could you describe the iteration mechanism in devism in detail?

The direct solver is using what can be called “Newton-Raphson” iterative method. Each solution is based on the previous solution. This is since such a method relies on a good initial guess in order to converge.

The second question is about the calculation of vector. The results shows that the Jouleheat in the calculation will get negative value in some area. But this is non-physical. The current and electricfield has same direction in theroy. Besides, the current also gets negative value. How to understand the vector value in the devsim? And the calculation between vectors？

My understanding is that there will be “cooling” when the current and electric field are in the opposite direction. There is potential for this to be a non-physical effect. For example in this paper:

We had to switch from using the Electric Field to the gradient of the electron and hole quasiFermi levels. Otherwise, there was suppressed current due to the velocity saturation term. You can see this if you compare the i-v curves in versions 2 and 3 of the paper.

The visualization of edge models complicated. It is possible to set them to be a `vector` or `scalar` type, but it has not worked out well in my experience.

My recommendation is to take an approach like this to convert the edge field to an element field for visualization.

``````element_from_edge_model(edge_model="EField", device=device, region=region)
element_model(device=device, region=region, name="Emag", equation="(EField_xˆ2 + EField_yˆ2)ˆ(0.5)")
``````

Hi @Juan
Thanks for your explation and help. Finally, I actully understand the reason of the negative Jouleheat source. I use another form of Heat soure which is mentioned in this paper:

Chowdhury T. Study of self-heating effects in GaN HEMTs[R]. Arizona State University, 2013.

I changed the code and I have upload in my github (I don’t know if i have right to do this. If it is wrong, please tell me).
https://github.com/ZhenglaiT/MOSFET-simulation
Besides, I am going to validate it with commercial software. Anyone interested, please discuss with me.

HI @ghost