1 /*$Id: s_tr_swp.cc,v 26.137 2010/04/10 02:37:05 al Exp $ -*- C++ -*-
2 * Copyright (C) 2001 Albert Davis
3 * Author: Albert Davis <aldavis@gnu.org>
4 *
5 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 3, or (at your option)
10 * any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 * 02110-1301, USA.
21 *------------------------------------------------------------------
22 * sweep time and simulate. output results.
23 * manage event queue
24 */
25 //testing=script 2007.11.22
26 #include "u_time_pair.h"
27 #include "u_sim_data.h"
28 #include "u_status.h"
29 #include "declare.h" /* gen */
30 #include "s_tr.h"
31 /*--------------------------------------------------------------------------*/
32 // void TRANSIENT::sweep(void);
33 // void TRANSIENT::first(void);
34 // bool TRANSIENT::next(void);
35 // void TRANSIENT::accept(void);
36 // void TRANSIENT::reject(void);
37 /*--------------------------------------------------------------------------*/
38 /*--------------------------------------------------------------------------*/
39 namespace TR {
40 static std::string step_cause[] = {
41 "impossible",
42 "user requested",
43 "event queue",
44 "command line \"skip\"",
45 "convergence failure, reducing (itl4)",
46 "slow convergence, holding (itl3)",
47 "truncation error",
48 "ambiguous event",
49 "limit growth",
50 "initial guess"
51 };
52 }
53 /*--------------------------------------------------------------------------*/
sweep()54 void TRANSIENT::sweep()
55 {
56 _sim->_phase = p_INIT_DC;
57 head(_tstart, _tstop, "Time");
58 _sim->_bypass_ok = false;
59 _sim->set_inc_mode_bad();
60
61 if (_cont) { // use the data from last time
62 _sim->_phase = p_RESTORE;
63 _sim->restore_voltages();
64 CARD_LIST::card_list.tr_restore();
65 }else{
66 _sim->clear_limit();
67 CARD_LIST::card_list.tr_begin();
68 }
69
70 first();
71 _sim->_genout = gen();
72
73 if (_sim->uic_now()) {
74 advance_time();
75 _sim->zero_voltages();
76 CARD_LIST::card_list.do_tr(); //evaluate_models
77 while (!_sim->_late_evalq.empty()) {itested(); //BUG// encapsulation violation
78 _sim->_late_evalq.front()->do_tr_last();
79 _sim->_late_evalq.pop_front();
80 }
81 _converged = true;
82 _sim->_loadq.clear(); // fake solve, clear the queue
83 //BUG// UIC needs further analysis.
84 }else{
85 _converged = solve_with_homotopy(OPT::DCBIAS,_trace);
86 if (!_converged) {
87 error(bWARNING, "did not converge\n");
88 }else{
89 }
90 }
91 review();
92 _accepted = true;
93 accept();
94
95 {
96 bool printnow = (_sim->_time0 == _tstart || _trace >= tALLTIME);
97 if (printnow) {
98 _sim->keep_voltages();
99 outdata(_sim->_time0);
100 }else{untested();
101 }
102 }
103
104 while (next()) {
105 _sim->_bypass_ok = false;
106 _sim->_phase = p_TRAN;
107 _sim->_genout = gen();
108 _converged = solve(OPT::TRHIGH,_trace);
109
110 _accepted = _converged && review();
111
112 if (_accepted) {
113 assert(_converged);
114 assert(_sim->_time0 <= _time_by_user_request);
115 accept();
116 if (step_cause() == scUSER) {
117 assert(up_order(_sim->_time0-_sim->_dtmin, _time_by_user_request, _sim->_time0+_sim->_dtmin));
118 ++_stepno;
119 _time_by_user_request += _tstep; /* advance user time */
120 }else{
121 }
122 assert(_sim->_time0 < _time_by_user_request);
123 }else{
124 reject();
125 assert(time1 < _time_by_user_request);
126 }
127 {
128 bool printnow =
129 (_trace >= tREJECTED)
130 || (_accepted && ((_trace >= tALLTIME)
131 || (step_cause() == scUSER && _sim->_time0+_sim->_dtmin > _tstart)));
132 if (printnow) {
133 _sim->keep_voltages();
134 outdata(_sim->_time0);
135 }else{
136 }
137 }
138
139 if (!_converged && OPT::quitconvfail) {untested();
140 outdata(_sim->_time0);
141 throw Exception("convergence failure, giving up");
142 }else{
143 }
144 }
145 }
146 /*--------------------------------------------------------------------------*/
set_step_cause(STEP_CAUSE C)147 void TRANSIENT::set_step_cause(STEP_CAUSE C)
148 {
149 switch (C) {
150 case scITER_A:untested();
151 case scADT:untested();
152 case scITER_R:
153 case scINITIAL:
154 case scSKIP:
155 case scTE:
156 case scAMBEVENT:
157 case scEVENTQ:
158 case scUSER:
159 ::status.control = C;
160 break;
161 case scNO_ADVANCE:untested();
162 case scZERO:untested();
163 case scSMALL:itested();
164 case scREJECT:
165 ::status.control += C;
166 break;
167 }
168 }
169 /*--------------------------------------------------------------------------*/
step_cause() const170 int TRANSIENT::step_cause()const
171 {
172 return ::status.control;
173 }
174 /*--------------------------------------------------------------------------*/
first()175 void TRANSIENT::first()
176 {
177 /* usually, _sim->_time0, time1 == 0, from setup */
178 assert(_sim->_time0 == time1);
179 assert(_sim->_time0 <= _tstart);
180 ::status.review.start();
181 _time_by_user_request = _sim->_time0 + _tstep; /* set next user step */
182 //_eq.Clear(); /* empty the queue */
183 while (!_sim->_eq.empty()) {
184 _sim->_eq.pop();
185 }
186 _stepno = 0;
187 set_step_cause(scUSER);
188 ++::status.hidden_steps;
189 ::status.review.stop();
190 }
191 /*--------------------------------------------------------------------------*/
192 #define check_consistency() { \
193 trace4("", __LINE__, newtime, almost_fixed_time, fixed_time); \
194 assert(almost_fixed_time <= fixed_time); \
195 assert(newtime <= fixed_time); \
196 /*assert(newtime == fixed_time || newtime <= fixed_time -_sim->_dtmin);*/ \
197 assert(newtime <= almost_fixed_time); \
198 /*assert(newtime == almost_fixed_time || newtime <= almost_fixed_time - _sim->_dtmin);*/ \
199 assert(newtime > time1); \
200 assert(newtime > reftime); \
201 assert(new_dt > 0.); \
202 assert(new_dt >= _sim->_dtmin); \
203 assert(newtime <= _time_by_user_request); \
204 /*assert(newtime == _time_by_user_request*/ \
205 /* || newtime < _time_by_user_request - _sim->_dtmin); */ \
206 }
207 #define check_consistency2() { \
208 assert(newtime > time1); \
209 assert(new_dt > 0.); \
210 assert(new_dt >= _sim->_dtmin); \
211 assert(newtime <= _time_by_user_request); \
212 /*assert(newtime == _time_by_user_request */ \
213 /* || newtime < _time_by_user_request - _sim->_dtmin);*/ \
214 }
215 /*--------------------------------------------------------------------------*/
216 /* next: go to next time step
217 * Set _sim->_time0 to the next time step, store the old one in time1.
218 * Try several methods. Take the one that gives the shortest step.
219 */
next()220 bool TRANSIENT::next()
221 {
222 ::status.review.start();
223
224 double old_dt = _sim->_time0 - time1;
225 assert(old_dt >= 0);
226
227 double newtime = NEVER;
228 double new_dt = NEVER;
229 STEP_CAUSE new_control = scNO_ADVANCE;
230
231 if (_sim->_time0 == time1) {
232 // initial step -- could be either t==0 or continue
233 // for the first time, just guess
234 // make it 100x smaller than expected
235 new_dt = std::max(_dtmax/100., _sim->_dtmin);
236 newtime = _sim->_time0 + new_dt;
237 new_control = scINITIAL;
238 }else if (!_converged) {
239 new_dt = old_dt / OPT::trstepshrink;
240 newtime = _time_by_iteration_count = time1 + new_dt;
241 new_control = scITER_R;
242 }else{
243 double reftime;
244 if (_accepted) {
245 reftime = _sim->_time0;
246 trace0("accepted");
247 }else{
248 reftime = time1;
249 trace0("rejected");
250 }
251 trace2("", step_cause(), old_dt);
252 trace3("", time1, _sim->_time0, reftime);
253
254 newtime = _time_by_user_request;
255 new_dt = newtime - reftime;
256 new_control = scUSER;
257 double fixed_time = newtime;
258 double almost_fixed_time = newtime;
259 check_consistency();
260
261 // event queue, events that absolutely will happen
262 // exact time. NOT ok to move or omit, even by _sim->_dtmin
263 // some action is associated with it.
264 if (!_sim->_eq.empty() && _sim->_eq.top() < newtime) {
265 newtime = _sim->_eq.top();
266 new_dt = newtime - reftime;
267 if (new_dt < _sim->_dtmin) {
268 //new_dt = _sim->_dtmin;
269 //newtime = reftime + new_dt;
270 }else{
271 }
272 new_control = scEVENTQ;
273 fixed_time = newtime;
274 almost_fixed_time = newtime;
275 check_consistency();
276 }else{
277 }
278
279 // device events that may not happen
280 // not sure of exact time. will be rescheduled if wrong.
281 // ok to move by _sim->_dtmin. time is not that accurate anyway.
282 if (_time_by_ambiguous_event < newtime - _sim->_dtmin) {
283 if (_time_by_ambiguous_event < time1 + 2*_sim->_dtmin) {untested();
284 double mintime = time1 + 2*_sim->_dtmin;
285 if (newtime - _sim->_dtmin < mintime) {untested();
286 newtime = mintime;
287 new_control = scAMBEVENT;
288 }else{untested();
289 }
290 }else{
291 newtime = _time_by_ambiguous_event;
292 new_control = scAMBEVENT;
293 }
294 new_dt = newtime - reftime;
295 almost_fixed_time = newtime;
296 check_consistency();
297 }else{
298 }
299
300 // device error estimates
301 if (_time_by_error_estimate < newtime - _sim->_dtmin) {
302 newtime = _time_by_error_estimate;
303 new_dt = newtime - reftime;
304 new_control = scTE;
305 check_consistency();
306 }else{
307 }
308
309 // skip parameter
310 if (new_dt > _dtmax) {
311 if (new_dt > _dtmax + _sim->_dtmin) {
312 new_control = scSKIP;
313 }else{
314 }
315 new_dt = _dtmax;
316 newtime = reftime + new_dt;
317 check_consistency();
318 }else{
319 }
320
321 // converged but with more iterations than we like
322 if ((new_dt > (old_dt + _sim->_dtmin) * OPT::trstephold)
323 && _sim->exceeds_iteration_limit(OPT::TRLOW)) {untested();
324 assert(_accepted);
325 new_dt = old_dt * OPT::trstephold;
326 newtime = reftime + new_dt;
327 new_control = scITER_A;
328 check_consistency();
329 }else{
330 }
331
332 // limit growth
333 if (_sim->analysis_is_tran_dynamic() && new_dt > old_dt * OPT::trstepgrow) {untested();
334 new_dt = old_dt * OPT::trstepgrow;
335 newtime = reftime + new_dt;
336 new_control = scADT;
337 check_consistency();
338 }else{
339 }
340
341 // quantize
342 if (newtime < almost_fixed_time) {
343 assert(new_dt >= 0);
344 if (newtime > reftime + old_dt*.8
345 && newtime < reftime + old_dt*1.5
346 && reftime + old_dt <= almost_fixed_time) {
347 // new_dt is close enough to old_dt.
348 // use old_dt, to avoid a step change.
349 new_dt = old_dt;
350 newtime = reftime + new_dt;
351 if (newtime > almost_fixed_time) {untested();
352 new_control = scAMBEVENT;
353 newtime = almost_fixed_time;
354 new_dt = newtime - reftime;
355 }else{
356 }
357 check_consistency();
358 }else{
359 // There will be a step change.
360 // Try to choose one that we will keep for a while.
361 // Choose new_dt to be in integer fraction of target_dt.
362 double target_dt = fixed_time - reftime;
363 assert(target_dt >= new_dt);
364 double steps = 1 + floor((target_dt - _sim->_dtmin) / new_dt);
365 assert(steps > 0);
366 new_dt = target_dt / steps;
367 newtime = reftime + new_dt;
368 check_consistency();
369 }
370 }else{
371 assert(newtime == almost_fixed_time);
372 }
373
374 // trap time step too small
375 if (!_accepted && new_dt < _sim->_dtmin) {untested();
376 new_dt = _sim->_dtmin;
377 newtime = reftime + new_dt;
378 new_control = scSMALL;
379 check_consistency();
380 }else{
381 }
382
383 // if all that makes it close to user_requested, make it official
384 if (up_order(newtime-_sim->_dtmin, _time_by_user_request, newtime+_sim->_dtmin)) {
385 //newtime = _time_by_user_request;
386 //new_dt = newtime - reftime;
387 new_control = scUSER;
388 check_consistency();
389 }else{
390 }
391 check_consistency();
392 assert(!_accepted || newtime > _sim->_time0);
393 assert(_accepted || newtime <= _sim->_time0);
394 }
395
396 set_step_cause(new_control);
397
398 /* got it, I think */
399
400 /* check to be sure */
401 if (newtime < time1 + _sim->_dtmin) {itested();
402 /* It's really bad. */
403 /* Reject the most recent step, back up as much as possible, */
404 /* and creep along */
405 assert(!_accepted);
406 assert(step_cause() < scREJECT);
407 assert(step_cause() >= 0);
408 error(bDANGER,"non-recoverable " + TR::step_cause[step_cause()] + "\n");
409 error(bDANGER, "newtime=%e rejectedtime=%e oldtime=%e using=%e\n",
410 newtime, _sim->_time0, time1, time1 + _sim->_dtmin);
411 newtime = time1 + _sim->_dtmin;
412 set_step_cause(scSMALL);
413 //check_consistency2();
414 throw Exception("tried everything, still doesn't work, giving up");
415 //}else if (newtime <= _sim->_time0 - _sim->_dtmin) {
416 }else if (newtime < _sim->_time0) {
417 /* Reject the most recent step. */
418 /* We have faith that it will work with a smaller time step. */
419 assert(!_accepted);
420 assert(newtime >= time1 + _sim->_dtmin);
421 error(bLOG, "backwards time step\n");
422 error(bLOG, "newtime=%e rejectedtime=%e oldtime=%e\n", newtime, _sim->_time0, time1);
423 set_step_cause(scREJECT);
424 _sim->mark_inc_mode_bad();
425 check_consistency2();
426 }else if (newtime < _sim->_time0 + _sim->_dtmin) {untested();
427 /* Another evaluation at the same time. */
428 /* Keep the most recent step, but creep along. */
429 assert(newtime > _sim->_time0 - _sim->_dtmin);
430 error(bDANGER, "zero time step\n");
431 error(bDANGER, "newtime=%e rejectedtime=%e oldtime=%e\n", newtime, _sim->_time0, time1);
432 if (_accepted) {untested();
433 time1 = _sim->_time0;
434 }else{untested();
435 assert(_converged);
436 }
437 check_consistency2();
438 newtime = _sim->_time0 + _sim->_dtmin;
439 if (newtime > _time_by_user_request) {untested();
440 newtime = _time_by_user_request;
441 set_step_cause(scUSER);
442 }else{untested();
443 }
444 set_step_cause(scZERO);
445 check_consistency2();
446 }else{
447 assert(_accepted);
448 assert(newtime >= _sim->_time0 + _sim->_dtmin);
449 /* All is OK. Moving on. */
450 /* Keep value of newtime */
451 time1 = _sim->_time0;
452 check_consistency2();
453 }
454 _sim->_time0 = newtime;
455
456 /* advance event queue (maybe) */
457 /* We already looked at it. Dump what's on top if we took it. */
458 while (!_sim->_eq.empty() && _sim->_eq.top() <= _sim->_time0) {
459 trace1("eq", _sim->_eq.top());
460 _sim->_eq.pop();
461 }
462 while (!_sim->_eq.empty() && _sim->_eq.top() < _sim->_time0 + _sim->_dtmin) {itested();
463 trace1("eq-extra", _sim->_eq.top());
464 _sim->_eq.pop();
465 }
466 //BUG// what if it is later rejected? It's lost!
467
468 check_consistency2();
469 ++::status.hidden_steps;
470 ++steps_total_;
471 ::status.review.stop();
472 trace0("next");
473 return (_sim->_time0 <= _tstop + _sim->_dtmin);
474 }
475 /*--------------------------------------------------------------------------*/
review()476 bool TRANSIENT::review()
477 {
478 ::status.review.start();
479 _sim->count_iterations(iTOTAL);
480
481 TIME_PAIR time_by = CARD_LIST::card_list.tr_review();
482 _time_by_error_estimate = time_by._error_estimate;
483
484 // limit minimum time step
485 // 2*_sim->_dtmin because _time[1] + _sim->_dtmin might be == _time[0].
486 if (time_by._event < time1 + 2*_sim->_dtmin) {
487 _time_by_ambiguous_event = time1 + 2*_sim->_dtmin;
488 }else{
489 _time_by_ambiguous_event = time_by._event;
490 }
491 // force advance when time too close to previous
492 if (std::abs(_time_by_ambiguous_event - _sim->_time0) < 2*_sim->_dtmin) {
493 _time_by_ambiguous_event = _sim->_time0 + 2*_sim->_dtmin;
494 }else{
495 }
496
497 if (time_by._error_estimate < time1 + 2*_sim->_dtmin) {
498 _time_by_error_estimate = time1 + 2*_sim->_dtmin;
499 }else{
500 _time_by_error_estimate = time_by._error_estimate;
501 }
502 if (std::abs(_time_by_error_estimate - _sim->_time0) < 1.1*_sim->_dtmin) {
503 _time_by_error_estimate = _sim->_time0 + 1.1*_sim->_dtmin;
504 }else{
505 }
506
507 ::status.review.stop();
508
509 return (_time_by_error_estimate > _sim->_time0 && _time_by_ambiguous_event > _sim->_time0);
510 }
511 /*--------------------------------------------------------------------------*/
accept()512 void TRANSIENT::accept()
513 {
514 ::status.accept.start();
515 _sim->set_limit();
516 if (OPT::traceload) {
517 while (!_sim->_acceptq.empty()) {
518 _sim->_acceptq.back()->tr_accept();
519 _sim->_acceptq.pop_back();
520 }
521 }else{itested();
522 _sim->_acceptq.clear();
523 CARD_LIST::card_list.tr_accept();
524 }
525 ++steps_accepted_;
526 ::status.accept.stop();
527 }
528 /*--------------------------------------------------------------------------*/
reject()529 void TRANSIENT::reject()
530 {
531 ::status.accept.start();
532 _sim->_acceptq.clear();
533 ++steps_rejected_;
534 ::status.accept.stop();
535 }
536 /*--------------------------------------------------------------------------*/
537 /*--------------------------------------------------------------------------*/
538