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