1 /**********
2 Copyright 1990 Regents of the University of California.  All rights reserved.
3 Author: 1985 Thomas L. Quarles
4 **********/
5 
6     /* CKTtrunc(ckt)
7      * this is a driver program to iterate through all the various
8      * truncation error functions provided for the circuit elements in the
9      * given circuit
10      */
11 
12 #include "ngspice/ngspice.h"
13 #include "ngspice/cktdefs.h"
14 #include "ngspice/smpdefs.h"
15 #include "ngspice/devdefs.h"
16 #include "ngspice/sperror.h"
17 
18 
19 int
CKTtrunc(CKTcircuit * ckt,double * timeStep)20 CKTtrunc(CKTcircuit *ckt, double *timeStep)
21 {
22 #ifndef NEWTRUNC
23     int i;
24     double timetemp;
25 #ifdef STEPDEBUG
26     double debugtemp;
27 #endif /* STEPDEBUG */
28     double startTime;
29     int error = OK;
30 
31     startTime = SPfrontEnd->IFseconds();
32 
33     timetemp = HUGE;
34     for (i=0;i<DEVmaxnum;i++) {
35         if (DEVices[i] && DEVices[i]->DEVtrunc && ckt->CKThead[i]) {
36 #ifdef STEPDEBUG
37             debugtemp = timetemp;
38 #endif /* STEPDEBUG */
39 	    error = DEVices[i]->DEVtrunc (ckt->CKThead[i], ckt, &timetemp);
40 	    if(error) {
41                 ckt->CKTstat->STATtranTruncTime += SPfrontEnd->IFseconds()
42                     - startTime;
43                 return(error);
44             }
45 #ifdef STEPDEBUG
46             if(debugtemp != timetemp) {
47                 printf("timestep cut by device type %s from %g to %g\n",
48                         DEVices[i]->DEVpublic.name, debugtemp, timetemp);
49             }
50 #endif /* STEPDEBUG */
51         }
52     }
53     *timeStep = MIN(2 * *timeStep,timetemp);
54 
55     ckt->CKTstat->STATtranTruncTime += SPfrontEnd->IFseconds() - startTime;
56     return(OK);
57 #else /* NEWTRUNC */
58     int i;
59     CKTnode *node;
60     double timetemp;
61     double tmp;
62     double diff;
63     double tol;
64     double startTime;
65     int size;
66 
67     startTime = SPfrontEnd->IFseconds();
68 
69     timetemp = HUGE;
70     size = SMPmatSize(ckt->CKTmatrix);
71 #ifdef STEPDEBUG
72     printf("at time %g, delta %g\n",ckt->CKTtime,ckt->CKTdeltaOld[0]);
73 #endif
74     node = ckt->CKTnodes;
75     switch(ckt->CKTintegrateMethod) {
76 
77     case TRAPEZOIDAL:
78         switch(ckt->CKTorder) {
79         case 1:
80             for(i=1;i<size;i++) {
81                 tol = MAX( fabs(ckt->CKTrhs[i]),fabs(ckt->CKTpred[i]))*
82                         ckt->CKTlteReltol+ckt->CKTlteAbstol;
83                 node = node->next;
84                 if(node->type!= SP_VOLTAGE) continue;
85                 diff = ckt->CKTrhs[i]-ckt->CKTpred[i];
86 #ifdef STEPDEBUG
87                 printf("%s: cor=%g, pred=%g ",node->name,
88                         ckt->CKTrhs[i],ckt->CKTpred[i]);
89 #endif
90                 if(diff != 0) {
91                     tmp = ckt->CKTtrtol * tol * 2 /diff;
92                     tmp = ckt->CKTdeltaOld[0]*sqrt(fabs(tmp));
93                     timetemp = MIN(timetemp,tmp);
94 #ifdef STEPDEBUG
95                     printf("tol = %g, diff = %g, h->%g\n",tol,diff,tmp);
96 #endif
97                 } else {
98 #ifdef STEPDEBUG
99                     printf("diff is 0\n");
100 #endif
101                 }
102             }
103             break;
104         case 2:
105             for(i=1;i<size;i++) {
106                 tol = MAX( fabs(ckt->CKTrhs[i]),fabs(ckt->CKTpred[i]))*
107                         ckt->CKTlteReltol+ckt->CKTlteAbstol;
108                 node = node->next;
109                 if(node->type!= SP_VOLTAGE) continue;
110                 diff = ckt->CKTrhs[i]-ckt->CKTpred[i];
111 #ifdef STEPDEBUG
112                 printf("%s: cor=%g, pred=%g ",node->name,ckt->CKTrhs[i],
113                         ckt->CKTpred[i]);
114 #endif
115                 if(diff != 0) {
116                     tmp = ckt->CKTdeltaOld[0]*ckt->CKTtrtol * tol * 3 *
117                             (ckt->CKTdeltaOld[0]+ckt->CKTdeltaOld[1])/diff;
118                     tmp = fabs(tmp);
119                     timetemp = MIN(timetemp,tmp);
120 #ifdef STEPDEBUG
121                     printf("tol = %g, diff = %g, h->%g\n",tol,diff,tmp);
122 #endif
123                 } else {
124 #ifdef STEPDEBUG
125                     printf("diff is 0\n");
126 #endif
127                 }
128             }
129             break;
130         default:
131             return(E_ORDER);
132         break;
133 
134         }
135     break;
136 
137     case GEAR: {
138         double delsum=0;
139         for(i=0;i<=ckt->CKTorder;i++) {
140             delsum += ckt->CKTdeltaOld[i];
141         }
142         for(i=1;i<size;i++) {
143             node = node->next;
144             if(node->type!= SP_VOLTAGE) continue;
145             tol = MAX( fabs(ckt->CKTrhs[i]),fabs(ckt->CKTpred[i]))*
146                     ckt->CKTlteReltol+ckt->CKTlteAbstol;
147             diff = (ckt->CKTrhs[i]-ckt->CKTpred[i]);
148 #ifdef STEPDEBUG
149             printf("%s: cor=%g, pred=%g ",node->name,ckt->CKTrhs[i],
150                     ckt->CKTpred[i]);
151 #endif
152             if(diff != 0) {
153                 tmp = tol*ckt->CKTtrtol*delsum/(diff*ckt->CKTdelta);
154                 tmp = fabs(tmp);
155                 switch(ckt->CKTorder) {
156                     case 0:
157                         break;
158                     case 1:
159                         tmp = sqrt(tmp);
160                         break;
161                     default:
162                         tmp = exp(log(tmp)/(ckt->CKTorder+1));
163                         break;
164                 }
165                 tmp *= ckt->CKTdelta;
166                 timetemp = MIN(timetemp,tmp);
167 #ifdef STEPDEBUG
168                 printf("tol = %g, diff = %g, h->%g\n",tol,diff,tmp);
169 #endif
170             } else {
171 #ifdef STEPDEBUG
172                 printf("diff is 0\n");
173 #endif
174             }
175         }
176     }
177     break;
178 
179     default:
180         return(E_METHOD);
181 
182     }
183     *timeStep = MIN(2 * *timeStep,timetemp);
184     ckt->CKTstat->STATtranTruncTime += SPfrontEnd->IFseconds() - startTime;
185     return(OK);
186 #endif /* NEWTRUNC */
187 }
188