twocont.c

来自「ngspice又一个电子CAD仿真软件代码.功能更全」· C语言 代码 · 共 1,053 行 · 第 1/3 页

/**********
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: twocont.c,v 1.3 2005/05/21 12:37:24 sjborley Exp $
**********/

#include "ngspice.h"
#include "numglobs.h"
#include "numenum.h"
#include "twomesh.h"
#include "twodev.h"
#include "bool.h"
#include "spmatrix.h"
#include "twoddefs.h"
#include "twodext.h"
#include "cidersupt.h"
#include "../../maths/misc/bernoull.h"

/*
 * Functions to setup and solve the continuity equations.
 * Both continuity equations are solved.
 * Separate functions are used for one continuity equation.
 */


/*
 * setup matrix pointers to Jacobian values and 
 * store direct pointers with the nodes
 */

void 
  TWO_jacBuild(TWOdevice *pDevice)
{
  char *matrix = pDevice->matrix;
  double *spGetElement();
  TWOelem *pElem;
  TWOnode *pNode;
  TWOchannel *pCh;
  int eIndex, nIndex;
  int nextIndex;			/* index of node to find next element */
  int psiEqn, nEqn, pEqn;		/* scratch for deref'd eqn numbers */
  int psiEqnTL = 0, nEqnTL = 0, pEqnTL = 0;
  int psiEqnTR = 0, nEqnTR = 0, pEqnTR = 0;
  int psiEqnBR = 0, nEqnBR = 0, pEqnBR = 0;
  int psiEqnBL = 0, nEqnBL = 0, pEqnBL = 0;
  int psiEqnInM = 0, psiEqnInP = 0;	/* scratch for deref'd surface eqns */
  int psiEqnOxM = 0, psiEqnOxP = 0;	/* M= more negative, P= more positive */
  
  for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
    pElem = pDevice->elements[ eIndex ];
    
    /* first the self terms */
    for ( nIndex = 0; nIndex <= 3; nIndex++ ) {
      pNode = pElem->pNodes[ nIndex ];
      /* get poisson-only pointer */
      psiEqn = pNode->psiEqn;
      pNode->fPsiPsi = spGetElement( matrix, psiEqn, psiEqn );
      
      if ( pElem->elemType == SEMICON ) {
	/* get continuity-coupling terms */
	nEqn = pNode->nEqn;
	pEqn = pNode->pEqn;
	/* pointers for additional terms */
	pNode->fPsiN = spGetElement( matrix, psiEqn, nEqn );
	pNode->fPsiP = spGetElement( matrix, psiEqn, pEqn );
	pNode->fNPsi = spGetElement( matrix, nEqn, psiEqn );
	pNode->fNN = spGetElement( matrix, nEqn, nEqn );
	pNode->fNP = spGetElement( matrix, nEqn, pEqn );
	pNode->fPPsi = spGetElement( matrix, pEqn, psiEqn );
	pNode->fPN = spGetElement( matrix, pEqn, nEqn );
	pNode->fPP = spGetElement( matrix, pEqn, pEqn );
      } else {
	nEqn = 0;
	pEqn = 0;
      }
      /* save equation indices */
      switch ( nIndex ) {
      case 0: /* TL Node */
	psiEqnTL = psiEqn;
	nEqnTL = nEqn;
	pEqnTL = pEqn;
	break;
      case 1: /* TR Node */
	psiEqnTR = psiEqn;
	nEqnTR = nEqn;
	pEqnTR = pEqn;
	break;
      case 2: /* BR Node */
	psiEqnBR = psiEqn;
	nEqnBR = nEqn;
	pEqnBR = pEqn;
	break;
      case 3: /* BL Node */
	psiEqnBL = psiEqn;
	nEqnBL = nEqn;
	pEqnBL = pEqn;
	break;
      default:
	break;
      }
    }
    
    /* now terms to couple to adjacent nodes */
    pNode = pElem->pTLNode;
    pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnTL, psiEqnTR );
    pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTL, psiEqnBL );
    if ( pElem->elemType == SEMICON ) {
      /* continuity equation pointers */
      pNode->fNPsiiP1    = spGetElement( matrix, nEqnTL, psiEqnTR );
      pNode->fNNiP1      = spGetElement( matrix, nEqnTL, nEqnTR );
      pNode->fNPsijP1    = spGetElement( matrix, nEqnTL, psiEqnBL );
      pNode->fNNjP1      = spGetElement( matrix, nEqnTL, nEqnBL );
      pNode->fPPsiiP1    = spGetElement( matrix, pEqnTL, psiEqnTR );
      pNode->fPPiP1      = spGetElement( matrix, pEqnTL, pEqnTR );
      pNode->fPPsijP1    = spGetElement( matrix, pEqnTL, psiEqnBL );
      pNode->fPPjP1      = spGetElement( matrix, pEqnTL, pEqnBL );
      /* Surface Mobility Model depends on diagonal node values */
      if ( MobDeriv && SurfaceMobility && pElem->channel ) {
	pNode->fNPsiiP1jP1 = spGetElement( matrix, nEqnTL, psiEqnBR );
	pNode->fNNiP1jP1   = spGetElement( matrix, nEqnTL, nEqnBR );
	pNode->fPPsiiP1jP1 = spGetElement( matrix, pEqnTL, psiEqnBR );
	pNode->fPPiP1jP1   = spGetElement( matrix, pEqnTL, pEqnBR );
      }
    }
    
    pNode = pElem->pTRNode;
    pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnTR, psiEqnTL );
    pNode->fPsiPsijP1 = spGetElement(matrix, psiEqnTR, psiEqnBR );
    if ( pElem->elemType == SEMICON ) {
      /* continuity equation pointers */
      pNode->fNPsiiM1    = spGetElement( matrix, nEqnTR, psiEqnTL );
      pNode->fNNiM1      = spGetElement( matrix, nEqnTR, nEqnTL );
      pNode->fNPsijP1    = spGetElement( matrix, nEqnTR, psiEqnBR );
      pNode->fNNjP1      = spGetElement( matrix, nEqnTR, nEqnBR );
      pNode->fPPsiiM1    = spGetElement( matrix, pEqnTR, psiEqnTL );
      pNode->fPPiM1      = spGetElement( matrix, pEqnTR, pEqnTL );
      pNode->fPPsijP1    = spGetElement( matrix, pEqnTR, psiEqnBR );
      pNode->fPPjP1      = spGetElement( matrix, pEqnTR, pEqnBR );
      /* Surface Mobility Model depends on diagonal node values */
      if ( MobDeriv && SurfaceMobility && pElem->channel ) {
	pNode->fNPsiiM1jP1 = spGetElement( matrix, nEqnTR, psiEqnBL );
	pNode->fNNiM1jP1   = spGetElement( matrix, nEqnTR, nEqnBL );
	pNode->fPPsiiM1jP1 = spGetElement( matrix, pEqnTR, psiEqnBL );
	pNode->fPPiM1jP1   = spGetElement( matrix, pEqnTR, pEqnBL );
      }
    }
    
    pNode = pElem->pBRNode;
    pNode->fPsiPsiiM1 = spGetElement(matrix, psiEqnBR, psiEqnBL );
    pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBR, psiEqnTR );
    if ( pElem->elemType == SEMICON ) {
      /* continuity equation pointers */
      pNode->fNPsiiM1    = spGetElement( matrix, nEqnBR, psiEqnBL );
      pNode->fNNiM1      = spGetElement( matrix, nEqnBR, nEqnBL );
      pNode->fNPsijM1    = spGetElement( matrix, nEqnBR, psiEqnTR );
      pNode->fNNjM1      = spGetElement( matrix, nEqnBR, nEqnTR );
      pNode->fPPsiiM1    = spGetElement( matrix, pEqnBR, psiEqnBL );
      pNode->fPPiM1      = spGetElement( matrix, pEqnBR, pEqnBL );
      pNode->fPPsijM1    = spGetElement( matrix, pEqnBR, psiEqnTR );
      pNode->fPPjM1      = spGetElement( matrix, pEqnBR, pEqnTR );
      /* Surface Mobility Model depends on diagonal node values */
      if ( MobDeriv && SurfaceMobility && pElem->channel ) {
	pNode->fNPsiiM1jM1 = spGetElement( matrix, nEqnBR, psiEqnTL );
	pNode->fNNiM1jM1   = spGetElement( matrix, nEqnBR, nEqnTL );
	pNode->fPPsiiM1jM1 = spGetElement( matrix, pEqnBR, psiEqnTL );
	pNode->fPPiM1jM1   = spGetElement( matrix, pEqnBR, pEqnTL );
      }
    }
    
    pNode = pElem->pBLNode;
    pNode->fPsiPsiiP1 = spGetElement(matrix, psiEqnBL, psiEqnBR );
    pNode->fPsiPsijM1 = spGetElement(matrix, psiEqnBL, psiEqnTL );
    if ( pElem->elemType == SEMICON ) {
      /* continuity equation pointers */
      pNode->fNPsiiP1    = spGetElement( matrix, nEqnBL, psiEqnBR );
      pNode->fNNiP1      = spGetElement( matrix, nEqnBL, nEqnBR );
      pNode->fNPsijM1    = spGetElement( matrix, nEqnBL, psiEqnTL );
      pNode->fNNjM1      = spGetElement( matrix, nEqnBL, nEqnTL );
      pNode->fPPsiiP1    = spGetElement( matrix, pEqnBL, psiEqnBR );
      pNode->fPPiP1      = spGetElement( matrix, pEqnBL, pEqnBR );
      pNode->fPPsijM1    = spGetElement( matrix, pEqnBL, psiEqnTL );
      pNode->fPPjM1      = spGetElement( matrix, pEqnBL, pEqnTL );
      /* Surface Mobility Model depends on diagonal node values */
      if ( MobDeriv && SurfaceMobility && pElem->channel ) {
	pNode->fNPsiiP1jM1 = spGetElement( matrix, nEqnBL, psiEqnTR );
	pNode->fNNiP1jM1   = spGetElement( matrix, nEqnBL, nEqnTR );
	pNode->fPPsiiP1jM1 = spGetElement( matrix, pEqnBL, psiEqnTR );
	pNode->fPPiP1jM1   = spGetElement( matrix, pEqnBL, pEqnTR );
      }
    }
  }
  /* 
   * Add terms for surface-field of inversion-layer mobility model.
   * Elements MUST be made from silicon for this to work.
   * No empty elements are allowed.
   * Don't need these pointers if SurfaceMobility isn't set.
   */
  if ( MobDeriv && SurfaceMobility ) {
    for ( pCh = pDevice->pChannel; pCh != NIL(TWOchannel);
	 pCh = pCh->next ) {
      pElem = pCh->pNElem;
      switch (pCh->type) {
      case 0:
	psiEqnInM = pElem->pBLNode->psiEqn;
	psiEqnInP = pElem->pBRNode->psiEqn;
	psiEqnOxM = pElem->pTLNode->psiEqn;
	psiEqnOxP = pElem->pTRNode->psiEqn;
	break;
      case 1:
	psiEqnInM = pElem->pTLNode->psiEqn;
	psiEqnInP = pElem->pBLNode->psiEqn;
	psiEqnOxM = pElem->pTRNode->psiEqn;
	psiEqnOxP = pElem->pBRNode->psiEqn;
	break;
      case 2:
	psiEqnInM = pElem->pTLNode->psiEqn;
	psiEqnInP = pElem->pTRNode->psiEqn;
	psiEqnOxM = pElem->pBLNode->psiEqn;
	psiEqnOxP = pElem->pBRNode->psiEqn;
	break;
      case 3:
	psiEqnInM = pElem->pTRNode->psiEqn;
	psiEqnInP = pElem->pBRNode->psiEqn;
	psiEqnOxM = pElem->pTLNode->psiEqn;
	psiEqnOxP = pElem->pBLNode->psiEqn;
	break;
      }
      pElem = pCh->pSeed;
      nextIndex = (pCh->type + 2)%4;
      while (pElem && pElem->channel == pCh->id) {
	for ( nIndex = 0; nIndex <= 3; nIndex++ ) {
	  pNode = pElem->pNodes[ nIndex ];
	  psiEqn = pNode->psiEqn;
	  nEqn = pNode->nEqn;
	  pEqn = pNode->pEqn;
	  if ( pCh->type % 2 == 0 ) { /* Vertical Slice */
	    if ( nIndex == 0 || nIndex == 3 ) { /* Left Side */
	      pNode->fNPsiIn   = spGetElement( matrix, nEqn, psiEqnInM );
	      pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP );
	      pNode->fNPsiOx   = spGetElement( matrix, nEqn, psiEqnOxM );
	      pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP );
	      pNode->fPPsiIn   = spGetElement( matrix, pEqn, psiEqnInM );
	      pNode->fPPsiInP1 = spGetElement( matrix, pEqn, psiEqnInP );
	      pNode->fPPsiOx   = spGetElement( matrix, pEqn, psiEqnOxM );
	      pNode->fPPsiOxP1 = spGetElement( matrix, pEqn, psiEqnOxP );
	    } else { /* Right Side */
	      pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM );
	      pNode->fNPsiIn   = spGetElement( matrix, nEqn, psiEqnInP );
	      pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM );
	      pNode->fNPsiOx   = spGetElement( matrix, nEqn, psiEqnOxP );
	      pNode->fPPsiInM1 = spGetElement( matrix, pEqn, psiEqnInM );
	      pNode->fPPsiIn   = spGetElement( matrix, pEqn, psiEqnInP );
	      pNode->fPPsiOxM1 = spGetElement( matrix, pEqn, psiEqnOxM );
	      pNode->fPPsiOx   = spGetElement( matrix, pEqn, psiEqnOxP );
	    }
	  } else { /* Horizontal Slice */
	    if ( nIndex <= 1 ) { /* Top Side */
	      pNode->fNPsiIn   = spGetElement( matrix, nEqn, psiEqnInM );
	      pNode->fNPsiInP1 = spGetElement( matrix, nEqn, psiEqnInP );
	      pNode->fNPsiOx   = spGetElement( matrix, nEqn, psiEqnOxM );
	      pNode->fNPsiOxP1 = spGetElement( matrix, nEqn, psiEqnOxP );
	      pNode->fPPsiIn   = spGetElement( matrix, pEqn, psiEqnInM );
	      pNode->fPPsiInP1 = spGetElement( matrix, pEqn, psiEqnInP );
	      pNode->fPPsiOx   = spGetElement( matrix, pEqn, psiEqnOxM );
	      pNode->fPPsiOxP1 = spGetElement( matrix, pEqn, psiEqnOxP );
	    } else { /* Bottom Side */
	      pNode->fNPsiInM1 = spGetElement( matrix, nEqn, psiEqnInM );
	      pNode->fNPsiIn   = spGetElement( matrix, nEqn, psiEqnInP );
	      pNode->fNPsiOxM1 = spGetElement( matrix, nEqn, psiEqnOxM );
	      pNode->fNPsiOx   = spGetElement( matrix, nEqn, psiEqnOxP );
	      pNode->fPPsiInM1 = spGetElement( matrix, pEqn, psiEqnInM );
	      pNode->fPPsiIn   = spGetElement( matrix, pEqn, psiEqnInP );
	      pNode->fPPsiOxM1 = spGetElement( matrix, pEqn, psiEqnOxM );
	      pNode->fPPsiOx   = spGetElement( matrix, pEqn, psiEqnOxP );
	    }
	  }
	} /* endfor nIndex */
	pElem = pElem->pElems[ nextIndex ];
      } /* endwhile pElem */
    } /* endfor pCh */
  } /* endif SurfaceMobility */
}


