onesolve.c

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

	  if (ABS(newPhiP - phiP) > tol) {
	    converged = FALSE;
	    goto done;
	  }
	}
      }
    }
  }
done:
  return (converged);
}

/*
 * See if the update to the solution is small enough. Returns TRUE if it is.
 */
BOOLEAN
ONEdeltaConverged(pDevice)
  ONEdevice *pDevice;
{
  int index;
  BOOLEAN converged = TRUE;
  double *solution = pDevice->dcSolution;
  double *delta = pDevice->dcDeltaSolution;
  double xOld, xNew, tol;

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

/*
 * See if the update to the solution is small enough and the solution is
 * physically reasonable. Returns TRUE if it is.
 */
BOOLEAN
ONEdeviceConverged(pDevice)
  ONEdevice *pDevice;
{
  int index, eIndex;
  BOOLEAN converged = TRUE;
  double *solution = pDevice->dcSolution;
  ONEnode *pNode;
  ONEelem *pElem;
  double startTime;

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

  /* CHECK CONVERGENCE */
  startTime = SPfrontEnd->IFseconds();
  if ((converged = ONEdeltaConverged(pDevice)) == TRUE) {
    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 (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);
}

/*
 * Load and factor the Jacobian so that it is consistent with the current
 * solution.
 */
void
ONEresetJacobian(pDevice)
  ONEdevice *pDevice;
{
  int error;

  ONE_jacLoad(pDevice);
  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
ONEstoreNeutralGuess(pDevice)
  ONEdevice *pDevice;
{
  int nIndex, eIndex;
  ONEelem *pElem;
  ONEnode *pNode;
  double nie, conc, absConc, refPsi, psi, ni, pi, sign;

  for (eIndex = 1; eIndex < pDevice->numNodes; eIndex++) {
    pElem = pDevice->elemArray[eIndex];
    refPsi = pElem->matlInfo->refPsi;
    if (pElem->elemType IS INSULATOR) {
      for (nIndex = 0; nIndex <= 1; nIndex++) {
	if (pElem->evalNodes[nIndex]) {
	  pNode = pElem->pNodes[nIndex];
	  if (pNode->nodeType IS CONTACT) {
	    /*
	     * A metal contact to insulating domain, so use work function
	     * difference as the value of psi.
	     */
	    pNode->psi = RefPsi - pNode->eaff;
	  } else {
	    pNode->psi = refPsi;
	  }
	}
      }
    }
    if (pElem->elemType IS SEMICON) {
      for (nIndex = 0; nIndex <= 1; nIndex++) {
	if (pElem->evalNodes[nIndex]) {
	  pNode = pElem->pNodes[nIndex];
	  nie = pNode->nie;
	  conc = pNode->netConc / nie;
	  psi = 0.0;
	  ni = nie;
	  pi = nie;
	  sign = SGN(conc);
	  absConc = ABS(conc);
	  if (conc ISNOT 0.0) {
	    psi = sign * log(0.5 * absConc + sqrt(1.0 + 0.25 * absConc * absConc));
	    ni = nie * exp(psi);
	    pi = nie * exp(-psi);
	    if (FreezeOut) {
	      /* Use Newton's Method to solve for psi. */
	      int ctr, maxiter = 10;
	      double rhs, deriv, fNa, fNd, fdNa, fdNd;
	      for (ctr = 0; ctr < maxiter; ctr++) {
		pNode->nConc = ni;
		pNode->pConc = pi;
		ONEQfreezeOut(pNode, &fNd, &fNa, &fdNd, &fdNa);
		rhs = pi - ni + pNode->nd * fNd - pNode->na * fNa;
		deriv = pi + ni - pNode->nd * fdNd + pNode->na * fdNa;
		psi += rhs / deriv;
		ni = nie * exp(psi);
		pi = nie * exp(-psi);
	      }
	    }
	  }
	  pNode->psi = refPsi + psi;
	  pNode->psi0 = pNode->psi;
	  pNode->vbe = refPsi;
	  pNode->nConc = ni;
	  pNode->pConc = pi;
	  /* Now store the initial guess in the dc solution vector. */
	  pDevice->dcSolution[pNode->poiEqn] = pNode->psi;
	}
      }
    }
  }
}

/*
 * Compute the device state at thermal equilibrium. This state is equal to
 * the solution of just Poisson's equation. The charge-neutral solution is
 * taken as an initial guess.
 */
void
ONEequilSolve(pDevice)
  ONEdevice *pDevice;
{
  BOOLEAN newSolver = FALSE;
  int error;
  int nIndex, eIndex;
  ONEelem *pElem;
  ONEnode *pNode;
  double startTime, setupTime, miscTime;

  setupTime = miscTime = 0.0;

  /* SETUP */
  startTime = SPfrontEnd->IFseconds();
  switch (pDevice->solverType) {
  case SLV_SMSIG:
  case SLV_BIAS:
    /* free up memory allocated for the bias solution */
    FREE(pDevice->dcSolution);
    FREE(pDevice->dcDeltaSolution);
    FREE(pDevice->copiedSolution);
    FREE(pDevice->rhs);
    FREE(pDevice->rhsImag);
    spDestroy(pDevice->matrix);
    /* FALLTHRU */
  case SLV_NONE:
    pDevice->poissonOnly = TRUE;
    pDevice->numEqns = pDevice->dimEquil - 1;
    ALLOC(pDevice->dcSolution, double, pDevice->dimEquil);
    ALLOC(pDevice->dcDeltaSolution, double, pDevice->dimEquil);
    ALLOC(pDevice->copiedSolution, double, pDevice->dimEquil);
    ALLOC(pDevice->rhs, double, pDevice->dimEquil);
    pDevice->matrix = spCreate(pDevice->numEqns, 0, &error);
    if (error IS spNO_MEMORY) {
      printf("ONEequilSolve: Out of Memory\n");
      exit(-1);
    }
    newSolver = TRUE;
    spSetReal(pDevice->matrix);
    ONEQjacBuild(pDevice);
    pDevice->numOrigEquil = spElementCount(pDevice->matrix);
    pDevice->numFillEquil = 0;
    /* FALLTHRU */
  case SLV_EQUIL:
    pDevice->solverType = SLV_EQUIL;
    break;
  default:
    fprintf(stderr, "Panic: Unknown solver type in equil solution.\n");
    exit(-1);
    break;
  }
  ONEstoreNeutralGuess(pDevice);
  setupTime += SPfrontEnd->IFseconds() - startTime;

  /* SOLVE */
  ONEdcSolve(pDevice, MaxIterations, newSolver, FALSE, (ONEtranInfo *) NULL);

  /* MISCELLANEOUS */
  startTime = SPfrontEnd->IFseconds();
  if (newSolver) {
    pDevice->numFillEquil = spFillinCount(pDevice->matrix);
  }
  if (pDevice->converged) {
    ONEQcommonTerms(pDevice);

    /* Save equilibrium potential. */
    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];
	  pNode->psi0 = pNode->psi;
	}
      }
    }
  } else {
    printf("ONEequilSolve: No Convergence\n");
  }
  miscTime += SPfrontEnd->IFseconds() - startTime;
  pDevice->pStats->setupTime[STAT_SETUP] += setupTime;
  pDevice->pStats->miscTime[STAT_SETUP] += miscTime;
}

