1 /**********
2 Copyright 1991 Regents of the University of California. All rights reserved.
3 Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
4 Author: 1991 David A. Gates, U. C. Berkeley CAD Group
5 **********/
6
7 #include "ngspice/ngspice.h"
8 #include "ngspice/numglobs.h"
9 #include "ngspice/numenum.h"
10 #include "ngspice/twomesh.h"
11 #include "ngspice/twodev.h"
12 #include "ngspice/spmatrix.h"
13 #include "twoddefs.h"
14 #include "twodext.h"
15
16
17 /* functions to setup and solve the 2D poisson equation */
18
19 void
TWOQjacBuild(TWOdevice * pDevice)20 TWOQjacBuild(TWOdevice *pDevice)
21 {
22 SMPmatrix *matrix = pDevice->matrix;
23 TWOelem *pElem;
24 TWOnode *pNode, *pNode1;
25 int eIndex, nIndex;
26
27 /* set up matrix pointers */
28 /* establish main diagonal first */
29 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
30 pElem = pDevice->elements[ eIndex ];
31 for ( nIndex = 0; nIndex <= 3; nIndex++ ) {
32 if ( pElem->evalNodes[ nIndex ] ) {
33 pNode = pElem->pNodes[ nIndex ];
34 pNode->fPsiPsi = spGetElement( matrix, pNode->poiEqn, pNode->poiEqn );
35 }
36 }
37 }
38 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
39 pElem = pDevice->elements[ eIndex ];
40
41 pNode = pElem->pTLNode;
42 pNode1 = pElem->pTRNode;
43 pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
44 pNode1 = pElem->pBLNode;
45 pNode->fPsiPsijP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
46 pNode = pElem->pTRNode;
47 pNode1 = pElem->pTLNode;
48 pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
49 pNode1 = pElem->pBRNode;
50 pNode->fPsiPsijP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
51 pNode = pElem->pBRNode;
52 pNode1 = pElem->pBLNode;
53 pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
54 pNode1 = pElem->pTRNode;
55 pNode->fPsiPsijM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
56 pNode = pElem->pBLNode;
57 pNode1 = pElem->pBRNode;
58 pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
59 pNode1 = pElem->pTLNode;
60 pNode->fPsiPsijM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
61 }
62 /*
63 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
64 pElem = pDevice->elements[ eIndex ];
65
66 pNode = pElem->pTLNode;
67 pNode->fPsiPsi = spGetElement( matrix, pNode->poiEqn, pNode->poiEqn );
68 pNode1 = pElem->pTRNode;
69 pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
70 pNode1 = pElem->pBLNode;
71 pNode->fPsiPsijP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
72 pNode = pElem->pTRNode;
73 pNode->fPsiPsi = spGetElement( matrix, pNode->poiEqn, pNode->poiEqn );
74 pNode1 = pElem->pTLNode;
75 pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
76 pNode1 = pElem->pBRNode;
77 pNode->fPsiPsijP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
78 pNode = pElem->pBRNode;
79 pNode->fPsiPsi = spGetElement( matrix, pNode->poiEqn, pNode->poiEqn );
80 pNode1 = pElem->pBLNode;
81 pNode->fPsiPsiiM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
82 pNode1 = pElem->pTRNode;
83 pNode->fPsiPsijM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
84 pNode = pElem->pBLNode;
85 pNode->fPsiPsi = spGetElement( matrix, pNode->poiEqn, pNode->poiEqn );
86 pNode1 = pElem->pBRNode;
87 pNode->fPsiPsiiP1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
88 pNode1 = pElem->pTLNode;
89 pNode->fPsiPsijM1 = spGetElement(matrix, pNode->poiEqn, pNode1->poiEqn );
90 }
91 */
92 }
93
94 void
TWOQsysLoad(TWOdevice * pDevice)95 TWOQsysLoad(TWOdevice *pDevice)
96 {
97 TWOelem *pElem;
98 TWOnode *pNode;
99 TWOedge *pHEdge, *pVEdge;
100 int index, eIndex;
101 double *pRhs = pDevice->rhs;
102 double dyOverDx, dxOverDy, dPsiT, dPsiB, dPsiL, dPsiR;
103
104 TWOQcommonTerms( pDevice );
105
106 /* zero the rhs vector */
107 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) {
108 pRhs[ index ] = 0.0;
109 }
110
111 /* zero the matrix */
112 spClear( pDevice->matrix );
113
114 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
115 pElem = pDevice->elements[ eIndex ];
116
117 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy;
118 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx;
119 dPsiT = pElem->pTopEdge->dPsi;
120 dPsiB = pElem->pBotEdge->dPsi;
121 dPsiL = pElem->pLeftEdge->dPsi;
122 dPsiR = pElem->pRightEdge->dPsi;
123
124 /* load for all i,j */
125 for ( index = 0; index <= 3; index++ ) {
126 pNode = pElem->pNodes[ index ];
127 if ( pNode->nodeType != CONTACT ) {
128 *(pNode->fPsiPsi) += dyOverDx + dxOverDy;
129 if ( index <= 1 ) {
130 pHEdge = pElem->pTopEdge;
131 } else {
132 pHEdge = pElem->pBotEdge;
133 }
134 if ( index == 0 || index == 3 ) {
135 pVEdge = pElem->pLeftEdge;
136 } else {
137 pVEdge = pElem->pRightEdge;
138 }
139 /* add surface state charges */
140 pRhs[ pNode->poiEqn ] += 0.5 * pElem->dx * pHEdge->qf;
141 pRhs[ pNode->poiEqn ] += 0.5 * pElem->dy * pVEdge->qf;
142 if ( pElem->elemType == SEMICON ) {
143 *(pNode->fPsiPsi) += 0.25 * pElem->dx * pElem->dy *
144 (pNode->nConc + pNode->pConc);
145 pRhs[ pNode->poiEqn ] += 0.25 * pElem->dx * pElem->dy *
146 (pNode->netConc + pNode->pConc - pNode->nConc);
147 }
148 }
149 }
150
151 pNode = pElem->pTLNode;
152 pRhs[ pNode->poiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL;
153 *(pNode->fPsiPsiiP1) -= dyOverDx;
154 *(pNode->fPsiPsijP1) -= dxOverDy;
155
156 pNode = pElem->pTRNode;
157 pRhs[ pNode->poiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR;
158 *(pNode->fPsiPsiiM1) -= dyOverDx;
159 *(pNode->fPsiPsijP1) -= dxOverDy;
160
161 pNode = pElem->pBRNode;
162 pRhs[ pNode->poiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR;
163 *(pNode->fPsiPsiiM1) -= dyOverDx;
164 *(pNode->fPsiPsijM1) -= dxOverDy;
165
166 pNode = pElem->pBLNode;
167 pRhs[ pNode->poiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL;
168 *(pNode->fPsiPsiiP1) -= dyOverDx;
169 *(pNode->fPsiPsijM1) -= dxOverDy;
170 }
171 }
172
173
174 void
TWOQrhsLoad(TWOdevice * pDevice)175 TWOQrhsLoad(TWOdevice *pDevice)
176 {
177 TWOelem *pElem;
178 TWOnode *pNode;
179 TWOedge *pHEdge, *pVEdge;
180 int index, eIndex;
181 double *pRhs = pDevice->rhs;
182 double dyOverDx, dxOverDy, dPsiT, dPsiB, dPsiL, dPsiR;
183
184 TWOQcommonTerms( pDevice );
185
186 /* zero the rhs vector */
187 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) {
188 pRhs[ index ] = 0.0;
189 }
190
191 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
192 pElem = pDevice->elements[ eIndex ];
193
194 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy;
195 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx;
196 dPsiT = pElem->pTopEdge->dPsi;
197 dPsiB = pElem->pBotEdge->dPsi;
198 dPsiL = pElem->pLeftEdge->dPsi;
199 dPsiR = pElem->pRightEdge->dPsi;
200
201 /* load nodal terms */
202 for ( index = 0; index <= 3; index++ ) {
203 pNode = pElem->pNodes[ index ];
204 if ( (pNode->nodeType != CONTACT) && (pElem->elemType == SEMICON) ) {
205 pRhs[ pNode->poiEqn ] += 0.25 * pElem->dx * pElem->dy *
206 (pNode->netConc + pNode->pConc - pNode->nConc);
207 }
208 if ( index <= 1 ) {
209 pHEdge = pElem->pTopEdge;
210 } else {
211 pHEdge = pElem->pBotEdge;
212 }
213 if ( index == 0 || index == 3 ) {
214 pVEdge = pElem->pLeftEdge;
215 } else {
216 pVEdge = pElem->pRightEdge;
217 }
218 /* add surface state charges */
219 pRhs[ pNode->poiEqn ] += 0.5 * pElem->dx * pHEdge->qf;
220 pRhs[ pNode->poiEqn ] += 0.5 * pElem->dy * pVEdge->qf;
221 }
222
223 /* load edge terms */
224 pNode = pElem->pTLNode;
225 pRhs[ pNode->poiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL;
226
227 pNode = pElem->pTRNode;
228 pRhs[ pNode->poiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR;
229
230 pNode = pElem->pBRNode;
231 pRhs[ pNode->poiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR;
232
233 pNode = pElem->pBLNode;
234 pRhs[ pNode->poiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL;
235 }
236 }
237
238 void
TWOQcommonTerms(TWOdevice * pDevice)239 TWOQcommonTerms(TWOdevice *pDevice)
240 {
241 TWOelem *pElem;
242 TWOedge *pEdge;
243 TWOnode *pNode;
244 int index, eIndex;
245 double psi1, psi2, refPsi;
246
247 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
248 pElem = pDevice->elements[ eIndex ];
249 refPsi = pElem->matlInfo->refPsi;
250 for ( index = 0; index <= 3; index++ ) {
251 if ( pElem->evalNodes[ index ] ) {
252 pNode = pElem->pNodes[ index ];
253 if ( pNode->nodeType != CONTACT ) {
254 pNode->psi = pDevice->dcSolution[ pNode->poiEqn ];
255 if ( pElem->elemType == SEMICON ) {
256 pNode->nConc = pNode->nie * exp( pNode->psi - refPsi );
257 pNode->pConc = pNode->nie * exp( - pNode->psi + refPsi );
258 }
259 }
260 }
261 if ( pElem->evalEdges[ index ] ) {
262 pEdge = pElem->pEdges[ index ];
263 pNode = pElem->pNodes[ index ];
264 if ( pNode->nodeType != CONTACT ) {
265 psi1 = pDevice->dcSolution[pNode->poiEqn];
266 } else {
267 psi1 = pNode->psi;
268 }
269 pNode = pElem->pNodes[ (index + 1) % 4 ];
270 if ( pNode->nodeType != CONTACT ) {
271 psi2 = pDevice->dcSolution[pNode->poiEqn];
272 } else {
273 psi2 = pNode->psi;
274 }
275 if ( index <= 1 ) {
276 pEdge->dPsi = psi2 - psi1;
277 } else {
278 pEdge->dPsi = psi1 - psi2;
279 }
280 }
281 }
282 }
283 }
284