twosolve.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

/**********
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
**********/

#include <math.h>
#include "numglobs.h"
#include "numenum.h"
#include "nummacs.h"
#include "twodev.h"
#include "twomesh.h"
#include "spmatrix.h"

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

/* External Declarations */
extern void TWOQjacBuild(), TWOQcommonTerms(), TWOQrhsLoad(), TWOQsysLoad();
extern void TWO_jacBuild(), TWO_commonTerms(), TWO_rhsLoad(), TWO_sysLoad();
extern void TWONjacBuild(), TWONcommonTerms(), TWONrhsLoad(), TWONsysLoad();
extern void TWOPjacBuild(), TWOPcommonTerms(), TWOPrhsLoad(), TWOPsysLoad();
extern BOOLEAN foundError();
extern double oneNorm(), l2Norm(), maxNorm(), predict(); 

/* Forward Declarations */
void TWOdcSolve();
BOOLEAN TWOdeltaConverged();
int TWOnewDelta();
void TWOstoreNeutralGuess();
void TWOstoreEquilibGuess();
void TWOstoreInitialGuess();
double TWOnuNorm();
void TWOstoreInitialGuess();
void TWOjacCheck();

/* functions to calculate the 2D solutions */
BOOLEAN TWOdcDebug = TRUE;
BOOLEAN TWOtranDebug = TRUE;
BOOLEAN TWOjacDebug = FALSE;

/* the iteration driving loop and convergence check */
void
TWOdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info)
  TWOdevice *pDevice;
  int iterationLimit;
  BOOLEAN newSolver;
  BOOLEAN tranAnalysis;
  TWOtranInfo *info;
{
  TWOnode *pNode;
  TWOelem *pElem;
  int size = pDevice->numEqns;
  int index, eIndex, error;
  int timesConverged = 0;
  BOOLEAN quitLoop;
  BOOLEAN debug;
  BOOLEAN negConc = FALSE;
  double *rhs = pDevice->rhs;
  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;

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

  quitLoop = FALSE;
  debug = (NOT tranAnalysis AND TWOdcDebug) OR (tranAnalysis AND TWOtranDebug);
  pDevice->iterationNumber = 0;
  pDevice->converged = FALSE;

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

    if ((NOT pDevice->poissonOnly) AND(iterationLimit > 0)
	AND(NOT tranAnalysis)) {
      TWOjacCheck(pDevice, tranAnalysis, info);
    }

    /* LOAD */
    startTime = SPfrontEnd->IFseconds();
    if (pDevice->poissonOnly) {
      TWOQsysLoad(pDevice);
    } else if (NOT OneCarrier) {
      TWO_sysLoad(pDevice, tranAnalysis, info);
    } else if (OneCarrier IS N_TYPE) {
      TWONsysLoad(pDevice, tranAnalysis, info);
    } else if (OneCarrier IS P_TYPE) {
      TWOPsysLoad(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);
      }
      pDevice->converged = FALSE;
      quitLoop = TRUE;
      continue;
    }

    /* 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 for DC bias solutions only. */
    if ((NOT pDevice->poissonOnly) AND (iterationLimit > 0)
	AND (NOT tranAnalysis) AND (pDevice->rhsNorm > 1e-1)) {
      error = TWOnewDelta(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 ISNOT 1) {
      pDevice->converged = TWOdeltaConverged(pDevice);
    }
    /* Check if the residual is smaller than abstol. */
    if (pDevice->converged AND (NOT pDevice->poissonOnly)
	AND (NOT tranAnalysis)) {
      if (NOT OneCarrier) {
	TWO_rhsLoad(pDevice, tranAnalysis, info);
      } else if (OneCarrier IS N_TYPE) {
	TWONrhsLoad(pDevice, tranAnalysis, info);
      } else if (OneCarrier IS P_TYPE) {
	TWOPrhsLoad(pDevice, tranAnalysis, info);
      }
      pDevice->rhsNorm = maxNorm(rhs, size);
      if (pDevice->rhsNorm > pDevice->abstol) {
	pDevice->converged = FALSE;
      }
      if ((++timesConverged >= 2)
	  AND (pDevice->rhsNorm < 1e3 * pDevice->abstol)) {
	pDevice->converged = TRUE;
      } else if (timesConverged >= 5) {
	pDevice->converged = FALSE;
	quitLoop = TRUE;
	continue;
      }
    } else if (pDevice->converged AND pDevice->poissonOnly) {
      TWOQrhsLoad(pDevice);
      pDevice->rhsNorm = maxNorm(rhs, size);
      if (pDevice->rhsNorm > pDevice->abstol) {
	pDevice->converged = FALSE;
      }
      if (++timesConverged >= 5) {
	pDevice->converged = TRUE;
      }
    }
    /* Check if the carrier concentrations are negative. */
    if (pDevice->converged AND (NOT pDevice->poissonOnly)) {
      /* Clear garbage entry since carrier-free elements reference it. */
      solution[0] = 0.0;
      for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
	pElem = pDevice->elements[eIndex];
	for (index = 0; index <= 3; 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;
	      }
	    }
	  }
	}
      }
      if (NOT pDevice->converged) {
	if (NOT OneCarrier) {
	  TWO_rhsLoad(pDevice, tranAnalysis, info);
	} else if (OneCarrier IS N_TYPE) {
	  TWONrhsLoad(pDevice, tranAnalysis, info);
	} else if (OneCarrier IS P_TYPE) {
	  TWOPrhsLoad(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 (NOT 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 (NOT pDevice->poissonOnly) {
	poissNorm = contNorm = 0.0;
	rhs[0] = 0.0;		/* Make sure garbage entry is clear. */
	for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
	  pElem = pDevice->elements[eIndex];
	  for (index = 0; index <= 3; 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 poisson, %11.4e A/um continuity\n",
	    poissNorm * EpsNorm * VNorm * 1e-4,
	    contNorm * JNorm * LNorm * 1e-4);
      } else {
	fprintf(stdout, "Residual: %11.4e C/um poisson\n",
	    pDevice->rhsNorm * EpsNorm * VNorm * 1e-4);
      }
    }
  }
}

BOOLEAN
TWOdeltaConverged(pDevice)
  TWOdevice *pDevice;
{
  /* This function returns a TRUE if converged else a FALSE. */
  int index;
  BOOLEAN converged = TRUE;
  double xOld, xNew, tol;

  for (index = 1; index <= pDevice->numEqns; index++) {
    xOld = pDevice->dcSolution[index];
    xNew = xOld + pDevice->dcDeltaSolution[index];
    tol = pDevice->abstol + pDevice->reltol * MAX(ABS(xOld), ABS(xNew));
    if (ABS(xOld - xNew) > tol) {
      converged = FALSE;
      break;
    }
  }
  return (converged);
}

BOOLEAN
TWOdeviceConverged(pDevice)
  TWOdevice *pDevice;
{
  int index, eIndex;
  BOOLEAN converged = TRUE;
  double *solution = pDevice->dcSolution;
  TWOnode *pNode;
  TWOelem *pElem;
  double startTime;

  /* If the update is sufficiently small, and the carrier concentrations
   * are all positive, then return TRUE, else return FALSE.
   */

  /* CHECK CONVERGENCE */
  startTime = SPfrontEnd->IFseconds();
  if ((converged = TWOdeltaConverged(pDevice)) == TRUE) {
    for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
      pElem = pDevice->elements[eIndex];
      for (index = 0; index <= 3; index++) {
	if (pElem->evalNodes[index]) {
	  pNode = pElem->pNodes[index];
	  if (pNode->nEqn != 0 AND solution[pNode->nEqn] < 0.0) {
	    converged = FALSE;
	    solution[pNode->nEqn] = pNode->nConc = 0.0;
	  }
	  if (pNode->pEqn != 0 AND solution[pNode->pEqn] < 0.0) {
	    converged = FALSE;
	    solution[pNode->pEqn] = pNode->pConc = 0.0;
	  }
	}
      }
    }
  }
  pDevice->pStats->checkTime[STAT_TRAN] += SPfrontEnd->IFseconds() - startTime;

  return (converged);
}