/*
 *  The Jacobian and Rhs are loaded by the following function.
 *  Inputs are the transient analysis flag and the transient
 *  information structure
 */

void 
  TWO_sysLoad(TWOdevice *pDevice, BOOLEAN tranAnalysis, TWOtranInfo *info)
{
  TWOelem *pElem;
  TWOnode *pNode;
  TWOedge *pHEdge, *pVEdge;
  TWOedge *pTEdge, *pBEdge, *pLEdge, *pREdge;
  TWOchannel *pCh;
  int index, eIndex;
  int nextIndex;			/* index of node to find next element */
  double *pRhs = pDevice->rhs;
  double dx, dy, dxdy, dyOverDx, dxOverDy;
  double ds;
  double dPsiT, dPsiB, dPsiL, dPsiR;
  double rhsN, rhsP;
  double generation, TWOavalanche();
  double nConc, pConc;
  double perTime = 0.0;
  
  /* first compute the currents and derivatives */
  TWO_commonTerms( pDevice, FALSE, tranAnalysis, info );
  
  /* find reciprocal timestep */
  if ( tranAnalysis ) {
    perTime = info->intCoeff[0];
  }
  
  /* zero the rhs vector */
  for ( index = 1 ; index <= pDevice->numEqns ; index++ ) {
    pRhs[ index ] = 0.0;
  }
  
  /* zero the matrix */
  spClear( pDevice->matrix );
  
  for ( eIndex = 1; eIndex <= pDevice->numElems; eIndex++ ) {
    pElem = pDevice->elements[ eIndex ];
    
    dx = 0.5 * pElem->dx;
    dy = 0.5 * pElem->dy;
    dxdy = dx * dy;
    dxOverDy = 0.5 * pElem->epsRel * pElem->dxOverDy;
    dyOverDx = 0.5 * pElem->epsRel * pElem->dyOverDx;
    
    pTEdge = pElem->pTopEdge;
    pBEdge = pElem->pBotEdge;
    pLEdge = pElem->pLeftEdge;
    pREdge = pElem->pRightEdge;
    dPsiT = pTEdge->dPsi;
    dPsiB = pBEdge->dPsi;
    dPsiL = pLEdge->dPsi;
    dPsiR = pREdge->dPsi;
    
    /* load for all i,j */
    for ( index = 0; index <= 3; index++ ) {
      pNode = pElem->pNodes[ index ];
      if ( pNode->nodeType != CONTACT ) {
	*(pNode->fPsiPsi) += dyOverDx + dxOverDy;
	if ( index <= 1 ) {
	  pHEdge = pTEdge;