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/bool.h" 13 #include "ngspice/spmatrix.h" 14 #include "twoddefs.h" 15 #include "twodext.h" 16 #include "ngspice/cidersupt.h" 17 #include "../../maths/misc/bernoull.h" 18 19 20 /* 21 * Functions to setup and solve the continuity equations. 22 * Both continuity equations are solved. 23 * Separate functions are used for one continuity equation. 24 */ 25 26 27 /* 28 * Setup matrix pointers to Jacobian entries and 29 * store direct pointers with the nodes. 30 */ 31 32 void TWONjacBuild(TWOdevice * pDevice)33 TWONjacBuild(TWOdevice *pDevice) 34 { 35 SMPmatrix *matrix = pDevice->matrix; 36 TWOelem *pElem; 37 TWOnode *pNode; 38 TWOchannel *pCh; 39 int eIndex, nIndex; 40 int nextIndex; /* index of node to find next element */ 41 int psiEqn, nEqn; /* scratch for deref'd eqn numbers */ 42 int psiEqnTL = 0, nEqnTL = 0; 43 int psiEqnTR = 0, nEqnTR = 0; 44 int psiEqnBR = 0, nEqnBR = 0; 45 int psiEqnBL = 0, nEqnBL = 0; 46 int psiEqnInM = 0, psiEqnInP = 0; /* scratch for deref'd surface eqns */ 47 int psiEqnOxM = 0, psiEqnOxP = 0; /* M= more negative, P= more positive */ 48 49 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 50 pElem = pDevice->elements[ eIndex ]; 51 52 /* first the self terms */ 53 for ( nIndex = 0; nIndex <= 3; nIndex++ ) { 54 pNode = pElem->pNodes[ nIndex ]; 55 /* get poisson-only pointer */ 56 psiEqn = pNode->psiEqn; 57 pNode->fPsiPsi = spGetElement( matrix, psiEqn, psiEqn ); 58 59 if ( pElem->elemType == SEMICON ) { 60 /* get continuity-coupling terms */ 61 nEqn = pNode->nEqn; 62 pNode->pEqn = 0; /* Throw pEqn number into garbage. */ 63 /* pointers for additional terms */ 64 pNode->fPsiN = spGetElement( matrix, psiEqn, nEqn ); 65 pNode->fNPsi = spGetElement( matrix, nEqn, psiEqn ); 66 pNode->fNN = spGetElement( matrix, nEqn, nEqn ); 67 } else { 68 nEqn = 0; 69 } 70 /* save equation indices */ 71 switch ( nIndex ) { 72 case 0: /* TL Node */ 73 psiEqnTL = psiEqn; 74 nEqnTL = nEqn; 75 break; 76 case 1: /* TR Node */ 77 psiEqnTR = psiEqn; 78 nEqnTR = nEqn; 79 break; 80 case 2: /* BR Node */ 81 psiEqnBR = psiEqn; 82 nEqnBR = nEqn; 83 break; 84 case 3: /* BL Node */ 85 psiEqnBL = psiEqn; 86 nEqnBL = nEqn; 87 break; 88 default: 89 break; 90 } 91 } 92 93 /* now terms to couple to adjacent nodes */ 94 pNode = pElem->pTLNode; 95 pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnTL, psiEqnTR ); 96 pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTL, psiEqnBL ); 97 if ( pElem->elemType == SEMICON ) { 98 /* continuity equation pointers */ 99 pNode->fNPsiiP1 = spGetElement( matrix, nEqnTL, psiEqnTR ); 100 pNode->fNNiP1 = spGetElement( matrix, nEqnTL, nEqnTR ); 101 pNode->fNPsijP1 = spGetElement( matrix, nEqnTL, psiEqnBL ); 102 pNode->fNNjP1 = spGetElement( matrix, nEqnTL, nEqnBL ); 103 /* Surface Mobility Model depends on diagonal node values */ 104 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 105 pNode->fNPsiiP1jP1 = spGetElement( matrix, nEqnTL, psiEqnBR ); 106 pNode->fNNiP1jP1 = spGetElement( matrix, nEqnTL, nEqnBR ); 107 } 108 } 109 110 pNode = pElem->pTRNode; 111 pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnTR, psiEqnTL ); 112 pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTR, psiEqnBR ); 113 if ( pElem->elemType == SEMICON ) { 114 /* continuity equation pointers */ 115 pNode->fNPsiiM1 = spGetElement( matrix, nEqnTR, psiEqnTL ); 116 pNode->fNNiM1 = spGetElement( matrix, nEqnTR, nEqnTL ); 117 pNode->fNPsijP1 = spGetElement( matrix, nEqnTR, psiEqnBR ); 118 pNode->fNNjP1 = spGetElement( matrix, nEqnTR, nEqnBR ); 119 /* Surface Mobility Model depends on diagonal node values */ 120 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 121 pNode->fNPsiiM1jP1 = spGetElement( matrix, nEqnTR, psiEqnBL ); 122 pNode->fNNiM1jP1 = spGetElement( matrix, nEqnTR, nEqnBL ); 123 } 124 } 125 126 pNode = pElem->pBRNode; 127 pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnBR, psiEqnBL ); 128 pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBR, psiEqnTR ); 129 if ( pElem->elemType == SEMICON ) { 130 /* continuity equation pointers */ 131 pNode->fNPsiiM1 = spGetElement( matrix, nEqnBR, psiEqnBL ); 132 pNode->fNNiM1 = spGetElement( matrix, nEqnBR, nEqnBL ); 133 pNode->fNPsijM1 = spGetElement( matrix, nEqnBR, psiEqnTR ); 134 pNode->fNNjM1 = spGetElement( matrix, nEqnBR, nEqnTR ); 135 /* Surface Mobility Model depends on diagonal node values */ 136 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 137 pNode->fNPsiiM1jM1 = spGetElement( matrix, nEqnBR, psiEqnTL ); 138 pNode->fNNiM1jM1 = spGetElement( matrix, nEqnBR, nEqnTL ); 139 } 140 } 141 142 pNode = pElem->pBLNode; 143 pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnBL, psiEqnBR ); 144 pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBL, psiEqnTL ); 145 if ( pElem->elemType == SEMICON ) { 146 /* continuity equation pointers */ 147 pNode->fNPsiiP1 = spGetElement( matrix, nEqnBL, psiEqnBR ); 148 pNode->fNNiP1 = spGetElement( matrix, nEqnBL, nEqnBR ); 149 pNode->fNPsijM1 = spGetElement( matrix, nEqnBL, psiEqnTL ); 150 pNode->fNNjM1 = spGetElement( matrix, nEqnBL, nEqnTL ); 151 /* Surface Mobility Model depends on diagonal node values */ 152 if ( MobDeriv && SurfaceMobility && pElem->channel ) { 153 pNode->fNPsiiP1jM1 = spGetElement( matrix, nEqnBL, psiEqnTR ); 154 pNode->fNNiP1jM1 = spGetElement( matrix, nEqnBL, nEqnTR ); 155 } 156 } 157 } 158 /* 159 * Add terms for surface-field of inversion-layer mobility model. 160 * Elements MUST be made from silicon for this to work. 161 * No empty elements are allowed. 162 * Don't need these pointers if SurfaceMobility isn't set. 163 */ 164 if ( MobDeriv && SurfaceMobility ) { 165 for ( pCh = pDevice->pChannel; pCh != NULL; 166 pCh = pCh->next ) { 167 pElem = pCh->pNElem; 168 switch (pCh->type) { 169 case 0: 170 psiEqnInM = pElem->pBLNode->psiEqn; 171 psiEqnInP = pElem->pBRNode->psiEqn; 172 psiEqnOxM = pElem->pTLNode->psiEqn; 173 psiEqnOxP = pElem->pTRNode->psiEqn; 174 break; 175 case 1: 176 psiEqnInM = pElem->pTLNode->psiEqn; 177 psiEqnInP = pElem->pBLNode->psiEqn; 178 psiEqnOxM = pElem->pTRNode->psiEqn; 179 psiEqnOxP = pElem->pBRNode->psiEqn; 180 break; 181 case 2: 182 psiEqnInM = pElem->pTLNode->psiEqn; 183 psiEqnInP = pElem->pTRNode->psiEqn; 184 psiEqnOxM = pElem->pBLNode->psiEqn; 185 psiEqnOxP = pElem->pBRNode->psiEqn; 186 break; 187 case 3: 188 psiEqnInM = pElem->pTRNode->psiEqn; 189 psiEqnInP = pElem->pBRNode->psiEqn; 190 psiEqnOxM = pElem->pTLNode->psiEqn; 191 psiEqnOxP = pElem->pBLNode->psiEqn; 192 break; 193 } 194 pElem = pCh->pSeed; 195 nextIndex = (pCh->type + 2)%4; 196 while (pElem && pElem->channel == pCh->id) { 197 for ( nIndex = 0; nIndex <= 3; nIndex++ ) { 198 pNode = pElem->pNodes[ nIndex ]; 199 psiEqn = pNode->psiEqn; 200 nEqn = pNode->nEqn; 201 if ( pCh->type % 2 == 0 ) { /* Vertical Slice */ 202 if ( nIndex == 0 || nIndex == 3 ) { /* Left Side */ 203 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInM ); 204 pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP ); 205 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxM ); 206 pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP ); 207 } else { /* Right Side */ 208 pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM ); 209 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInP ); 210 pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM ); 211 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxP ); 212 } 213 } else { /* Horizontal Slice */ 214 if ( nIndex == 0 || nIndex == 3 ) { /* Left (Top?) Side : bug 483 */ 215 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInM ); 216 pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP ); 217 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxM ); 218 pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP ); 219 } else { /* Bottom Side */ 220 pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM ); 221 pNode->fNPsiIn = spGetElement( matrix, nEqn, psiEqnInP ); 222 pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM ); 223 pNode->fNPsiOx = spGetElement( matrix, nEqn, psiEqnOxP ); 224 } 225 } 226 } /* endfor nIndex */ 227 pElem = pElem->pElems[ nextIndex ]; 228 } /* endwhile pElem */ 229 } /* endfor pCh */ 230 } /* endif SurfaceMobility */ 231 } 232 233 234 /* 235 * The Jacobian and Rhs are loaded by the following function. 236 * Inputs are the transient analysis flag and the transient 237 * information structure 238 */ 239 240 void TWONsysLoad(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)241 TWONsysLoad(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info) 242 { 243 TWOelem *pElem; 244 TWOnode *pNode; 245 TWOedge *pHEdge, *pVEdge; 246 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 247 TWOchannel *pCh; 248 int index, eIndex; 249 int nextIndex; /* index of node to find next element */ 250 double *pRhs = pDevice->rhs; 251 double dx, dy, dxdy, dyOverDx, dxOverDy; 252 double ds; 253 double dPsiT, dPsiB, dPsiL, dPsiR; 254 double rhsN; 255 double nConc, pConc; 256 double perTime = 0.0; 257 258 /* first compute the currents and derivatives */ 259 TWONcommonTerms( pDevice, FALSE, tranAnalysis, info ); 260 261 /* find reciprocal timestep */ 262 if ( tranAnalysis ) { 263 perTime = info->intCoeff[0]; 264 } 265 266 /* zero the rhs vector */ 267 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) { 268 pRhs[ index ] = 0.0; 269 } 270 271 /* zero the matrix */ 272 spClear( pDevice->matrix ); 273 274 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 275 pElem = pDevice->elements[ eIndex ]; 276 277 dx = 0.5 * pElem->dx; 278 dy = 0.5 * pElem->dy; 279 dxdy = dx * dy; 280 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 281 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 282 283 pTEdge = pElem->pTopEdge; 284 pBEdge = pElem->pBotEdge; 285 pLEdge = pElem->pLeftEdge; 286 pREdge = pElem->pRightEdge; 287 dPsiT = pTEdge->dPsi; 288 dPsiB = pBEdge->dPsi; 289 dPsiL = pLEdge->dPsi; 290 dPsiR = pREdge->dPsi; 291 292 /* load for all i,j */ 293 for ( index = 0; index <= 3; index++ ) { 294 pNode = pElem->pNodes[ index ]; 295 if ( pNode->nodeType != CONTACT ) { 296 *(pNode->fPsiPsi) += dyOverDx + dxOverDy; 297 if ( index <= 1 ) { 298 pHEdge = pTEdge; 299 } else { 300 pHEdge = pBEdge; 301 } 302 if ( index == 0 || index == 3 ) { 303 pVEdge = pLEdge; 304 } else { 305 pVEdge = pREdge; 306 } 307 /* Add surface state charges. */ 308 pRhs[ pNode->psiEqn ] += dx * pHEdge->qf; 309 pRhs[ pNode->psiEqn ] += dy * pVEdge->qf; 310 if ( pElem->elemType == SEMICON ) { 311 nConc = pDevice->devState0 [pNode->nodeN]; 312 pConc = pDevice->devState0 [pNode->nodeP]; 313 314 *(pNode->fPsiN) += dxdy; 315 *(pNode->fPsiPsi) += dxdy * pConc; 316 *(pNode->fNPsi) -= dy * pHEdge->dJnDpsiP1 + dx * pVEdge->dJnDpsiP1; 317 pRhs[ pNode->psiEqn ] += dxdy * (pNode->netConc + pConc - nConc); 318 319 /* Handle generation terms */ 320 *(pNode->fNN) -= dxdy * pNode->dUdN; 321 *(pNode->fNPsi) += dxdy * pNode->dUdP * pConc; 322 rhsN = - dxdy * pNode->uNet; 323 pRhs[ pNode->nEqn ] -= rhsN; 324 325 /* Handle dXdT continuity terms */ 326 if ( tranAnalysis ) { 327 *(pNode->fNN) -= dxdy * perTime; 328 pRhs[ pNode->nEqn ] += dxdy * pNode->dNdT; 329 } 330 } 331 } 332 } 333 334 /* Handle neighbor and edge dependent terms */ 335 pNode = pElem->pTLNode; 336 if ( pNode->nodeType != CONTACT ) { 337 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL; 338 *(pNode->fPsiPsiiP1) -= dyOverDx; 339 *(pNode->fPsiPsijP1) -= dxOverDy; 340 if ( pElem->elemType == SEMICON ) { 341 pRhs[ pNode->nEqn ] -= dy * pTEdge->jn + dx * pLEdge->jn; 342 *(pNode->fNN) += dy * pTEdge->dJnDn + dx * pLEdge->dJnDn; 343 *(pNode->fNPsiiP1) += dy * pTEdge->dJnDpsiP1; 344 *(pNode->fNNiP1) += dy * pTEdge->dJnDnP1; 345 *(pNode->fNPsijP1) += dx * pLEdge->dJnDpsiP1; 346 *(pNode->fNNjP1) += dx * pLEdge->dJnDnP1; 347 } 348 } 349 pNode = pElem->pTRNode; 350 if ( pNode->nodeType != CONTACT ) { 351 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR; 352 *(pNode->fPsiPsiiM1) -= dyOverDx; 353 *(pNode->fPsiPsijP1) -= dxOverDy; 354 if ( pElem->elemType == SEMICON ) { 355 pRhs[ pNode->nEqn ] -= -dy * pTEdge->jn + dx * pREdge->jn; 356 *(pNode->fNN) += -dy * pTEdge->dJnDnP1 + dx * pREdge->dJnDn; 357 *(pNode->fNPsiiM1) += dy * pTEdge->dJnDpsiP1; 358 *(pNode->fNNiM1) -= dy * pTEdge->dJnDn; 359 *(pNode->fNPsijP1) += dx * pREdge->dJnDpsiP1; 360 *(pNode->fNNjP1) += dx * pREdge->dJnDnP1; 361 } 362 } 363 pNode = pElem->pBRNode; 364 if ( pNode->nodeType != CONTACT ) { 365 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR; 366 *(pNode->fPsiPsiiM1) -= dyOverDx; 367 *(pNode->fPsiPsijM1) -= dxOverDy; 368 if ( pElem->elemType == SEMICON ) { 369 pRhs[ pNode->nEqn ] -= -dy * pBEdge->jn - dx * pREdge->jn; 370 *(pNode->fNN) += -dy * pBEdge->dJnDnP1 - dx * pREdge->dJnDnP1; 371 *(pNode->fNPsiiM1) += dy * pBEdge->dJnDpsiP1; 372 *(pNode->fNNiM1) -= dy * pBEdge->dJnDn; 373 *(pNode->fNPsijM1) += dx * pREdge->dJnDpsiP1; 374 *(pNode->fNNjM1) -= dx * pREdge->dJnDn; 375 } 376 } 377 pNode = pElem->pBLNode; 378 if ( pNode->nodeType != CONTACT ) { 379 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL; 380 *(pNode->fPsiPsiiP1) -= dyOverDx; 381 *(pNode->fPsiPsijM1) -= dxOverDy; 382 if ( pElem->elemType == SEMICON ) { 383 pRhs[ pNode->nEqn ] -= dy * pBEdge->jn - dx * pLEdge->jn; 384 *(pNode->fNN) += dy * pBEdge->dJnDn - dx * pLEdge->dJnDnP1; 385 *(pNode->fNPsiiP1) += dy * pBEdge->dJnDpsiP1; 386 *(pNode->fNNiP1) += dy * pBEdge->dJnDnP1; 387 *(pNode->fNPsijM1) += dx * pLEdge->dJnDpsiP1; 388 *(pNode->fNNjM1) -= dx * pLEdge->dJnDn; 389 } 390 } 391 } 392 393 /* Calculate the Inversion-Layer Mobility Dependent Terms in Jac. */ 394 if ( MobDeriv && SurfaceMobility ) { 395 for ( pCh = pDevice->pChannel; pCh != NULL; 396 pCh = pCh->next ) { 397 /* Find effective height of oxide element at interface. */ 398 if ( pCh->type%2 == 0 ) { /* Vertical slice */ 399 ds = pCh->pNElem->dy / pCh->pNElem->epsRel; 400 } else { /* Horizontal slice */ 401 ds = pCh->pNElem->dx / pCh->pNElem->epsRel; 402 } 403 pElem = pCh->pSeed; 404 nextIndex = (pCh->type + 2)%4; 405 while (pElem && pElem->channel == pCh->id) { 406 TWONmobDeriv( pElem, pCh->type, ds ); 407 pElem = pElem->pElems[ nextIndex ]; 408 } 409 } /* endfor pCh != NULL */ 410 } /* endif MobDeriv and SurfaceMobility */ 411 } 412 413 414 /* 415 * This function used only for direct method ac analysis. 416 * Used to load only the dc Jacobian matrix. Rhs is unaffected 417 */ 418 419 void TWONjacLoad(TWOdevice * pDevice)420 TWONjacLoad(TWOdevice *pDevice) 421 { 422 TWOelem *pElem; 423 TWOnode *pNode; 424 TWOedge *pHEdge, *pVEdge; 425 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 426 TWOchannel *pCh; 427 int index, eIndex; 428 int nextIndex; /* index of node to find next element */ 429 double dx, dy, dxdy, dyOverDx, dxOverDy; 430 double ds; 431 double pConc; 432 433 /* first compute the currents and derivatives */ 434 TWONcommonTerms( pDevice, FALSE, FALSE, NULL ); 435 436 /* zero the matrix */ 437 spClear( pDevice->matrix ); 438 439 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 440 pElem = pDevice->elements[ eIndex ]; 441 dx = 0.5 * pElem->dx; 442 dy = 0.5 * pElem->dy; 443 dxdy = dx * dy; 444 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 445 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 446 447 pTEdge = pElem->pTopEdge; 448 pBEdge = pElem->pBotEdge; 449 pLEdge = pElem->pLeftEdge; 450 pREdge = pElem->pRightEdge; 451 452 /* load for all i,j */ 453 for ( index = 0; index <= 3; index++ ) { 454 pNode = pElem->pNodes[ index ]; 455 if ( pNode->nodeType != CONTACT ) { 456 *(pNode->fPsiPsi) += dyOverDx + dxOverDy; 457 if ( pElem->elemType == SEMICON ) { 458 if ( index <= 1 ) { 459 pHEdge = pTEdge; 460 } else { 461 pHEdge = pBEdge; 462 } 463 if ( index == 0 || index == 3 ) { 464 pVEdge = pLEdge; 465 } else { 466 pVEdge = pREdge; 467 } 468 pConc = pDevice->devState0 [pNode->nodeP]; 469 *(pNode->fPsiN) += dxdy; 470 *(pNode->fPsiPsi) += dxdy * pConc; 471 *(pNode->fNPsi) -= dy * pHEdge->dJnDpsiP1 + dx * pVEdge->dJnDpsiP1; 472 473 /* Handle generation terms */ 474 *(pNode->fNN) -= dxdy * pNode->dUdN; 475 *(pNode->fNPsi) += dxdy * pNode->dUdP * pConc; 476 } 477 } 478 } 479 480 /* Handle neighbor and edge dependent terms */ 481 pNode = pElem->pTLNode; 482 if ( pNode->nodeType != CONTACT ) { 483 *(pNode->fPsiPsiiP1) -= dyOverDx; 484 *(pNode->fPsiPsijP1) -= dxOverDy; 485 if ( pElem->elemType == SEMICON ) { 486 *(pNode->fNN) += dy * pTEdge->dJnDn + dx * pLEdge->dJnDn; 487 *(pNode->fNPsiiP1) += dy * pTEdge->dJnDpsiP1; 488 *(pNode->fNNiP1) += dy * pTEdge->dJnDnP1; 489 *(pNode->fNPsijP1) += dx * pLEdge->dJnDpsiP1; 490 *(pNode->fNNjP1) += dx * pLEdge->dJnDnP1; 491 } 492 } 493 pNode = pElem->pTRNode; 494 if ( pNode->nodeType != CONTACT ) { 495 *(pNode->fPsiPsiiM1) -= dyOverDx; 496 *(pNode->fPsiPsijP1) -= dxOverDy; 497 if ( pElem->elemType == SEMICON ) { 498 *(pNode->fNN) += -dy * pTEdge->dJnDnP1 + dx * pREdge->dJnDn; 499 *(pNode->fNPsiiM1) += dy * pTEdge->dJnDpsiP1; 500 *(pNode->fNNiM1) -= dy * pTEdge->dJnDn; 501 *(pNode->fNPsijP1) += dx * pREdge->dJnDpsiP1; 502 *(pNode->fNNjP1) += dx * pREdge->dJnDnP1; 503 } 504 } 505 pNode = pElem->pBRNode; 506 if ( pNode->nodeType != CONTACT ) { 507 *(pNode->fPsiPsiiM1) -= dyOverDx; 508 *(pNode->fPsiPsijM1) -= dxOverDy; 509 if ( pElem->elemType == SEMICON ) { 510 *(pNode->fNN) += -dy * pBEdge->dJnDnP1 - dx * pREdge->dJnDnP1; 511 *(pNode->fNPsiiM1) += dy * pBEdge->dJnDpsiP1; 512 *(pNode->fNNiM1) -= dy * pBEdge->dJnDn; 513 *(pNode->fNPsijM1) += dx * pREdge->dJnDpsiP1; 514 *(pNode->fNNjM1) -= dx * pREdge->dJnDn; 515 } 516 } 517 pNode = pElem->pBLNode; 518 if ( pNode->nodeType != CONTACT ) { 519 *(pNode->fPsiPsiiP1) -= dyOverDx; 520 *(pNode->fPsiPsijM1) -= dxOverDy; 521 if ( pElem->elemType == SEMICON ) { 522 *(pNode->fNN) += dy * pBEdge->dJnDn - dx * pLEdge->dJnDnP1; 523 *(pNode->fNPsiiP1) += dy * pBEdge->dJnDpsiP1; 524 *(pNode->fNNiP1) += dy * pBEdge->dJnDnP1; 525 *(pNode->fNPsijM1) += dx * pLEdge->dJnDpsiP1; 526 *(pNode->fNNjM1) -= dx * pLEdge->dJnDn; 527 } 528 } 529 } 530 531 /* Calculate the Inversion-Layer Mobility Dependent Terms in Jac. */ 532 if ( MobDeriv && SurfaceMobility ) { 533 for ( pCh = pDevice->pChannel; pCh != NULL; 534 pCh = pCh->next ) { 535 /* Find effective height of oxide element at interface. */ 536 if ( pCh->type%2 == 0 ) { /* Vertical slice */ 537 ds = pCh->pNElem->dy / pCh->pNElem->epsRel; 538 } else { /* Horizontal slice */ 539 ds = pCh->pNElem->dx / pCh->pNElem->epsRel; 540 } 541 pElem = pCh->pSeed; 542 nextIndex = (pCh->type + 2)%4; 543 while (pElem && pElem->channel == pCh->id) { 544 TWONmobDeriv( pElem, pCh->type, ds ); 545 pElem = pElem->pElems[ nextIndex ]; 546 } 547 } /* endfor pCh != NULL */ 548 } /* endif MobDeriv and SurfaceMobility */ 549 } 550 551 /* load only the Rhs vector */ 552 void TWONrhsLoad(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)553 TWONrhsLoad(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info) 554 { 555 TWOelem *pElem; 556 TWOnode *pNode; 557 TWOedge *pHEdge, *pVEdge; 558 TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge; 559 int index, eIndex; 560 double *pRhs = pDevice->rhs; 561 double dx, dy, dxdy, dyOverDx, dxOverDy; 562 double dPsiT, dPsiB, dPsiL, dPsiR; 563 double rhsN; 564 double nConc, pConc; 565 double perTime; 566 567 /* first compute the currents */ 568 TWONcommonTerms( pDevice, TRUE, tranAnalysis, info ); 569 570 /* find reciprocal timestep */ 571 if ( tranAnalysis ) { 572 perTime = info->intCoeff[0]; 573 } 574 575 /* zero the rhs vector */ 576 for ( index = 1 ; index <= pDevice->numEqns ; index++ ) { 577 pRhs[ index ] = 0.0; 578 } 579 580 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 581 pElem = pDevice->elements[ eIndex ]; 582 583 dx = 0.5 * pElem->dx; 584 dy = 0.5 * pElem->dy; 585 dxdy = dx * dy; 586 dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy; 587 dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx; 588 589 pTEdge = pElem->pTopEdge; 590 pBEdge = pElem->pBotEdge; 591 pLEdge = pElem->pLeftEdge; 592 pREdge = pElem->pRightEdge; 593 dPsiT = pTEdge->dPsi; 594 dPsiB = pBEdge->dPsi; 595 dPsiL = pLEdge->dPsi; 596 dPsiR = pREdge->dPsi; 597 598 /* load for all i,j */ 599 for ( index = 0; index <= 3; index++ ) { 600 pNode = pElem->pNodes[ index ]; 601 if ( pNode->nodeType != CONTACT ) { 602 if ( index <= 1 ) { 603 pHEdge = pTEdge; 604 } else { 605 pHEdge = pBEdge; 606 } 607 if ( index == 0 || index == 3 ) { 608 pVEdge = pLEdge; 609 } else { 610 pVEdge = pREdge; 611 } 612 /* Add surface state charges. */ 613 pRhs[ pNode->psiEqn ] += dx * pHEdge->qf; 614 pRhs[ pNode->psiEqn ] += dy * pVEdge->qf; 615 if ( pElem->elemType == SEMICON ) { 616 nConc = pDevice->devState0 [pNode->nodeN]; 617 pConc = pDevice->devState0 [pNode->nodeP]; 618 pRhs[ pNode->psiEqn ] += dxdy * (pNode->netConc + pConc - nConc); 619 620 /* Handle generation terms */ 621 rhsN = - dxdy * pNode->uNet; 622 pRhs[ pNode->nEqn ] -= rhsN; 623 624 /* Handle dXdT continuity terms */ 625 if ( tranAnalysis ) { 626 pRhs[ pNode->nEqn ] += dxdy * pNode->dNdT; 627 } 628 } 629 } 630 } 631 632 /* Handle neighbor and edge dependent terms */ 633 pNode = pElem->pTLNode; 634 if ( pNode->nodeType != CONTACT ) { 635 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiT - dxOverDy * dPsiL; 636 if ( pElem->elemType == SEMICON ) { 637 pRhs[ pNode->nEqn ] -= dy * pTEdge->jn + dx * pLEdge->jn; 638 } 639 } 640 pNode = pElem->pTRNode; 641 if ( pNode->nodeType != CONTACT ) { 642 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiT - dxOverDy * dPsiR; 643 if ( pElem->elemType == SEMICON ) { 644 pRhs[ pNode->nEqn ] -= -dy * pTEdge->jn + dx * pREdge->jn; 645 } 646 } 647 pNode = pElem->pBRNode; 648 if ( pNode->nodeType != CONTACT ) { 649 pRhs[ pNode->psiEqn ] -= dyOverDx * dPsiB + dxOverDy * dPsiR; 650 if ( pElem->elemType == SEMICON ) { 651 pRhs[ pNode->nEqn ] -= -dy * pBEdge->jn - dx * pREdge->jn; 652 } 653 } 654 pNode = pElem->pBLNode; 655 if ( pNode->nodeType != CONTACT ) { 656 pRhs[ pNode->psiEqn ] -= -dyOverDx * dPsiB + dxOverDy * dPsiL; 657 if ( pElem->elemType == SEMICON ) { 658 pRhs[ pNode->nEqn ] -= dy * pBEdge->jn - dx * pLEdge->jn; 659 } 660 } 661 } 662 } 663 664 /* 665 * computation of current densities, recombination rates, 666 * mobilities and their derivatives 667 */ 668 void TWONcommonTerms(TWOdevice * pDevice,BOOLEAN currentOnly,BOOLEAN tranAnalysis,TWOtranInfo * info)669 TWONcommonTerms(TWOdevice *pDevice, BOOLEAN currentOnly, 670 BOOLEAN tranAnalysis, TWOtranInfo *info) 671 { 672 TWOelem *pElem; 673 TWOedge *pEdge; 674 TWOnode *pNode; 675 int index, eIndex; 676 int nextIndex; /* index of node to find next element */ 677 double psi1, psi2, refPsi, nC, nP1; 678 double dPsiN; 679 double bPsiN, dbPsiN, bMPsiN, dbMPsiN; 680 double muN, dMuN, rDx, rDy; 681 double psi, nConc = 0.0, pConc = 0.0; 682 double cnAug, cpAug; 683 double eSurf = 0.0; /* For channel mobilities */ 684 double qInt = 0.0; 685 TWOchannel *pCh; 686 687 /* evaluate all node (including recombination) and edge quantities */ 688 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 689 pElem = pDevice->elements[ eIndex ]; 690 refPsi = pElem->matlInfo->refPsi; 691 cnAug = pElem->matlInfo->cAug[ELEC]; 692 cpAug = pElem->matlInfo->cAug[HOLE]; 693 for ( index = 0; index <= 3; index++ ) { 694 if ( pElem->evalNodes[ index ] ) { 695 pNode = pElem->pNodes[ index ]; 696 if ( pNode->nodeType != CONTACT ) { 697 psi = pDevice->dcSolution[ pNode->psiEqn ]; 698 if ( pElem->elemType == SEMICON ) { 699 nConc = pDevice->dcSolution[ pNode->nEqn ]; 700 pConc = pNode->nie * exp( - psi + refPsi ); 701 if ( Srh ) { 702 recomb(nConc, pConc, 703 pNode->tn, pNode->tp, cnAug, cpAug, pNode->nie, 704 &pNode->uNet, &pNode->dUdN, &pNode->dUdP); 705 } else { 706 pNode->uNet = 0.0; 707 pNode->dUdN = 0.0; 708 pNode->dUdP = 0.0; 709 } 710 } 711 } else { 712 /* a contact node */ 713 psi = pNode->psi; 714 if ( pElem->elemType == SEMICON ) { 715 nConc = pNode->nConc; 716 pConc = pNode->pConc; 717 } 718 } 719 720 /* store info in the state tables */ 721 pDevice->devState0 [pNode->nodePsi] = psi; 722 if ( pElem->elemType == SEMICON ) { 723 pDevice->devState0 [pNode->nodeN] = nConc; 724 pDevice->devState0 [pNode->nodeP] = pConc; 725 if ( tranAnalysis && pNode->nodeType != CONTACT ) { 726 pNode->dNdT = integrate( pDevice->devStates, info, pNode->nodeN ); 727 } 728 } 729 } 730 } 731 for ( index = 0; index <= 3; index++ ) { 732 if ( pElem->evalEdges[ index ] ) { 733 pEdge = pElem->pEdges[ index ]; 734 pNode = pElem->pNodes[ index ]; 735 if ( pNode->nodeType != CONTACT ) { 736 psi1 = pDevice->dcSolution[pNode->psiEqn]; 737 } else { 738 psi1 = pNode->psi; 739 } 740 pNode = pElem->pNodes[ (index + 1) % 4 ]; 741 if ( pNode->nodeType != CONTACT ) { 742 psi2 = pDevice->dcSolution[pNode->psiEqn]; 743 } else { 744 psi2 = pNode->psi; 745 } 746 if ( index <= 1 ) { 747 pEdge->dPsi = psi2 - psi1; 748 } else { 749 pEdge->dPsi = psi1 - psi2; 750 } 751 pDevice->devState0 [pEdge->edgeDpsi] = pEdge->dPsi; 752 753 if ( pElem->elemType == SEMICON ) { 754 /* Calculate weighted driving forces - wdfn & wdfp for the edge */ 755 dPsiN = pEdge->dPsi + pEdge->dCBand; 756 bernoulli( dPsiN, &bPsiN, &dbPsiN, 757 &bMPsiN, &dbMPsiN, !currentOnly ); 758 if ( index <= 1 ) { 759 nC = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeN); 760 nP1 = *(pDevice->devState0 + pElem->pNodes[ index+1 ]->nodeN); 761 } else { 762 nC = *(pDevice->devState0 + pElem->pNodes[(index+1)%4]->nodeN); 763 nP1 = *(pDevice->devState0 + pElem->pNodes[ index ]->nodeN); 764 } 765 pEdge->wdfn = bPsiN * nP1 - bMPsiN * nC; 766 pEdge->jn = 0.0; 767 if ( !currentOnly ) { 768 pEdge->dWnDpsiP1 = dbPsiN * nP1 - dbMPsiN * nC; 769 pEdge->dWnDn = - bMPsiN; 770 pEdge->dWnDnP1 = bPsiN; 771 pEdge->dJnDpsiP1 = 0.0; 772 pEdge->dJnDn = 0.0; 773 pEdge->dJnDnP1 = 0.0; 774 } 775 } 776 } 777 } 778 } 779 780 /* DAG: calculate mobilities for channel elems */ 781 if ( SurfaceMobility ) { 782 for ( pCh = pDevice->pChannel; 783 pCh != NULL; pCh = pCh->next ) { 784 pElem = pCh->pNElem; 785 switch (pCh->type) { 786 case 0: 787 eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi ) 788 * pElem->epsRel / pElem->dy; 789 qInt = 0.5 * pElem->pBotEdge->qf; 790 break; 791 case 1: 792 eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi ) 793 * pElem->epsRel / pElem->dx; 794 qInt = 0.5 * pElem->pLeftEdge->qf; 795 break; 796 case 2: 797 eSurf = - 0.5 * (pElem->pLeftEdge->dPsi + pElem->pRightEdge->dPsi ) 798 * pElem->epsRel / pElem->dy; 799 qInt = 0.5 * pElem->pTopEdge->qf; 800 break; 801 case 3: 802 eSurf = - 0.5 * (pElem->pTopEdge->dPsi + pElem->pBotEdge->dPsi ) 803 * pElem->epsRel / pElem->dx; 804 qInt = 0.5 * pElem->pRightEdge->qf; 805 break; 806 } 807 eSurf += qInt; 808 pElem = pCh->pSeed; 809 nextIndex = (pCh->type + 2)%4; 810 while (pElem && pElem->channel == pCh->id) { 811 TWONmobility( pElem, eSurf ); 812 pElem = pElem->pElems[ nextIndex ]; 813 } 814 } /* endfor pCH != NULL */ 815 } /* endif SurfaceMobility */ 816 817 /* calculate the current densities assuming mobility value depend on Ex,Ey*/ 818 for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) { 819 pElem = pDevice->elements[ eIndex ]; 820 821 rDx = 1.0 / pElem->dx; 822 rDy = 1.0 / pElem->dy; 823 for ( index = 0; index <= 3; index++ ) { 824 pEdge = pElem->pEdges[ index ]; 825 /* calculate conductive currents */ 826 if ( pElem->elemType == SEMICON ) { 827 /* get mobility for this edge */ 828 if ( !pElem->channel ) { 829 /* Calculate mobility for non-channel elements */ 830 muN = pElem->mun0; 831 dMuN = 0.0; 832 dPsiN = pEdge->dPsi + pEdge->dCBand; 833 if ( index%2 == 0 ) { 834 MOBfieldDep( pElem->matlInfo, ELEC, - dPsiN * rDx, &muN, &dMuN ); 835 } else { 836 MOBfieldDep( pElem->matlInfo, ELEC, - dPsiN * rDy, &muN, &dMuN ); 837 } 838 } else { 839 /* Retrieve previously calculated value. */ 840 muN = pElem->mun; 841 dMuN = 0.0; 842 } 843 switch ( index ) { 844 case 0: 845 muN *= pEdge->kPos * rDx; 846 dMuN *= pEdge->kPos * rDx * rDx; 847 break; 848 case 1: 849 muN *= pEdge->kNeg * rDy; 850 dMuN *= pEdge->kNeg * rDy * rDy; 851 break; 852 case 2: 853 muN *= pEdge->kNeg * rDx; 854 dMuN *= pEdge->kNeg * rDx * rDx; 855 break; 856 case 3: 857 muN *= pEdge->kPos * rDy; 858 dMuN *= pEdge->kPos * rDy * rDy; 859 break; 860 } 861 /* Now that the mobility for this edge is known, do current */ 862 pEdge->jn += muN * pEdge->wdfn; 863 if ( !currentOnly ) { 864 pEdge->dJnDpsiP1 += muN * pEdge->dWnDpsiP1; 865 pEdge->dJnDn += muN * pEdge->dWnDn; 866 pEdge->dJnDnP1 += muN * pEdge->dWnDnP1; 867 if ( MobDeriv && (!pElem->channel) ) { 868 pEdge->dJnDpsiP1 -= dMuN * pEdge->wdfn; 869 } 870 } 871 } 872 /* calculate displacement current only once */ 873 if ( pElem->evalEdges[ index ] ) { 874 if ( tranAnalysis ) { 875 if ( index == 0 || index == 2 ) { 876 /* horizontal edges */ 877 pEdge->jd = -integrate(pDevice->devStates, info, 878 pEdge->edgeDpsi) * rDx; 879 } else { 880 /* vertical edges */ 881 pEdge->jd = -integrate(pDevice->devStates, info, 882 pEdge->edgeDpsi) * rDy; 883 } 884 } 885 } 886 } 887 } 888 } 889