Dean Huang
233 lines
7.3 KiB
Python
233 lines
7.3 KiB
Python
# Copyright 2013 DEVSIM LLC
|
|
#
|
|
# SPDX-License-Identifier: Apache-2.0
|
|
|
|
# The purpose is to verify our triangle element field calculation.
|
|
# It is based on Laux's weighting scheme
|
|
# @article{Laux:1985,
|
|
# author = {Laux, Steven E. and Byrnes, Robert G.},
|
|
# title = {Semiconductor device simulation using generalized mobility models},
|
|
# journal = {IBM J. Res. Dev.},
|
|
# issue_date = {May 1985},
|
|
# volume = {29},
|
|
# number = {3},
|
|
# month = may,
|
|
# year = {1985},
|
|
# issn = {0018-8646},
|
|
# pages = {289--301},
|
|
# numpages = {13},
|
|
# url = {http://dx.doi.org/10.1147/rd.293.0289},
|
|
# doi = {10.1147/rd.293.0289},
|
|
# acmid = {1012099},
|
|
# publisher = {IBM Corp.},
|
|
# address = {Riverton, NJ, USA},
|
|
# }
|
|
|
|
import sys
|
|
|
|
try:
|
|
import numpy
|
|
import numpy.linalg
|
|
except ImportError:
|
|
print(
|
|
"numpy is not available with your installation and is required for this program to run."
|
|
)
|
|
sys.exit(-1)
|
|
|
|
|
|
from devsim import (
|
|
edge_from_node_model,
|
|
element_from_edge_model,
|
|
element_from_node_model,
|
|
element_model,
|
|
get_edge_model_values,
|
|
get_element_model_values,
|
|
load_devices,
|
|
)
|
|
|
|
|
|
def calculateValues(scalar_efield, eecouple, sx, sy):
|
|
"""
|
|
calculates values, using the exact algorithm used internally to the software
|
|
"""
|
|
_row0 = numpy.array([0, 0, 1])
|
|
_row1 = numpy.array([1, 2, 2])
|
|
ans = numpy.zeros((3, 1))
|
|
M = numpy.zeros((2, 2))
|
|
B = numpy.zeros((2, 1))
|
|
e_this = [0] * 2
|
|
e_opp = [0] * 2
|
|
efx = [0] * len(scalar_efield)
|
|
efy = [0] * len(scalar_efield)
|
|
for ei in range(len(scalar_efield) // 3):
|
|
for i in range(3):
|
|
e_base = 3 * ei + i
|
|
k = 0
|
|
for j in range(3):
|
|
if i == j:
|
|
continue
|
|
e_this[k] = 3 * ei + j
|
|
k = k + 1
|
|
|
|
e_opp[0] = e_this[1]
|
|
e_opp[1] = e_this[0]
|
|
|
|
M[0][0] = sx[e_base]
|
|
M[0][1] = sy[e_base]
|
|
B[0] = scalar_efield[e_base]
|
|
|
|
ex = 0
|
|
ey = 0
|
|
den = 0
|
|
for j in range(2):
|
|
M[1][0] = sx[e_this[j]]
|
|
M[1][1] = sy[e_this[j]]
|
|
B[1] = scalar_efield[e_this[j]]
|
|
ans = numpy.linalg.solve(M, B)
|
|
ex = ex + ans[0] * eecouple[e_this[j]]
|
|
ey = ey + ans[1] * eecouple[e_this[j]]
|
|
den = den + eecouple[e_this[j]]
|
|
efx[e_base] = ex / den
|
|
efy[e_base] = ey / den
|
|
return (efx, efy)
|
|
|
|
|
|
load_devices(file="trimesh3.msh")
|
|
device = "MyDevice"
|
|
region = "MyRegion"
|
|
# trying to verify Laux current waiting implementation
|
|
# these are the scalar values
|
|
element_model(
|
|
device=device, region=region, name="scalar_efield", equation="ElectricField"
|
|
)
|
|
element_model(device=device, region=region, name="sx", equation="unitx")
|
|
element_model(device=device, region=region, name="sy", equation="unity")
|
|
|
|
scalar_efield = get_element_model_values(
|
|
device=device, region=region, name="scalar_efield"
|
|
)
|
|
sx = get_element_model_values(device=device, region=region, name="sx")
|
|
sy = get_element_model_values(device=device, region=region, name="sy")
|
|
|
|
element_from_edge_model(device=device, region=region, edge_model="ElectricField")
|
|
efieldx = get_element_model_values(device=device, region=region, name="ElectricField_x")
|
|
efieldy = get_element_model_values(device=device, region=region, name="ElectricField_y")
|
|
eecouple = get_element_model_values(
|
|
device=device, region=region, name="ElementEdgeCouple"
|
|
)
|
|
|
|
|
|
(ex, ey) = calculateValues(scalar_efield, eecouple, sx, sy)
|
|
for i in range(len(ex)):
|
|
print("%g\t%g\t%g\t%g" % (ex[i][0], ey[i][0], efieldx[i], efieldy[i]))
|
|
|
|
#
|
|
# now verify derivatives
|
|
#
|
|
edge_from_node_model(node_model="node_index", device=device, region=region)
|
|
edge_nodes = []
|
|
potential_diff = []
|
|
nedge = 0
|
|
for i in range(2):
|
|
tmp = get_edge_model_values(device=device, region=region, name="node_index@n%d" % i)
|
|
if not edge_nodes:
|
|
nedge = len(tmp)
|
|
edge_nodes = [0] * nedge
|
|
potential_diff = [0] * nedge
|
|
for j in range(nedge):
|
|
edge_nodes[j] = [0] * 2
|
|
potential_diff[j] = [0] * 2
|
|
for j in range(nedge):
|
|
edge_nodes[j][i] = tmp[j]
|
|
tmp = get_edge_model_values(
|
|
device=device, region=region, name="ElectricField:Potential@n%d" % i
|
|
)
|
|
for j in range(nedge):
|
|
potential_diff[j][i] = tmp[j]
|
|
|
|
|
|
element_model(
|
|
device=device, region=region, name="scalar_edge_index", equation="edge_index"
|
|
)
|
|
scalar_edge_index = get_element_model_values(
|
|
device=device, region=region, name="scalar_edge_index"
|
|
)
|
|
scalar_edge_index = [int(x) for x in scalar_edge_index]
|
|
element_from_node_model(node_model="node_index", device=device, region=region)
|
|
|
|
nelem = len(efieldx)
|
|
element_node_indexes = [0] * nelem
|
|
ddata_x = [0] * nelem
|
|
ddata_y = [0] * nelem
|
|
cdata_x = [0] * nelem
|
|
cdata_y = [0] * nelem
|
|
for i in range(nelem):
|
|
element_node_indexes[i] = [0] * 3
|
|
ddata_x[i] = [0] * 3
|
|
ddata_y[i] = [0] * 3
|
|
cdata_x[i] = [0] * 3
|
|
cdata_y[i] = [0] * 3
|
|
for i in range(3):
|
|
tmp = get_element_model_values(
|
|
device=device, region=region, name="node_index@en%d" % i
|
|
)
|
|
for j in range(nelem):
|
|
element_node_indexes[j][i] = int(tmp[j])
|
|
|
|
element_from_edge_model(
|
|
device=device, region=region, edge_model="ElectricField", derivative="Potential"
|
|
)
|
|
for i in range(3):
|
|
tmp = get_element_model_values(
|
|
device=device, region=region, name="ElectricField_x:Potential@en%d" % i
|
|
)
|
|
for j in range(nelem):
|
|
cdata_x[j][i] = tmp[j]
|
|
for i in range(3):
|
|
tmp = get_element_model_values(
|
|
device=device, region=region, name="ElectricField_y:Potential@en%d" % i
|
|
)
|
|
for j in range(nelem):
|
|
cdata_y[j][i] = tmp[j]
|
|
|
|
diff_field = [0] * nelem
|
|
dest_map = [0] * nelem # maps global node of element to this particular edge
|
|
for dnode in range(3): # this is the node we are taking the derivative with respect to
|
|
# it is the global node index and will be figure out shortly
|
|
for i in range(nelem // 3):
|
|
bi = 3 * i # this is the base element edge(0)
|
|
dindex = element_node_indexes[bi][
|
|
dnode
|
|
] # this is the node we are taking the derivative with respect to
|
|
for j in range(3): # this is the edge of interest
|
|
ei = bi + j # this is the edge we are looking at
|
|
sei = scalar_edge_index[ei]
|
|
if dindex in edge_nodes[sei]:
|
|
pos = edge_nodes[sei].index(dindex)
|
|
diff_field[ei] = potential_diff[sei][pos]
|
|
else:
|
|
diff_field[ei] = 0.0
|
|
dest_map[ei] = element_node_indexes[ei].index(dindex)
|
|
(ex, ey) = calculateValues(diff_field, eecouple, sx, sy)
|
|
for i in range(nelem):
|
|
ddata_x[i][dest_map[i]] = ex[i]
|
|
ddata_y[i][dest_map[i]] = ey[i]
|
|
# these are the node indexes w.r.t. the first edge on the triangle
|
|
# ddata_x[i][dnode] = ex[i]
|
|
# ddata_y[i][dnode] = ey[i]
|
|
|
|
# print element_node_indexes
|
|
|
|
for i in range(3):
|
|
print("index %d" % i)
|
|
for j in range(nelem):
|
|
print(
|
|
"%g\t%g\t%g\t%g"
|
|
% (ddata_x[j][i][0], ddata_y[j][i][0], cdata_x[j][i], cdata_y[j][i])
|
|
)
|
|
|
|
# print element_node_indexes
|
|
# print potential_diff
|
|
# print get_edge_model_values(device=device, region=region, name="ElectricField:Potential@n0")
|
|
# print get_edge_model_values(device=device, region=region, name="ElectricField:Potential@n1")
|