/*
 * Compute the device state under an applied bias. The equilibrium solution
 * is taken as an initial guess the first time this is called.
 */
void
ONEbiasSolve(pDevice, iterationLimit, tranAnalysis, info)
  ONEdevice *pDevice;
  int iterationLimit;
  BOOLEAN tranAnalysis;
  ONEtranInfo *info;
{
  BOOLEAN newSolver = FALSE;
  int error;
  int index, eIndex;
  double *solution;
  ONEelem *pElem;
  ONEnode *pNode;
  double startTime, setupTime, miscTime;

  setupTime = miscTime = 0.0;

  /* SETUP */
  startTime = SPfrontEnd->IFseconds();
  switch (pDevice->solverType) {
  case SLV_EQUIL:
    /* Free up the vectors allocated in the equilibrium solution. */
    FREE(pDevice->dcSolution);
    FREE(pDevice->dcDeltaSolution);
    FREE(pDevice->copiedSolution);
    FREE(pDevice->rhs);
    spDestroy(pDevice->matrix);
    /* FALLTHRU */
  case SLV_NONE:
    pDevice->poissonOnly = FALSE;
    pDevice->numEqns = pDevice->dimBias - 1;
    ALLOC(pDevice->dcSolution, double, pDevice->dimBias);
    ALLOC(pDevice->dcDeltaSolution, double, pDevice->dimBias);
    ALLOC(pDevice->copiedSolution, double, pDevice->dimBias);
    ALLOC(pDevice->rhs, double, pDevice->dimBias);
    ALLOC(pDevice->rhsImag, double, pDevice->dimBias);
    pDevice->matrix = spCreate(pDevice->numEqns, 1, &error);
    if (error IS spNO_MEMORY) {
      printf("ONEbiasSolve: Out of Memory\n");
      exit(-1);
    }
    newSolver = TRUE;
    ONE_jacBuild(pDevice);
    pDevice->numOrigBias = spElementCount(pDevice->matrix);
    pDevice->numFillBias = 0;
    ONEstoreInitialGuess(pDevice);
    /* FALLTHRU */
  case SLV_SMSIG:
    spSetReal(pDevice->matrix);
    /* FALLTHRU */
  case SLV_BIAS:
    pDevice->solverType = SLV_BIAS;
    break;
  default:
    fprintf(stderr, "Panic: Unknown solver type in bias solution.\n");
    exit(-1);
    break;
  }
  setupTime += SPfrontEnd->IFseconds() - startTime;

  /* SOLVE */
  ONEdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info);

  /* MISCELLANEOUS */
  startTime = SPfrontEnd->IFseconds();
  if (newSolver) {
    pDevice->numFillBias = spFillinCount(pDevice->matrix);
  }
  solution = pDevice->dcSolution;
  if ((NOT pDevice->converged) AND iterationLimit > 1) {
    printf("ONEbiasSolve: No Convergence\n");
  } else if (pDevice->converged) {
    /* Update the nodal quantities. */
    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 (pNode->psiEqn != 0) {
	    pNode->psi = solution[pNode->psiEqn];
	  }
	  if (pNode->nEqn != 0) {
	    pNode->nConc = solution[pNode->nEqn];
	  }
	  if (pNode->pEqn != 0) {
	    pNode->pConc = solution[pNode->pEqn];
	  }
	}
      }
    }
    /* Update the current terms. */
    ONE_commonTerms(pDevice, FALSE, tranAnalysis, info);
  } else if (iterationLimit <= 1) {
    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 (pNode->nodeType ISNOT CONTACT) {
	    pNode->psi = solution[pNode->psiEqn];