onesolve.c

ngspice又一个电子CAD仿真软件代码.功能更全

/**********
Copyright 1991 Regents of the University of California.  All rights reserved.
Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
Author:	1991 David A. Gates, U. C. Berkeley CAD Group
$Id: onesolve.c,v 1.3 2005/05/21 12:37:23 sjborley Exp $
**********/

/*
 * Functions needed to calculate solutions for 1D devices.
 */

#include "ngspice.h"
#include "numglobs.h"
#include "numenum.h"
#include "onedev.h"
#include "onemesh.h"
#include "spmatrix.h"
#include "bool.h"
#include "macros.h"
#include "onedext.h"
#include "oneddefs.h"
#include "cidersupt.h"
#include "../../maths/misc/norm.h"

#include "ifsim.h"
extern IFfrontEnd *SPfrontEnd;



/* Forward Declarations */
void ONEdcSolve(ONEdevice *, int, BOOLEAN, BOOLEAN, ONEtranInfo *);
BOOLEAN ONEpsiDeltaConverged(ONEdevice *, int *);
BOOLEAN ONEdeltaConverged(ONEdevice *);
int ONEnewDelta(ONEdevice *, BOOLEAN , ONEtranInfo *);
void ONEstoreNeutralGuess(ONEdevice *);
void ONEstoreEquilibGuess(ONEdevice *);
void ONEstoreInitialGuess(ONEdevice *pDevice);
double ONEnuNorm(ONEdevice *);
void ONEjacCheck(ONEdevice *, BOOLEAN , ONEtranInfo *);


