twoproj.c

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

		solution[ pNode->pEqn ] = guessNewConc( pConc, delP );
	      }
	      else {
		solution[ pNode->pEqn ] = newP;
	      }
	    }
	  }
	}
      }
    }
  }
  
  if ( ABS( delVgb ) > MIN_DELV ) {
    
    incVgb = pDevice->dcDeltaSolution;
    storeNewRhs( pDevice, pGContact );
    spSolve( pDevice->matrix, pDevice->rhs, incVgb, NIL(spREAL), NIL(spREAL));
    
    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->nodeType != CONTACT ) {
	    delPsi = incVgb[ pNode->psiEqn ] * delVgb;
	    solution[ pNode->psiEqn ] += delPsi;
	    if ( pElem->elemType == SEMICON
		&& (!OneCarrier || (OneCarrier == N_TYPE)) ) {
	      delN = incVgb[ pNode->nEqn ] * delVgb;
	      nConc = solution[ pNode->nEqn ];
	      newN = nConc + delN;
	      if ( newN <= 0.0 ) {
		solution[ pNode->nEqn ] = guessNewConc( nConc, delN );
	      }
	      else {
		solution[ pNode->nEqn ] = newN;
	      }
	    }
	    if ( pElem->elemType == SEMICON
		&& (!OneCarrier || (OneCarrier == P_TYPE)) ) {
	      delP = incVgb[ pNode->pEqn ] * delVgb;
	      pConc = solution[ pNode->pEqn ];
	      newP = pConc + delP;
	      if ( newP <= 0.0 ) {
		solution[ pNode->pEqn ] = guessNewConc( pConc, delP );
	      }
	      else {
		solution[ pNode->pEqn ] = newP;
	      }
	    }
	  }
	}
      }
    }
  }
}

/* functions to update the solution for the full-LU and
   modified two-level Newton methods
   */

void NUMD2update(TWOdevice *pDevice, double delV, BOOLEAN updateBoundary)
{
  TWOnode *pNode;
  TWOelem *pElem;
  int index, eIndex, numContactNodes;
  TWOcontact *pContact = pDevice->pLastContact;
  double delPsi, delN, delP, *incVpn, *solution = pDevice->dcSolution;
  
  delV = -delV / VNorm;
  if ( updateBoundary ) {
    /* update the boundary condition for the last contact */
    numContactNodes = pContact->numNodes;
    for ( index = 0; index < numContactNodes; index++ ) {
      pNode = pContact->pNodes[ index ];
      pNode->psi += delV;
    }
  }
  
  /* the equations have been solved while computing the conductances */
  /* solution is in dcDeltaSolution, so use it */

  incVpn = pDevice->dcDeltaSolution;
  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->nodeType != CONTACT ) {
	  delPsi = incVpn[ pNode->psiEqn ] * delV;
	  solution[ pNode->psiEqn ] = pNode->psi + delPsi;
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == N_TYPE)) ) {
	    delN = incVpn[ pNode->nEqn ] * delV;
	    solution[ pNode->nEqn ] = pNode->nConc + delN;
	  }
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == P_TYPE)) ) {
	    delP = incVpn[ pNode->pEqn ] * delV;
	    solution[ pNode->pEqn ] = pNode->pConc + delP;
	  }
	}
      }
    }
  }
}


