twosolve.c

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

	  } else {
	    pNode->psi = refPsi;
	  }
	}
      }
    }
    if (pElem->elemType IS SEMICON) {
      for (nIndex = 0; nIndex <= 3; nIndex++) {
	if (pElem->evalNodes[nIndex]) {
	  pNode = pElem->pNodes[nIndex];
	  /* silicon nodes */
	  nie = pNode->nie;
	  conc = pNode->netConc / pNode->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);
	  }
	  pNode->psi = refPsi + psi;
	  pNode->nConc = ni;
	  pNode->pConc = pi;
	  /* store the initial guess in the dc solution vector */
	  if (pNode->nodeType ISNOT CONTACT) {
	    pDevice->dcSolution[pNode->poiEqn] = pNode->psi;
	  }
	}
      }
    }
  }
}

/* computing the equilibrium solution; solution of Poisson's eqn */
/* the solution is used as an initial starting point for bias conditions */

void
TWOequilSolve(pDevice)
  TWOdevice *pDevice;
{
  BOOLEAN newSolver = FALSE;
  int error;
  int nIndex, eIndex;
  TWOelem *pElem;
  TWOnode *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);
  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("TWOequilSolve: Out of Memory\n");
      exit(-1);
    }
    newSolver = TRUE;
    spSetReal(pDevice->matrix);
    TWOQjacBuild(pDevice);
    pDevice->numOrigEquil = spElementCount(pDevice->matrix);
    pDevice->numFillEquil = 0;
  case SLV_EQUIL:
    pDevice->solverType = SLV_EQUIL;
    break;
  default:
    fprintf(stderr, "Panic: Unknown solver type in equil solution.\n");
    exit(-1);
    break;
  }
  TWOstoreNeutralGuess(pDevice);
  setupTime += SPfrontEnd->IFseconds() - startTime;

  /* SOLVE */
  TWOdcSolve(pDevice, MaxIterations, newSolver, FALSE, (TWOtranInfo *) NULL);

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

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

/* compute the solution under an applied bias */
/* the equilibrium solution is taken as an initial guess */
void
TWObiasSolve(pDevice, iterationLimit, tranAnalysis, info)
  TWOdevice *pDevice;
  int iterationLimit;
  BOOLEAN tranAnalysis;
  TWOtranInfo *info;
{
  BOOLEAN newSolver = FALSE;
  int error;
  int index, eIndex;
  TWOelem *pElem;
  TWOnode *pNode;
  double refPsi;
  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);
  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("TWObiasSolve: Out of Memory\n");
      exit(-1);
    }
    newSolver = TRUE;
    if (NOT OneCarrier) {
      TWO_jacBuild(pDevice);
    } else if (OneCarrier IS N_TYPE) {
      TWONjacBuild(pDevice);
    } else if (OneCarrier IS P_TYPE) {
      TWOPjacBuild(pDevice);
    }
    pDevice->numOrigBias = spElementCount(pDevice->matrix);
    pDevice->numFillBias = 0;
    TWOstoreInitialGuess(pDevice);
  case SLV_SMSIG:
    spSetReal(pDevice->matrix);
  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 */
  TWOdcSolve(pDevice, iterationLimit, newSolver, tranAnalysis, info);

  /* MISCELLANEOUS */
  startTime = SPfrontEnd->IFseconds();
  if (newSolver) {
    pDevice->numFillBias = spFillinCount(pDevice->matrix);
  }
  if ((NOT pDevice->converged) AND iterationLimit > 1) {
    printf("TWObiasSolve: No Convergence\n");
  } else if (pDevice->converged) {
    /* update the nodal quantities */
    for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
      pElem = pDevice->elements[eIndex];
      refPsi = pElem->matlInfo->refPsi;
      for (index = 0; index <= 3; index++) {
	if (pElem->evalNodes[index]) {
	  pNode = pElem->pNodes[index];
	  if (pNode->nodeType ISNOT CONTACT) {
	    pNode->psi = pDevice->dcSolution[pNode->psiEqn];
	    if (pElem->elemType IS SEMICON) {
	      if (NOT OneCarrier) {
		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
	      } else if (OneCarrier IS N_TYPE) {
		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
		pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
	      } else if (OneCarrier IS P_TYPE) {
		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
		pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
	      }
	    }
	  }
	}
      }
    }

    /* update the current terms */
    if (NOT OneCarrier) {
      TWO_commonTerms(pDevice, FALSE, tranAnalysis, info);
    } else if (OneCarrier IS N_TYPE) {
      TWONcommonTerms(pDevice, FALSE, tranAnalysis, info);
    } else if (OneCarrier IS P_TYPE) {
      TWOPcommonTerms(pDevice, FALSE, tranAnalysis, info);
    }
  } else if (iterationLimit <= 1) {
    for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
      pElem = pDevice->elements[eIndex];
      refPsi = pElem->matlInfo->refPsi;
      for (index = 0; index <= 3; index++) {
	if (pElem->evalNodes[index]) {
	  pNode = pElem->pNodes[index];
	  if (pNode->nodeType ISNOT CONTACT) {
	    pNode->psi = pDevice->dcSolution[pNode->psiEqn];
	    *(pDevice->devState0 + pNode->nodePsi) = pNode->psi;
	    if (pElem->elemType IS SEMICON) {
	      if (NOT OneCarrier) {
		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
	      } else if (OneCarrier IS N_TYPE) {
		pNode->nConc = pDevice->dcSolution[pNode->nEqn];
		pNode->pConc = pNode->nie * exp(-pNode->psi + refPsi);
	      } else if (OneCarrier IS P_TYPE) {
		pNode->pConc = pDevice->dcSolution[pNode->pEqn];
		pNode->nConc = pNode->nie * exp(pNode->psi - refPsi);
	      }
	      *(pDevice->devState0 + pNode->nodeN) = pNode->nConc;
	      *(pDevice->devState0 + pNode->nodeP) = pNode->pConc;
	    }
	  }
	}
      }
    }
  }
  miscTime += SPfrontEnd->IFseconds() - startTime;
  if (tranAnalysis) {
    pDevice->pStats->setupTime[STAT_TRAN] += setupTime;
    pDevice->pStats->miscTime[STAT_TRAN] += miscTime;
  } else {
    pDevice->pStats->setupTime[STAT_DC] += setupTime;
    pDevice->pStats->miscTime[STAT_DC] += miscTime;
  }
}

