ngspice/src/spicelib/analysis/cktop.c

394 lines
15 KiB
C

/**********
Copyright 1990 Regents of the University of California. All rights reserved.
Author: 1985 Thomas L. Quarles
Modified: 2000 AlansFixes
**********/
#include "ngspice.h"
#include "cktdefs.h"
#include "devdefs.h"
#include "sperror.h"
int
CKTop(CKTcircuit *ckt, long int firstmode, long int continuemode, int iterlim)
{
int converged;
int i;
CKTnode *n;
double raise, ConvFact, NumNodes;
double *OldRhsOld, *OldCKTstate0;
int iters;
ckt->CKTmode = firstmode;
if(!ckt->CKTnoOpIter) {
converged = NIiter(ckt,iterlim);
} else {
converged = 1; /* the 'go directly to gmin stepping' option */
}
if(converged != 0) {
/* no convergence on the first try, so we do something else */
/* first, check if we should try gmin stepping */
/* note that no path out of this code allows ckt->CKTdiagGmin to be
* anything but 0.000000000
*/
if(ckt->CKTnumGminSteps ==1) {
double OldGmin, gtarget, factor;
int success, failed;
ckt->CKTmode = firstmode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"trying dynamic Gmin stepping",(IFuid *)NULL);
NumNodes=0;
for (n = ckt->CKTnodes; n; n = n->next) {
NumNodes++;
};
OldRhsOld=(double *)MALLOC((NumNodes+1)*sizeof(double));
OldCKTstate0=(double *)
MALLOC((ckt->CKTnumStates+1)*sizeof(double));
for (n = ckt->CKTnodes; n; n = n->next) {
*(ckt->CKTrhsOld+n->number)=0;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(ckt->CKTstate0+i) = 0;
};
factor = ckt->CKTgminFactor;
OldGmin = 1e-2;
ckt->CKTdiagGmin = OldGmin / factor;
gtarget = MAX(ckt->CKTgmin,ckt->CKTgshunt);
success = failed = 0;
while ( (!success) && (!failed) ) {
fprintf(stderr, "\rTrying gmin = %12.4E ", ckt->CKTdiagGmin);
ckt->CKTnoncon =1;
iters = ckt->CKTstat->STATnumIter;
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
iters = (ckt->CKTstat->STATnumIter)-iters;
if(converged == 0) {
ckt->CKTmode=continuemode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful Gmin step",(IFuid *)NULL);
if (ckt->CKTdiagGmin <= gtarget) {
success = 1;
} else {
i=0;
for (n = ckt->CKTnodes; n; n = n->next) {
OldRhsOld[i]=*(ckt->CKTrhsOld+n->number);
i++;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(OldCKTstate0+i) = *(ckt->CKTstate0+i);
};
if (iters <= (ckt->CKTdcTrcvMaxIter/4)) {
factor *= sqrt(factor);
if (factor > ckt->CKTgminFactor)
factor = ckt->CKTgminFactor;
};
if (iters > (3*ckt->CKTdcTrcvMaxIter/4)) {
factor = sqrt(factor);
};
OldGmin = ckt->CKTdiagGmin;
if ((ckt->CKTdiagGmin) < (factor*gtarget)) {
factor = ckt->CKTdiagGmin / gtarget;
ckt->CKTdiagGmin = gtarget;
} else {
ckt->CKTdiagGmin /= factor;
};
};
} else {
if (factor < 1.00005) {
failed = 1;
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"last gmin step failed",(IFuid *)NULL);
} else {
factor=sqrt(sqrt(factor));
ckt->CKTdiagGmin = OldGmin / factor;
i=0;
for (n = ckt->CKTnodes; n; n = n->next) {
*(ckt->CKTrhsOld+n->number)=OldRhsOld[i];
i++;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(ckt->CKTstate0+i) = *(OldCKTstate0+i);
};
};
}
}
ckt->CKTdiagGmin=ckt->CKTgshunt;
FREE(OldRhsOld);
FREE(OldCKTstate0);
converged = NIiter(ckt,iterlim);
if (converged!=0) {
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"dynamic gmin stepping failed",(IFuid *)NULL);
} else {
(*(SPfrontEnd->IFerror))(ERR_INFO,
"Dynamic gmin stepping completed",(IFuid *)NULL);
return(0);
};
} else if(ckt->CKTnumGminSteps >1) {
ckt->CKTmode = firstmode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"starting Gmin stepping",(IFuid *)NULL);
if (ckt->CKTgshunt==0) {
ckt->CKTdiagGmin = ckt->CKTgmin;
} else {
ckt->CKTdiagGmin = ckt->CKTgshunt;
};
for(i=0;i<ckt->CKTnumGminSteps;i++) {
ckt->CKTdiagGmin *= ckt->CKTgminFactor;
}
for(i=0;i<=ckt->CKTnumGminSteps;i++) {
fprintf(stderr, "Trying gmin = %12.4E ", ckt->CKTdiagGmin);
ckt->CKTnoncon =1;
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
if(converged != 0) {
ckt->CKTdiagGmin = ckt->CKTgshunt;
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"Gmin step failed",(IFuid *)NULL);
break;
}
ckt->CKTdiagGmin /= ckt->CKTgminFactor;
ckt->CKTmode=continuemode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful Gmin step",(IFuid *)NULL);
}
ckt->CKTdiagGmin = ckt->CKTgshunt;
converged = NIiter(ckt,iterlim);
if(converged == 0) {
(*(SPfrontEnd->IFerror))(ERR_INFO,
"Gmin stepping completed",(IFuid *)NULL);
return(0);
}
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"Gmin stepping failed",(IFuid *)NULL);
}
/* now, we'll try source stepping - we scale the sources
* to 0, converge, then start stepping them up until they
* are at their normal values
*
* note that no path out of this code allows ckt->CKTsrcFact to be
* anything but 1.000000000
*/
if(ckt->CKTnumSrcSteps >=1) {
ckt->CKTmode = firstmode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"starting source stepping",(IFuid *)NULL);
if(ckt->CKTnumSrcSteps==1) {
ckt->CKTsrcFact=0; raise=0.001; ConvFact=0;
NumNodes=0;
for (n = ckt->CKTnodes; n; n = n->next) {
NumNodes++;
};
OldRhsOld=(double *)MALLOC((NumNodes+1)*sizeof(double));
OldCKTstate0=(double *)
MALLOC((ckt->CKTnumStates+1)*sizeof(double));
for (n = ckt->CKTnodes; n; n = n->next) {
*(ckt->CKTrhsOld+n->number)=0;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(ckt->CKTstate0+i) = 0;
};
/* First, try a straight solution with all sources at zero */
fprintf(stderr, "\rSupplies reduced to %8.4f%% ",
ckt->CKTsrcFact*100);
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
/* If this doesn't work, try gmin stepping as well for the first solution */
if(converged != 0) {
fprintf(stderr, "\n");
if (ckt->CKTgshunt<=0) {
ckt->CKTdiagGmin = ckt->CKTgmin;
} else {
ckt->CKTdiagGmin = ckt->CKTgshunt;
};
for(i=0;i<10;i++) {
ckt->CKTdiagGmin *= 10;
}
for(i=0;i<=10;i++) {
fprintf(stderr, "Trying gmin = %12.4E ",
ckt->CKTdiagGmin);
ckt->CKTnoncon =1;
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
if(converged != 0) {
ckt->CKTdiagGmin = ckt->CKTgshunt;
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"Gmin step failed",(IFuid *)NULL);
break;
}
ckt->CKTdiagGmin /= 10;
ckt->CKTmode=continuemode;
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful Gmin step",(IFuid *)NULL);
}
ckt->CKTdiagGmin = ckt->CKTgshunt;
};
/* If we've got convergence, then try stepping up the sources */
if(converged == 0) {
i=0;
for (n = ckt->CKTnodes; n; n = n->next) {
OldRhsOld[i]=*(ckt->CKTrhsOld+n->number);
i++;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(OldCKTstate0+i) = *(ckt->CKTstate0+i);
};
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful source step",(IFuid *)NULL);
ckt->CKTsrcFact=ConvFact+raise;
};
if(converged == 0) do {
fprintf(stderr, "\rSupplies reduced to %8.4f%% ",
ckt->CKTsrcFact*100);
iters = ckt->CKTstat->STATnumIter;
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
iters = (ckt->CKTstat->STATnumIter)-iters;
ckt->CKTmode = continuemode;
if (converged == 0) {
ConvFact=ckt->CKTsrcFact;
i=0;
for (n = ckt->CKTnodes; n; n = n->next) {
OldRhsOld[i]=*(ckt->CKTrhsOld+n->number);
i++;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(OldCKTstate0+i) = *(ckt->CKTstate0+i);
};
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful source step",(IFuid *)NULL);
ckt->CKTsrcFact=ConvFact+raise;
if (iters <= (ckt->CKTdcTrcvMaxIter/4)) {
raise=raise*1.5;
};
if (iters > (3*ckt->CKTdcTrcvMaxIter/4)) {
raise=raise*0.5;
};
/* if (raise>0.01) raise=0.01; */
} else {
if ((ckt->CKTsrcFact-ConvFact)<1e-8) break;
raise=raise/10;
if (raise>0.01) raise=0.01;
ckt->CKTsrcFact=ConvFact;
i=0;
for (n = ckt->CKTnodes; n; n = n->next) {
*(ckt->CKTrhsOld+n->number)=OldRhsOld[i];
i++;
};
for(i=0;i<ckt->CKTnumStates;i++) {
*(ckt->CKTstate0+i) = *(OldCKTstate0+i);
};
};
if ((ckt->CKTsrcFact)>1) ckt->CKTsrcFact=1;
} while ((raise>=1e-7) && (ConvFact<1));
FREE(OldRhsOld);
FREE(OldCKTstate0);
ckt->CKTsrcFact = 1;
if (ConvFact!=1) {
ckt->CKTsrcFact = 1;
ckt->CKTcurrentAnalysis = DOING_TRAN;
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"source stepping failed",(IFuid *)NULL);
return(E_ITERLIM);
} else {
(*(SPfrontEnd->IFerror))(ERR_INFO,
"Source stepping completed",(IFuid *)NULL);
return(0);
};
} else {
for(i=0;i<=ckt->CKTnumSrcSteps;i++) {
ckt->CKTsrcFact = ((double)i)/((double)ckt->CKTnumSrcSteps);
converged = NIiter(ckt,ckt->CKTdcTrcvMaxIter);
ckt->CKTmode = continuemode;
if(converged != 0) {
ckt->CKTsrcFact = 1;
ckt->CKTcurrentAnalysis = DOING_TRAN;
(*(SPfrontEnd->IFerror))(ERR_WARNING,
"source stepping failed",(IFuid *)NULL);
return(converged);
}
(*(SPfrontEnd->IFerror))(ERR_INFO,
"One successful source step",(IFuid *)NULL);
}
(*(SPfrontEnd->IFerror))(ERR_INFO,
"Source stepping completed",(IFuid *)NULL);
ckt->CKTsrcFact = 1;
return(0);
};
} else {
return(converged);
}
}
return(0);
}
/* CKTconvTest(ckt)
* this is a driver program to iterate through all the various
* convTest functions provided for the circuit elements in the
* given circuit
*/
int
CKTconvTest(CKTcircuit *ckt)
{
extern SPICEdev **DEVices;
int i;
int error = OK;
#ifdef PARALLEL_ARCH
int ibuf[2];
long type = MT_CONV, length = 2;
#endif /* PARALLEL_ARCH */
for (i=0;i<DEVmaxnum;i++) {
if (((*DEVices[i]).DEVconvTest != NULL) && (ckt->CKThead[i] != NULL)) {
error = (*((*DEVices[i]).DEVconvTest))(ckt->CKThead[i],ckt);
}
#ifdef PARALLEL_ARCH
if (error || ckt->CKTnoncon) goto combine;
#else
if (error) return(error);
if (ckt->CKTnoncon) {
/* printf("convTest: device %s failed\n",
(*DEVices[i]).DEVpublic.name); */
return(OK);
}
#endif /* PARALLEL_ARCH */
}
#ifdef PARALLEL_ARCH
combine:
/* See if any of the DEVconvTest functions bailed. If not, proceed. */
ibuf[0] = error;
ibuf[1] = ckt->CKTnoncon;
IGOP_( &type, ibuf, &length, "+" );
ckt->CKTnoncon = ibuf[1];
if ( ibuf[0] != error ) {
error = E_MULTIERR;
}
return (error);
#else
return(OK);
#endif /* PARALLEL_ARCH */
}