Opening MSH files

Hello,

Using the gmsh_diode3d.py example code as reference, I was trying to load in a MSH file (which I have made sure runs fine in Gmsh).
However, I keep getting “ERROR: MeshFormat 4.1 0 8 not supported”. Using Gmsh, I tried changing the MeshFormat version to 208, 308, and 2.208 but I keep getting a similar error. I noticed that all the MSH files used in the DevSim examples are of MeshFormat 2.108, but I haven’t been able to make a file of that format.
Does anyone know how to deal with this issue?
Any help would be appreciated!

Hello @katzir
Thanks for trying devsim. Mesh formats of 2.1 0 8 and 2.2 0 8 are supported. Please check the number you are seeing in the file.

When I run:

gmsh -2 -format 'msh2' 3dblock.geo

with the latest gmsh 4.10.3 from https://gmsh.info the resulting .msh file is version “2.2 0 8”.

Hello @Juan
Thank you for the quick reply! And you’re right, I just realized that when I’m using the 2.2 0 8 version the error message is different, sorry about that.
When I try creating a mesh using version 2.2 0 8, I initially get an error in the $Elements section saying “ERROR: physical number 0 must be positive”. According the the Gmsh manual on section 4.3.1 where it describes the syntax for this version, the only part where I have 0’s in the Elements is in the “elementary” position, which in my case starts like this:
$Elements
6014
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
The file runs fine in Gmsh, so do you know why it doesn’t like the 0’s here?

Each element needs to be in a physical group. In your .geo file, you should have lines like:

Physical Surface("top") = {20};
Physical Surface("bot") = {16};
Physical Volume("bulk") = {26};

which maps the appropriate surfaces to contacts and interfaces, and volumes to regions. In this mode, Gmsh should automatically remove elements that are not in any group. From your example, an element type of 15 is a point, and would automatically be filtered out.

Okay that makes sense now, thank you so much for the help!

Hello @Juan
I have a new .msh file, but I am having trouble loading it into DevSim. I have made sure that every single volume entity belongs in a Physical Group, and I have also added some Plane Surfaces into groups since I will need them later.
Since I am having trouble with reading the file, I simplified my code for now to this (kind of follows the examples of example/gmsh_diode3d.py and testing/mesh3.py):

  load_device(file="gmshDiode2.msh")
  write_devices(file="gmshDiode_write.msh", type="devsim_data")

and also tried this:

  create_gmsh_mesh(mesh="mesh1", file="gmshDiode2.msh")
  write_devices(file="gmshDiode_write.msh", type="devsim_data")

The first block of code gives me the following error:
error
and the second block of code gives the same error but at “line: 5: syntax error” instead.

I compared my file to ‘gmsh_diode3d.msh’ (the file loaded in example/diode_common.py) but I can’t find anything different in the formatting (I changed both to MeshFormat 2.2 0 8 when comparing). For reference, this is how my .msh file starts:

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
5
2 45 “Left Side”
2 46 “Right Side”
3 1 “Bottom V”
3 2 “Cylinder V”
3 6 “Top V”
$EndPhysicalNames
$Nodes
3256
1 0 0 0
2 0 8e-06 0
3 8e-06 8e-06 0

