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 "../../maths/misc/norm.h"
8 #include "ngspice/bool.h"
9 #include "ngspice/cidersupt.h"
10 #include "ngspice/cpextern.h"
11 #include "ngspice/ifsim.h"
12 #include "ngspice/macros.h"
13 #include "ngspice/ngspice.h"
14 #include "ngspice/numenum.h"
15 #include "ngspice/numglobs.h"
16 #include "ngspice/spmatrix.h"
17 #include "ngspice/twodev.h"
18 #include "ngspice/twomesh.h"
19 #include "twoddefs.h"
20 #include "twodext.h"
21 extern IFfrontEnd *SPfrontEnd;
22 
23 
24 /* functions to calculate the 2D solutions */
25 
26 
27 /* the iteration driving loop and convergence check */
28 void
TWOdcSolve(TWOdevice * pDevice,int iterationLimit,BOOLEAN newSolver,BOOLEAN tranAnalysis,TWOtranInfo * info)29 TWOdcSolve(TWOdevice *pDevice, int iterationLimit, BOOLEAN newSolver,
30            BOOLEAN tranAnalysis, TWOtranInfo *info)
31 {
32   TWOnode *pNode;
33   TWOelem *pElem;
34   int size = pDevice->numEqns;
35   int index, eIndex, error;
36   int timesConverged = 0;
37   BOOLEAN quitLoop;
38   BOOLEAN debug;
39   BOOLEAN negConc = FALSE;
40   double *rhs = pDevice->rhs;
41   double *solution = pDevice->dcSolution;
42   double *delta = pDevice->dcDeltaSolution;
43   double poissNorm, contNorm;
44   double startTime, totalStartTime;
45   double totalTime, loadTime, factorTime, solveTime, updateTime, checkTime;
46   double orderTime = 0.0;
47 
48   totalTime = loadTime = factorTime = solveTime = updateTime = checkTime = 0.0;
49   totalStartTime = SPfrontEnd->IFseconds();
50 
51   quitLoop = FALSE;
52   debug = (!tranAnalysis && TWOdcDebug) || (tranAnalysis && TWOtranDebug);
53   pDevice->iterationNumber = 0;
54   pDevice->converged = FALSE;
55 
56   if (debug) {
57     if (pDevice->poissonOnly) {
58       fprintf(stdout, "Equilibrium Solution:\n");
59     } else {
60       fprintf(stdout, "Bias Solution:\n");
61     }
62     fprintf(stdout, "Iteration  RHS Norm\n");
63   }
64   while (!(pDevice->converged || (pDevice->iterationNumber > iterationLimit)
65 	  || quitLoop)) {
66     pDevice->iterationNumber++;
67 
68     if ((!pDevice->poissonOnly) && (iterationLimit > 0)
69 	&&(!tranAnalysis)) {
70       TWOjacCheck(pDevice, tranAnalysis, info);
71     }
72 
73     /* LOAD */
74     startTime = SPfrontEnd->IFseconds();
75     if (pDevice->poissonOnly) {
76       TWOQsysLoad(pDevice);
77     } else if (!OneCarrier) {
78       TWO_sysLoad(pDevice, tranAnalysis, info);
79     } else if (OneCarrier == N_TYPE) {
80       TWONsysLoad(pDevice, tranAnalysis, info);
81     } else if (OneCarrier == P_TYPE) {
82       TWOPsysLoad(pDevice, tranAnalysis, info);
83     }
84     pDevice->rhsNorm = maxNorm(rhs, size);
85     loadTime += SPfrontEnd->IFseconds() - startTime;
86     if (debug) {
87       fprintf(stdout, "%7d   %11.4e%s\n",
88 	  pDevice->iterationNumber - 1, pDevice->rhsNorm,
89 	  negConc ? "   negative conc encountered" : "");
90       negConc = FALSE;
91     }
92 
93     /* FACTOR */
94     startTime = SPfrontEnd->IFseconds();
95     error = spFactor(pDevice->matrix);
96     factorTime += SPfrontEnd->IFseconds() - startTime;
97     if (newSolver) {
98         if (pDevice->iterationNumber == 1) {
99             orderTime = factorTime;
100         }
101         else if (pDevice->iterationNumber == 2) {
102             orderTime -= factorTime - orderTime;
103             factorTime -= orderTime;
104             if (pDevice->poissonOnly) {
105                 pDevice->pStats->orderTime[STAT_SETUP] += orderTime;
106             }
107             else {
108                 pDevice->pStats->orderTime[STAT_DC] += orderTime;
109             }
110             /* After first two iterations, no special handling for a
111              * new solver */
112             newSolver = FALSE;
113         } /* end of case of iteratin 2 */
114     } /* end of special processing for a new solver */
115 
116     if (foundError(error)) {
117       if (error == spSINGULAR) {
118 	int badRow, badCol;
119 	spWhereSingular(pDevice->matrix, &badRow, &badCol);
120 	printf("*****  singular at (%d,%d)\n", badRow, badCol);
121       }
122       pDevice->converged = FALSE;
123       quitLoop = TRUE;
124       continue;
125     }
126 
127     /* SOLVE */
128     startTime = SPfrontEnd->IFseconds();
129     spSolve(pDevice->matrix, rhs, delta, NULL, NULL);
130     solveTime += SPfrontEnd->IFseconds() - startTime;
131 
132     /* UPDATE */
133     startTime = SPfrontEnd->IFseconds();
134     /* Use norm reducing Newton method for DC bias solutions only. */
135     if ((!pDevice->poissonOnly) && (iterationLimit > 0)
136 	&& (!tranAnalysis) && (pDevice->rhsNorm > 1e-1)) {
137       error = TWOnewDelta(pDevice, tranAnalysis, info);
138       if (error) {
139 	pDevice->converged = FALSE;
140 	quitLoop = TRUE;
141 	updateTime += SPfrontEnd->IFseconds() - startTime;
142 	continue;
143       }
144     }
145     for (index = 1; index <= size; index++) {
146       solution[index] += delta[index];
147     }
148     updateTime += SPfrontEnd->IFseconds() - startTime;
149 
150     /* CHECK CONVERGENCE */
151     startTime = SPfrontEnd->IFseconds();
152     /* Check if updates have gotten sufficiently small. */
153     if (pDevice->iterationNumber != 1) {
154       pDevice->converged = TWOdeltaConverged(pDevice);
155     }
156     /* Check if the residual is smaller than abstol. */
157     if (pDevice->converged && (!pDevice->poissonOnly)
158 	&& (!tranAnalysis)) {
159       if (!OneCarrier) {
160 	TWO_rhsLoad(pDevice, tranAnalysis, info);
161       } else if (OneCarrier == N_TYPE) {
162 	TWONrhsLoad(pDevice, tranAnalysis, info);
163       } else if (OneCarrier == P_TYPE) {
164 	TWOPrhsLoad(pDevice, tranAnalysis, info);
165       }
166       pDevice->rhsNorm = maxNorm(rhs, size);
167       if (pDevice->rhsNorm > pDevice->abstol) {
168 	pDevice->converged = FALSE;
169       }
170       if ((++timesConverged >= 2)
171 	  && (pDevice->rhsNorm < 1e3 * pDevice->abstol)) {
172 	pDevice->converged = TRUE;
173       } else if (timesConverged >= 5) {
174 	pDevice->converged = FALSE;
175 	quitLoop = TRUE;
176 	continue;
177       }
178     } else if (pDevice->converged && pDevice->poissonOnly) {
179       TWOQrhsLoad(pDevice);
180       pDevice->rhsNorm = maxNorm(rhs, size);
181       if (pDevice->rhsNorm > pDevice->abstol) {
182 	pDevice->converged = FALSE;
183       }
184       if (++timesConverged >= 5) {
185 	pDevice->converged = TRUE;
186       }
187     }
188     /* Check if the carrier concentrations are negative. */
189     if (pDevice->converged && (!pDevice->poissonOnly)) {
190       /* Clear garbage entry since carrier-free elements reference it. */
191       solution[0] = 0.0;
192       for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
193 	pElem = pDevice->elements[eIndex];
194 	for (index = 0; index <= 3; index++) {
195 	  if (pElem->evalNodes[index]) {
196 	    pNode = pElem->pNodes[index];
197 	    if (solution[pNode->nEqn] < 0.0) {
198 	      pDevice->converged = FALSE;
199 	      negConc = TRUE;
200 	      if (tranAnalysis) {
201 		quitLoop = TRUE;
202 	      } else {
203 		solution[pNode->nEqn] = 0.0;
204 	      }
205 	    }
206 	    if (solution[pNode->pEqn] < 0.0) {
207 	      pDevice->converged = FALSE;
208 	      negConc = TRUE;
209 	      if (tranAnalysis) {
210 		quitLoop = TRUE;
211 	      } else {
212 		solution[pNode->pEqn] = 0.0;
213 	      }
214 	    }
215 	  }
216 	}
217       }
218       if (!pDevice->converged) {
219 	if (!OneCarrier) {
220 	  TWO_rhsLoad(pDevice, tranAnalysis, info);
221 	} else if (OneCarrier == N_TYPE) {
222 	  TWONrhsLoad(pDevice, tranAnalysis, info);
223 	} else if (OneCarrier == P_TYPE) {
224 	  TWOPrhsLoad(pDevice, tranAnalysis, info);
225 	}
226 	pDevice->rhsNorm = maxNorm(rhs, size);
227       }
228     }
229     checkTime += SPfrontEnd->IFseconds() - startTime;
230   }
231   totalTime += SPfrontEnd->IFseconds() - totalStartTime;
232 
233   if (tranAnalysis) {
234     pDevice->pStats->loadTime[STAT_TRAN] += loadTime;
235     pDevice->pStats->factorTime[STAT_TRAN] += factorTime;
236     pDevice->pStats->solveTime[STAT_TRAN] += solveTime;
237     pDevice->pStats->updateTime[STAT_TRAN] += updateTime;
238     pDevice->pStats->checkTime[STAT_TRAN] += checkTime;
239     pDevice->pStats->numIters[STAT_TRAN] += pDevice->iterationNumber;
240   } else if (pDevice->poissonOnly) {
241     pDevice->pStats->loadTime[STAT_SETUP] += loadTime;
242     pDevice->pStats->factorTime[STAT_SETUP] += factorTime;
243     pDevice->pStats->solveTime[STAT_SETUP] += solveTime;
244     pDevice->pStats->updateTime[STAT_SETUP] += updateTime;
245     pDevice->pStats->checkTime[STAT_SETUP] += checkTime;
246     pDevice->pStats->numIters[STAT_SETUP] += pDevice->iterationNumber;
247   } else {
248     pDevice->pStats->loadTime[STAT_DC] += loadTime;
249     pDevice->pStats->factorTime[STAT_DC] += factorTime;
250     pDevice->pStats->solveTime[STAT_DC] += solveTime;
251     pDevice->pStats->updateTime[STAT_DC] += updateTime;
252     pDevice->pStats->checkTime[STAT_DC] += checkTime;
253     pDevice->pStats->numIters[STAT_DC] += pDevice->iterationNumber;
254   }
255 
256   if (debug) {
257     if (!tranAnalysis) {
258       pDevice->rhsNorm = maxNorm(rhs, size);
259       fprintf(stdout, "%7d   %11.4e%s\n",
260 	  pDevice->iterationNumber, pDevice->rhsNorm,
261 	  negConc ? "   negative conc in solution" : "");
262     }
263     if (pDevice->converged) {
264       if (!pDevice->poissonOnly) {
265 	poissNorm = contNorm = 0.0;
266 	rhs[0] = 0.0;		/* Make sure garbage entry is clear. */
267 	for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
268 	  pElem = pDevice->elements[eIndex];
269 	  for (index = 0; index <= 3; index++) {
270 	    if (pElem->evalNodes[index]) {
271 	      pNode = pElem->pNodes[index];
272 	      poissNorm = MAX(poissNorm,ABS(rhs[pNode->psiEqn]));
273 	      contNorm = MAX(contNorm,ABS(rhs[pNode->nEqn]));
274 	      contNorm = MAX(contNorm,ABS(rhs[pNode->pEqn]));
275 	    }
276 	  }
277 	}
278 	fprintf(stdout,
279 	    "Residual: %11.4e C/um poisson, %11.4e A/um continuity\n",
280 	    poissNorm * EpsNorm * VNorm * 1e-4,
281 	    contNorm * JNorm * LNorm * 1e-4);
282       } else {
283 	fprintf(stdout, "Residual: %11.4e C/um poisson\n",
284 	    pDevice->rhsNorm * EpsNorm * VNorm * 1e-4);
285       }
286     }
287   }
288 }
289 
290 BOOLEAN
TWOdeltaConverged(TWOdevice * pDevice)291 TWOdeltaConverged(TWOdevice *pDevice)
292 {
293   /* This function returns a TRUE if converged else a FALSE. */
294   int index;
295   BOOLEAN converged = TRUE;
296   double xOld, xNew, tol;
297 
298   for (index = 1; index <= pDevice->numEqns; index++) {
299     xOld = pDevice->dcSolution[index];
300     xNew = xOld + pDevice->dcDeltaSolution[index];
301     tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
302     if (ABS(xOld - xNew) > tol) {
303       converged = FALSE;
304       break;
305     }
306   }
307   return (converged);
308 }
309 
310 BOOLEAN
TWOdeviceConverged(TWOdevice * pDevice)311 TWOdeviceConverged(TWOdevice *pDevice)
312 {
313   int index, eIndex;
314   BOOLEAN converged = TRUE;
315   double *solution = pDevice->dcSolution;
316   TWOnode *pNode;
317   TWOelem *pElem;
318   double startTime;
319 
320   /* If the update is sufficiently small, and the carrier concentrations
321    * are all positive, then return TRUE, else return FALSE.
322    */
323 
324   /* CHECK CONVERGENCE */
325   startTime = SPfrontEnd->IFseconds();
326   if ((converged = TWOdeltaConverged(pDevice)) == TRUE) {
327     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
328       pElem = pDevice->elements[eIndex];
329       for (index = 0; index <= 3; index++) {
330 	if (pElem->evalNodes[index]) {
331 	  pNode = pElem->pNodes[index];
332 	  if (pNode->nEqn != 0 && solution[pNode->nEqn] < 0.0) {
333 	    converged = FALSE;
334 	    solution[pNode->nEqn] = pNode->nConc = 0.0;
335 	  }
336 	  if (pNode->pEqn != 0 && solution[pNode->pEqn] < 0.0) {
337 	    converged = FALSE;
338 	    solution[pNode->pEqn] = pNode->pConc = 0.0;
339 	  }
340 	}
341       }
342     }
343   }
344   pDevice->pStats->checkTime[STAT_TRAN] += SPfrontEnd->IFseconds() - startTime;
345 
346   return (converged);
347 }
348 
349 void
TWOresetJacobian(TWOdevice * pDevice)350 TWOresetJacobian(TWOdevice *pDevice)
351 {
352   int error;
353 
354 
355   if (!OneCarrier) {
356     TWO_jacLoad(pDevice);
357   } else if (OneCarrier == N_TYPE) {
358     TWONjacLoad(pDevice);
359   } else if (OneCarrier == P_TYPE) {
360     TWOPjacLoad(pDevice);
361   } else {
362     printf("TWOresetJacobian: unknown carrier type\n");
363     exit(-1);
364   }
365   error = spFactor(pDevice->matrix);
366   if (foundError(error)) {
367     exit(-1);
368   }
369 }
370 
371 /*
372  * Compute the device state assuming charge neutrality exists everywhere in
373  * the device.  In insulators, where there is no charge, assign the potential
374  * at half the insulator band gap (refPsi).
375  */
376 void
TWOstoreNeutralGuess(TWOdevice * pDevice)377 TWOstoreNeutralGuess(TWOdevice *pDevice)
378 {
379   int nIndex, eIndex;
380   TWOelem *pElem;
381   TWOnode *pNode;
382   double nie, conc, absConc, refPsi, psi, ni, pi, sign;
383 
384   /* assign the initial guess for Poisson's equation */
385   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
386     pElem = pDevice->elements[eIndex];
387     refPsi = pElem->matlInfo->refPsi;
388     if (pElem->elemType == INSULATOR) {
389       for (nIndex = 0; nIndex <= 3; nIndex++) {
390 	if (pElem->evalNodes[nIndex]) {
391 	  pNode = pElem->pNodes[nIndex];
392 	  if (pNode->nodeType == CONTACT) {
393 	    /*
394 	     * a metal contact to insulator domain so use work function
395 	     * difference as the value of psi
396 	     */
397 	    pNode->psi = RefPsi - pNode->eaff;
398 	  } else {
399 	    pNode->psi = refPsi;
400 	  }
401 	}
402       }
403     }
404     if (pElem->elemType == SEMICON) {
405       for (nIndex = 0; nIndex <= 3; nIndex++) {
406 	if (pElem->evalNodes[nIndex]) {
407 	  pNode = pElem->pNodes[nIndex];
408 	  /* silicon nodes */
409 	  nie = pNode->nie;
410 	  conc = pNode->netConc / pNode->nie;
411 	  psi = 0.0;
412 	  ni = nie;
413 	  pi = nie;
414 	  sign = SGN(conc);
415 	  absConc = ABS(conc);
416 	  if (conc != 0.0) {
417 	    psi = sign * log(0.5 * absConc
418 		+ sqrt(1.0 + 0.25 * absConc * absConc));
419 	    ni = nie * exp(psi);
420 	    pi = nie * exp(-psi);
421 	  }
422 	  pNode->psi = refPsi + psi;
423 	  pNode->nConc = ni;
424 	  pNode->pConc = pi;
425 	  /* store the initial guess in the dc solution vector */
426 	  if (pNode->nodeType != CONTACT) {
427 	    pDevice->dcSolution[pNode->poiEqn] = pNode->psi;
428 	  }
429 	}
430       }
431     }
432   }
433 }
434 
435 /* computing the equilibrium solution; solution of Poisson's eqn */
436 /* the solution is used as an initial starting point for bias conditions */
TWOequilSolve(TWOdevice * pDevice)437 int TWOequilSolve(TWOdevice *pDevice)
438 {
439     BOOLEAN newSolver = FALSE;
440     int error;
441     int nIndex, eIndex;
442     TWOelem *pElem;
443     TWOnode *pNode;
444     double startTime, setupTime, miscTime;
445 
446     setupTime = miscTime = 0.0;
447 
448     /* SETUP */
449     startTime = SPfrontEnd->IFseconds();
450 
451     /* Set up pDevice to compute the equilibrium solution. If the solver
452      * is for bias, the arrays must be freed and allocated to the correct
453      * sizes for an equilibrium solution; if it is a new solver, they are
454      * only allocated; and if already an equilibrium solve, nothing
455      * needs to be done */
456     /* FALLTHROUGH added to suppress GCC warning due to
457      * -Wimplicit-fallthrough flag */
458     switch (pDevice->solverType) {
459         case SLV_SMSIG:
460         case SLV_BIAS:
461             /* Free memory allocated for the bias solution */
462             FREE(pDevice->dcSolution);
463             FREE(pDevice->dcDeltaSolution);
464             FREE(pDevice->copiedSolution);
465             FREE(pDevice->rhs);
466             FREE(pDevice->rhsImag);
467             spDestroy(pDevice->matrix);
468             /* FALLTHROUGH */
469         case SLV_NONE: {
470             /* Allocate memory needed for an equilibrium solution */
471             const int n_dim = pDevice->dimEquil;
472             const int n_eqn = n_dim - 1;
473             pDevice->poissonOnly = TRUE;
474             pDevice->numEqns = n_eqn;
475             XCALLOC(pDevice->dcSolution, double, n_dim);
476             XCALLOC(pDevice->dcDeltaSolution, double, n_dim);
477             XCALLOC(pDevice->copiedSolution, double, n_dim);
478             XCALLOC(pDevice->rhs, double, n_dim);
479             pDevice->matrix = spCreate(n_eqn, 0, &error);
480             if (error == spNO_MEMORY) {
481                 (void) fprintf(cp_err, "TWOequilSolve: Out of Memory\n");
482                 return E_NOMEM;
483             }
484             newSolver = TRUE;
485             spSetReal(pDevice->matrix); /* set to a real matrix */
486             TWOQjacBuild(pDevice);
487             pDevice->numOrigEquil = spElementCount(pDevice->matrix);
488             pDevice->numFillEquil = 0;
489             pDevice->solverType = SLV_EQUIL;
490             break;
491         }
492         case SLV_EQUIL: /* Nothing to do if already equilibrium solver */
493             break;
494         default: /* Invalid data */
495             fprintf(stderr, "Panic: Unknown solver type in equil solution.\n");
496             return E_PANIC;
497     } /* end of switch over solve type */
498     TWOstoreNeutralGuess(pDevice);
499     setupTime += SPfrontEnd->IFseconds() - startTime;
500 
501   /* SOLVE */
502   TWOdcSolve(pDevice, MaxIterations, newSolver, FALSE, NULL);
503 
504   /* MISCELLANEOUS */
505   startTime = SPfrontEnd->IFseconds();
506   if (newSolver) {
507     pDevice->numFillEquil = spFillinCount(pDevice->matrix);
508   }
509   if (pDevice->converged) {
510     TWOQcommonTerms(pDevice);
511 
512     /* save equilibrium potential */
513     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
514       pElem = pDevice->elements[eIndex];
515       for (nIndex = 0; nIndex <= 3; nIndex++) {
516 	if (pElem->evalNodes[nIndex]) {
517 	  pNode = pElem->pNodes[nIndex];
518 	  pNode->psi0 = pNode->psi;
519 	}
520       }
521     }
522   } else {
523     printf("TWOequilSolve: No Convergence\n");
524   }
525   miscTime += SPfrontEnd->IFseconds() - startTime;
526   pDevice->pStats->setupTime[STAT_SETUP] += setupTime;
527   pDevice->pStats->miscTime[STAT_SETUP] += miscTime;
528 
529   return 0;
530 }
531 
532 /* compute the solution under an applied bias */
533 /* the equilibrium solution is taken as an initial guess */
534 void
TWObiasSolve(TWOdevice * pDevice,int iterationLimit,BOOLEAN tranAnalysis,TWOtranInfo * info)535 TWObiasSolve(TWOdevice *pDevice, int iterationLimit, BOOLEAN tranAnalysis,
536              TWOtranInfo *info)
537 {
538   BOOLEAN newSolver = FALSE;
539   int error;
540   int index, eIndex;
541   TWOelem *pElem;
542   TWOnode *pNode;
543   double refPsi;
544   double startTime, setupTime, miscTime;
545 
546   setupTime = miscTime = 0.0;
547 
548   /* SETUP */
549   startTime = SPfrontEnd->IFseconds();
550     /* Set up for solving for bias */
551     /* FALLTHROUGH added to suppress GCC warning due to
552      * -Wimplicit-fallthrough flag */
553     switch (pDevice->solverType) {
554     case SLV_EQUIL:
555         /* free up the vectors allocated in the equilibrium solution */
556         FREE(pDevice->dcSolution);
557         FREE(pDevice->dcDeltaSolution);
558         FREE(pDevice->copiedSolution);
559         FREE(pDevice->rhs);
560         spDestroy(pDevice->matrix);
561         /* FALLTHROUGH */
562   case SLV_NONE:
563     /* Set up for bias */
564     pDevice->poissonOnly = FALSE;
565     pDevice->numEqns = pDevice->dimBias - 1;
566     XCALLOC(pDevice->dcSolution, double, pDevice->dimBias);
567     XCALLOC(pDevice->dcDeltaSolution, double, pDevice->dimBias);
568     XCALLOC(pDevice->copiedSolution, double, pDevice->dimBias);
569     XCALLOC(pDevice->rhs, double, pDevice->dimBias);
570     XCALLOC(pDevice->rhsImag, double, pDevice->dimBias);
571     pDevice->matrix = spCreate(pDevice->numEqns, 1, &error);
572     if (error == spNO_MEMORY) {
573       printf("TWObiasSolve: Out of Memory\n");
574       exit(-1);
575     }
576     newSolver = TRUE;
577     if (!OneCarrier) {
578       TWO_jacBuild(pDevice);
579     } else if (OneCarrier == N_TYPE) {
580       TWONjacBuild(pDevice);
581     } else if (OneCarrier == P_TYPE) {
582       TWOPjacBuild(pDevice);
583     }
584     pDevice->numOrigBias = spElementCount(pDevice->matrix);
585     pDevice->numFillBias = 0;
586     TWOstoreInitialGuess(pDevice);
587         /* FALLTHROUGH */
588   case SLV_SMSIG:
589     spSetReal(pDevice->matrix);
590     pDevice->solverType = SLV_BIAS;
591     break;
592   case SLV_BIAS:
593     break;
594   default:
595     fprintf(stderr, "Panic: Unknown solver type in bias solution.\n");
596     exit(-1);
597     break;
598   }
599   setupTime += SPfrontEnd->IFseconds() - startTime;
600 
601   /* SOLVE */
602   TWOdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info);
603 
604   /* MISCELLANEOUS */
605   startTime = SPfrontEnd->IFseconds();
606   if (newSolver) {
607     pDevice->numFillBias = spFillinCount(pDevice->matrix);
608   }
609   if ((!pDevice->converged) && iterationLimit > 1) {
610     printf("TWObiasSolve: No Convergence\n");
611   } else if (pDevice->converged) {
612     /* update the nodal quantities */
613     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
614       pElem = pDevice->elements[eIndex];
615       refPsi = pElem->matlInfo->refPsi;
616       for (index = 0; index <= 3; index++) {
617 	if (pElem->evalNodes[index]) {
618 	  pNode = pElem->pNodes[index];
619 	  if (pNode->nodeType != CONTACT) {
620 	    pNode->psi = pDevice->dcSolution[pNode->psiEqn];
621 	    if (pElem->elemType == SEMICON) {
622 	      if (!OneCarrier) {
623 		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
624 		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
625 	      } else if (OneCarrier == N_TYPE) {
626 		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
627 		pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
628 	      } else if (OneCarrier == P_TYPE) {
629 		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
630 		pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
631 	      }
632 	    }
633 	  }
634 	}
635       }
636     }
637 
638     /* update the current terms */
639     if (!OneCarrier) {
640       TWO_commonTerms(pDevice, FALSE, tranAnalysis, info);
641     } else if (OneCarrier == N_TYPE) {
642       TWONcommonTerms(pDevice, FALSE, tranAnalysis, info);
643     } else if (OneCarrier == P_TYPE) {
644       TWOPcommonTerms(pDevice, FALSE, tranAnalysis, info);
645     }
646   } else if (iterationLimit <= 1) {
647     for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
648       pElem = pDevice->elements[eIndex];
649       refPsi = pElem->matlInfo->refPsi;
650       for (index = 0; index <= 3; index++) {
651 	if (pElem->evalNodes[index]) {
652 	  pNode = pElem->pNodes[index];
653 	  if (pNode->nodeType != CONTACT) {
654 	    pNode->psi = pDevice->dcSolution[pNode->psiEqn];
655 	    pDevice->devState0 [pNode->nodePsi] = pNode->psi;
656 	    if (pElem->elemType == SEMICON) {
657 	      if (!OneCarrier) {
658 		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
659 		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
660 	      } else if (OneCarrier == N_TYPE) {
661 		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
662 		pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
663 	      } else if (OneCarrier == P_TYPE) {
664 		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
665 		pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
666 	      }
667 	      pDevice->devState0 [pNode->nodeN] = pNode->nConc;
668 	      pDevice->devState0 [pNode->nodeP] = pNode->pConc;
669 	    }
670 	  }
671 	}
672       }
673     }
674   }
675   miscTime += SPfrontEnd->IFseconds() - startTime;
676   if (tranAnalysis) {
677     pDevice->pStats->setupTime[STAT_TRAN] += setupTime;
678     pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
679   } else {
680     pDevice->pStats->setupTime[STAT_DC] += setupTime;
681     pDevice->pStats->miscTime[STAT_DC] += miscTime;
682   }
683 }
684 
685 /* Copy the device's equilibrium state into the solution vector. */
686 void
TWOstoreEquilibGuess(TWOdevice * pDevice)687 TWOstoreEquilibGuess(TWOdevice *pDevice)
688 {
689   int nIndex, eIndex;
690   double *solution = pDevice->dcSolution;
691   double refPsi;
692   TWOelem *pElem;
693   TWOnode *pNode;
694 
695   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
696     pElem = pDevice->elements[eIndex];
697     refPsi = pElem->matlInfo->refPsi;
698     for (nIndex = 0; nIndex <= 3; nIndex++) {
699       if (pElem->evalNodes[nIndex]) {
700 	pNode = pElem->pNodes[nIndex];
701 	if (pNode->nodeType != CONTACT) {
702 	  pDevice->dcSolution[pNode->psiEqn] = pNode->psi0;
703 	  if (pElem->elemType == SEMICON) {
704 	    if (!OneCarrier) {
705 	      solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
706 	      solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
707 	    } else if (OneCarrier == N_TYPE) {
708 	      solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
709 	    } else if (OneCarrier == P_TYPE) {
710 	      solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
711 	    }
712 	  }
713 	}
714       }
715     }
716   }
717 }
718 
719 /* Copy the device's internal state into the solution vector. */
720 void
TWOstoreInitialGuess(TWOdevice * pDevice)721 TWOstoreInitialGuess(TWOdevice *pDevice)
722 {
723   int nIndex, eIndex;
724   double *solution = pDevice->dcSolution;
725   TWOelem *pElem;
726   TWOnode *pNode;
727 
728   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
729     pElem = pDevice->elements[eIndex];
730     for (nIndex = 0; nIndex <= 3; nIndex++) {
731       if (pElem->evalNodes[nIndex]) {
732 	pNode = pElem->pNodes[nIndex];
733 	if (pNode->nodeType != CONTACT) {
734 	  solution[pNode->psiEqn] = pNode->psi;
735 	  if (pElem->elemType == SEMICON) {
736 	    if (!OneCarrier) {
737 	      solution[pNode->nEqn] = pNode->nConc;
738 	      solution[pNode->pEqn] = pNode->pConc;
739 	    } else if (OneCarrier == N_TYPE) {
740 	      solution[pNode->nEqn] = pNode->nConc;
741 	    } else if (OneCarrier == P_TYPE) {
742 	      solution[pNode->pEqn] = pNode->pConc;
743 	    }
744 	  }
745 	}
746       }
747     }
748   }
749 }
750 
751 
752 /*
753  * function computeNewDelta computes an acceptable delta
754  */
755 
756 void
oldTWOnewDelta(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)757 oldTWOnewDelta(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
758 {
759   int index;
760   double newNorm, fib, lambda, fibn, fibp;
761   BOOLEAN acceptable = FALSE;
762   lambda = 1.0;
763   fibn = 1.0;
764   fibp = 1.0;
765 
766   /*
767    * copy the contents of dcSolution into copiedSolution and modify
768    * dcSolution by adding the deltaSolution
769    */
770 
771   for (index = 1; index <= pDevice->numEqns; ++index) {
772     pDevice->copiedSolution[index] = pDevice->dcSolution[index];
773     pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index];
774   }
775 
776   /* compute L2 norm of the deltaX vector */
777   pDevice->rhsNorm = l2Norm(pDevice->dcDeltaSolution, pDevice->numEqns);
778 
779   if (pDevice->poissonOnly) {
780     TWOQrhsLoad(pDevice);
781   } else if (!OneCarrier) {
782     TWO_rhsLoad(pDevice, tranAnalysis, info);
783   } else if (OneCarrier == N_TYPE) {
784     TWONrhsLoad(pDevice, tranAnalysis, info);
785   } else if (OneCarrier == P_TYPE) {
786     TWOPrhsLoad(pDevice, tranAnalysis, info);
787   }
788   newNorm = TWOnuNorm(pDevice);
789   if (newNorm <= pDevice->rhsNorm) {
790     acceptable = TRUE;
791   } else {
792     /* chop the step size so that deltax is acceptable */
793     for (; !acceptable;) {
794       fib = fibp;
795       fibp = fibn;
796       fibn += fib;
797       lambda *= (fibp / fibn);
798 
799       for (index = 1; index <= pDevice->numEqns; index++) {
800 	pDevice->dcSolution[index] = pDevice->copiedSolution[index]
801 	    + lambda * pDevice->dcDeltaSolution[index];
802       }
803       if (pDevice->poissonOnly) {
804 	TWOQrhsLoad(pDevice);
805       } else if (!OneCarrier) {
806 	TWO_rhsLoad(pDevice, tranAnalysis, info);
807       } else if (OneCarrier == N_TYPE) {
808 	TWONrhsLoad(pDevice, tranAnalysis, info);
809       } else if (OneCarrier == P_TYPE) {
810 	TWOPrhsLoad(pDevice, tranAnalysis, info);
811       }
812       newNorm = TWOnuNorm(pDevice);
813       if (newNorm <= pDevice->rhsNorm) {
814 	acceptable = TRUE;
815       }
816     }
817   }
818   /* restore the previous dcSolution and the new deltaSolution */
819   pDevice->rhsNorm = newNorm;
820   for (index = 1; index <= pDevice->numEqns; index++) {
821     pDevice->dcSolution[index] = pDevice->copiedSolution[index];
822     pDevice->dcDeltaSolution[index] *= lambda;
823   }
824 }
825 
826 
827 
828 int
TWOnewDelta(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)829 TWOnewDelta(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
830 {
831   int index, iterNum = 0;
832   double newNorm;
833   double fib, lambda, fibn, fibp;
834   BOOLEAN acceptable = FALSE, error = FALSE;
835 
836   lambda = 1.0;
837   fibn = 1.0;
838   fibp = 1.0;
839 
840   /*
841    * copy the contents of dcSolution into copiedSolution and modify
842    * dcSolution by adding the deltaSolution
843    */
844 
845   for (index = 1; index <= pDevice->numEqns; ++index) {
846     pDevice->copiedSolution[index] = pDevice->dcSolution[index];
847     pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index];
848   }
849 
850   if (pDevice->poissonOnly) {
851     TWOQrhsLoad(pDevice);
852   } else if (!OneCarrier) {
853     TWO_rhsLoad(pDevice, tranAnalysis, info);
854   } else if (OneCarrier == N_TYPE) {
855     TWONrhsLoad(pDevice, tranAnalysis, info);
856   } else if (OneCarrier == P_TYPE) {
857     TWOPrhsLoad(pDevice, tranAnalysis, info);
858   }
859   newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
860   if (pDevice->rhsNorm <= pDevice->abstol) {
861     lambda = 0.0;
862     newNorm = pDevice->rhsNorm;
863   } else if (newNorm < pDevice->rhsNorm) {
864     acceptable = TRUE;
865   } else {
866     /* chop the step size so that deltax is acceptable */
867     if (TWOdcDebug) {
868       fprintf(stdout, "          %11.4e  %11.4e\n",
869 	  newNorm, lambda);
870     }
871     while (!acceptable) {
872       iterNum++;
873 
874       if (iterNum > NORM_RED_MAXITERS) {
875 	error = TRUE;
876 	lambda = 0.0;
877 	/* Don't break out until after we've reset the device. */
878       }
879       fib = fibp;
880       fibp = fibn;
881       fibn += fib;
882       lambda *= (fibp / fibn);
883 
884       for (index = 1; index <= pDevice->numEqns; index++) {
885 	pDevice->dcSolution[index] = pDevice->copiedSolution[index]
886 	    + lambda * pDevice->dcDeltaSolution[index];
887       }
888       if (pDevice->poissonOnly) {
889 	TWOQrhsLoad(pDevice);
890       } else if (!OneCarrier) {
891 	TWO_rhsLoad(pDevice, tranAnalysis, info);
892       } else if (OneCarrier == N_TYPE) {
893 	TWONrhsLoad(pDevice, tranAnalysis, info);
894       } else if (OneCarrier == P_TYPE) {
895 	TWOPrhsLoad(pDevice, tranAnalysis, info);
896       }
897       newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
898       if (error) {
899 	break;
900       }
901       if (newNorm <= pDevice->rhsNorm) {
902 	acceptable = TRUE;
903       }
904       if (TWOdcDebug) {
905 	fprintf(stdout, "          %11.4e  %11.4e\n",
906 	    newNorm, lambda);
907       }
908     }
909   }
910   /* restore the previous dcSolution and the new deltaSolution */
911   pDevice->rhsNorm = newNorm;
912   for (index = 1; index <= pDevice->numEqns; index++) {
913     pDevice->dcSolution[index] = pDevice->copiedSolution[index];
914     pDevice->dcDeltaSolution[index] *= lambda;
915   }
916   return (error);
917 }
918 
919 
920 /* Predict the values of the internal variables at the next timepoint. */
921 /* Needed for Predictor-Corrector LTE estimation */
922 void
TWOpredict(TWOdevice * pDevice,TWOtranInfo * info)923 TWOpredict(TWOdevice *pDevice, TWOtranInfo *info)
924 {
925   int nIndex, eIndex;
926   TWOnode *pNode;
927   TWOelem *pElem;
928   double startTime, miscTime = 0.0;
929 
930   /* TRANSIENT MISC */
931   startTime = SPfrontEnd->IFseconds();
932   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
933     pElem = pDevice->elements[eIndex];
934     for (nIndex = 0; nIndex <= 3; nIndex++) {
935       if (pElem->evalNodes[nIndex]) {
936 	pNode = pElem->pNodes[nIndex];
937 	pNode->psi = pDevice->devState1 [pNode->nodePsi];
938 	if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
939 	  if (!OneCarrier) {
940 	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
941 	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
942 	  } else if (OneCarrier == N_TYPE) {
943 	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
944 	    pNode->pPred = pDevice->devState1 [pNode->nodeP];
945 	  } else if (OneCarrier == P_TYPE) {
946 	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
947 	    pNode->nPred = pDevice->devState1 [pNode->nodeN];
948 	  }
949 	  pNode->nConc = pNode->nPred;
950 	  pNode->pConc = pNode->pPred;
951 	}
952       }
953     }
954   }
955   miscTime += SPfrontEnd->IFseconds() - startTime;
956   pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
957 }
958 
959 
960 /* Estimate the device's overall truncation error. */
961 double
TWOtrunc(TWOdevice * pDevice,TWOtranInfo * info,double delta)962 TWOtrunc(TWOdevice *pDevice, TWOtranInfo *info, double delta)
963 {
964   int nIndex, eIndex;
965   TWOelem *pElem;
966   TWOnode *pNode;
967   double tolN, tolP, lte, relError, temp;
968   double lteCoeff = info->lteCoeff;
969   double mult = 10.0;
970   double reltol;
971   double startTime, lteTime = 0.0;
972 
973   /* TRANSIENT LTE */
974   startTime = SPfrontEnd->IFseconds();
975 
976   relError = 0.0;
977   reltol = pDevice->reltol * mult;
978 
979   /* need to get the predictor for the current order */
980   /* the scheme currently used is very dumb. need to fix later */
981   computePredCoeff(info->method, info->order, info->predCoeff, info->delta);
982 
983   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
984     pElem = pDevice->elements[eIndex];
985     for (nIndex = 0; nIndex <= 3; nIndex++) {
986       if (pElem->evalNodes[nIndex]) {
987 	pNode = pElem->pNodes[nIndex];
988 	if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
989 	  if (!OneCarrier) {
990 	    tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
991 	    tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
992 	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
993 	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
994 	    lte = lteCoeff * (pNode->nConc - pNode->nPred);
995 	    temp = lte / tolN;
996 	    relError += temp * temp;
997 	    lte = lteCoeff * (pNode->pConc - pNode->pPred);
998 	    temp = lte / tolP;
999 	    relError += temp * temp;
1000 	  } else if (OneCarrier == N_TYPE) {
1001 	    tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
1002 	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
1003 	    lte = lteCoeff * (pNode->nConc - pNode->nPred);
1004 	    temp = lte / tolN;
1005 	    relError += temp * temp;
1006 	  } else if (OneCarrier == P_TYPE) {
1007 	    tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
1008 	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
1009 	    lte = lteCoeff * (pNode->pConc - pNode->pPred);
1010 	    temp = lte / tolP;
1011 	    relError += temp * temp;
1012 	  }
1013 	}
1014       }
1015     }
1016   }
1017 
1018   relError = MAX(pDevice->abstol, relError);	/* make sure it is non zero */
1019 
1020   /* the total relative error has been calculated norm-2 squared */
1021   relError = sqrt(relError / pDevice->numEqns);
1022 
1023   /* depending on the order of the integration method compute new delta */
1024   temp = delta / pow(relError, 1.0 / (info->order + 1));
1025 
1026   lteTime += SPfrontEnd->IFseconds() - startTime;
1027   pDevice->pStats->lteTime += lteTime;
1028 
1029   /* return the new delta (stored as temp) */
1030   return (temp);
1031 }
1032 
1033 /* Save info from state table into the internal state. */
1034 void
TWOsaveState(TWOdevice * pDevice)1035 TWOsaveState(TWOdevice *pDevice)
1036 {
1037   int nIndex, eIndex;
1038   TWOnode *pNode;
1039   TWOelem *pElem;
1040 
1041   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1042     pElem = pDevice->elements[eIndex];
1043     for (nIndex = 0; nIndex <= 3; nIndex++) {
1044       if (pElem->evalNodes[nIndex]) {
1045 	pNode = pElem->pNodes[nIndex];
1046 	pNode->psi = pDevice->devState1 [pNode->nodePsi];
1047 	if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
1048 	  pNode->nConc = pDevice->devState1 [pNode->nodeN];
1049 	  pNode->pConc = pDevice->devState1 [pNode->nodeP];
1050 	}
1051       }
1052     }
1053   }
1054 }
1055 
1056 /*
1057  * A function that checks convergence based on the convergence of the quasi
1058  * Fermi levels. This should be better since we are looking at potentials in
1059  * all (psi, phin, phip)
1060  */
1061 BOOLEAN
TWOpsiDeltaConverged(TWOdevice * pDevice)1062 TWOpsiDeltaConverged(TWOdevice *pDevice)
1063 {
1064   int index, nIndex, eIndex;
1065   TWOnode *pNode;
1066   TWOelem *pElem;
1067   double xOld, xNew, xDelta, tol;
1068   double psi, newPsi, nConc, pConc, newN, newP;
1069   double phiN, phiP, newPhiN, newPhiP;
1070   BOOLEAN converged = TRUE;
1071 
1072   /* equilibrium solution */
1073   if (pDevice->poissonOnly) {
1074     for (index = 1; index <= pDevice->numEqns; index++) {
1075       xOld = pDevice->dcSolution[index];
1076       xDelta = pDevice->dcDeltaSolution[index];
1077       xNew = xOld + xDelta;
1078       tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
1079       if (ABS(xDelta) > tol) {
1080 	converged = FALSE;
1081 	goto done;
1082       }
1083     }
1084     return (converged);
1085   }
1086   /* bias solution. check convergence on psi, phin, phip */
1087   for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
1088     pElem = pDevice->elements[eIndex];
1089     for (nIndex = 0; nIndex <= 3; nIndex++) {
1090       if (pElem->evalNodes[nIndex]) {
1091 	pNode = pElem->pNodes[nIndex];
1092 	if (pNode->nodeType != CONTACT) {
1093 	  /* check convergence on psi */
1094 	  xOld = pDevice->dcSolution[pNode->psiEqn];
1095 	  xDelta = pDevice->dcDeltaSolution[pNode->psiEqn];
1096 	  xNew = xOld + xDelta;
1097 	  tol = pDevice->abstol +
1098 	      pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
1099 	  if (ABS(xDelta) > tol) {
1100 	    converged = FALSE;
1101 	    goto done;
1102 	  }
1103 	}
1104 	/* check convergence on phin and phip */
1105 	if ((pElem->elemType == SEMICON) && (pNode->nodeType != CONTACT)) {
1106 	  psi = pDevice->dcSolution[pNode->psiEqn];
1107 	  nConc = pDevice->dcSolution[pNode->nEqn];
1108 	  pConc = pDevice->dcSolution[pNode->pEqn];
1109 	  newPsi = psi + pDevice->dcDeltaSolution[pNode->psiEqn];
1110 	  newN = nConc + pDevice->dcDeltaSolution[pNode->nEqn];
1111 	  newP = pConc + pDevice->dcDeltaSolution[pNode->pEqn];
1112 	  phiN = psi - log(nConc / pNode->nie);
1113 	  phiP = psi + log(pConc / pNode->nie);
1114 	  newPhiN = newPsi - log(newN / pNode->nie);
1115 	  newPhiP = newPsi + log(newP / pNode->nie);
1116 	  tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiN), ABS(newPhiN));
1117 	  if (ABS(newPhiN - phiN) > tol) {
1118 	    converged = FALSE;
1119 	    goto done;
1120 	  }
1121 	  tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiP), ABS(newPhiP));
1122 	  if (ABS(newPhiP - phiP) > tol) {
1123 	    converged = FALSE;
1124 	    goto done;
1125 	  }
1126 	}
1127       }
1128     }
1129   }
1130 done:
1131   return (converged);
1132 }
1133 
1134 
1135 /* Function to compute Nu norm of the rhs vector. */
1136 /* nu-Norm calculation based upon work done at Stanford. */
1137 double
TWOnuNorm(TWOdevice * pDevice)1138 TWOnuNorm(TWOdevice *pDevice)
1139 {
1140   double norm = 0.0;
1141   double temp;
1142   int index;
1143 
1144   /* the LU Decomposed matrix is available. use it to calculate x */
1145 
1146   spSolve(pDevice->matrix, pDevice->rhs, pDevice->rhsImag,
1147       NULL, NULL);
1148 
1149   /* the solution is in the rhsImag vector */
1150   /* compute L2-norm of the rhsImag vector */
1151 
1152   for (index = 1; index <= pDevice->numEqns; index++) {
1153     temp = pDevice->rhsImag[index];
1154     norm += temp * temp;
1155   }
1156   norm = sqrt(norm);
1157   return (norm);
1158 }
1159 
1160 /*
1161  * Check for numerical errors in the Jacobian.  Useful debugging aid when new
1162  * models are being incorporated.
1163  */
1164 void
TWOjacCheck(TWOdevice * pDevice,BOOLEAN tranAnalysis,TWOtranInfo * info)1165 TWOjacCheck(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
1166 {
1167   int index, rIndex;
1168   double del, diff, tol, *dptr;
1169 
1170   if (TWOjacDebug) {
1171     if (!OneCarrier) {
1172       TWO_sysLoad(pDevice, tranAnalysis, info);
1173     } else if (OneCarrier == N_TYPE) {
1174       TWONsysLoad(pDevice, tranAnalysis, info);
1175     } else if (OneCarrier == P_TYPE) {
1176       TWOPsysLoad(pDevice, tranAnalysis, info);
1177     }
1178     /*
1179      * spPrint( pDevice->matrix, 0, 1, 1 );
1180      */
1181     pDevice->rhsNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
1182     for (index = 1; index <= pDevice->numEqns; index++) {
1183       if (1e3 * ABS(pDevice->rhs[index]) > pDevice->rhsNorm) {
1184 	fprintf(stderr, "eqn %d: res %11.4e, norm %11.4e\n",
1185 	    index, pDevice->rhs[index], pDevice->rhsNorm);
1186       }
1187     }
1188     for (index = 1; index <= pDevice->numEqns; index++) {
1189       pDevice->rhsImag[index] = pDevice->rhs[index];
1190     }
1191     for (index = 1; index <= pDevice->numEqns; index++) {
1192       pDevice->copiedSolution[index] = pDevice->dcSolution[index];
1193       del = 1e-4 * pDevice->abstol + 1e-6 * ABS(pDevice->dcSolution[index]);
1194       pDevice->dcSolution[index] += del;
1195       if (!OneCarrier) {
1196 	TWO_rhsLoad(pDevice, tranAnalysis, info);
1197       } else if (OneCarrier == N_TYPE) {
1198 	TWONrhsLoad(pDevice, tranAnalysis, info);
1199       } else if (OneCarrier == P_TYPE) {
1200 	TWOPrhsLoad(pDevice, tranAnalysis, info);
1201       }
1202       pDevice->dcSolution[index] = pDevice->copiedSolution[index];
1203       for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) {
1204 	diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del;
1205 	dptr = spFindElement(pDevice->matrix, rIndex, index);
1206 	if (dptr != NULL) {
1207 	  tol = (1e-4 * pDevice->abstol) + (1e-2 * MAX(ABS(diff), ABS(*dptr)));
1208 	  if ((diff != 0.0) && (ABS(diff - *dptr) > tol)) {
1209 	    fprintf(stderr, "Mismatch[%d][%d]: FD = %11.4e, AJ = %11.4e\n\t FD-AJ = %11.4e vs. %11.4e\n",
1210 		rIndex, index, diff, *dptr,
1211 		ABS(diff - *dptr), tol);
1212 	  }
1213 	} else {
1214 	  if (diff != 0.0) {
1215 	    fprintf(stderr, "Missing [%d][%d]: FD = %11.4e, AJ = 0.0\n",
1216 		rIndex, index, diff);
1217 	  }
1218 	}
1219       }
1220     }
1221   }
1222 }
1223