/* The iteration driving loop and convergence check */
void
ONEdcSolve(ONEdevice *pDevice, int iterationLimit, BOOLEAN newSolver, 
           BOOLEAN tranAnalysis, ONEtranInfo *info)
{
  ONEnode *pNode;
  ONEelem *pElem;
  int index, eIndex, error;
  int timesConverged = 0, negConc = FALSE;
  int size = pDevice->numEqns;
  BOOLEAN quitLoop;
  BOOLEAN debug = FALSE;
  double *rhs = pDevice->rhs;
/*  double *intermediate = pDevice->copiedSolution; */
  double *solution = pDevice->dcSolution;
  double *delta = pDevice->dcDeltaSolution;
  double poissNorm, contNorm;
  double startTime, totalStartTime;
  double totalTime, loadTime, factorTime, solveTime, updateTime, checkTime;
  double orderTime = 0.0;

  quitLoop = FALSE;
  debug =   (!tranAnalysis && ONEdcDebug) 
          ||(tranAnalysis && ONEtranDebug); 
  pDevice->iterationNumber = 0;
  pDevice->converged = FALSE;


  totalTime = loadTime = factorTime = solveTime = updateTime = checkTime = 0.0;
  totalStartTime = SPfrontEnd->IFseconds();

  if (debug) {
    if (pDevice->poissonOnly) {
      fprintf(stdout, "Equilibrium Solution:\n");
    } else {
      fprintf(stdout, "Bias Solution:\n");
    }
    fprintf(stdout, "Iteration  RHS Norm\n");
  }
  while (! (pDevice->converged 
         || pDevice->iterationNumber > iterationLimit
	 || quitLoop)) {
    pDevice->iterationNumber++;

    if ((!pDevice->poissonOnly) && (iterationLimit > 0)
	&&(!tranAnalysis)) {
      ONEjacCheck(pDevice, tranAnalysis, info);
    }
    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    if (pDevice->poissonOnly) {
      ONEQsysLoad(pDevice);
    } else {
      ONE_sysLoad(pDevice, tranAnalysis, info);
    }
    pDevice->rhsNorm = maxNorm(rhs, size);
    loadTime += SPfrontEnd->IFseconds() - startTime;
    if (debug) {
      fprintf(stdout, "%7d   %11.4e%s\n",
	  pDevice->iterationNumber - 1, pDevice->rhsNorm,
	  negConc ? "   negative conc encountered" : "");
      negConc = FALSE;
    }
    /* FACTOR */
    startTime = SPfrontEnd->IFseconds();
    error = spFactor(pDevice->matrix);
    factorTime += SPfrontEnd->IFseconds() - startTime;
    if (newSolver) {
      if (pDevice->iterationNumber == 1) {
	orderTime = factorTime;
      } else if (pDevice->iterationNumber == 2) {
	orderTime -= factorTime - orderTime;
	factorTime -= orderTime;
	if (pDevice->poissonOnly) {
	  pDevice->pStats->orderTime[STAT_SETUP] += orderTime;
	} else {
	  pDevice->pStats->orderTime[STAT_DC] += orderTime;
	}
	newSolver = FALSE;
      }
    }
    if (foundError(error)) {
      if (error == spSINGULAR) {
	int badRow, badCol;
	spWhereSingular(pDevice->matrix, &badRow, &badCol);
	printf("*****  singular at (%d,%d)\n", badRow, badCol);
      }
      /* Should probably try to recover now, but we'll punt instead. */
      exit(-1);
    }
    /* SOLVE */
    startTime = SPfrontEnd->IFseconds();
    spSolve(pDevice->matrix, rhs, delta, NIL(spREAL), NIL(spREAL));
    solveTime += SPfrontEnd->IFseconds() - startTime;

    /* UPDATE */
    startTime = SPfrontEnd->IFseconds();
    /*
     * Use norm reducing Newton method only for DC bias solutions. Since norm
     * reducing can get trapped by numerical errors, turn it off when we are
     * somewhat close to the solution.
     */
    if ((!pDevice->poissonOnly) && (iterationLimit > 0)
	&& (!tranAnalysis) && (pDevice->rhsNorm > 1e-6)) {
      error = ONEnewDelta(pDevice, tranAnalysis, info);
      if (error) {
	pDevice->converged = FALSE;
	quitLoop = TRUE;
	updateTime += SPfrontEnd->IFseconds() - startTime;
	continue;
      }
    }
    for (index = 1; index <= size; index++) {
      solution[index] += delta[index];
    }
    updateTime += SPfrontEnd->IFseconds() - startTime;

    /* CHECK CONVERGENCE */
    startTime = SPfrontEnd->IFseconds();
    /* Check if updates have gotten sufficiently small. */
    if (pDevice->iterationNumber != 1) {
      /*
       * pDevice->converged = ONEdeltaConverged(pDevice);
       */
      pDevice->converged = ONEpsiDeltaConverged(pDevice, &negConc);
    }
    /* Check if the rhs residual is smaller than abstol. */
    if (pDevice->converged &&(!pDevice->poissonOnly)
	&&(!tranAnalysis)) {
      ONE_rhsLoad(pDevice, tranAnalysis, info);
      pDevice->rhsNorm = maxNorm(rhs, size);
      if (pDevice->rhsNorm > pDevice->abstol) {
	pDevice->converged = FALSE;
      }
      if ((++timesConverged >= 2)
	  &&(pDevice->rhsNorm < 1e3 * pDevice->abstol)) {
	pDevice->converged = TRUE;
      } else if (timesConverged >= 5) {
	pDevice->converged = FALSE;
	quitLoop = TRUE;
      }
    } else if (pDevice->converged && pDevice->poissonOnly) {
      ONEQrhsLoad(pDevice);
      pDevice->rhsNorm = maxNorm(rhs, size);
      if (pDevice->rhsNorm > pDevice->abstol) {
	pDevice->converged = FALSE;
      }
      if (++timesConverged >= 5) {
	pDevice->converged = TRUE;
      }
    }
    /* Check if any of the carrier concentrations are negative. */
    if (pDevice->converged &&(!pDevice->poissonOnly)) {
      /* Clear garbage entry since carrier-free elements reference it */
      solution[0] = 0.0;

      for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
	pElem = pDevice->elemArray[eIndex];
	for (index = 0; index <= 1; index++) {
	  if (pElem->evalNodes[index]) {
	    pNode = pElem->pNodes[index];
	    if (solution[pNode->nEqn] < 0.0) {
	      pDevice->converged = FALSE;
	      negConc = TRUE;
	      if (tranAnalysis) {
		quitLoop = TRUE;
	      } else {
		solution[pNode->nEqn] = 0.0;
	      }
	    }
	    if (solution[pNode->pEqn] < 0.0) {
	      pDevice->converged = FALSE;
	      negConc = TRUE;
	      if (tranAnalysis) {
		quitLoop = TRUE;
	      } else {
		solution[pNode->pEqn] = 0.0;
	      }
	    }
	  }
	}
      }
      /* Set to a consistent state if negative conc was encountered. */
      if (!pDevice->converged) {
	ONE_rhsLoad(pDevice, tranAnalysis, info);
	pDevice->rhsNorm = maxNorm(rhs, size);
      }
    }
    checkTime += SPfrontEnd->IFseconds() - startTime;
  }
  totalTime += SPfrontEnd->IFseconds() - totalStartTime;

  if (tranAnalysis) {
    pDevice->pStats->loadTime[STAT_TRAN] += loadTime;
    pDevice->pStats->factorTime[STAT_TRAN] += factorTime;
    pDevice->pStats->solveTime[STAT_TRAN] += solveTime;
    pDevice->pStats->updateTime[STAT_TRAN] += updateTime;
    pDevice->pStats->checkTime[STAT_TRAN] += checkTime;
    pDevice->pStats->numIters[STAT_TRAN] += pDevice->iterationNumber;
  } else if (pDevice->poissonOnly) {
    pDevice->pStats->loadTime[STAT_SETUP] += loadTime;
    pDevice->pStats->factorTime[STAT_SETUP] += factorTime;
    pDevice->pStats->solveTime[STAT_SETUP] += solveTime;
    pDevice->pStats->updateTime[STAT_SETUP] += updateTime;
    pDevice->pStats->checkTime[STAT_SETUP] += checkTime;
    pDevice->pStats->numIters[STAT_SETUP] += pDevice->iterationNumber;
  } else {
    pDevice->pStats->loadTime[STAT_DC] += loadTime;
    pDevice->pStats->factorTime[STAT_DC] += factorTime;
    pDevice->pStats->solveTime[STAT_DC] += solveTime;
    pDevice->pStats->updateTime[STAT_DC] += updateTime;
    pDevice->pStats->checkTime[STAT_DC] += checkTime;
    pDevice->pStats->numIters[STAT_DC] += pDevice->iterationNumber;
  }

  if (debug) {
    if (!tranAnalysis) {
      pDevice->rhsNorm = maxNorm(rhs, size);
      fprintf(stdout, "%7d   %11.4e%s\n",
	pDevice->iterationNumber, pDevice->rhsNorm,
	negConc ? "   negative conc in solution" : "");
    }
    if (pDevice->converged) {
      if (!pDevice->poissonOnly) {
	rhs[0] = 0.0;
	poissNorm = contNorm = 0.0;
	for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
	  pElem = pDevice->elemArray[eIndex];
	  for (index = 0; index <= 1; index++) {
	    if (pElem->evalNodes[index]) {
	      pNode = pElem->pNodes[index];
	      poissNorm = MAX(poissNorm, ABS(rhs[pNode->psiEqn]));
	      contNorm = MAX(contNorm, ABS(rhs[pNode->nEqn]));
	      contNorm = MAX(contNorm, ABS(rhs[pNode->pEqn]));
	    }
	  }
	}
	fprintf(stdout,
	    "Residual: %11.4e C/um^2 poisson, %11.4e A/um^2 continuity\n",
	    poissNorm * EpsNorm * ENorm * 1e-8,
	    contNorm * JNorm * 1e-8);
      } else {
	fprintf(stdout, "Residual: %11.4e C/um^2 poisson\n",
	    pDevice->rhsNorm * EpsNorm * ENorm * 1e-8);
      }
    }
  }
}