void NBJT2update(TWOdevice *pDevice, double delVce, double delVbe, 
                 BOOLEAN updateBoundary)
{
  TWOnode *pNode;
  TWOelem *pElem;
  int index, eIndex, numContactNodes;
  TWOcontact *pColContact = pDevice->pFirstContact;
  TWOcontact *pBaseContact = pDevice->pFirstContact->next;
  double delPsi, delN, delP, *incVce, *incVbe, *solution = pDevice->dcSolution;
  
  /* normalize the voltages for calculations */
  
  if ( delVce != 0.0 ) {
    delVce = delVce / VNorm;
    if ( updateBoundary ) {
      numContactNodes = pColContact->numNodes;
      for ( index = 0; index < numContactNodes; index++ ) {
	pNode = pColContact->pNodes[ index ];
	pNode->psi += delVce;
      }
    }
  }
  if ( delVbe != 0.0 ) {
    delVbe = delVbe / VNorm;
    if ( updateBoundary ) {
      numContactNodes = pBaseContact->numNodes;
      for ( index = 0; index < numContactNodes; index++ ) {
	pNode = pBaseContact->pNodes[ index ];
	pNode->psi += delVbe;
      }
    }
  }
  
  /* use solution from computeConductance to do update */
  /* set incVce = dcDeltaSolution; incVbe = copiedSolution */
  
  incVce = pDevice->dcDeltaSolution;
  incVbe = pDevice->copiedSolution;
  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->nodeType != CONTACT ) {
	  delPsi = (incVce[ pNode->psiEqn ] * delVce
		    + incVbe[ pNode->psiEqn ] * delVbe);
	  solution[ pNode->psiEqn ] = pNode->psi + delPsi;
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == N_TYPE)) ) {
	    delN = (incVce[ pNode->nEqn ] * delVce
		    + incVbe[ pNode->nEqn ] * delVbe);
	    solution[ pNode->nEqn ] = pNode->nConc + delN;
	  }
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == P_TYPE)) ) {
	    delP = (incVce[ pNode->pEqn ] * delVce
		    + incVbe[ pNode->pEqn ] * delVbe);
	    solution[ pNode->pEqn ] = pNode->pConc + delP;
	  }
	}
      }
    }
  }
}


void NUMOSupdate(TWOdevice *pDevice, double delVdb, double delVsb, 
                 double delVgb, BOOLEAN updateBoundary)
{
  TWOnode *pNode;
  TWOelem *pElem;
  int index, eIndex, numContactNodes;
  TWOcontact *pDContact = pDevice->pFirstContact;
  TWOcontact *pGContact = pDevice->pFirstContact->next;
  TWOcontact *pSContact = pDevice->pFirstContact->next->next;
  double delPsi, delN, delP;
  double *incVdb, *incVsb, *incVgb, *solution = pDevice->dcSolution;
  
  /* normalize the voltages for calculations */
  if ( delVdb != 0.0 ) {
    delVdb = delVdb / VNorm;
    if ( updateBoundary ) {
      numContactNodes = pDContact->numNodes;
      for ( index = 0; index < numContactNodes; index++ ) {
	pNode = pDContact->pNodes[ index ];
	pNode->psi += delVdb;
      }
    }
  }
  if ( delVsb != 0.0 ) {
    delVsb = delVsb / VNorm;
    if ( updateBoundary ) {
      numContactNodes = pSContact->numNodes;
      for ( index = 0; index < numContactNodes; index++ ) {
	pNode = pSContact->pNodes[ index ];
	pNode->psi += delVsb;
      }
    }
  }
  if ( delVgb != 0.0 ) {
    delVgb = delVgb / VNorm;
    if ( updateBoundary ) {
      numContactNodes = pGContact->numNodes;
      for ( index = 0; index < numContactNodes; index++ ) {
	pNode = pGContact->pNodes[ index ];
	pNode->psi += delVgb;
      }
    }
  }
  
  /* use solution from computeConductance to do update */
  
  incVdb = pDevice->dcDeltaSolution;
  incVsb = pDevice->copiedSolution;
  incVgb = pDevice->rhsImag;
  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->nodeType != CONTACT ) {
	  delPsi = (incVdb[ pNode->psiEqn ] * delVdb
		    + incVsb[ pNode->psiEqn ] * delVsb
		    + incVgb[ pNode->psiEqn ] * delVgb);
	  solution[ pNode->psiEqn ] = pNode->psi + delPsi;
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == N_TYPE)) ) {
	    delN = (incVdb[ pNode->nEqn ] * delVdb
		    + incVsb[ pNode->nEqn ] * delVsb
		    + incVgb[ pNode->nEqn ] * delVgb);
	    solution[ pNode->nEqn ] = pNode->nConc + delN;
	  }
	  if ( pElem->elemType == SEMICON
	      && (!OneCarrier || (OneCarrier == P_TYPE)) ) {
	    delP = (incVdb[ pNode->pEqn ] * delVdb
		    + incVsb[ pNode->pEqn ] * delVsb
		    + incVgb[ pNode->pEqn ] * delVgb);
	    solution[ pNode->pEqn ] = pNode->pConc + delP;
	  }
	}
      }
    }
  }
}


