Dear Juan, I want to know How to use the function custom_equation in time mode? Could you can give me an example?
In the document, it states that A value of True should be used for assembling bulk equations, and a value of False should be used for assembling contact and interface boundary conditions. I want to know why?
If I use this function in the setting of interface condition in transient mode, How can I get the time step and whether the original equations at the interface nodes are substituted by the equations which I set in the function custom_equation? Thank you!
As mentioned in:
https://devsim.net/models.html#custom-matrix-assembly
The time entry is provided in the call back function when requested from the simulator. The quantities should be the argument to the d/dt operation. For example:
\frac{\partial n}{\partial t}
you would provide to the simulator n \cdot \text{NodeVolume} since this needs to be integrated by volume.
Contacts
A bulk equation is assembled for each node in a region, even if it is at a contact node. If a contact equation is specified for the same nodes, the bulk equation entries are erased and replaced. For a custom equation, True means that the custom matrix entries corresponding to the contact nodes will be erased. A value of False means that the custom matrix entries will be kept and added to the terms specified in the contact_equation command.
For example, the Poisson equation is assembled at each bulk node. When assembled into the matrix, these entries are deleted for the rows corresponding to the contact nodes. The replacement contact equation would enforce the external applied bias to the device.
Interface equation for continuous mode
The bulk equation in the second region will be assembled into the same row as the bulk equation for the first region. The corresponding row in the second region will be assembled with interface_node_model specified for the interface_equation.
When True the custom matrix entry corresponding to the nodes in the second region will be permutated to the same row as the corresponding nodes in the first region. When False the custom matrix entry for the nodes in the second region will be added to the interface node model.
Interface equation for fluxterm mode
The bulk equation corresponding to the interface node in the second region will not be permutated. So the True and False values will not result in any change to the simulation matrix.
Interface equation for hybrid mode
This would be more complicated and not recommended
Since the custom_equation is a Python callback, it can track the time step in a variable captured outside the function. In this example, you can change the value for time before each solve command.
time = 0.0
def callback_function(w, t):
if t == 'TIME':
x = f(time)
.
.
.
If you give me some more information about what types of equations you wish to solve at the interface nodes, I may be able to provide the sketch of an example.
For a contact_equation, you can set:
contact_equation node_model=""
and the equation row corresponding to the contact will have no entries.
For an interface_equation, you can set:
interface_equation type="continuous" interface_model=""
and the row corresponding to the interface nodes in the second region will be empty.
If you want to replace the bulk equations on both sides of the interface, that would be tricky, but can still be done by specifying a contact equation for both regions for the same coincident nodes.
A very simple example for circuit equations is presented in testing/equation1.py.
For testing purposes, you can debug using the following command to get the sparse matrix out of the simulator.
https://devsim.net/CommandReference.html#devsim.get_matrix_and_rhs
Thanks for your reply!
The discrete format of different transient methods are different, for example bfd1, bdf2, tr, etc. Where do I define the format dn/dt ? Could you give me the script including the time entry about how to define dn/dt?
In fact, I want to know whether dn/dt needs to be discretized by myself according to the discrete format and give the corresponding matrix and rhs? When I debug my code, I have found that there are two callbacks at each iteration in the transient mode in which the first is DC and the twice callback is TIME, I want to know why?
Thank you very much for your reply！！！
If you are assembling dn/dt into a bulk equation, the RHS entry would be ElectronDensity*NodeVolume into the row corresponding to the node. The MATRIX entry would be the derivative with respect to the ElectronDensity and would be NodeVolume in the row and column corresponding to the node. To get the equation number, you could use:
https://devsim.net/CommandReference.html#devsim.get_equation_numbers
which would return the list of RHS row numbers for each node in the region.
If you need the node indexes corresponding to a certain interface or contact, you can use:
https://devsim.net/CommandReference.html#devsim.get_element_node_list
to get the corresponding node indexes in each region.
Is dn/dt the term you are looking for?
Volume integration would need to be applied by scaling the nodal quantity with NodeVolume.
Edge flux integration would need to be applied by scaling the edge quantity with EdgeCouple. This is described in:
https://doi.org/10.1109/TED.2021.3094776
and is handled automatically when models are specified in the equation command. Since you are adding directly to the matrix, the integrations would need to be applied in your code.
The reason for the separate dc and time assembly is that the transient solve will apply the time derivative operation itself. Please see chapter 7 of:
https://github.com/devsim/devsim_misc/blob/main/devsim_docs/TCADdocs.pdf
for a brief description on how the simulator applies the time integration for the different methods. Since an initial solution is required, the transient_dc solution step is required first. The easiest time method to use is bdf1. Please see testing/transient_rc.py and testing/transient_circ.py for simple rc circuit examples.
@fendou1997 if your time-derivative model is fairly simple, it should be possible to extend the simulator to support it with the standard equation commands. Please let me know what your model looks like.
That is I want to know, thank you very much!
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Aug 14 08:04:55 2021
@author: zhaochengliu
"""
import math
import numpy
from devsim import *
#devsim.set_parameter(name="debug_level", value="verbose")
device="onedimention"
create_1d_mesh (mesh=device)
add_1d_mesh_line(mesh=device, pos=0, ps=1, tag='top')
add_1d_mesh_line(mesh=device, pos=1, ps=1, tag='mid')
add_1d_mesh_line(mesh=device, pos=2, ps=1, tag='bot')
add_1d_contact (mesh=device, name='top', tag='top', material='metal')
add_1d_contact (mesh=device, name='bot', tag='bot', material='metal')
add_1d_interface(mesh=device, name="interface_insulator", tag='mid')
add_1d_region(mesh=device, material="insulator_1", region="insulator1", tag1='top', tag2='mid')
add_1d_region(mesh=device, material="insulator_2", region="insulator2", tag1='mid', tag2='bot')
finalize_mesh (mesh=device)
create_device (mesh=device, device=device)
set_parameter(device=device, region="insulator1", name="Permitivity", value=1)
set_parameter(device=device, region="insulator2", name="Permitivity", value=2)
def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)
def CreateInterfaceModel(device, interface, model, expression):
result=interface_model(device=device, interface=interface, name=model, equation=expression)
def CreateSiliconOxideInterface(device, interface):
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")
print(model_name)
def CreateContinuousInterfaceModel(device, interface, variable):
mname = "continuous{0}".format(variable)
meq = "{0}@r0 - {0}@r1".format(variable)
mname0 = "{0}:{1}@r0".format(mname, variable)
mname1 = "{0}:{1}@r1".format(mname, variable)
CreateInterfaceModel(device, interface, mname, meq)
CreateInterfaceModel(device, interface, mname0, "1")
CreateInterfaceModel(device, interface, mname1, "-1")
return mname
for region in ("insulator1", "insulator2"):
node_solution(device=device, region=region, name="Potential")
edge_from_node_model(device=device, region=region, node_model="Potential")
edge_model(device=device, region=region, name="ElectricField",
equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")
edge_model(device=device, region=region, name="ElectricField:Potential@n0",
equation="EdgeInverseLength")
edge_model(device=device, region=region, name="ElectricField:Potential@n1",
equation="-EdgeInverseLength")
edge_model(device=device, region=region, name="DField",
equation="Permitivity*ElectricField")
edge_model(device=device, region=region, name="DField:Potential@n0",
equation="diff(Permitivity*ElectricField, Potential@n0)")
edge_model(device=device, region=region, name="DField:Potential@n1",
equation="-DField:Potential@n0")
for region in ("insulator1", "insulator2"):
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="DField", variable_update="default")
CreateSiliconOxideInterface(device, "interface_insulator")
for c in ("top", "bot"):
contact_node_model(device=device, contact=c, name="%s_bc" % c,
equation="Potential - %s_bias" % c)
contact_node_model(device=device, contact=c, name="%s_bc:Potential" % c,
equation="1")
contact_equation(device=device, contact=c, name="PotentialEquation",
node_model="%s_bc" % c, edge_charge_model="DField")
set_parameter(device=device, name="top_bias", value=0.0)
set_parameter(device=device, name="bot_bias", value=0.0)
solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=30)
class InterfaceBC:
def __init__(self, device, region1, region2, interface, variable):
self.device = device
self.region1 = region1
self.region2 = region2
self.interface = interface
self.variable = variable
self.interface_nodes1 = None
self.interface_nodes2 = None
self.node_volume1 = None
self.node_volume2 = None
self.interface_bias_name = None
self.equation_numbers1= None
self.equation_numbers2= None
def initialize(self):
self.equation_numbers1=get_equation_numbers(device=self.device, region=self.region1, variable=self.variable)
interface_nodes1 = []
for e in get_element_node_list(device=self.device, region=self.region1, interface=self.interface):
interface_nodes1.extend(e)
self.interface_nodes1 = sorted(set(interface_nodes1))
self.node_volume1 = get_node_model_values(device=self.device, region=self.region1, name="NodeVolume")
self.equation_numbers2=get_equation_numbers(device=self.device, region=self.region2, variable=self.variable)
interface_nodes2 = []
for e in get_element_node_list(device=self.device, region=self.region2, interface=self.interface):
interface_nodes2.extend(e)
self.interface_nodes2 = sorted(set(interface_nodes2))
self.node_volume2 = get_node_model_values(device=self.device, region=self.region2, name="NodeVolume")
print("----")
print(interface_nodes1)
print(interface_nodes2)
def assemble(self, what, timemode):
'''
This assumes that constant bias conditions
'''
if not self.interface_nodes1:
self.initialize()
rcv=[]
rv=[]
if timemode != "DC":
print("123")
if what != "MATRIXONLY":
var = get_node_model_values(device=self.device, region=self.region1, name=self.variable)
for c in self.interface_nodes1:
r = self.equation_numbers1[c]
v = (var[c] - 1) * self.node_volume1[c]
rv.extend([r+1, 0])
if what != "RHS":
# this is the derivative
for c in self.interface_nodes1:
r = self.equation_numbers1[c]
v = self.node_volume1[c]
rcv.extend([r+1, r+1, 0])
rcv.extend([r+1, r, 0])
print("--------------------------")
print(rv)
print(rcv)
# always return False for boundary conditions
return rcv, rv, False
def getAssembleFunc(self):
return lambda what, timemode: self.assemble(what, timemode)
# how return function operating on one class
f=InterfaceBC(device=device, region1="insulator1", region2="insulator2", interface="interface_insulator", variable="Potential")
custom_equation(name="test1", procedure=f.getAssembleFunc())
### deactivate the existing contact model
node_model(device=device, region=region, name="continuousPotential", equation="0")
node_model(device=device, region=region, name="continuousPotential:Potential@r0", equation="0")
node_model(device=device, region=region, name="continuousPotential:Potential@r1", equation="0")
####
#### Ramp the bias to 0.5 Volts
####
v = 1
set_parameter(device=device, name="top_bias", value=v)
solve(type="dc", absolute_error=1e10, relative_error=1e-10, maximum_iterations=30)
write_devices(file="data.devsim", type="devsim_data")
This is the one dimensional example only contain 3 points which has 1 interface. I found that the calculated result doesn’t change when I change the rcv in the above script. However, the result will change if I change rv. By using the function get_matrix_and_rhs, I have found that the matrix is indeed changed but the result doesn’t change
If the function custom_eqaution has some error or my script has some bugs. What’s more, the matrix I return in the custom_eqaution is the equation solved in the solver? Waiting for your response.
Hi @fendou1997
Thank you for sharing your example.
The node_model command in your original script does not replace interface models. To disable the interface models, you would need to use these commands, instead of node_model.
interface_model(device=device, interface="interface_insulator", name="continuousPotential", equation="0")
interface_model(device=device, interface="interface_insulator", name="continuousPotential:Potential@r0", equation="0")
interface_model(device=device, interface="interface_insulator", name="continuousPotential:Potential@r1", equation="0")
When debugging for a capacitor, you can debug using a function like this:
def print_charge():
for c in ("top", "bot"):
print("charge %s %g" % (c, get_contact_charge(device=device, contact=c, equation="PotentialEquation")))
With the interface models, you will see this after the final solve.
charge top 0.666667
charge bot -0.666667
If you change the equations to “0” using the interface_model command above, you will see:
charge top 0
charge bot 0
If you are trying to get the same result, you should be able to use this for debugging.
To find the node pairs at each interface. You can use the coordinate_index model. You also don’t necessarily need to integrate the interface node model with NodeVolume.
In initialize:
coordinate_indexes = [
get_node_model_values(device=self.device, region=self.region1, name="coordinate_index"),
get_node_model_values(device=self.device, region=self.region2, name="coordinate_index")
]
cimap = {}
for i in self.interface_nodes1:
cimap[coordinate_indexes[0][i]] = [i, -1]
for i in self.interface_nodes2:
cimap[coordinate_indexes[1][i]][1] = i
self.node_pairs = []
for _, v in cimap.items():
self.node_pairs.append(v)
And for the equation assembly:
if what != "MATRIXONLY":
v0 = get_node_model_values(device=self.device, region=self.region1, name=self.variable)
v1 = get_node_model_values(device=self.device, region=self.region2, name=self.variable)
for n0, n1 in self.node_pairs:
r = self.equation_numbers2[n1]
v = v0[n0] - v1[n1]
rv.extend([r, v])
if what != "RHS":
for n0, n1 in self.node_pairs:
c0 = self.equation_numbers1[n0]
c1 = self.equation_numbers2[n1]
r = c1
rcv.extend([r, c0, 1.0])
rcv.extend([r, c1, -1.0])
Seem to give consistent results, including in the AbsError in the solve:
Replacing Interface Node Model continuousPotential in interface interface_insulator of material
Replacing Interface Node Model continuousPotential:Potential@r0 in interface interface_insulator of material
Replacing Interface Node Model continuousPotential:Potential@r1 in interface interface_insulator of material
number of equations 4
----
[1]
[0]
--------------------------
[2, 0.0]
[2, 1, 1.0, 2, 2, -1.0]
Iteration: 0
Device: "onedimention" RelError: 2.00000e+00 AbsError: 1.00000e+00
Region: "insulator1" RelError: 1.00000e+00 AbsError: 1.00000e+00
Equation: "PotentialEquation" RelError: 1.00000e+00 AbsError: 1.00000e+00
Region: "insulator2" RelError: 1.00000e+00 AbsError: 3.33333e-01
Equation: "PotentialEquation" RelError: 1.00000e+00 AbsError: 3.33333e-01
--------------------------
[2, 0.0]
[2, 1, 1.0, 2, 2, -1.0]
Iteration: 1
Device: "onedimention" RelError: 2.22045e-16 AbsError: 3.70074e-17
Region: "insulator1" RelError: 1.11022e-16 AbsError: 3.70074e-17
Equation: "PotentialEquation" RelError: 1.11022e-16 AbsError: 3.70074e-17
Region: "insulator2" RelError: 1.11022e-16 AbsError: 3.70074e-17
Equation: "PotentialEquation" RelError: 1.11022e-16 AbsError: 3.70074e-17
charge top 0.666667
charge bot -0.666667
Thanks for your apply. But if I want to modify the interface equation, what should I do? The following is my code which is not convergent.
class InterfaceBC:
def __init__(self, device, region1, region2, interface, variable):
self.device = device
self.region1 = region1
self.region2 = region2
self.interface = interface
self.variable = variable
self.interface_nodes1 = None
self.interface_nodes2 = None
self.node_volume1 = None
self.node_volume2 = None
self.interface_bias_name = None
self.equation_numbers1= None
self.equation_numbers2= None
def initialize(self):
self.equation_numbers1=get_equation_numbers(device=self.device, region=self.region1, variable=self.variable)
interface_nodes1 = []
for e in get_element_node_list(device=self.device, region=self.region1, interface=self.interface):
interface_nodes1.extend(e)
self.interface_nodes1 = sorted(set(interface_nodes1))
self.node_volume1 = get_node_model_values(device=self.device, region=self.region1, name="NodeVolume")
self.equation_numbers2=get_equation_numbers(device=self.device, region=self.region2, variable=self.variable)
interface_nodes2 = []
for e in get_element_node_list(device=self.device, region=self.region2, interface=self.interface):
interface_nodes2.extend(e)
self.interface_nodes2 = sorted(set(interface_nodes2))
self.node_volume2 = get_node_model_values(device=self.device, region=self.region2, name="NodeVolume")
print("----")
print(interface_nodes1)
print(interface_nodes2)
def assemble(self, what, timemode):
'''
This assumes that constant bias conditions
'''
if not self.interface_nodes1:
self.initialize()
rcv=[]
rv=[]
if timemode != "DC":
print("123")
if what == "MATRIXANDRHS":
rv.extend([2, 0.0])
rcv.extend([2, 1, 1.0])
rcv.extend([2, 2, -1.0])
rv.extend([1, 0])
rcv.extend([1, 1, 0.5])
rcv.extend([1, 0, -0.5])
print("--------------------------")
print(rv)
print(rcv)
# always return False for boundary conditions
return rcv, rv, True
def getAssembleFunc(self):
return lambda what, timemode: self.assemble(what, timemode)
# how return function operating on one class
f=InterfaceBC(device=device, region1="insulator1", region2="insulator2", interface="interface_insulator", variable="Potential")
custom_equation(name="test1", procedure=f.getAssembleFunc())
delete_interface_equation(device=device, interface="interface_insulator",name="PotentialEquation")
### deactivate the existing contact model
interface_model(device=device, interface="interface_insulator", name="continuousPotential", equation="0")
interface_model(device=device, interface="interface_insulator", name="continuousPotential:Potential@r0", equation="0")
interface_model(device=device, interface="interface_insulator", name="continuousPotential:Potential@r1", equation="0")
####
#### Ramp the bias to 0.5 Volts
####
v = 1
set_parameter(device=device, name="top_bias", value=v)
solve(type="dc", absolute_error=1e-4, relative_error=1e-4, maximum_iterations=30)
print_charge()
write_devices(file="data.devsim", type="devsim_data")
This line results in no permutation of the bulk equations.
delete_interface_equation(device=device, interface="interface_insulator",name="PotentialEquation")
Please remove the line. The interface equation command is needed so that bulk equations in row 2 are permutated into row 1. Row 2 is now empty, and can be replaced with your equation.
Please also note:
rv.extend([2, 0.0])
should be replaced with a value based on the solution variables in order to drive convergence to a correct solution.