(Just in case, I made sure that these Physical Groups also show up in the .geo file, like this:

Physical Surface(“Left Side”) = {21, 142};
Physical Surface(“Right Side”) = {23, 150};
Physical Volume(“Bottom V”) = {26};
Physical Volume(“Cylinder V”) = {27, 28, 29, 30};
Physical Volume(“Top V”) = {31};

Is there some issue that I am not seeing? Or could this perhaps be an error caused by something else completely and it’s just showing up as a syntax error in the console?

(PS: I also tried making a msh file where all entities were in physical groups, but that didn’t seem to make a difference)

Hi @katzir

Sorry for the confusion about the “.msh” extension, because it is used as either a devsim format file or a gmsh format file. If you are trying to read a gmsh format, you should use create_gmsh_mesh. You also need to map the physical groups to contacts, interfaces, and regions. From diode_common.py:

def Create3DGmshMesh(device, region):
    #this reads in the gmsh format
    create_gmsh_mesh (mesh="diode3d", file="gmsh_diode3d.msh")
    add_gmsh_region  (mesh="diode3d", gmsh_name="Bulk",    region=region, material="Silicon")
    add_gmsh_contact (mesh="diode3d", gmsh_name="Base",    region=region, material="metal", name="top")
    add_gmsh_contact (mesh="diode3d", gmsh_name="Emitter", region=region, material="metal", name="bot")
    finalize_mesh    (mesh="diode3d")
    create_device    (mesh="diode3d", device=device)

Please note that a device is not created until finalize_mesh and create_device are called. The write_devices command is for writing out a file in the devsim format, and can contain simulation data, as well as equations for the models.

When physical groups are defined, Gmsh should automatically prune out entities that are not associated with any region.

Looking at GmshScanner.l in the devsim source code, I can see that the parser is not expecting spaces in the physical names, so I would expect that to be the source of your issue on line 5.

<PHYSICALNAMES_SEC>\"([A-Za-z]|{UTF8})([A-Za-z_0-9:@]|{UTF8})*\" {

Please suggest how you would want this issue to be resolved.
At the very least, the simulator can provide a more specific error. Supporting spaces in the names may also be possible.

Hello @Juan
Thank you for the clarification. Since I didn’t have too many physical groups, it was easy to just change the names so they don’t have spaces and now it works, thank you!
This doesn’t seem like too big of an issue to resolve, so I guess simply adding something about naming physical groups in the documentation would be nice.

I also have a follow-up question to this. The way I understand it (please correct me if I’m wrong), in the example of gmsh_diode3d.py, only one region is created since there was only one physical group containing a volume, the “Bulk”, and then the doping within the volume is based on an equation, using node_model().

In my case, I have 3 different physical groups of volumes. When put together, they make a rectangular volume, as shown in the picture (normally there are more diodes next to each other and their volumes become part of the physical groups mentioned in the earlier post. I just simplified this image to try and make it clearer what I’m trying to do):

The “CylinderV” (the volumes of the 4 cylinders shown in the picture) would act as P-junctions and the rest of the volumes as N-junctions. So they would be made of the same material but need to have different doping.
Is there a way to do that without basing the doping on an equation, as done in gmsh_diode3d.py?
Also, since I have these different volumes within the big structure that I want to keep track of, how exactly would I need to create them, meaning would they all be created as regions? Or maybe add a physical group with the total volume to make a region that, and then somehow keep track of the smaller components?

Hi @katzir

For step doping like that, it may be better to have a separate region for each. Creating a Physical Surface for all of the interfaces is tedious, but you can use a script like this:
https://github.com/devsim/devsim_misc/tree/main/gmsh

One issue I ran into with Gmsh and 3d meshes was the poor quality of some of the elements. In this paper I discuss some of the issues I ran into:
https://doi.org/10.36227/techrxiv.14129081.v3

I just want to you to be aware of this when working with complicated structures. I have some ideas how to fix this, which we can discuss this further if you encounter any issues.

Thank you for the resources and the heads-up!
I’ll take a look at those and try to work it out. I’ll definitely let you know if I run into any other issues.

Hi @Juan

So I tried using the script you mentioned, but I get an error when I input my MSH file (I also double checked that I can actually run the example in the link you sent, just to make sure I’m running this right):


The furthest back that I could trace the error from the functions called is in the run() function of test_convert.py (line 384). The way I understand it, it seems like yaml_map is populated in line 391 with the new interface ‘N_junctions_P_junctions’, which is then assigned to interface_names (in line 422), which is then passed into get_interface_map (line 424-425), which is the function where the error occurs (line 258). Since ‘N_junctions’ and ‘P_junctions’ are my only volume physical groups, I think it makes sense that this interface was created, but I’m not sure how to go about dealing with this.
Do you know why this is happening?

Also, since last time, I change my Physical Groups a bit, which now are (from .geo file):

Physical Surface(“Top”) = {27, 51, 74, 97, . . .};
Physical Surface(“Bottom”) = {1};
Physical Volume(“N_junctions”) = {1, 38};
Physical Volume(“P_junctions”) = {2, 3, 4, 5, . . .};

Using the following picture as reference, ‘P_junctions’ are the cylinder volumes, ‘N_junctions’ are the rest of the volumes, ‘Bottom’ is the blue plane, and ‘Top’ are the orange planes which are made up of the circles and the plane around them.

When using DevSim, I also tried having different regions for each of the volume Physical Groups, but ran into another problem. I made ‘Bottom’ and ‘Top’ to be contacts in the region created by the ‘N_junctions’, but I get an error that many of the coordinates in the ‘Top’ contact are not on that region. I get a similar error when I try to make the contact using the region created by the ‘P_junctions’.
Do you know if this is related to why I’m getting the error when I try to use test_convert.py? Or is this a completely new problem, and if so, do you know how I can fix it?

HI @katzir

I think the first issue you present is a bug. The script may not be handling the case where you have one physical group made up of multiple elementary ids for surfaces.

For the second issue, the script is not handling the case where you have one plane over multiple regions. You may be able to cut holes in the plane. Alternatively, I can make the script flexible to allow you to create a contact that ignores coordinates not on the region of interest.

It should be noted that having contacts touching each other or near interface may lead to convergence problems in the simulation. In this case, you may want to add surface rings between the plane and the surface circles for the p regions.

Would you be willing to share your geo file with me? I will not redistribute it without your permission. If possible please send it to info@devsim.com.

Hi @Juan

I just sent an email with my geo file.

Regarding the second issue, for the purposes of the project I’m working on, unfortunately it wouldn’t make sense to make holes in the plane of the contact or add something between the P-regions and the contact.

If you could make the script flexible with the contacts like you suggested, I would really appreciate it! If that were to be adjusted however, similar to gmsh_diode3d.py, if I were to make the material of the contacts to be metal and pass a current through the Top one, would it still affect both the P and N regions, since you said it would “ignore coordinates not on the region of interest”?
Because essentially this is my goal: to have a block of silicone of P and N junctions with metal contacts at the Top and Bottom, where there is a current passing through the Top contact, affecting both the P and N junctions.

Hi @katzir

I apologize for the delay. There was a bug in the script implementation, please try the newest version of:
https://github.com/devsim/devsim_misc/blob/main/gmsh/test_convert.py
The script will then support finding interfaces when there are multiple volumes in the same region.

As for your other issue concerning contacts, there are 2 options:
Option 1:
You can separate the contacts out manually:

Physical Surface("Top_p") = {27, 51, 74, 97};
Physical Surface("Top_n") = {120};

so that the contacts are separated on a per volume basis.

Option 2:

Make both “Top_volume” and “Bottom_volume” into volumes sitting above and below the device. Then create a yaml file containing lines like this:

# contact_regions are used to create contacts at their interface with other regions
# they can be removed from the final mesh by setting "remove: True"
contact_regions:
  - contact: Top_volume
    remove: True
  - contact: Bottom_volume
    remove: True

Of course, the yaml file is useful in that these volumes will be removed from the final mesh. In addition, these names have a higher priority, so there will not be any elements erased if they conflict with any of the generated interfaces. Please note that did not try this with your example, but I know it to work with the ring example.

Option 3:

Fixing this in the simulator would take a lot of time, and I am currently not able to devote resources to working on this. However, code contributions are welcome.

As an additional note, it is possible to drive multiple contacts with the same circuit node.

For example, you could create a circuit node “Top”, which is driving “Top_n” and “Top_p”.

Please see examples/diode/ssac_diode.py for an example of driving a diode from an external voltage source. This was to drive a diode with a small-signal voltage, and then get both the real and imaginary components of the small-signal current.

Hello @Juan
Thank you so much for your help!
I’m not sure when I’ll be able to try out these options, but when I do I will let you know if something doesn’t work or if I figure out something helpful.

Hi @Juan

Since last time, I ended up simplifying my mesh model a lot in order to try and get that to work before attempting the more complicated model again, and now it’s kind of similar to the one you mentioned here: https://github.com/devsim/devsim_misc/blob/main/gmsh/test.msh
I “corrected” the MSH file by running test_convert.py on test.msh and loaded the new file into test.py.
Based on the code in https://github.com/devsim/devsim_misc/blob/main/gmsh/test.py and the “corrected” MSH file, I then added the following code:

write_devices (file="test_out.msh")
diode_common.SetParameters(device=device, region="bulk")
diode_common.SetParameters(device=device, region="disk")
diode_common.SetParameters(device=device, region="nlayer")

# Same carrier number, reference at z=0.5e-5
node_model(device=device, region="bulk", name="Acceptors", equation="1.0e18*step(0.5e-5-z);")
node_model(device=device, region="bulk", name="Donors",    equation="1.0e18*step(z-0.5e-5);")
node_model(device=device, region="bulk", name="NetDoping", equation="Donors-Acceptors;")

node_model(device=device, region="disk", name="Acceptors", equation="1.0e18*step(0.5e-5-z);")
node_model(device=device, region="disk", name="Donors",    equation="1.0e18*step(z-0.5e-5);")
node_model(device=device, region="disk", name="NetDoping", equation="Donors-Acceptors;")

node_model(device=device, region="nlayer", name="Acceptors", equation="1.0e18*step(0.5e-5-z);")
node_model(device=device, region="nlayer", name="Donors",    equation="1.0e18*step(z-0.5e-5);")
node_model(device=device, region="nlayer", name="NetDoping", equation="Donors-Acceptors;")

diode_common.InitialSolution(device, "bulk")
diode_common.InitialSolution(device, "disk")
diode_common.InitialSolution(device, "nlayer")

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

diode_common.DriftDiffusionInitialSolution(device, "bulk")
diode_common.DriftDiffusionInitialSolution(device, "disk")
diode_common.DriftDiffusionInitialSolution(device, "nlayer")

solve(type="dc", absolute_error=1e15, relative_error=1e-3, maximum_iterations=50)

# Comment-out this while-loop for noCurrent case
v = 0.0
while v < 0.51:
    set_parameter(device=device, name=GetContactBiasName("top"), value=v)
    solve(type="dc", absolute_error=1e5, relative_error=1e-3, maximum_iterations=30)
    PrintCurrents(device, "top")
    PrintCurrents(device, "bot")
    v += 0.1

write_devices(file="test_withCurrent", type="vtk")

Please correct me if my logic for the above code is wrong.
Assuming it is correct, I realized a few things when running this:

  1. When the doping is defined using z = 0.5e-5 as reference (which lies in the ‘disk’ region), whether there is the same number of carriers in both regions or not, changes are only shown in the ‘disk’ region after passing a current through the top contact.

  2. Similarly, when the doping is defined using let’s say z = -0.5e-5 as reference (which lies in the ‘bulk’ region), whether there is the same number of carriers in both regions or not, changes are only shown in the ‘bulk’ region after passing a current through the top contact.

  3. When the doping is defined using z = 0 as reference (which lies on ‘disk_bulk’, the plane/interface between the ‘disk’ and ‘bulk’ regions) and the number of carriers is different in the two regions, changes are only shown in the ‘disk’ region. This is how I changed the doping equations for this (I haven’t written down the equations again in all my points, but they follow the same logic):

# Different carrier numbers, reference at z=0
node_model(device=device, region="bulk", name="Acceptors", equation="1.0e17*step(0-z);")
node_model(device=device, region="bulk", name="Donors",    equation="1.0e14*step(z-0);")
node_model(device=device, region="bulk", name="NetDoping", equation="Donors-Acceptors;")
# same Acceptor/Donor equations for the rest of the regions
  1. Finally, when the doping is defined using z = 0 as reference and the number of carriers is the same in all the regions, changes are shown in the whole model.

My question is, is there a way to have different doping in the two regions and, after passing a current, see the changes/effects in both regions?

I used VisIt for visualizing these results and making the above conclusions, using the file created after calling write_devices().
Since the numbers that are printed out change when I run test.py, I’m open to the possibility that perhaps I need to change the type of the file where the data is saved or maybe the way I display the data in VisIt. But I come to similar conclusions when I use test.py and when I run my own file, which is making me think that maybe the issue is related to the interfaces or with the way I’m setting up the doping equations (I checked that the Acceptors/Holes and Donors/Electrons are displayed in the correct regions of the model before moving on to checking other information).

As reference, here are some of the results I get in VisIt, for first piece of code I showed, corresponding to my 1st point (the change in the noCurrent vs withCurrent cases can be seen in the ‘disk’ region, where the range of Potential values changes):


Update: For testing purposes, I tried creating the whole model as one volume and changing the doping based on the z-coordinate. This gave me results that make sense (in VisIt), so I’m starting to lean towards the case where my issue may be related to the interfaces.

Your update is correct in that the issue is related to interfaces. Please see these functions in the python_packages/simple_physics module:

The CreateSiliconOxideInterface function can be used to create a continuous potential condition for the potential only condition.

The CreateSiliconSiliconInterface command can be used to create a continuous potential, electron, and hole condition across an interface.

There are examples in the code for the first function, but not the second.

Also in the tests, I can see that testing/res2.py and testing/res3.py uses

test_common.SetupContinuousElectronsAtInterface(device, interface)
test_common.SetupElectronSRVAtInterface(device, interface)

which enforce either continuous electron density or a surface recombination velocity limited electron flow across a boundary.

For future use, it is possible to use these functions to programmatically setup these conditions:

# use on device
devsim.get_interface_list
# use on interface
devsim.get_region_list
# use on region
devsim.get_material

Also, for debugging purposes you can use the asinh and log math functions to create a model showing signed doping quantities logarithmically. VisIt also has similar capabilities for creating new models.

asinh(NetDoping)/log(10)

Hi @Juan

Thank you for the list of functions and sorry for the late reply.
I’ve been trying to get these functions to work, but I haven’t had much luck with them.
For easier reference, I tried making this work in the test.py file that I mentioned last time:
https://github.com/devsim/devsim_misc/blob/main/gmsh/test.py

For my case, I think CreateSiliconSiliconInterface would be more helpful, but I couldn’t get it to work.
So instead I tried to start with using CreateSiliconOxideInterface, for which I found an example in examples\mobility\gmsh_mos2d_kla.py.
At first, I tried simply calling CreateSiliconOxideInterface right after setting the Net Doping equations, but when the solve() function was called it gave me an error that the equation for continuousPotential evaluates to an invalid.
So I modified it to try and follow the structure of gmsh_mos2d_kla.py. With this, I was able to get the Potential equation, but then when the solve() function is called I get a DEVSIM FATAL error, and this is printed out right before the error:

Device: mydevice Region: nlayer contact: bot Contact Equation: ElectronContinuityEquation Circuit node “bot_bias” does not exist

I tried to find where the ElectronContinuityEquation is set up but I haven’t been able to make any helpful changes.
As of right now, this is what I have in the file test.py up to the point where the error occurs:
(For now, I avoided calling functions that call multiple functions in them to be able to see better what is going in on)

from devsim import *
from devsim.python_packages.simple_physics import *
import diode_common

mesh="mymesh"
device="mydevice"

create_gmsh_mesh(mesh=mesh, file='newTest.msh')
add_gmsh_region  (mesh=mesh, gmsh_name="bulk",    region="bulk", material="Silicon")
add_gmsh_region  (mesh=mesh, gmsh_name="disk",    region="disk", material="Silicon")
add_gmsh_region  (mesh=mesh, gmsh_name="nlayer",    region="nlayer", material="Silicon")
add_gmsh_contact (mesh=mesh, gmsh_name="top",    region="disk", material="metal", name="top")
add_gmsh_contact (mesh=mesh, gmsh_name="bot",    region="nlayer", material="metal", name="bot")
add_gmsh_interface (mesh=mesh, gmsh_name="disk_bulk", region0="disk", region1="bulk", name="disk_bulk")
add_gmsh_interface (mesh=mesh, gmsh_name="bulk_nlayer", region0="bulk", region1="nlayer", name="bulk_nlayer")
finalize_mesh    (mesh=mesh)
create_device    (mesh=mesh, device=device)

write_devices (file="test_out.msh")

# copied for other file
def SetupContinuousPotentialAtInterface(device, interface):
    # type continuous means that regular equations in both regions are swapped into the primary region
    interface_model(device=device, interface=interface, name="continuousPotential", equation="Potential@r0-Potential@r1")
    interface_model(device=device, interface=interface, name="continuousPotential:Potential@r0", equation= "1")
    interface_model(device=device, interface=interface, name="continuousPotential:Potential@r1", equation="-1")
    interface_equation(device=device, interface=interface, name="PotentialEquation", variable_name="Potential", interface_model="continuousPotential", type="continuous")

# copied for other file
def SetupContinuousElectronsAtInterface(device, interface):
    # type continuous means that regular equations in both regions are swapped into the primary region
    interface_model(device=device, interface=interface, name="continuousElectrons", equation="Electrons@r0-Electrons@r1")
    interface_model(device=device, interface=interface, name="continuousElectrons:Electrons@r0", equation= "1")
    interface_model(device=device, interface=interface, name="continuousElectrons:Electrons@r1", equation="-1")
    interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", variable_name="Electrons", interface_model="continuousElectrons", type="continuous")

# copied for other file
def DriftDiffusionInitialSolution(device, region, circuit_contacts=None):
    #### drift diffusion solution variables
    CreateSolution(device, region, "Electrons")
    CreateSolution(device, region, "Holes")

    #### create initial guess from dc only solution
    set_node_values(device=device, region=region, name="Electrons", init_from="IntrinsicElectrons")
    set_node_values(device=device, region=region, name="Holes",     init_from="IntrinsicHoles")

    ###
    ### Set up equations
    ###
    CreateSiliconDriftDiffusion(device, region)
    for i in get_contact_list(device=device):
        if circuit_contacts and i in circuit_contacts:
            CreateSiliconDriftDiffusionAtContact(device, region, i, True)
        else:
            CreateSiliconDriftDiffusionAtContact(device, region, i)


diode_common.SetParameters(device=device, region="bulk")
diode_common.SetParameters(device=device, region="disk")
diode_common.SetParameters(device=device, region="nlayer")

# Same carrier number, reference at z=0.5e-5
node_model(device=device, region="bulk", name="Acceptors", equation="1.0e18*step(0-z);")
node_model(device=device, region="bulk", name="Donors",    equation="1.0e18*step(z-0);")
node_model(device=device, region="bulk", name="NetDoping", equation="Donors-Acceptors;")

node_model(device=device, region="disk", name="Acceptors", equation="1.0e18*step(0-z);")
node_model(device=device, region="disk", name="Donors",    equation="1.0e18*step(z-0);")
node_model(device=device, region="disk", name="NetDoping", equation="Donors-Acceptors;")

node_model(device=device, region="nlayer", name="Acceptors", equation="1.0e18*step(0-z);")
node_model(device=device, region="nlayer", name="Donors",    equation="1.0e18*step(z-0);")
node_model(device=device, region="nlayer", name="NetDoping", equation="Donors-Acceptors;")

CreateSolution(device=device, region="bulk", name = "Potential")
CreateSolution(device=device, region="disk", name = "Potential")
CreateSolution(device=device, region="nlayer", name = "Potential")
# Having these next 6 lines or not, doesn't seem to make a difference
CreateSolution(device=device, region="bulk", name = "Electrons")
CreateSolution(device=device, region="disk", name = "Electrons")
CreateSolution(device=device, region="nlayer", name = "Electrons")
CreateSolution(device=device, region="bulk", name = "Holes")
CreateSolution(device=device, region="disk", name = "Holes")
CreateSolution(device=device, region="nlayer", name = "Holes")

SetSiliconParameters(device=device, region="bulk", T=300)
SetSiliconParameters(device=device, region="disk", T=300)
SetSiliconParameters(device=device, region="nlayer", T=300)

CreateSiliconPotentialOnly(device=device, region="bulk")
CreateSiliconPotentialOnly(device=device, region="disk")
CreateSiliconPotentialOnly(device=device, region="nlayer")

SetOxideParameters(device=device, region="bulk", T=300)
SetOxideParameters(device=device, region="disk", T=300)
SetOxideParameters(device=device, region="nlayer", T=300)

CreateOxidePotentialOnly(device=device, region="bulk")
CreateOxidePotentialOnly(device=device, region="disk")
CreateOxidePotentialOnly(device=device, region="nlayer")

CreateSiliconPotentialOnlyContact(device=device, region="disk", contact="top")
CreateSiliconPotentialOnlyContact(device=device, region="nlayer", contact="bot")

# If next 2 lines are commented out, error changes to:
# Value for "bot_bias" not available
#set_parameter(device=device, region="disk", name=GetContactBiasName(contact="top"), value=0.0)
#set_parameter(device=device, region="nlayer", name=GetContactBiasName(contact="bot"), value=0.0)

CreateSiliconOxideInterface(device=device, interface="disk_bulk")
CreateSiliconOxideInterface(device=device, interface="bulk_nlayer")
###########################################################

# Error occurs in this line
solve(type="dc", absolute_error=1.0, relative_error=1e-5, maximum_iterations=30)

(I also tried using SetupContinuousElectronsAtInterface but there seems to be a similar issue with getting the ElectronContinuityEquation)

Do you know what I’m missing or what I’m doing wrong?
If possible, could you explain how to correctly use CreateSiliconSiliconInterface in the case of test.py?

@katzir
Thanks for providing your script. The “invalid” error is the result of there not being a parameter by the name of “bot_bias”. This parameter must be available before the first solve command. It is generated by with these commands:

set_parameter(device=device, region="disk", name=GetContactBiasName(contact="top"), value=0.0)
set_parameter(device=device, region="nlayer", name=GetContactBiasName(contact="bot"), value=0.0)

where the GetContactBiasName is a simple function to name the parameter related to the contact.

In doing this the simulation converges in 2 iterations, since you are solving oxide equations where the doping is not accounted for.

By commenting out these lines:

#SetOxideParameters(device=device, region="bulk", T=300)
#SetOxideParameters(device=device, region="disk", T=300)
#SetOxideParameters(device=device, region="nlayer", T=300)

#CreateOxidePotentialOnly(device=device, region="bulk")
#CreateOxidePotentialOnly(device=device, region="disk")
#CreateOxidePotentialOnly(device=device, region="nlayer")

you are then solving a Silicon system where only the Poisson equation is being solved at equilibrium. The CreateSiliconOxideInterface command is equivalent to enforcing the continuous potential across the boundary. Please think of it as being equivalent to a function named CreateSiliconSiliconPotentialOnlyInterface.

The solution for this problem is then the initial condition for the full drift diffusion equation. Which can then be copied for each region using:

diode_common.DriftDiffusionInitialSolution(device, region)

But you will need to add the interfaces equations using CreateSiliconSiliconInterface, before issuing the solve command for the initial drift diffusion simulation.

1 Like