/* Copy the device's equilibrium state into the solution vector. */
void
TWOstoreEquilibGuess(pDevice)
  TWOdevice *pDevice;
{
  int nIndex, eIndex;
  double *solution = pDevice->dcSolution;
  double refPsi;
  TWOelem *pElem;
  TWOnode *pNode;

  for (eIndex = 1; eIndex <= pDevice->numElems; eIndex++) {
    pElem = pDevice->elements[eIndex];
    refPsi = pElem->matlInfo->refPsi;
    for (nIndex = 0; nIndex <= 3; nIndex++) {
      if (pElem->evalNodes[nIndex]) {
	pNode = pElem->pNodes[nIndex];
	if (pNode->nodeType ISNOT CONTACT) {
	  pDevice->dcSolution[pNode->psiEqn] = pNode->psi0;
	  if (pElem->elemType IS SEMICON) {
	    if (NOT OneCarrier) {
	      solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
	      solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
	    } else if (OneCarrier IS N_TYPE) {
	      solution[pNode->nEqn] = pNode->nie * exp(pNode->psi0 - refPsi);
	    } else if (OneCarrier IS P_TYPE) {
	      solution[pNode->pEqn] = pNode->nie * exp(-pNode->psi0 + refPsi);
	    }
	  }
	}
      }
    }
  }
}

/* Copy the device's internal state into the solution vector. */
void
TWOstoreInitialGuess(pDevice)
  TWOdevice *pDevice;
{
  int nIndex, eIndex;
  double *solution = pDevice->dcSolution;
  TWOelem *pElem;
  TWOnode *pNode;

  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) {
	  solution[pNode->psiEqn] = pNode->psi;
	  if (pElem->elemType IS SEMICON) {
	    if (NOT OneCarrier) {
	      solution[pNode->nEqn] = pNode->nConc;
	      solution[pNode->pEqn] = pNode->pConc;
	    } else if (OneCarrier IS N_TYPE) {
	      solution[pNode->nEqn] = pNode->nConc;
	    } else if (OneCarrier IS P_TYPE) {
	      solution[pNode->pEqn] = pNode->pConc;
	    }
	  }
	}
      }
    }
  }
}


/*
 * function computeNewDelta computes an acceptable delta
 */

void
oldTWOnewDelta(pDevice, tranAnalysis, info)
  TWOdevice *pDevice;
  BOOLEAN tranAnalysis;
  TWOtranInfo *info;
{
  int index;
  double newNorm, fib, lambda, fibn, fibp;
  double TWOnuNorm(), l2Norm(), maxNorm();
  BOOLEAN acceptable = 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];
  }

  /* compute L2 norm of the deltaX vector */
  pDevice->rhsNorm = l2Norm(pDevice->dcDeltaSolution, pDevice->numEqns);

  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 = TWOnuNorm(pDevice);
  if (newNorm <= pDevice->rhsNorm) {
    acceptable = TRUE;
  } else {
    /* chop the step size so that deltax is acceptable */
    for (; NOT acceptable;) {
      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 = TWOnuNorm(pDevice);
      if (newNorm <= pDevice->rhsNorm) {
	acceptable = TRUE;
      }
    }
  }
  /* 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;
  }