1 /**********
2 Author: 2010-05 Stefano Perticaroli ``spertica''
3 First Review: 2012-02 Francesco Lannutti and Stefano Perticaroli ``spertica''
4 Second Review: 2012-10 Stefano Perticaroli ``spertica'' and Francesco Lannutti
5 **********/
6
7 /* Include files for the PSS analysis */
8 #include "ngspice/ngspice.h"
9 #include "ngspice/cktdefs.h"
10 #include "cktaccept.h"
11 #include "ngspice/pssdefs.h"
12 #include "ngspice/sperror.h"
13 #include "ngspice/fteext.h"
14
15 #ifdef XSPICE
16 /* gtri - add - wbk - Add headers */
17 #include "ngspice/miftypes.h"
18
19 #include "ngspice/evt.h"
20 #include "ngspice/enh.h"
21 #include "ngspice/mif.h"
22 #include "ngspice/evtproto.h"
23 #include "ngspice/ipctiein.h"
24 /* gtri - end - wbk - Add headers */
25 #endif
26
27 #ifdef CLUSTER
28 #include "ngspice/cluster.h"
29 #endif
30
31
32 #define INIT_STATS() \
33 do { \
34 startTime = SPfrontEnd->IFseconds(); \
35 startIters = ckt->CKTstat->STATnumIter; \
36 startdTime = ckt->CKTstat->STATdecompTime; \
37 startsTime = ckt->CKTstat->STATsolveTime; \
38 startlTime = ckt->CKTstat->STATloadTime; \
39 startkTime = ckt->CKTstat->STATsyncTime; \
40 } while(0)
41
42 #define UPDATE_STATS(analysis) \
43 do { \
44 ckt->CKTcurrentAnalysis = analysis; \
45 ckt->CKTstat->STATtranTime += SPfrontEnd->IFseconds() - startTime; \
46 ckt->CKTstat->STATtranIter += ckt->CKTstat->STATnumIter - startIters; \
47 ckt->CKTstat->STATtranDecompTime += ckt->CKTstat->STATdecompTime - startdTime; \
48 ckt->CKTstat->STATtranSolveTime += ckt->CKTstat->STATsolveTime - startsTime; \
49 ckt->CKTstat->STATtranLoadTime += ckt->CKTstat->STATloadTime - startlTime; \
50 ckt->CKTstat->STATtranSyncTime += ckt->CKTstat->STATsyncTime - startkTime; \
51 } while(0)
52
53
54 /* Define some useful macro */
55 #define HISTORY 1024
56 #define GF_LAST 313
57
58
59 static int
60 DFT(long int, int, double *, double *, double *, double, double *, double *, double *, double *, double *);
61
62
63 int
DCpss(CKTcircuit * ckt,int restart)64 DCpss(CKTcircuit *ckt,
65 int restart) /* forced restart flag */
66 {
67 PSSan *job = (PSSan *) ckt->CKTcurJob;
68
69 int i;
70 double olddelta;
71 double delta;
72 double newdelta;
73 double *temp;
74 double startdTime;
75 double startsTime;
76 double startlTime;
77 double startkTime;
78 double startTime;
79 int startIters;
80 int converged;
81 int firsttime;
82 int error;
83 int save_order;
84 long save_mode;
85 IFuid timeUid;
86 IFuid *nameList;
87 int numNames;
88 double maxstepsize = 0.0;
89
90 int ltra_num;
91 CKTnode *node;
92 #ifdef XSPICE
93 /* gtri - add - wbk - 12/19/90 - Add IPC stuff */
94 Ipc_Boolean_t ipc_firsttime = IPC_TRUE;
95 Ipc_Boolean_t ipc_secondtime = IPC_FALSE;
96 Ipc_Boolean_t ipc_delta_cut = IPC_FALSE;
97 double ipc_last_time = 0.0;
98 double ipc_last_delta = 0.0;
99 /* gtri - end - wbk - 12/19/90 - Add IPC stuff */
100 #endif
101 #ifdef CLUSTER
102 int redostep;
103 #endif
104
105 /* New variables */
106 int j, oscnNode;
107 IFuid freqUid;
108
109 enum {STABILIZATION, SHOOTING, PSS} pss_state = STABILIZATION;
110
111 double err = 0, predsum = 0 ;
112 double time_temp = 0, gf_history [HISTORY], rr_history [HISTORY], predsum_history [HISTORY], nextstep ;
113 int msize, shooting_cycle_counter = 0;
114 double *RHS_copy_se, *RHS_copy_der, *RHS_derivative, *pred, err_0 = HUGE_VAL ;
115 double time_err_min_1 = 0, time_err_min_0 = 0, err_min_0 = HUGE_VAL, err_min_1 = 0 ;
116 double err_1 = 0, err_max = HUGE_VAL ;
117 int pss_points_cycle = 0, dynamic_test = 0 ;
118 double gf_last_0 = HUGE_VAL, gf_last_1 = GF_LAST ;
119 double thd = 0 ;
120 double *psstimes, *pssvalues;
121 double *RHS_max, *RHS_min, *err_conv ;
122
123 /* Francesco Lannutti's MOD */
124 /* Stuff needed by frequency estimation reiteration, based on the DFT result */
125 int position;
126 double max_freq;
127
128
129 /* Print some useful information */
130 fprintf (stderr, "Periodic Steady State Analysis Started\n\n") ;
131 fprintf (stderr, "PSS Guessed Frequency %g\n", ckt->CKTguessedFreq) ;
132 fprintf (stderr, "PSS Points %ld\n", ckt->CKTpsspoints) ;
133 fprintf (stderr, "PSS Harmonics number %d\n", ckt->CKTharms) ;
134 fprintf (stderr, "PSS Steady Coefficient %g\n", ckt->CKTsteady_coeff) ;
135 fprintf (stderr, "PSS sc_iter %d\n", ckt->CKTsc_iter) ;
136 fprintf (stderr, "PSS Stabilization Time %g\n", ckt->CKTstabTime) ;
137
138
139 oscnNode = job->PSSoscNode->number ;
140
141
142 /* Variables and memory initialization */
143
144 for (i = 0 ; i < HISTORY ; i++)
145 {
146 rr_history [i] = 0.0 ;
147 gf_history [i] = 0.0 ;
148 }
149
150 msize = SMPmatSize (ckt->CKTmatrix) ;
151 RHS_copy_se = TMALLOC (double, msize) ; /* Set the current RHS reference for next Shooting Evaluation */
152 RHS_copy_der = TMALLOC (double, msize) ; /* Used to compute current Derivative */
153 RHS_derivative = TMALLOC (double, msize) ;
154 pred = TMALLOC (double, msize) ;
155 RHS_max = TMALLOC (double, msize) ;
156 RHS_min = TMALLOC (double, msize) ;
157 err_conv = TMALLOC (double, msize) ;
158
159 for (i = 0 ; i < msize ; i++)
160 {
161 RHS_copy_se [i] = 0.0 ;
162 RHS_copy_der [i] = 0.0 ;
163 RHS_derivative [i] = 0.0 ;
164 pred [i] = 0.0 ;
165 }
166
167 psstimes = TMALLOC (double, ckt->CKTpsspoints + 1) ;
168 pssvalues = TMALLOC (double, msize * (ckt->CKTpsspoints + 1)) ;
169
170 for (i = 0 ; i < ckt->CKTpsspoints + 1 ; i++)
171 psstimes [i] = 0.0 ;
172
173 for (i = 0 ; i < msize * (ckt->CKTpsspoints + 1) ; i++)
174 pssvalues [i] = 0.0 ;
175
176 /* Delta timestep and circuit time setup */
177 delta = ckt->CKTstep ;
178 ckt->CKTtime = ckt->CKTinitTime ;
179 ckt->CKTfinalTime = ckt->CKTstabTime ;
180
181 /* Starting PSS Algorithm, based on Transient Analysis */
182 if(restart || ckt->CKTtime == 0) {
183 delta = MIN (1 / ckt->CKTguessedFreq / 100, ckt->CKTstep) ;
184
185 #ifdef STEPDEBUG
186 fprintf (stderr, "delta = %g finalTime/200: %g CKTstep: %g\n", delta, ckt->CKTfinalTime / 200, ckt->CKTstep) ;
187 #endif
188 /* begin LTRA code addition */
189 if (ckt->CKTtimePoints != NULL)
190 FREE(ckt->CKTtimePoints);
191
192 if (ckt->CKTstep >= ckt->CKTmaxStep)
193 maxstepsize = ckt->CKTstep;
194 else
195 maxstepsize = ckt->CKTmaxStep;
196
197 ckt->CKTsizeIncr = 10;
198 ckt->CKTtimeIndex = -1; /* before the DC soln has been stored */
199 ckt->CKTtimeListSize = (int)(1 / ckt->CKTguessedFreq / maxstepsize + 0.5);
200 ltra_num = CKTtypelook("LTRA");
201 if (ltra_num >= 0 && ckt->CKThead[ltra_num] != NULL)
202 ckt->CKTtimePoints = TMALLOC(double, ckt->CKTtimeListSize);
203 /* end LTRA code addition */
204
205 /* Breakpoints initialization */
206 if(ckt->CKTbreaks) FREE(ckt->CKTbreaks);
207 ckt->CKTbreaks = TMALLOC(double, 2);
208 if(ckt->CKTbreaks == NULL) return(E_NOMEM);
209 ckt->CKTbreaks[0] = 0;
210 ckt->CKTbreaks[1] = ckt->CKTfinalTime;
211 ckt->CKTbreakSize = 2;
212
213 #ifdef XSPICE
214 /* gtri - begin - wbk - 12/19/90 - Modify setting of CKTminBreak */
215 /* if (ckt->CKTminBreak == 0)
216 ckt->CKTminBreak = ckt->CKTmaxStep * 5e-5 ; */
217
218 /* Set to 10 times delmin for ATESSE 1 compatibity */
219 if(ckt->CKTminBreak==0) ckt->CKTminBreak = 10.0 * ckt->CKTdelmin;
220 /* gtri - end - wbk - 12/19/90 - Modify setting of CKTminBreak */
221 #else
222 /* Minimum Breakpoint Setup */
223 if(ckt->CKTminBreak==0) ckt->CKTminBreak=ckt->CKTmaxStep*5e-5;
224 #endif
225
226 #ifdef XSPICE
227 /* gtri - add - wbk - 12/19/90 - Add IPC stuff and set anal_init and anal_type */
228 /* Tell the beginPlot routine what mode we're in */
229 g_ipc.anal_type = IPC_ANAL_TRAN;
230
231 /* Tell the code models what mode we're in */
232 g_mif_info.circuit.anal_type = MIF_DC;
233
234 g_mif_info.circuit.anal_init = MIF_TRUE;
235 /* gtri - end - wbk */
236 #endif
237
238 /* Time Domain plot start and prepared to be filled in later */
239 error = CKTnames(ckt,&numNames,&nameList);
240 if(error) return(error);
241 SPfrontEnd->IFnewUid (ckt, &timeUid, NULL, "time", UID_OTHER, NULL);
242 error = SPfrontEnd->OUTpBeginPlot (ckt, ckt->CKTcurJob,
243 "Time Domain Periodic Steady State Analysis",
244 timeUid, IF_REAL,
245 numNames, nameList, IF_REAL,
246 &(job->PSSplot_td));
247 tfree(nameList);
248 if(error) return(error);
249
250 /* Time initialization for Transient Analysis */
251 ckt->CKTtime = 0;
252 ckt->CKTdelta = 0;
253 ckt->CKTbreak = 1;
254 firsttime = 1;
255 save_mode = (ckt->CKTmode&MODEUIC) | MODETRANOP | MODEINITJCT;
256 save_order = ckt->CKTorder;
257
258 #ifdef XSPICE
259 /* gtri - begin - wbk - set a breakpoint at end of supply ramping time */
260 /* must do this after CKTtime set to 0 above */
261 if(ckt->enh->ramp.ramptime > 0.0)
262 CKTsetBreak(ckt, ckt->enh->ramp.ramptime);
263 /* gtri - end - wbk - set a breakpoint at end of supply ramping time */
264
265 /* gtri - begin - wbk - Call EVTop if event-driven instances exist */
266 if(ckt->evt->counts.num_insts != 0) {
267 /* use new DCOP algorithm */
268 converged = EVTop(ckt,
269 (ckt->CKTmode & MODEUIC) | MODETRANOP | MODEINITJCT,
270 (ckt->CKTmode & MODEUIC) | MODETRANOP | MODEINITFLOAT,
271 ckt->CKTdcMaxIter,
272 MIF_TRUE);
273 EVTdump(ckt, IPC_ANAL_DCOP, 0.0);
274
275 EVTop_save(ckt, MIF_FALSE, 0.0);
276
277 /* gtri - end - wbk - Call EVTop if event-driven instances exist */
278 } else
279 #endif
280
281 /* Looking for a working Operating Point */
282 converged = CKTop(ckt,
283 (ckt->CKTmode & MODEUIC) | MODETRANOP | MODEINITJCT,
284 (ckt->CKTmode & MODEUIC) | MODETRANOP | MODEINITFLOAT,
285 ckt->CKTdcMaxIter);
286
287 #ifdef STEPDEBUG
288 if(converged != 0) {
289 fprintf(stdout,"\nTransient solution failed -\n");
290 CKTncDump(ckt);
291 fprintf(stdout,"\n");
292 fflush(stdout);
293 } else if (!ft_noacctprint && !ft_noinitprint) {
294 fprintf(stdout,"\nInitial Transient Solution\n");
295 fprintf(stdout,"--------------------------\n\n");
296 fprintf(stdout,"%-30s %15s\n", "Node", "Voltage");
297 fprintf(stdout,"%-30s %15s\n", "----", "-------");
298 for(node=ckt->CKTnodes->next;node;node=node->next) {
299 if (strstr(node->name, "#branch") || !strchr(node->name, '#'))
300 fprintf(stdout,"%-30s %15g\n", node->name,
301 ckt->CKTrhsOld[node->number]);
302 }
303 fprintf(stdout,"\n");
304 fflush(stdout);
305 }
306 #endif
307
308 /* If no convergence reached - NO valid Operating Point */
309 if(converged != 0) return(converged);
310 #ifdef XSPICE
311 /* gtri - add - wbk - 12/19/90 - Add IPC stuff */
312
313 /* Send the operating point results for Mspice compatibility */
314 if(g_ipc.enabled) {
315 ipc_send_dcop_prefix();
316 CKTdump(ckt, 0.0, job->PSSplot_td);
317 ipc_send_dcop_suffix();
318 }
319
320 /* gtri - end - wbk */
321
322 /* gtri - add - wbk - 12/19/90 - set anal_init and anal_type */
323
324 g_mif_info.circuit.anal_init = MIF_TRUE;
325
326 /* Tell the code models what mode we're in */
327 g_mif_info.circuit.anal_type = MIF_TRAN;
328
329 /* gtri - end - wbk */
330
331 /* gtri - begin - wbk - Add Breakpoint stuff */
332
333 /* Initialize the temporary breakpoint variables to infinity */
334 g_mif_info.breakpoint.current = HUGE_VAL;
335 g_mif_info.breakpoint.last = HUGE_VAL;
336
337 /* gtri - end - wbk - Add Breakpoint stuff */
338 #endif
339 ckt->CKTstat->STATtimePts ++;
340
341 /* Setting Integration Order to Backward Euler */
342 ckt->CKTorder = 1;
343
344 /* Copying the maxStep to every deltaOld */
345 for(i=0;i<7;i++) {
346 ckt->CKTdeltaOld[i]=ckt->CKTmaxStep;
347 }
348
349 /* Setting DELTA */
350 ckt->CKTdelta = delta;
351 #ifdef STEPDEBUG
352 fprintf (stderr, "delta initialized to %g\n", ckt->CKTdelta);
353 #endif
354
355 ckt->CKTsaveDelta = ckt->CKTfinalTime/50;
356
357 ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODETRAN | MODEINITTRAN;
358 /* Changing Circuit MODE */
359 /* modeinittran set here */
360 ckt->CKTag[0]=ckt->CKTag[1]=0;
361
362 /* State0 copied into State1 - DEPRECATED LEGACY function - to be replaced with memmove() */
363 memcpy(ckt->CKTstate1, ckt->CKTstate0,
364 (size_t) ckt->CKTnumStates * sizeof(double));
365
366 /* Statistics Initialization using a macro at the beginning of this code */
367 INIT_STATS();
368 #ifdef CLUSTER
369 CLUsetup(ckt);
370 #endif
371 } else {
372 /* saj As traninit resets CKTmode */
373 ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODETRAN | MODEINITPRED;
374 /* saj */
375 INIT_STATS();
376 if(ckt->CKTminBreak==0) ckt->CKTminBreak=ckt->CKTmaxStep*5e-5;
377 firsttime=0;
378 /* To get rawfile working saj*/
379 error = SPfrontEnd->OUTpBeginPlot (NULL, NULL,
380 NULL,
381 NULL, 0,
382 666, NULL, 666,
383 &(job->PSSplot_td));
384 if(error) {
385 fprintf(stderr, "Couldn't relink rawfile\n");
386 return error;
387 }
388 /* end saj*/
389
390 /* Skip nextTime if it isn't the firsttime! :) */
391 goto resume;
392 }
393
394 /* 650 */
395 nextTime:
396
397 /* begin LTRA code addition */
398 if (ckt->CKTtimePoints) {
399 ckt->CKTtimeIndex++;
400 if (ckt->CKTtimeIndex >= ckt->CKTtimeListSize) {
401 /* need more space */
402 int need;
403 if (pss_state == STABILIZATION)
404 need = (int) ceil((ckt->CKTstabTime - ckt->CKTtime) / maxstepsize ) ;
405 else
406 need = (int) ceil((time_temp + 1 / ckt->CKTguessedFreq - ckt->CKTtime) / maxstepsize) ;
407
408 if (need < ckt->CKTsizeIncr)
409 need = ckt->CKTsizeIncr;
410 ckt->CKTtimeListSize += need;
411 ckt->CKTtimePoints = TREALLOC(double, ckt->CKTtimePoints, ckt->CKTtimeListSize);
412 ckt->CKTsizeIncr = (int) ceil(1.4 * ckt->CKTsizeIncr);
413 }
414 ckt->CKTtimePoints[ckt->CKTtimeIndex] = ckt->CKTtime;
415 }
416 /* end LTRA code addition */
417
418 /* Check for the timepoint acceptance */
419 error = CKTaccept(ckt);
420 /* check if current breakpoint is outdated; if so, clear */
421 if (ckt->CKTtime > ckt->CKTbreaks[0]) CKTclrBreak(ckt);
422
423 /*
424 * Breakpoint handling scheme:
425 * When a timepoint t is accepted (by CKTaccept), clear all previous
426 * breakpoints, because they will never be needed again.
427 *
428 * t may itself be a breakpoint, or indistinguishably close. DON'T
429 * clear t itself; recognise it as a breakpoint and act accordingly
430 *
431 * if t is not a breakpoint, limit the timestep so that the next
432 * breakpoint is not crossed
433 */
434
435 #ifdef STEPDEBUG
436 fprintf (stderr, "Delta %g accepted at time %g (finaltime: %g)\n", ckt->CKTdelta, ckt->CKTtime, ckt->CKTfinalTime) ;
437 fflush(stdout);
438 #endif /* STEPDEBUG */
439 ckt->CKTstat->STATaccepted ++;
440 ckt->CKTbreak = 0;
441 /* XXX Error will cause single process to bail. */
442 if(error) {
443 UPDATE_STATS(DOING_TRAN);
444 return(error);
445 }
446 #ifdef XSPICE
447 /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */
448
449 if ((g_ipc.enabled) || wantevtdata) {
450
451 if (pss_state == PSS)
452 {
453 /* Send event-driven results */
454 EVTdump(ckt, IPC_ANAL_TRAN, 0.0);
455
456 /* Then follow with analog results... */
457
458 /* Test to see if delta was cut by a breakpoint, */
459 /* a non-convergence, or a too large truncation error */
460 if(ipc_firsttime)
461 ipc_delta_cut = IPC_FALSE;
462 else if(ckt->CKTtime < (ipc_last_time + (0.999 * ipc_last_delta)))
463 ipc_delta_cut = IPC_TRUE;
464 else
465 ipc_delta_cut = IPC_FALSE;
466
467 /* Record the data required to check for delta cuts */
468 ipc_last_time = ckt->CKTtime;
469 ipc_last_delta = MIN(ckt->CKTdelta, ckt->CKTmaxStep);
470
471 /* Send results data if time since last dump is greater */
472 /* than 'mintime', or if first or second timepoints, */
473 /* or if delta was cut */
474 if( (ckt->CKTtime >= (g_ipc.mintime + g_ipc.last_time)) ||
475 ipc_firsttime || ipc_secondtime || ipc_delta_cut ) {
476
477 ipc_send_data_prefix(ckt->CKTtime);
478 CKTdump(ckt, ckt->CKTtime, job->PSSplot_td);
479 ipc_send_data_suffix();
480
481 if(ipc_firsttime) {
482 ipc_firsttime = IPC_FALSE;
483 ipc_secondtime = IPC_TRUE;
484 } else if(ipc_secondtime) {
485 ipc_secondtime = IPC_FALSE;
486 }
487
488 g_ipc.last_time = ckt->CKTtime;
489 }
490 }
491 } else
492 /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */
493 #endif
494 #ifdef CLUSTER
495 if (pss_state == PSS)
496 CLUoutput(ckt);
497 #endif
498
499 if (pss_state == PSS)
500 {
501 nextstep = time_temp + 1 / ckt->CKTguessedFreq * ((double)(pss_points_cycle) / (double)ckt->CKTpsspoints) ;
502
503 /* If in_pss, store data for Time Domain Plot and gather ordered data for FFT computing */
504 if ((AlmostEqualUlps (ckt->CKTtime, nextstep, 10)) || (ckt->CKTtime > time_temp + 1 / ckt->CKTguessedFreq))
505 {
506
507 #ifdef STEPDEBUG
508 fprintf (stderr, "IN_PSS: time point accepted in evolution for FFT calculations.\n") ;
509 fprintf (stderr, "Circuit time %1.15g, final time %1.15g, point index %d and total requested points %ld\n",
510 ckt->CKTtime, nextstep, pss_points_cycle, ckt->CKTpsspoints) ;
511 #endif
512
513 CKTdump (ckt, ckt->CKTtime, job->PSSplot_td) ;
514
515 /* Store the Time Base for the DFT */
516 psstimes [pss_points_cycle] = ckt->CKTtime ;
517
518 /* Store values for the FFT calculation */
519 for (i = 1 ; i <= msize ; i++)
520 pssvalues [i - 1 + pss_points_cycle * msize] = ckt->CKTrhsOld [i] ;
521
522 /* Update PSS counter cycle, used to stop the entire algorithm */
523 pss_points_cycle++ ;
524
525 /* Set the next BreakPoint for PSS */
526 CKTsetBreak (ckt, time_temp + (1 / ckt->CKTguessedFreq) * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ;
527
528 #ifdef STEPDEBUG
529 fprintf (stderr, "Next breakpoint set in: %1.15g\n", time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ;
530 #endif
531
532 } else {
533 /* Algo can enter here but should do nothing */
534
535 #ifdef STEPDEBUG
536 fprintf (stderr, "IN_PSS: time point accepted in evolution but dropped for FFT calculations\n") ;
537 #endif
538
539 }
540 }
541
542 #ifdef XSPICE
543 /* gtri - begin - wbk - Update event queues/data for accepted timepoint */
544 /* Note: this must be done AFTER sending results to SI so it can't */
545 /* go next to CKTaccept() above */
546 if(ckt->evt->counts.num_insts > 0)
547 EVTaccept(ckt, ckt->CKTtime);
548 /* gtri - end - wbk - Update event queues/data for accepted timepoint */
549 #endif
550 ckt->CKTstat->STAToldIter = ckt->CKTstat->STATnumIter;
551
552 /* ***********************************/
553 /* ******* SHOOTING CODE BLOCK *******/
554 /* ***********************************/
555 switch(pss_state) {
556
557 case STABILIZATION:
558 {
559 /* Test if stabTime has been reached */
560 if (AlmostEqualUlps (ckt->CKTtime, ckt->CKTstabTime, 100))
561 {
562 time_temp = ckt->CKTtime ;
563
564 /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */
565 ckt->CKTfinalTime = time_temp + 2 / ckt->CKTguessedFreq ;
566 fprintf (stderr, "Exiting from stabilization\n") ;
567 fprintf (stderr, "Time of first shooting evaluation will be %1.10g\n", time_temp + 1 / ckt->CKTguessedFreq) ;
568
569 /* Next time is no more in stabilization - Unset the flag */
570 pss_state = SHOOTING;
571
572 /* Save the RHS_copy_der as the NEW CKTrhsOld */
573 for (i = 1 ; i <= msize ; i++)
574 RHS_copy_der [i - 1] = ckt->CKTrhsOld [i] ;
575
576 /* Print RHS on exiting from stabilization */
577 fprintf (stderr, "RHS on exiting from stabilization: ") ;
578 for (i = 1 ; i <= msize ; i++)
579 {
580 RHS_copy_se [i - 1] = ckt->CKTrhsOld [i] ;
581 fprintf (stderr, "%-15g ", RHS_copy_se [i - 1]) ;
582 }
583 fprintf (stderr, "\n") ;
584
585 /* RHS_max and RHS_min initialization - HUGE_VAL is the maximum machine error */
586 for (i = 0 ; i < msize ; i++)
587 {
588 RHS_max [i] = -HUGE_VAL ;
589 RHS_min [i] = HUGE_VAL ;
590 }
591 }
592 }
593 break;
594
595 case SHOOTING:
596 {
597 double offset, interval, nextBreak ;
598 /* Calculation of error norms of RHS solution of every accepted nextTime */
599 err = 0 ;
600 for (i = 0 ; i < msize ; i++)
601 {
602 /* Save max per node or branch of every estimated period */
603 if (RHS_max [i] < ckt->CKTrhsOld [i + 1])
604 RHS_max [i] = ckt->CKTrhsOld [i + 1] ;
605
606 /* Save min per node or branch of every estimated period */
607 if (RHS_min [i] > ckt->CKTrhsOld [i + 1])
608 RHS_min [i] = ckt->CKTrhsOld [i + 1] ;
609
610 /* CKTrhsOld is the last CORRECT value of RHS */
611 err_conv [i] = ckt->CKTrhsOld [i + 1] - RHS_copy_se [i] ;
612 err += err_conv [i] * err_conv [i] ;
613
614 /* Compute and store derivative */
615 RHS_derivative [i] = (ckt->CKTrhsOld [i + 1] - RHS_copy_der [i]) / ckt->CKTdelta ;
616
617 /* Save the RHS_copy_der as the NEW CKTrhsOld */
618 RHS_copy_der [i] = ckt->CKTrhsOld [i + 1] ;
619
620 #ifdef STEPDEBUG
621 fprintf (stderr, "Pred is so high or so low! Diff is: %g\n", err_conv [i]) ;
622 #endif
623
624 }
625 err = sqrt (err) ;
626
627 /* Start frequency estimation */
628 if ((err < err_0) && (ckt->CKTtime >= time_temp + 0.5 / ckt->CKTguessedFreq)) /* far enough from time temp... */
629 {
630 if (err < err_min_0)
631 {
632 err_min_1 = err_min_0 ; /* store previous minimum of RHS vector error */
633 err_min_0 = err ; /* store minimum of RHS vector error */
634 time_err_min_1 = time_err_min_0 ; /* store previous minimum of RHS vector error time */
635 time_err_min_0 = ckt->CKTtime ; /* store minimum of RHS vector error time */
636 }
637 }
638 err_0 = err ;
639
640 if ((err > err_1) && (ckt->CKTtime >= time_temp + 0.1 / ckt->CKTguessedFreq)) /* far enough from time temp... */
641 {
642 if (err > err_max)
643 err_max = err ; /* store maximum of RHS vector error */
644 }
645 err_1 = err ;
646
647
648 /* *************************************** */
649 /* ********** Breakpoint update ********** */
650 /* *************************************** */
651
652 /* Force the tran analysis to evaluate requested breakpoints. Breakpoints are even more closer as
653 the next occurence of guessed period is approaching. La lunga notte dei robot viventi... */
654
655 /* double offset, interval, nextBreak ;
656 int i ;
657 */
658 if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.995) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq)))
659 {
660 offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.995 ;
661 interval = (1 / ckt->CKTguessedFreq) * (1 - 0.995) * (ckt->CKTsteady_coeff / 10) ;
662 i = (int)((ckt->CKTtime - offset) / interval) ;
663 nextBreak = offset + (i + 1) * interval ;
664 CKTsetBreak (ckt, nextBreak) ;
665 }
666 else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.8) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.995))
667 {
668 offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.8 ;
669 interval = (1 / ckt->CKTguessedFreq) * (0.995 - 0.8) * (ckt->CKTsteady_coeff / 5) ;
670 i = (int)((ckt->CKTtime - offset) / interval) ;
671 nextBreak = offset + (i + 1) * interval ;
672 CKTsetBreak (ckt, nextBreak) ;
673 }
674 else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.5) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.8))
675 {
676 offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.5 ;
677 interval = (1 / ckt->CKTguessedFreq) * (0.8 - 0.5) * (ckt->CKTsteady_coeff / 3) ;
678 i = (int)((ckt->CKTtime - offset) / interval) ;
679 nextBreak = offset + (i + 1) * interval ;
680 CKTsetBreak (ckt, nextBreak) ;
681 }
682 else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.1) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.5))
683 {
684 offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.1 ;
685 interval = (1 / ckt->CKTguessedFreq) * (0.5 - 0.1) * (ckt->CKTsteady_coeff / 2) ;
686 i = (int)((ckt->CKTtime - offset) / interval) ;
687 nextBreak = offset + (i + 1) * interval ;
688 CKTsetBreak (ckt, nextBreak) ;
689 }
690 else if ((ckt->CKTtime > time_temp) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.1))
691 {
692 offset = time_temp ;
693 interval = (1 / ckt->CKTguessedFreq) * (0.1) * (ckt->CKTsteady_coeff) ;
694 i = (int)((ckt->CKTtime - offset) / interval) ;
695 nextBreak = offset + (i + 1) * interval ;
696 CKTsetBreak (ckt, nextBreak) ;
697 } else {
698 fprintf (stderr, "Strange behavior\n\n") ;
699 fprintf (stderr, "CKTtime: %g\ntime_temp: %g\n\n", ckt->CKTtime, time_temp) ;
700 }
701
702 /* *************************************** */
703 /* ******* END Breakpoint update ********* */
704 /* *************************************** */
705
706
707 /* If evolution is near shooting... */
708 if ((AlmostEqualUlps (ckt->CKTtime, time_temp + 1 / ckt->CKTguessedFreq, 10)) || (ckt->CKTtime > time_temp + 1 / ckt->CKTguessedFreq))
709 {
710 int excessive_err_nodes = 0 ;
711
712 /* Calculation of error norms of RHS solution of every accepted nextTime */
713 predsum = 0 ;
714 for (i = 0 ; i < msize ; i++)
715 {
716 /* Pitagora ha sempre ragione!!! :))) */
717 /* pred is treated as FREQUENCY to avoid numerical overflow when derivative is close to ZERO */
718 pred [i] = RHS_derivative [i] / err_conv [i] ;
719
720 #ifdef STEPDEBUG
721 fprintf (stderr, "Pred is so high or so low! Diff is: %g\n", err_conv [i]) ;
722 #endif
723
724 if ((fabs (pred [i]) > 1.0e6 * ckt->CKTguessedFreq) || (err_conv [i] == 0))
725 {
726 if (pred [i] > 0)
727 pred [i] = 1.0e6 * ckt->CKTguessedFreq ;
728 else
729 pred [i] = -1.0e6 * ckt->CKTguessedFreq ;
730 }
731
732 predsum += pred [i] ;
733
734 #ifdef STEPDEBUG
735 fprintf (stderr, "Predsum in time before to be divided by dynamic_test has value %g\n", 1 / predsum) ;
736 fprintf (stderr, "Current Diff: %g, Derivative: %g, Frequency Projection: %g\n", err_conv [i], RHS_derivative [i], pred [i]) ;
737 #endif
738
739 }
740
741 // int excessive_err_nodes = 0 ;
742
743 if (shooting_cycle_counter == 0)
744 {
745 /* If first time in shooting we warn about that ! */
746 fprintf (stderr, "In shooting...\n") ;
747 }
748
749 //#ifdef STEPDEBUG
750 /* For debugging purpose */
751 fprintf (stderr, "\n----------------\n") ;
752 fprintf (stderr, "Shooting cycle iteration number: %3d ||", shooting_cycle_counter) ;
753
754 if (shooting_cycle_counter > 0)
755 fprintf (stderr, " rr: %g || predsum: %g\n", rr_history [shooting_cycle_counter - 1], 1 / predsum) ;
756 else
757 fprintf (stderr, " rr: %g || predsum: %g\n", 0.0, 1 / predsum) ;
758
759 // fprintf (stderr, "Print of dynamically consistent nodes voltages or branches currents:\n") ;
760 /* --------------------- */
761 //#endif
762
763 for (i = 0, node = ckt->CKTnodes->next ; node ; i++, node = node->next)
764 {
765 /* Voltage Node */
766 if (!strchr (node->name, '#'))
767 {
768 if (fabs (err_conv [i]) > (fabs (RHS_max [i] - RHS_min [i]) * ckt->CKTreltol + ckt->CKTvoltTol) *
769 ckt->CKTtrtol * ckt->CKTsteady_coeff)
770 {
771 excessive_err_nodes++ ;
772 }
773
774 /* If the dynamic is below 10uV, it's dropped */
775 if (fabs (RHS_max [i] - RHS_min [i]) > 10 * 1e-6)
776 {
777 dynamic_test++ ; /* test on voltage dynamic consistence */
778 }
779
780 /* Current Node */
781 } else {
782 if (fabs (err_conv [i]) > (fabs (RHS_max [i] - RHS_min [i]) * ckt->CKTreltol + ckt->CKTabstol) *
783 ckt->CKTtrtol * ckt->CKTsteady_coeff)
784 {
785 excessive_err_nodes++ ;
786 }
787
788 /* If the dynamic is below 10nA, it's dropped */
789 if (fabs (RHS_max [i] - RHS_min [i]) > 10 * 1e-9)
790 {
791 dynamic_test++ ; /* test on current dynamic consistence */
792 }
793 }
794 }
795
796 if (dynamic_test == 0)
797 {
798 /* Test for dynamic existence */
799 fprintf (stderr, "No detectable dynamic on voltages nodes or currents branches. PSS analysis aborted\n") ;
800
801 /* Terminates plot in Time Domain and frees the allocated memory */
802 SPfrontEnd->OUTendPlot (job->PSSplot_td) ;
803 FREE (RHS_copy_se) ;
804 FREE (RHS_copy_der) ;
805 FREE (RHS_max) ;
806 FREE (RHS_min) ;
807 FREE (err_conv) ;
808 FREE (psstimes) ;
809 FREE (pssvalues) ;
810 return (E_PANIC) ; /* to be corrected with definition of new error macro in iferrmsg.h */
811 }
812 else if ((time_err_min_0 - time_temp) < 0)
813 {
814 /* Something has gone wrong... */
815 fprintf (stderr, "Cannot find a minimum for error vector in estimated period. Try to adjust tstab! PSS analysis aborted\n") ;
816
817 /* Terminates plot in Time Domain and frees the allocated memory */
818 SPfrontEnd->OUTendPlot (job->PSSplot_td) ;
819 FREE (RHS_copy_se) ;
820 FREE (RHS_copy_der) ;
821 FREE (RHS_max) ;
822 FREE (RHS_min) ;
823 FREE (err_conv) ;
824 FREE (psstimes) ;
825 FREE (pssvalues) ;
826 return (E_PANIC) ; /* to be corrected with definition of new error macro in iferrmsg.h */
827 }
828
829 //#ifdef STEPDEBUG
830 // fprintf (stderr, "Global Convergence Error reference: %g, Time Projection: %g.\n",
831 // err_conv_ref / dynamic_test, predsum) ;
832 //#endif
833
834 /* Take the mean value of time prediction trough the dynamic test variable - predsum becomes TIME */
835 predsum = 1 / (predsum * dynamic_test) ;
836
837 /* Store the predsum history as absolute value */
838 predsum_history [shooting_cycle_counter] = fabs (predsum) ;
839
840 /***********************************/
841 /*** FREQUENCY ESTIMATION UPDATE ***/
842 /***********************************/
843 if ((err_min_0 == err) || (err_min_0 == HUGE_VAL))
844 {
845 /* Enters here if guessed frequency is higher than the 'real' value */
846 ckt->CKTguessedFreq = 1 / (1 / ckt->CKTguessedFreq + fabs (predsum)) ;
847
848 #ifdef STEPDEBUG
849 fprintf (stderr, "Frequency DOWN: est per %g, err min %g, err min 1 %g, err max %g, err %g\n",
850 time_err_min_0 - time_temp, err_min_0, err_min_1, err_max, err) ;
851 #endif
852
853 } else {
854 /* Enters here if guessed frequency is lower than the 'real' value */
855 ckt->CKTguessedFreq = 1 / (time_err_min_0 - time_temp) ;
856
857 #ifdef STEPDEBUG
858 fprintf (stderr, "Frequency UP: est per %g, err min %g, err min 1 %g, err max %g, err %g\n",
859 time_err_min_0 - time_temp, err_min_0, err_min_1, err_max, err) ;
860 #endif
861
862 }
863
864 /* Temporary variables to store previous occurrence of guessed frequency */
865 gf_last_1 = gf_last_0 ;
866 gf_last_0 = ckt->CKTguessedFreq ;
867
868 /* Next evaluation of shooting will be updated time (time_temp) summed to updated guessed period */
869 time_temp = ckt->CKTtime ;
870
871 /* Store error history */
872 rr_history [shooting_cycle_counter] = err ;
873 gf_history [shooting_cycle_counter] = ckt->CKTguessedFreq ;
874 shooting_cycle_counter++ ;
875
876 fprintf (stderr, "Updated guessed frequency: %1.10lg .\n", ckt->CKTguessedFreq) ;
877 fprintf (stderr, "Next shooting evaluation time is %1.10g and current time is %1.10g.\n",
878 time_temp + 1 / ckt->CKTguessedFreq, ckt->CKTtime) ;
879
880 /* Restore maximum and minimum error for next search */
881 err_min_0 = HUGE_VAL ;
882 err_max = -HUGE_VAL ;
883 err_0 = HUGE_VAL ;
884 err_1 = -HUGE_VAL ;
885 dynamic_test = 0 ;
886
887 /* Reset actual RHS reference for next shooting evaluation */
888 for (i = 1 ; i <= msize ; i++)
889 RHS_copy_se [i - 1] = ckt->CKTrhsOld [i] ;
890
891 #ifdef STEPDEBUG
892 fprintf (stderr, "RHS on new shooting cycle: ") ;
893 for (i = 0 ; i < msize ; i++)
894 fprintf (stderr, "%-15g ", RHS_copy_se [i]) ;
895 fprintf (stderr, "\n") ;
896 #endif
897
898 for (i = 0 ; i < msize ; i++)
899 {
900 /* Reset max and min per node or branch on every shooting cycle */
901 RHS_max [i] = -HUGE_VAL ;
902 RHS_min [i] = HUGE_VAL ;
903 }
904
905 fprintf (stderr, "----------------\n\n") ;
906
907 /* Shooting Exit Condition */
908 if ((shooting_cycle_counter > ckt->CKTsc_iter) || (excessive_err_nodes == 0))
909 {
910 int k ;
911 double minimum ;
912
913 pss_state = PSS ;
914
915 #ifdef STEPDEBUG
916 fprintf (stderr, "\nFrequency estimation (FE) and RHS period residual (PR) evolution\n") ;
917 #endif
918
919 // minimum = rr_history [0] ;
920 minimum = predsum_history [0] ;
921 k = 0 ;
922 for (i = 0 ; i < shooting_cycle_counter ; i++)
923 {
924 /* Print some statistics */
925 fprintf (stderr, "%-3d -> FE: %-15.10g || RR: %15.10g", i, gf_history [i], rr_history [i]) ;
926
927 /* Take the minimum residual iteration */
928 // if (minimum > rr_history [i])
929 if (minimum > predsum_history [i])
930 {
931 // minimum = rr_history [i] ;
932 minimum = predsum_history [i] ;
933 k = i ;
934 }
935 fprintf (stderr, " || predsum/dynamic_test: %15.10g || minimum: %15.10g\n", predsum_history [i], minimum) ;
936 }
937
938 if (excessive_err_nodes == 0) /* SHOOTING has converged */
939 ckt->CKTguessedFreq = gf_history [shooting_cycle_counter - 1] ;
940 else
941 ckt->CKTguessedFreq = gf_history [k] ;
942
943 /* Save the current Time */
944 time_temp = ckt->CKTtime ;
945
946 /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */
947 ckt->CKTfinalTime = time_temp + 1 / ckt->CKTguessedFreq ;
948
949 /* Dump the first PSS point for the FFT */
950 CKTdump (ckt, ckt->CKTtime, job->PSSplot_td) ;
951 psstimes [pss_points_cycle] = ckt->CKTtime ;
952 for (i = 1 ; i <= msize ; i++)
953 pssvalues [i - 1 + pss_points_cycle * msize] = ckt->CKTrhsOld [i] ;
954
955 /* Update the PSS points counter and set the next Breakpoint */
956 pss_points_cycle++ ;
957 CKTsetBreak (ckt, time_temp + (1 / ckt->CKTguessedFreq) * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ;
958
959 if (excessive_err_nodes == 0)
960 fprintf (stderr, "\nConvergence reached. Final circuit time is %1.10g seconds (iteration n° %d) and predicted fundamental frequency is %15.10g Hz\n", ckt->CKTtime, shooting_cycle_counter - 1, ckt->CKTguessedFreq) ;
961 else
962 fprintf (stderr, "\nConvergence not reached. However the most near convergence iteration has predicted (iteration %d) a fundamental frequency of %15.10g Hz\n", k, ckt->CKTguessedFreq) ;
963
964 #ifdef STEPDEBUG
965 fprintf (stderr, "time_temp %g\n", time_temp) ;
966 fprintf (stderr, "IN_PSS: FIRST time point accepted in evolution for FFT calculations\n") ;
967 fprintf (stderr, "Circuit time %1.15g, final time %1.15g, point index %d and total requested points %ld\n",
968 ckt->CKTtime, time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints),
969 pss_points_cycle, ckt->CKTpsspoints) ;
970 fprintf (stderr, "Next breakpoint set in: %1.15g\n",
971 time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ;
972 #endif
973
974 } else {
975 /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */
976 ckt->CKTfinalTime = time_temp + 1 / ckt->CKTguessedFreq ;
977
978 /* Set next the breakpoint */
979 CKTsetBreak (ckt, time_temp + 1 / ckt->CKTguessedFreq) ;
980 }
981 }
982 }
983 break;
984
985 case PSS:
986 {
987 /* The algorithm enters here when in_pss is set */
988
989 #ifdef STEPDEBUG
990 fprintf (stderr, "ttemp %1.15g, final_time %1.15g, current_time %1.15g\n", time_temp, time_temp + 1 / ckt->CKTguessedFreq, ckt->CKTtime) ;
991 #endif
992
993 if ((pss_points_cycle == ckt->CKTpsspoints + 1) || (ckt->CKTtime > ckt->CKTfinalTime))
994 {
995 double *pssfreqs = TMALLOC (double, ckt->CKTharms);
996 double *pssmags = TMALLOC (double, ckt->CKTharms);
997 double *pssphases = TMALLOC (double, ckt->CKTharms);
998 double *pssnmags = TMALLOC (double, ckt->CKTharms);
999 double *pssnphases = TMALLOC (double, ckt->CKTharms);
1000 double *pssValues = TMALLOC (double, ckt->CKTpsspoints + 1);
1001 double *pssResults = TMALLOC (double, msize * ckt->CKTharms);
1002
1003 /* End plot in Time Domain */
1004 SPfrontEnd->OUTendPlot (job->PSSplot_td) ;
1005
1006 /* Frequency Plot Creation */
1007 error = CKTnames (ckt, &numNames, &nameList) ;
1008 if (error)
1009 return (error) ;
1010 SPfrontEnd->IFnewUid (ckt, &freqUid, NULL, "frequency", UID_OTHER, NULL) ;
1011 error = SPfrontEnd->OUTpBeginPlot (ckt, ckt->CKTcurJob,
1012 "Frequency Domain Periodic Steady State Analysis",
1013 freqUid, IF_REAL,
1014 numNames, nameList, IF_REAL,
1015 &(job->PSSplot_fd)) ;
1016 tfree (nameList) ;
1017 SPfrontEnd->OUTattributes (job->PSSplot_fd, NULL, PLOT_COMB, NULL) ;
1018
1019 /* ******************** */
1020 /* Starting DFT on data */
1021 /* ******************** */
1022 for (i = 0 ; i < msize ; i++)
1023 {
1024 for (j = 0 ; j < ckt->CKTpsspoints ; j++)
1025 pssValues [j] = pssvalues [j * msize + i] ;
1026
1027 DFT (ckt->CKTpsspoints, ckt->CKTharms, &thd, psstimes, pssValues, ckt->CKTguessedFreq,
1028 pssfreqs, pssmags, pssphases, pssnmags, pssnphases) ;
1029
1030 for (j = 0 ; j < ckt->CKTharms ; j++)
1031 pssResults [j * msize + i] = pssmags [j] ;
1032 }
1033
1034 for (j = 0 ; j < ckt->CKTharms ; j++)
1035 {
1036 for (i = 0 ; i < msize ; i++)
1037 ckt->CKTrhsOld [i + 1] = pssResults [j * msize + i] ;
1038
1039 CKTdump (ckt, pssfreqs [j], job->PSSplot_fd) ;
1040 }
1041 /* ****************** */
1042 /* Ending DFT on data */
1043 /* ****************** */
1044
1045 /* Terminates plot in Frequency Domain and frees the allocated memory */
1046 SPfrontEnd->OUTendPlot (job->PSSplot_fd) ;
1047
1048
1049
1050 /* Francesco Lannutti's MOD */
1051
1052 /* Verify the frequency found */
1053 max_freq = pssResults [msize] ; /* max_freq = pssResults [1 * msize + 0] ; */
1054 position = 1 ;
1055 for (j = 1 ; j < ckt->CKTharms ; j++)
1056 {
1057 for (i = 0 ; i < msize ; i++)
1058 {
1059 if (max_freq < pssResults [j * msize + i])
1060 {
1061 max_freq = pssResults [j * msize + i] ;
1062 position = j ;
1063 }
1064 }
1065 }
1066
1067 if (pssfreqs [position] != ckt->CKTguessedFreq)
1068 {
1069 ckt->CKTguessedFreq = pssfreqs [position] ;
1070 fprintf (stderr, "The predicted fundamental frequency is incorrect.\nRelaunching the analysis...\n\n") ;
1071 fprintf (stderr, "The new guessed fundamental frequency is: %.6g\n\n", ckt->CKTguessedFreq) ;
1072 DCpss (ckt, 1) ;
1073 }
1074 /****************************/
1075
1076
1077 FREE (pssResults) ;
1078 FREE (pssValues) ;
1079 FREE (pssnphases) ;
1080 FREE (pssnmags) ;
1081 FREE (pssphases) ;
1082 FREE (pssmags) ;
1083 FREE (pssfreqs) ;
1084
1085 FREE (RHS_copy_se) ;
1086 FREE (RHS_copy_der) ;
1087 FREE (RHS_max) ;
1088 FREE (RHS_min) ;
1089 FREE (err_conv) ;
1090 FREE (psstimes) ;
1091 FREE (pssvalues) ;
1092 return (OK) ;
1093 }
1094 }
1095 break;
1096
1097 } /* switch(pss_state) */
1098
1099 /* ********************************** */
1100 /* **** END SHOOTING CODE BLOCK ***** */
1101 /* ********************************** */
1102
1103 if(SPfrontEnd->IFpauseTest()) {
1104 /* user requested pause... */
1105 UPDATE_STATS(DOING_TRAN);
1106 return(E_PAUSE);
1107 }
1108 /* RESUME */
1109 resume:
1110 #ifdef STEPDEBUG
1111 if( (ckt->CKTdelta <= ckt->CKTfinalTime/50) &&
1112 (ckt->CKTdelta <= ckt->CKTmaxStep)) {
1113 ;
1114 } else {
1115 if(ckt->CKTfinalTime/50<ckt->CKTmaxStep) {
1116 fprintf (stderr, "limited by Tstop/50\n");
1117 } else {
1118 fprintf (stderr, "limited by Tmax == %g\n", ckt->CKTmaxStep);
1119 }
1120 }
1121 #endif
1122 #ifdef HAS_PROGREP
1123 if (ckt->CKTtime == 0.)
1124 SetAnalyse( "tran init", 0);
1125 else if ((pss_state != PSS) && (shooting_cycle_counter > 0))
1126 SetAnalyse("shooting", shooting_cycle_counter) ;
1127 else
1128 SetAnalyse( "tran", (int)((ckt->CKTtime * 1000.) / ckt->CKTfinalTime));
1129 #endif
1130 ckt->CKTdelta =
1131 MIN(ckt->CKTdelta,ckt->CKTmaxStep);
1132 #ifdef XSPICE
1133 /* gtri - begin - wbk - Cut integration order if first timepoint after breakpoint */
1134 /* if(ckt->CKTtime == g_mif_info.breakpoint.last) */
1135 if ( AlmostEqualUlps( ckt->CKTtime, g_mif_info.breakpoint.last, 100 ) )
1136 ckt->CKTorder = 1;
1137 /* gtri - end - wbk - Cut integration order if first timepoint after breakpoint */
1138
1139 #endif
1140
1141 /* are we at a breakpoint, or indistinguishably close? */
1142 /* if ((ckt->CKTtime == ckt->CKTbreaks[0]) || (ckt->CKTbreaks[0] - */
1143 if (ckt->CKTbreaks [0] - ckt->CKTtime <= ckt->CKTdelmin)
1144 {
1145 /*if ( AlmostEqualUlps( ckt->CKTtime, ckt->CKTbreaks[0], 100 ) || (ckt->CKTbreaks[0] -
1146 * (ckt->CKTtime) <= ckt->CKTdelmin)) {*/
1147 /* first timepoint after a breakpoint - cut integration order */
1148 /* and limit timestep to .1 times minimum of time to next breakpoint,
1149 * and previous timestep
1150 */
1151 ckt->CKTorder = 1;
1152 #ifdef STEPDEBUG
1153 if( (ckt->CKTdelta > .1*ckt->CKTsaveDelta) ||
1154 (ckt->CKTdelta > .1*(ckt->CKTbreaks[1] - ckt->CKTbreaks[0])) ) {
1155 if(ckt->CKTsaveDelta < (ckt->CKTbreaks[1] - ckt->CKTbreaks[0])) {
1156 fprintf (stderr, "limited by pre-breakpoint delta (saveDelta: %1.10g, nxt_breakpt: %1.10g, curr_breakpt: %1.10g and CKTtime: %1.10g\n",
1157 ckt->CKTsaveDelta, ckt->CKTbreaks [1], ckt->CKTbreaks [0], ckt->CKTtime) ;
1158 } else {
1159 fprintf (stderr, "limited by next breakpoint\n") ;
1160 fprintf (stderr, "(saveDelta: %1.10g, Delta: %1.10g, CKTtime: %1.10g and delmin: %1.10g\n",
1161 ckt->CKTsaveDelta, ckt->CKTdelta, ckt->CKTtime, ckt->CKTdelmin) ;
1162 }
1163 }
1164 #endif
1165
1166 if (ckt->CKTbreaks [1] - ckt->CKTbreaks [0] == 0)
1167 ckt->CKTdelta = ckt->CKTdelmin ;
1168 else
1169 ckt->CKTdelta = MIN (ckt->CKTdelta, .1 * MIN (ckt->CKTsaveDelta,
1170 ckt->CKTbreaks[1] - ckt->CKTbreaks[0]));
1171
1172 if(firsttime) {
1173 ckt->CKTdelta /= 10;
1174 #ifdef STEPDEBUG
1175 fprintf(stderr, "delta cut for initial timepoint\n");
1176 #endif
1177 }
1178
1179 #ifndef XSPICE
1180 /* don't want to get below delmin for no reason */
1181 ckt->CKTdelta = MAX(ckt->CKTdelta, ckt->CKTdelmin*2.0);
1182 #endif
1183
1184 }
1185
1186 #ifndef XSPICE
1187 else if(ckt->CKTtime + ckt->CKTdelta >= ckt->CKTbreaks[0]) {
1188 ckt->CKTsaveDelta = ckt->CKTdelta;
1189 ckt->CKTdelta = ckt->CKTbreaks[0] - ckt->CKTtime;
1190 /* fprintf (stderr, "delta cut to %g to hit breakpoint\n" ,ckt->CKTdelta) ; */
1191 fflush(stdout);
1192 ckt->CKTbreak = 1; /* why? the current pt. is not a bkpt. */
1193 }
1194 #endif /* !XSPICE */
1195
1196
1197 #ifdef XSPICE
1198 /* gtri - begin - wbk - Add Breakpoint stuff */
1199
1200 if(ckt->CKTtime + ckt->CKTdelta >= g_mif_info.breakpoint.current) {
1201 /* If next time > temporary breakpoint, force it to the breakpoint */
1202 /* And mark that timestep was set by temporary breakpoint */
1203 ckt->CKTsaveDelta = ckt->CKTdelta;
1204 ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime;
1205 g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta;
1206 } else {
1207 /* Else, mark that timestep was not set by temporary breakpoint */
1208 g_mif_info.breakpoint.last = HUGE_VAL;
1209 }
1210
1211 /* gtri - end - wbk - Add Breakpoint stuff */
1212
1213 /* gtri - begin - wbk - Modify Breakpoint stuff */
1214 /* Throw out any permanent breakpoint times <= current time */
1215 for (;;) {
1216 #ifdef STEPDEBUG
1217 fprintf (stderr, " brk_pt: %g ckt_time: %g ckt_min_break: %g\n", ckt->CKTbreaks [0], ckt->CKTtime, ckt->CKTminBreak) ;
1218 #endif
1219 if(AlmostEqualUlps(ckt->CKTbreaks[0], ckt->CKTtime, 100) ||
1220 ckt->CKTbreaks[0] <= ckt->CKTtime + ckt->CKTminBreak) {
1221 #ifdef STEPDEBUG
1222 fprintf (stderr, "throwing out permanent breakpoint times <= current time (brk pt: %g)\n", ckt->CKTbreaks [0]) ;
1223 fprintf (stderr, "ckt_time: %g ckt_min_break: %g\n", ckt->CKTtime, ckt->CKTminBreak) ;
1224 #endif
1225 CKTclrBreak(ckt);
1226 } else {
1227 break;
1228 }
1229 }
1230 /* Force the breakpoint if appropriate */
1231 if(ckt->CKTtime + ckt->CKTdelta > ckt->CKTbreaks[0]) {
1232 ckt->CKTbreak = 1;
1233 ckt->CKTsaveDelta = ckt->CKTdelta;
1234 ckt->CKTdelta = ckt->CKTbreaks[0] - ckt->CKTtime;
1235 }
1236
1237 /* gtri - end - wbk - Modify Breakpoint stuff */
1238
1239 /* gtri - begin - wbk - Do event solution */
1240
1241 if(ckt->evt->counts.num_insts > 0) {
1242
1243 /* if time = 0 and op_alternate was specified as false during */
1244 /* dcop analysis, call any changed instances to let them */
1245 /* post their outputs with their associated delays */
1246 if((ckt->CKTtime == 0.0) && (! ckt->evt->options.op_alternate))
1247 EVTiter(ckt);
1248
1249 /* while there are events on the queue with event time <= next */
1250 /* projected analog time, process them */
1251 while((g_mif_info.circuit.evt_step = EVTnext_time(ckt))
1252 <= (ckt->CKTtime + ckt->CKTdelta)) {
1253
1254 /* Initialize temp analog bkpt to infinity */
1255 g_mif_info.breakpoint.current = HUGE_VAL;
1256
1257 /* Pull items off queue and process them */
1258 EVTdequeue(ckt, g_mif_info.circuit.evt_step);
1259 EVTiter(ckt);
1260
1261 /* If any instances have forced an earlier */
1262 /* next analog time, cut the delta */
1263 if(ckt->CKTbreaks[0] < g_mif_info.breakpoint.current)
1264 if(ckt->CKTbreaks[0] > ckt->CKTtime + ckt->CKTminBreak)
1265 g_mif_info.breakpoint.current = ckt->CKTbreaks[0];
1266 if(g_mif_info.breakpoint.current < ckt->CKTtime + ckt->CKTdelta) {
1267 /* Breakpoint must be > last accepted timepoint */
1268 /* and >= current event time */
1269 if(g_mif_info.breakpoint.current > ckt->CKTtime + ckt->CKTminBreak &&
1270 g_mif_info.breakpoint.current >= g_mif_info.circuit.evt_step) {
1271 ckt->CKTsaveDelta = ckt->CKTdelta;
1272 ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime;
1273 g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta;
1274 }
1275 }
1276
1277 } /* end while next event time <= next analog time */
1278 } /* end if there are event instances */
1279
1280 /* gtri - end - wbk - Do event solution */
1281 #else
1282
1283 #ifdef CLUSTER
1284 if(!CLUsync(ckt->CKTtime,&ckt->CKTdelta,0)) {
1285 fprintf (stderr, "Sync error!\n");
1286 exit(0);
1287 }
1288 #endif /* CLUSTER */
1289
1290 #endif
1291
1292 /* What is that??? */
1293 for(i=5; i>=0; i--)
1294 ckt->CKTdeltaOld[i+1] = ckt->CKTdeltaOld[i];
1295 ckt->CKTdeltaOld[0] = ckt->CKTdelta;
1296
1297 temp = ckt->CKTstates[ckt->CKTmaxOrder+1];
1298 for(i=ckt->CKTmaxOrder;i>=0;i--) {
1299 ckt->CKTstates[i+1] = ckt->CKTstates[i];
1300 }
1301 ckt->CKTstates[0] = temp;
1302
1303 /* 600 */
1304 for (;;) {
1305 #ifdef CLUSTER
1306 redostep = 1;
1307 #endif
1308 #ifdef XSPICE
1309 /* gtri - add - wbk - 4/17/91 - Fix Berkeley bug */
1310 /* This is needed here to allow CAPask to output currents */
1311 /* during Transient analysis. A grep for CKTcurrentAnalysis */
1312 /* indicates that it should not hurt anything else ... */
1313
1314 ckt->CKTcurrentAnalysis = DOING_TRAN;
1315
1316 /* gtri - end - wbk - 4/17/91 - Fix Berkeley bug */
1317 #endif
1318 olddelta=ckt->CKTdelta;
1319 /* time abort? */
1320 ckt->CKTtime += ckt->CKTdelta;
1321 #ifdef CLUSTER
1322 CLUinput(ckt);
1323 #endif
1324 ckt->CKTdeltaOld[0]=ckt->CKTdelta;
1325 NIcomCof(ckt);
1326 #ifdef PREDICTOR
1327 error = NIpred(ckt);
1328 #endif /* PREDICTOR */
1329 save_mode = ckt->CKTmode;
1330 save_order = ckt->CKTorder;
1331 #ifdef XSPICE
1332 /* gtri - begin - wbk - Add Breakpoint stuff */
1333
1334 /* Initialize temporary breakpoint to infinity */
1335 g_mif_info.breakpoint.current = HUGE_VAL;
1336
1337 /* gtri - end - wbk - Add Breakpoint stuff */
1338
1339
1340 /* gtri - begin - wbk - add convergence problem reporting flags */
1341 /* delta is forced to equal delmin on last attempt near line 650 */
1342 if(ckt->CKTdelta <= ckt->CKTdelmin)
1343 ckt->enh->conv_debug.last_NIiter_call = MIF_TRUE;
1344 else
1345 ckt->enh->conv_debug.last_NIiter_call = MIF_FALSE;
1346 /* gtri - begin - wbk - add convergence problem reporting flags */
1347
1348
1349 /* gtri - begin - wbk - Call all hybrids */
1350
1351 /* gtri - begin - wbk - Set evt_step */
1352
1353 if(ckt->evt->counts.num_insts > 0) {
1354 g_mif_info.circuit.evt_step = ckt->CKTtime;
1355 }
1356 /* gtri - end - wbk - Set evt_step */
1357 #endif
1358
1359 converged = NIiter(ckt,ckt->CKTtranMaxIter);
1360
1361 #ifdef XSPICE
1362 if(ckt->evt->counts.num_insts > 0) {
1363 g_mif_info.circuit.evt_step = ckt->CKTtime;
1364 EVTcall_hybrids(ckt);
1365 }
1366 /* gtri - end - wbk - Call all hybrids */
1367
1368 #endif
1369 ckt->CKTstat->STATtimePts ++;
1370 ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITPRED;
1371 if(firsttime) {
1372 memcpy(ckt->CKTstate2, ckt->CKTstate1,
1373 (size_t) ckt->CKTnumStates * sizeof(double));
1374 memcpy(ckt->CKTstate3, ckt->CKTstate1,
1375 (size_t) ckt->CKTnumStates * sizeof(double));
1376 }
1377 /* txl, cpl addition */
1378 if (converged == 1111) {
1379 return(converged);
1380 }
1381
1382 #ifdef STEPDEBUG
1383 if (pss_state == PSS)
1384 fprintf (stderr, "pss_state: %d, converged: %d\n", pss_state, converged) ;
1385 #endif
1386 if(converged != 0) {
1387 #ifndef CLUSTER
1388 ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta;
1389 ckt->CKTstat->STATrejected ++;
1390 #endif
1391 ckt->CKTdelta = ckt->CKTdelta/8;
1392 #ifdef STEPDEBUG
1393 fprintf (stderr, "delta cut to %g for non-convergence\n", ckt->CKTdelta) ;
1394 fflush(stdout);
1395 #endif
1396 if(firsttime) {
1397 ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODETRAN | MODEINITTRAN;
1398 }
1399 ckt->CKTorder = 1;
1400
1401 #ifdef XSPICE
1402 /* gtri - begin - wbk - Add Breakpoint stuff */
1403
1404 /* Force backup if temporary breakpoint is < current time */
1405 } else if(g_mif_info.breakpoint.current < ckt->CKTtime) {
1406 ckt->CKTsaveDelta = ckt->CKTdelta;
1407 ckt->CKTtime -= ckt->CKTdelta;
1408 ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime;
1409 g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta;
1410
1411 if(firsttime) {
1412 ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITTRAN;
1413 }
1414 ckt->CKTorder = 1;
1415
1416 /* gtri - end - wbk - Add Breakpoint stuff */
1417 #endif
1418
1419 } else {
1420 if (firsttime) {
1421 firsttime = 0;
1422 #ifndef CLUSTER
1423 goto nextTime; /* no check on
1424 * first time point
1425 */
1426 #else
1427 redostep = 0;
1428 goto chkStep;
1429 #endif
1430 }
1431 newdelta = ckt->CKTdelta;
1432 error = CKTtrunc(ckt,&newdelta);
1433 if(error) {
1434 UPDATE_STATS(DOING_TRAN);
1435 return(error);
1436 }
1437 if (newdelta > .9 * ckt->CKTdelta) {
1438 if ((ckt->CKTorder == 1) && (ckt->CKTmaxOrder > 1)) { /* don't rise the order for backward Euler */
1439 newdelta = ckt->CKTdelta;
1440 ckt->CKTorder = 2;
1441 error = CKTtrunc(ckt, &newdelta);
1442 if (error) {
1443 UPDATE_STATS(DOING_TRAN);
1444 return(error);
1445 }
1446 if (newdelta <= 1.05 * ckt->CKTdelta) {
1447 ckt->CKTorder = 1;
1448 }
1449 }
1450 /* time point OK - 630 */
1451 ckt->CKTdelta = newdelta;
1452 #ifdef NDEV
1453 if (!ft_norefprint) {
1454 /* show a time process indicator, by Gong Ding, gdiso@ustc.edu */
1455 if (ckt->CKTtime / ckt->CKTfinalTime * 100 < 10.0)
1456 fprintf(stderr, "%%%3.2lf\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100);
1457 else if (ckt->CKTtime / ckt->CKTfinalTime * 100 < 100.0)
1458 fprintf(stderr, "%%%4.2lf\b\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100);
1459 else
1460 fprintf(stderr, "%%%5.2lf\b\b\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100);
1461 fflush(stdout);
1462 }
1463 #endif
1464
1465 #ifdef STEPDEBUG
1466 fprintf (stderr, "delta set to truncation error result: %g. Point accepted at CKTtime: %g\n", ckt->CKTdelta, ckt->CKTtime) ;
1467 fflush(stdout);
1468 #endif
1469
1470 #ifndef CLUSTER
1471 /* go to 650 - trapezoidal */
1472 goto nextTime;
1473 #else
1474 redostep = 0;
1475 goto chkStep;
1476 #endif
1477 } else { /* newdelta <= .9 * ckt->CKTdelta */
1478 #ifndef CLUSTER
1479 ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta;
1480 ckt->CKTstat->STATrejected ++;
1481 #endif
1482 ckt->CKTdelta = newdelta;
1483 #ifdef STEPDEBUG
1484 fprintf (stderr, "delta set to truncation error result:point rejected\n") ;
1485 #endif
1486 }
1487 }
1488
1489 if (ckt->CKTdelta <= ckt->CKTdelmin) {
1490 if (olddelta > ckt->CKTdelmin) {
1491 ckt->CKTdelta = ckt->CKTdelmin;
1492 #ifdef STEPDEBUG
1493 fprintf (stderr, "delta at delmin\n");
1494 #endif
1495 } else {
1496 UPDATE_STATS(DOING_TRAN);
1497 errMsg = CKTtrouble(ckt, "Timestep too small");
1498 return(E_TIMESTEP);
1499 }
1500 }
1501 #ifdef XSPICE
1502 /* gtri - begin - wbk - Do event backup */
1503
1504 if(ckt->evt->counts.num_insts > 0)
1505 EVTbackup(ckt, ckt->CKTtime + ckt->CKTdelta);
1506
1507 /* gtri - end - wbk - Do event backup */
1508 #endif
1509 #ifdef CLUSTER
1510 chkStep:
1511 if(CLUsync(ckt->CKTtime,&ckt->CKTdelta,redostep)){
1512 goto nextTime;
1513 } else {
1514 ckt->CKTtime -= olddelta;
1515 ckt->CKTstat->STATrejected ++;
1516 }
1517 #endif
1518 }
1519 /* NOTREACHED */
1520 }
1521
1522 static int
DFT(long int ndata,int numFreq,double * thd,double * Time,double * Value,double FundFreq,double * Freq,double * Mag,double * Phase,double * nMag,double * nPhase)1523 DFT
1524 (
1525 long int ndata, /* number of entries in the Time and Value arrays */
1526 int numFreq, /* number of harmonics to calculate */
1527 double *thd, /* total harmonic distortion (percent) to be returned */
1528 double *Time, /* times at which the voltage/current values were measured */
1529 double *Value, /* voltage or current vector whose transform is desired */
1530 double FundFreq, /* the fundamental frequency of the analysis */
1531 double *Freq, /* the frequency value of the various harmonics */
1532 double *Mag, /* the Magnitude of the fourier transform */
1533 double *Phase, /* the Phase of the fourier transform */
1534 double *nMag, /* the normalized magnitude of the transform: nMag (fund) = 1 */
1535 double *nPhase /* the normalized phase of the transform: Nphase (fund) = 0 */
1536 )
1537 {
1538 /* Note: we can consider these as a set of arrays. The sizes are:
1539 * Time [ndata], Value [ndata], Freq [numFreq], Mag [numfreq],
1540 * Phase [numfreq], nMag [numfreq], nPhase [numfreq]
1541 *
1542 * The arrays must all be allocated by the caller.
1543 * The Time and Value array must be reasonably distributed over at
1544 * least one full period of the fundamental Frequency for the
1545 * fourier transform to be useful. The function will take the
1546 * last period of the frequency as data for the transform.
1547 *
1548 * We are assuming that the caller has provided exactly one period
1549 * of the fundamental frequency. */
1550 int i, j;
1551 double tmp;
1552
1553 NG_IGNORE (Time);
1554
1555 /* clear output/computation arrays */
1556
1557 for (i = 0; i < numFreq; i++) {
1558 Mag [i] = 0;
1559 Phase [i] = 0;
1560 }
1561
1562 for (i = 0; i < ndata; i++) {
1563 for (j = 0; j < numFreq; j++) {
1564 Mag [j] += (Value [i] * sin (j * 2.0 * M_PI * i / ((double)ndata)));
1565 Phase [j] += (Value [i] * cos (j * 2.0 * M_PI * i / ((double)ndata)));
1566 }
1567 }
1568
1569 Mag [0] = Phase [0] / (double)ndata;
1570 Phase [0] = 0;
1571 nMag [0] = 0;
1572 nPhase [0] = 0;
1573 Freq [0] = 0;
1574 *thd = 0;
1575
1576 for (i = 1; i < numFreq; i++) {
1577 tmp = Mag [i] * 2.0 / (double)ndata;
1578 Phase [i] *= 2.0 / (double)ndata;
1579 Freq [i] = i * FundFreq;
1580 Mag [i] = hypot (tmp, Phase [i]);
1581 Phase [i] = atan2 (Phase [i], tmp) * 180.0 / M_PI;
1582 nMag [i] = Mag [i] / Mag [1];
1583 nPhase [i] = Phase [i] - Phase [1];
1584 if (i > 1)
1585 *thd += nMag [i] * nMag [i];
1586 }
1587
1588 *thd = 100 * sqrt (*thd);
1589 return (OK);
1590 }
1591