/*
 * A function that checks convergence based on the convergence of the quasi
 * Fermi levels. In theory, this should work better than the one currently
 * being used since we are always looking at potentials: (psi, phin, phip).
 */
BOOLEAN
ONEpsiDeltaConverged(ONEdevice *pDevice, int *pNegConc)
{
  int index, nIndex, eIndex;
  ONEnode *pNode;
  ONEelem *pElem;
  double xOld, xNew, xDelta, tol;
  double psi, newPsi, nConc, pConc, newN, newP;
  double phiN, phiP, newPhiN, newPhiP;
  BOOLEAN converged = TRUE;

  /* Equilibrium solution. */
  if (pDevice->poissonOnly) {
    for (index = 1; index <= pDevice->numEqns; index++) {
      xOld = pDevice->dcSolution[index];
      xDelta = pDevice->dcDeltaSolution[index];
      xNew = xOld + xDelta;
      tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
      if (ABS(xDelta) > tol) {
	converged = FALSE;
	goto done;
      }
    }
    return (converged); 

  }
  /* Bias solution. Check convergence on psi, phin, phip. */
  for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
    pElem = pDevice->elemArray[eIndex];
    for (nIndex = 0; nIndex <= 1; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	if (pNode->nodeType != CONTACT) {
	  /* check convergence on psi */
	  xOld = pDevice->dcSolution[pNode->psiEqn];
	  xDelta = pDevice->dcDeltaSolution[pNode->psiEqn];
	  xNew = xOld + xDelta;
	  tol = pDevice->abstol +
	      pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
	  if (ABS(xDelta) > tol) {
	    converged = FALSE;
	    goto done;
	  }
	}
	/* check convergence on phin and phip */
	if (pElem->elemType == SEMICON && pNode->nodeType != CONTACT) {
	  psi = pDevice->dcSolution[pNode->psiEqn];
	  nConc = pDevice->dcSolution[pNode->nEqn];
	  pConc = pDevice->dcSolution[pNode->pEqn];
	  newPsi = psi + pDevice->dcDeltaSolution[pNode->psiEqn];
	  newN = nConc + pDevice->dcDeltaSolution[pNode->nEqn];
	  newP = pConc + pDevice->dcDeltaSolution[pNode->pEqn];
	  if (newN <= 0.0 || newP <= 0.0) {
	    *pNegConc = TRUE;
	    converged = FALSE;
	    goto done;
	  }
	  phiN = psi - log(nConc / pNode->nie);
	  phiP = psi + log(pConc / pNode->nie);
	  newPhiN = newPsi - log(newN / pNode->nie);
	  newPhiP = newPsi + log(newP / pNode->nie);
	  tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiN), ABS(newPhiN));
	  if (ABS(newPhiN - phiN) > tol) {
	    converged = FALSE;
	    goto done;
	  }