twosolve.c

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

}


#define NORM_RED_MAXITERS 10

int
TWOnewDelta(pDevice, tranAnalysis, info)
  TWOdevice *pDevice;
  BOOLEAN tranAnalysis;
  TWOtranInfo *info;
{
  int index, iterNum = 0;
  double newNorm, origNorm;
  double fib, lambda, fibn, fibp, maxNorm();
  BOOLEAN acceptable = FALSE, error = FALSE;

  lambda = 1.0;
  fibn = 1.0;
  fibp = 1.0;

  /*
   * copy the contents of dcSolution into copiedSolution and modify
   * dcSolution by adding the deltaSolution
   */

  for (index = 1; index <= pDevice->numEqns; ++index) {
    pDevice->copiedSolution[index] = pDevice->dcSolution[index];
    pDevice->dcSolution[index] += pDevice->dcDeltaSolution[index];
  }

  if (pDevice->poissonOnly) {
    TWOQrhsLoad(pDevice);
  } else 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);
  }
  newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
  if (pDevice->rhsNorm <= pDevice->abstol) {
    lambda = 0.0;
    newNorm = pDevice->rhsNorm;
  } else if (newNorm < pDevice->rhsNorm) {
    acceptable = TRUE;
  } else {
    /* chop the step size so that deltax is acceptable */
    if (TWOdcDebug) {
      fprintf(stdout, "          %11.4e  %11.4e\n",
	  newNorm, lambda);
    }
    while (NOT acceptable) {
      iterNum++;

      if (iterNum > NORM_RED_MAXITERS) {
	error = TRUE;
	lambda = 0.0;
	/* Don't break out until after we've reset the device. */
      }
      fib = fibp;
      fibp = fibn;
      fibn += fib;
      lambda *= (fibp / fibn);

      for (index = 1; index <= pDevice->numEqns; index++) {
	pDevice->dcSolution[index] = pDevice->copiedSolution[index]
	    + lambda * pDevice->dcDeltaSolution[index];
      }
      if (pDevice->poissonOnly) {
	TWOQrhsLoad(pDevice);
      } else 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);
      }
      newNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
      if (error) {
	break;
      }
      if (newNorm <= pDevice->rhsNorm) {
	acceptable = TRUE;
      }
      if (TWOdcDebug) {
	fprintf(stdout, "          %11.4e  %11.4e\n",
	    newNorm, lambda);
      }
    }
  }
  /* restore the previous dcSolution and the new deltaSolution */
  pDevice->rhsNorm = newNorm;
  for (index = 1; index <= pDevice->numEqns; index++) {
    pDevice->dcSolution[index] = pDevice->copiedSolution[index];
    pDevice->dcDeltaSolution[index] *= lambda;
  }
  return (error);
}


/* Predict the values of the internal variables at the next timepoint. */
/* Needed for Predictor-Corrector LTE estimation */
void
TWOpredict(pDevice, info)
  TWOdevice *pDevice;
  TWOtranInfo *info;
{
  int nIndex, eIndex;
  TWOnode *pNode;
  TWOelem *pElem;
  double startTime, miscTime = 0.0;

  /* TRANSIENT MISC */
  startTime = SPfrontEnd->IFseconds();
  for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    for (nIndex = 0; nIndex <= 3; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	pNode->psi = *(pDevice->devState1 + pNode->nodePsi);
	if (pElem->elemType IS SEMICON AND pNode->nodeType ISNOT CONTACT) {
	  if (NOT OneCarrier) {
	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
	  } else if (OneCarrier IS N_TYPE) {
	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
	    pNode->pPred = *(pDevice->devState1 + pNode->nodeP);
	  } else if (OneCarrier IS P_TYPE) {
	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
	    pNode->nPred = *(pDevice->devState1 + pNode->nodeN);
	  }
	  pNode->nConc = pNode->nPred;
	  pNode->pConc = pNode->pPred;
	}
      }
    }
  }
  miscTime += SPfrontEnd->IFseconds() - startTime;
  pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
}


/* Estimate the device's overall truncation error. */
double
TWOtrunc(pDevice, info, delta)
  TWOdevice *pDevice;
  TWOtranInfo *info;
  double delta;
{
  int nIndex, eIndex;
  TWOelem *pElem;
  TWOnode *pNode;
  double tolN, tolP, lte, relError, temp, relLTE;
  double lteCoeff = info->lteCoeff;
  double mult = 10.0;
  double reltol;
  double startTime, lteTime = 0.0;

  /* TRANSIENT LTE */
  startTime = SPfrontEnd->IFseconds();

  relError = 0.0;
  reltol = pDevice->reltol * mult;

  /* need to get the predictor for the current order */
  /* the scheme currently used is very dumb. need to fix later */
  computePredCoeff(info->method, info->order, info->predCoeff, info->delta);

  for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    for (nIndex = 0; nIndex <= 3; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	if (pElem->elemType IS SEMICON AND pNode->nodeType ISNOT CONTACT) {
	  if (NOT OneCarrier) {
	    tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
	    tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
	    lte = lteCoeff * (pNode->nConc - pNode->nPred);
	    temp = lte / tolN;
	    relError += temp * temp;
	    lte = lteCoeff * (pNode->pConc - pNode->pPred);
	    temp = lte / tolP;
	    relError += temp * temp;
	  } else if (OneCarrier IS N_TYPE) {
	    tolN = pDevice->abstol + reltol * ABS(pNode->nConc);
	    pNode->nPred = predict(pDevice->devStates, info, pNode->nodeN);
	    lte = lteCoeff * (pNode->nConc - pNode->nPred);
	    temp = lte / tolN;
	    relError += temp * temp;
	  } else if (OneCarrier IS P_TYPE) {
	    tolP = pDevice->abstol + reltol * ABS(pNode->pConc);
	    pNode->pPred = predict(pDevice->devStates, info, pNode->nodeP);
	    lte = lteCoeff * (pNode->pConc - pNode->pPred);
	    temp = lte / tolP;
	    relError += temp * temp;
	  }
	}
      }
    }
  }

  relError = MAX(pDevice->abstol, relError);	/* make sure it is non zero */

  /* the total relative error has been calculated norm-2 squared */
  relError = sqrt(relError / pDevice->numEqns);

  /* depending on the order of the integration method compute new delta */
  temp = delta / pow(relError, 1.0 / (info->order + 1));

  lteTime += SPfrontEnd->IFseconds() - startTime;
  pDevice->pStats->lteTime += lteTime;

  /* return the new delta (stored as temp) */
  return (temp);
}