void storeNewRhs(TWOdevice *pDevice, TWOcontact *pContact)
{
  int index, i, numContactNodes;
  TWOelem *pElem;
  TWOnode *pNode, *pHNode = NULL, *pVNode = NULL;
  TWOedge *pHEdge = NULL, *pVEdge = NULL;
  double *rhs = pDevice->rhs;
  
  /* zero the rhs before loading in the new rhs */
  for ( index = 1; index <= pDevice->numEqns ; index++ ) {
    rhs[ index ] = 0.0;
  }
  
  numContactNodes = pContact->numNodes;
  for ( index = 0; index < numContactNodes; index++ ) {
    pNode = pContact->pNodes[ index ];
    for ( i = 0; i <= 3; i++ ) {
      pElem = pNode->pElems[ i ];
      if ( pElem != NIL(TWOelem)) {
	/* found an element to which this node belongs */
	switch ( i ) {
	case 0:
	  /* the TL element of this node */
	  pHNode = pElem->pBLNode;
	  pVNode = pElem->pTRNode;
	  pHEdge = pElem->pBotEdge;
	  pVEdge = pElem->pRightEdge;
	  break;
	case 1:
	  /* the TR element */
	  pHNode = pElem->pBRNode;
	  pVNode = pElem->pTLNode;
	  pHEdge = pElem->pBotEdge;
	  pVEdge = pElem->pLeftEdge;
	  break;
	case 2:
	  /* the BR element */
	  pHNode = pElem->pTRNode;
	  pVNode = pElem->pBLNode;
	  pHEdge = pElem->pTopEdge;
	  pVEdge = pElem->pLeftEdge;
	  break;
	case 3:
	  /* the BL element */
	  pHNode = pElem->pTLNode;
	  pVNode = pElem->pBRNode;
	  pHEdge = pElem->pTopEdge;
	  pVEdge = pElem->pRightEdge;
	  break;
	default:
	  printf( "storeNewRhs: shouldn't be here\n");
	  break;
	}
	if ( pHNode->nodeType != CONTACT ) {
	  /* contribution to the x nodes */
	  rhs[ pHNode->psiEqn ] += pElem->epsRel * 0.5 * pElem->dyOverDx;
	  if ( pElem->elemType == SEMICON ) {
	    if ( !OneCarrier ) {
	      rhs[ pHNode->nEqn ] -= 0.5 * pElem->dy * pHEdge->dJnDpsiP1;
	      rhs[ pHNode->pEqn ] -= 0.5 * pElem->dy * pHEdge->dJpDpsiP1;
	    } else if ( OneCarrier == N_TYPE ) {
	      rhs[ pHNode->nEqn ] -= 0.5 * pElem->dy * pHEdge->dJnDpsiP1;
	    } else if ( OneCarrier == P_TYPE ) {
	      rhs[ pHNode->pEqn ] -= 0.5 * pElem->dy * pHEdge->dJpDpsiP1;
	    }
	  }
	}
	if ( pVNode->nodeType != CONTACT ) {
	  /* contribution to the y nodes */
	  rhs[ pVNode->psiEqn ] += pElem->epsRel * 0.5 * pElem->dxOverDy;
	  if ( pElem->elemType == SEMICON ) {
	    if ( !OneCarrier ) {
	      rhs[ pVNode->nEqn ] -= 0.5 * pElem->dx * pVEdge->dJnDpsiP1;
	      rhs[ pVNode->pEqn ] -= 0.5 * pElem->dx * pVEdge->dJpDpsiP1;
	    } else if ( OneCarrier == N_TYPE ) {
	      rhs[ pVNode->nEqn ] -= 0.5 * pElem->dx * pVEdge->dJnDpsiP1;
	    } else if ( OneCarrier == P_TYPE ) {
	      rhs[ pVNode->pEqn ] -= 0.5 * pElem->dx * pVEdge->dJpDpsiP1;
	    }
	  }
	}
      }
    }
  }
}