void
TWOresetJacobian(pDevice)
  TWOdevice *pDevice;
{
  void TWO_jacLoad(), TWONjacLoad(), TWOPjacLoad();
  int error;
  BOOLEAN foundError();

  if (NOT OneCarrier) {
    TWO_jacLoad(pDevice);
  } else if (OneCarrier IS N_TYPE) {
    TWONjacLoad(pDevice);
  } else if (OneCarrier IS P_TYPE) {
    TWOPjacLoad(pDevice);
  } else {
    printf("TWOresetJacobian: unknown carrier type\n");
    exit(-1);
  }
  error = spFactor(pDevice->matrix);
  if (foundError(error)) {
    exit(-1);
  }
}

/*
 * Compute the device state assuming charge neutrality exists everywhere in
 * the device.  In insulators, where there is no charge, assign the potential
 * at half the insulator band gap (refPsi).
 */
void
TWOstoreNeutralGuess(pDevice)
  TWOdevice *pDevice;
{
  int nIndex, eIndex;
  TWOelem *pElem;
  TWOnode *pNode;
  double nie, conc, absConc, refPsi, psi, ni, pi, sign;

  /* assign the initial guess for Poisson's equation */
  for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    refPsi = pElem->matlInfo->refPsi;
    if (pElem->elemType IS INSULATOR) {
      for (nIndex = 0; nIndex <= 3; nIndex++) {
	if (pElem->evalNodes[nIndex]) {
	  pNode = pElem->pNodes[nIndex];
	  if (pNode->nodeType IS CONTACT) {
	    /*
	     * a metal contact to insulator domain so use work function
	     * difference as the value of psi
	     */
	    pNode->psi = RefPsi - pNode->eaff;