/* Save info from state table into the internal state. */
void
TWOsaveState(pDevice)
  TWOdevice *pDevice;
{
  int nIndex, eIndex;
  TWOnode *pNode;
  TWOelem *pElem;

  for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    for (nIndex = 0; nIndex <= 3; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	pNode->psi = *(pDevice->devState1 + pNode->nodePsi);
	if (pElem->elemType IS SEMICON AND pNode->nodeType ISNOT CONTACT) {
	  pNode->nConc = *(pDevice->devState1 + pNode->nodeN);
	  pNode->pConc = *(pDevice->devState1 + pNode->nodeP);
	}
      }
    }
  }
}

/*
 * A function that checks convergence based on the convergence of the quasi
 * Fermi levels. This should be better since we are looking at potentials in
 * all (psi, phin, phip)
 */
BOOLEAN
TWOpsiDeltaConverged(pDevice)
  TWOdevice *pDevice;
{
  int index, nIndex, eIndex;
  TWOnode *pNode;
  TWOelem *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->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    for (nIndex = 0; nIndex <= 3; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	if (pNode->nodeType ISNOT 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 IS SEMICON AND pNode->nodeType ISNOT 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];
	  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;
	  }
	  tol = pDevice->abstol + pDevice->reltol * MAX(ABS(phiP), ABS(newPhiP));
	  if (ABS(newPhiP - phiP) > tol) {
	    converged = FALSE;
	    goto done;
	  }
	}
      }
    }
  }
done:
  return (converged);
}


/* Function to compute Nu norm of the rhs vector. */
/* nu-Norm calculation based upon work done at Stanford. */
double
TWOnuNorm(pDevice)
  TWOdevice *pDevice;
{
  double norm = 0.0;
  double temp;
  int index;

  /* the LU Decomposed matrix is available. use it to calculate x */

  spSolve(pDevice->matrix, pDevice->rhs, pDevice->rhsImag,
      NIL(spREAL), NIL(spREAL));

  /* the solution is in the rhsImag vector */
  /* compute L2-norm of the rhsImag vector */

  for (index = 1; index <= pDevice->numEqns; index++) {
    temp = pDevice->rhsImag[index];
    norm += temp * temp;
  }
  norm = sqrt(norm);
  return (norm);
}

/*
 * Check for numerical errors in the Jacobian.  Useful debugging aid when new
 * models are being incorporated.
 */
void
TWOjacCheck(pDevice, tranAnalysis, info)
  TWOdevice *pDevice;
  BOOLEAN tranAnalysis;
  TWOtranInfo *info;
{
  int index, rIndex;
  double del, diff, tol, *dptr;
  double *spFindElement();

  if (TWOjacDebug) {
    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);
    }
    /*
     * spPrint( pDevice->matrix, 0, 1, 1 );
     */
    pDevice->rhsNorm = maxNorm(pDevice->rhs, pDevice->numEqns);
    for (index = 1; index <= pDevice->numEqns; index++) {
      if (1e3 * ABS(pDevice->rhs[index]) > pDevice->rhsNorm) {
	fprintf(stderr, "eqn %d: res %11.4e, norm %11.4e\n",
	    index, pDevice->rhs[index], pDevice->rhsNorm);
      }
    }
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->rhsImag[index] = pDevice->rhs[index];
    }
    for (index = 1; index <= pDevice->numEqns; index++) {
      pDevice->copiedSolution[index] = pDevice->dcSolution[index];
      del = 1e-4 * pDevice->abstol + 1e-6 * ABS(pDevice->dcSolution[index]);
      pDevice->dcSolution[index] += del;
      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->dcSolution[index] = pDevice->copiedSolution[index];
      for (rIndex = 1; rIndex <= pDevice->numEqns; rIndex++) {
	diff = (pDevice->rhsImag[rIndex] - pDevice->rhs[rIndex]) / del;
	dptr = spFindElement(pDevice->matrix, rIndex, index);
	if (dptr ISNOT NIL(double)) {
	  tol = (1e-4 * pDevice->abstol) + (1e-2 * MAX(ABS(diff), ABS(*dptr)));
	  if ((diff != 0.0) && (ABS(diff - *dptr) > tol)) {
	    fprintf(stderr, "Mismatch[%d][%d]: FD = %11.4e, AJ = %11.4e\n\t FD-AJ = %11.4e vs. %11.4e\n",
		rIndex, index, diff, *dptr,
		ABS(diff - *dptr), tol);
	  }
	} else {
	  if (diff != 0.0) {
	    fprintf(stderr, "Missing [%d][%d]: FD = %11.4e, AJ = 0.0\n",
		rIndex, index, diff);
	  }
	}
      }
    }
  }
}