src/s_tr_swp.cc

   1 /*$Id: s_tr_swp.cc 2016/09/20 al $ -*- 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 #include "u_time_pair.h"
  26 #include "u_sim_data.h"
  27 #include "u_status.h"
  28 #include "e_card.h" // debugging...
  29 #include "declare.h"    /* gen */
  30 #include "s_tr.h"
  31 #include "e_adp.h" // hack. see below
  32 using namespace std;
  33 //#define ALT_CQ // alternatively clear queue (experimental)
  34 /*--------------------------------------------------------------------------*/
  35 //      void    TRANSIENT::sweep(void);
  36 //      void    TRANSIENT::first(void);
  37 //      bool    TRANSIENT::next(void);
  38 //      void    TRANSIENT::accept(void);
  39 //      void    TRANSIENT::reject(void);
  40 /*--------------------------------------------------------------------------*/
  41 /*--------------------------------------------------------------------------*/
  42 namespace TR {
  43   static std::string step_cause[] = {
  44     "impossible",
  45     "user requested",
  46     "event queue",
  47     "command line \"skip\"",
  48     "convergence failure, reducing (itl4)",
  49     "slow convergence, holding (itl3)",
  50     "truncation error",
  51     "ambiguous event",
  52     "limit growth",
  53     "initial guess"
  54   };
  55 }
  56 /*--------------------------------------------------------------------------*/
  57 void TRANSIENT::sweep()
  58 {
  59   trace2("TRANSIENT::sweep", _cont, _cont_dc);
  60   _sim->_phase = p_INIT_DC;
  61   head(_tstart, _tstop, "Time");
  62   _sim->_bypass_ok = false;
  63   _sim->set_inc_mode_bad();
  64
  65   if ( _print_only ) {
  66     _sim->_phase = p_RESTORE;
  67     _sim->restore_voltages();
  68     CARD_LIST::card_list.tr_restore();
  69       outdata(_sim->_time0);
  70       return;
  71
  72   } else if (_inside_tt && _cont_dc) {
  73     _sim->restore_voltages(); // required by some tr_begins.
  74                               // might be possible to move
  75                               // node dependency into tr_accept
  76     _sim->_cont = true; // keep coil from overwriting current
  77     _sim->_phase = p_RESTORE; // short cut differentiate
  78     CARD_LIST::card_list.tr_restore();
  79     CARD_LIST::card_list.do_tr();
  80     _sim->_cont = false;
  81     _sim->clear_limit();
  82     advance_time();
  83     _sim->reset_iteration_counter(iSTEP);
  84     CARD_LIST::card_list.do_tr();
  85   } else if (_cont_dc) {
  86     // continue from DC point.
  87     _sim->restore_voltages(); // required by some tr_begins.
  88                               // might be possible to move
  89                               // node dependency into tr_accept
  90     if(_inside_tt){ untested();
  91     }else if(_edge_detect & edYES){
  92       _sim->keep_voltages(true);
  93     }
  94     _sim->_cont = true; // keep coil from overwriting current
  95     _sim->_phase = p_RESTORE; // short cut differentiate
  96     CARD_LIST::card_list.tr_begin();
  97     _sim->_cont = false;
  98     _sim->clear_limit();
  99     advance_time();
 100     _sim->reset_iteration_counter(iSTEP);
 101     CARD_LIST::card_list.do_tr();
 102
 103   }else if (_inside_tt) {
 104     trace0("TRANSIENT::sweep inside tt");
 105     assert(  _sim->_mode == s_TTT );
 106
 107     _sim->restore_voltages();
 108     _sim->_phase = p_RESTORE;
 109 //     CARD_LIST::card_list.tr_restore(); no. tt_advance takes care of it.
 110     CARD_LIST::card_list.do_tr();
 111
 112   }else if (_cont) {
 113     // use the data from last time
 114     _sim->_phase = p_RESTORE;
 115     _sim->restore_voltages();
 116     CARD_LIST::card_list.tr_restore();
 117     if(_edge_detect & edYES){
 118       _sim->keep_voltages(true);
 119     }else{
 120     }
 121   }else{
 122     _sim->clear_limit();
 123     CARD_LIST::card_list.tr_begin();
 124   }
 125
 126   first();
 127   _sim->_genout = gen();
 128
 129   // assert (_sim->_loadq.empty());
 130   if ( /*_sim->more_uic_now() */ 0 ) {
 131     trace0("TRAN: more UIC now ");
 132     advance_time();
 133     _sim->zero_voltages();
 134     CARD_LIST::card_list.do_tr();
 135     while (!_sim->_late_evalq.empty()) {itested(); //BUG// encapsulation violation
 136       _sim->_late_evalq.front()->do_tr_last();
 137       _sim->_late_evalq.pop_front();
 138     }
 139 //    _converged = true;
 140     _converged = solve_with_homotopy(OPT::DCBIAS,_trace);
 141
 142
 143   } else if (_sim->uic_now() || _inside_tt ) {
 144     trace3("TRANSIENT::sweep uic_now solve", _time1, _sim->_time0, _inside_tt);
 145     advance_time();
 146     trace0("TRANSIENT::sweep advanced");
 147     if (!_inside_tt) {
 148       _sim->zero_voltages(); // ?
 149     }
 150     CARD_LIST::card_list.do_tr();    //evaluate_models
 151     while (!_sim->_late_evalq.empty()) {untested(); //BUG// encapsulation violation
 152       _sim->_late_evalq.front()->do_tr_last();
 153       _sim->_late_evalq.pop_front();
 154     }
 155     _converged = true;
 156   } else if (_cont_dc) {
 157     // continue from DC point.
 158     // advance_time();
 159     assert(_sim->_late_evalq.empty());
 160     while (!_sim->_late_evalq.empty()) {itested(); //BUG// encapsulation violation
 161       _sim->_late_evalq.front()->do_tr_last();
 162       _sim->_late_evalq.pop_front();
 163     }
 164     _converged = true;
 165   } else if (_cont_dc && _sim->analysis_is_static()) { untested();
 166   } else {
 167     _converged = solve_with_homotopy(OPT::DCBIAS,_trace);
 168     if (!_converged) {
 169       error(bWARNING, "did not converge\n");
 170     }else{
 171     }
 172   }
 173   trace0("TRANSIENT::sweep review...");
 174   review();
 175   trace0("TRANSIENT::sweep reviewed");
 176   _accepted = true;
 177   accept();
 178   trace0("TRANSIENT::sweep accepted");
 179
 180   {
 181     bool printnow = (_sim->_time0 == _tstart || _trace >= tALLTIME);
 182     if (printnow) {
 183       _sim->keep_voltages();
 184       outdata(_sim->_time0);
 185     }else{
 186       ++::status.hidden_steps;
 187     }
 188   }
 189
 190   assert(_tstrobe >=OPT::dtmin ); // == wont work because of CAUSE
 191                                 // we do only increase _time_by_user_request if
 192                                 // CAUSE == user.
 193                                 // the second step is always caused by initial guess...
 194                                 // BUG?
 195   trace3("TRANSIENT::sweep entering loop", (STEP_CAUSE)step_cause(), _sim->_time0, _sim->last_time());
 196   _edge_break = false;
 197   _dt_by_edge1 = -NEVER;
 198   _dt_by_edge0 = -NEVER;
 199   while (next()) {
 200     trace3("TRANSIENT::sweep loop ... ", (STEP_CAUSE)step_cause(), _sim->_time0, _sim->last_time());
 201     _sim->_bypass_ok = false;
 202     _sim->_phase = p_TRAN;
 203     _sim->_genout = gen();
 204     _converged = solve(OPT::TRHIGH,_trace);
 205
 206     _accepted = _converged && review();
 207
 208     if (_accepted) {
 209       trace1("TRANSIENT::sweep accepted", (STEP_CAUSE)step_cause());
 210       assert(_converged);
 211       accept();
 212       if (step_cause() == scUSER) {
 213         trace3("up order", _sim->_time0-_sim->_dtmin, _time_by_user_request, _sim->_time0+_sim->_dtmin);
 214         assert( up_order(_sim->_time0-_sim->_dtmin, _time_by_user_request, _sim->_time0+_sim->_dtmin)
 215             || _sim->_time0 == _tstop );
 216         ++(_sim->_stepno);
 217         trace1("TRANSIENT::sweep delivered req", _time_by_user_request);
 218         if (_time_by_user_request<_tstop){
 219           _time_by_user_request += _tstrobe;    /* advance user time */
 220           // _time_by_user_request = min(_time_by_user_request, (double)_tstop+_sim->_dtmin);
 221         } else {
 222           _time_by_user_request += _tstrobe;    /* advance user time */
 223         }
 224       }else{
 225       }
 226       assert(_sim->_time0 <= _time_by_user_request);
 227     }else{
 228       trace2("TRANSIENT::sweep NOT accepted", _sim->_time0, _trace);
 229       reject();
 230       assert(_time1 < _time_by_user_request);
 231     }
 232     {
 233       bool printnow =
 234         (_trace >= tREJECTED)
 235         || (_accepted && (_trace >= tALLTIME
 236                           || step_cause() == scUSER
 237                           || (!_tstrobe.has_hard_value() && _sim->_time0+_sim->_dtmin > _tstart)));
 238       if (printnow) {
 239         _sim->keep_voltages();
 240         trace2("TRANSIENT::sweep" ,_sim->last_time(), (double)_tstop);
 241         assert(_sim->last_time() < _tstop+_sim->_dtmin || _sim->_freezetime);
 242
 243         outdata(_sim->_time0);
 244         CKT_BASE::tr_behaviour_del = 0;
 245         CKT_BASE::tr_behaviour_rel = 0;
 246       }else{
 247         ++::status.hidden_steps;
 248       }
 249     }
 250
 251     if (!_converged && OPT::quitconvfail) { untested();
 252       outdata(_sim->_time0);
 253       throw Exception("convergence failure, giving up");
 254     }else{
 255     }
 256   }
 257
 258 //  _sim->_time0 = _sim->last_time(); // hack?
 259 //                                    // this is required because keep_voltage
 260 //                                    // calls set last_time to time0
 261
 262   for (uint_t ii = _sim->_lu.size(); ii > 0; --ii) {
 263     _sim->_nstat[ii].set_discont(disNONE);
 264   }
 265   assert(!_sim->_nstat[0].discont());
 266 }
 267 /*--------------------------------------------------------------------------*/
 268 void TRANSIENT::set_step_cause(STEP_CAUSE C)
 269 {
 270   switch (C) {
 271   case scITER_A:untested();
 272   case scADT:untested();
 273   case scUSER:
 274   case scEVENTQ:
 275   case scSKIP:
 276   case scITER_R:
 277   case scTE:
 278   case scAMBEVENT:
 279   case scINITIAL:
 280     ::status.control = C;
 281     break;
 282   case scNO_ADVANCE:untested();
 283   case scZERO:untested();
 284   case scSMALL:
 285   case scREJECT:
 286     ::status.control += C;
 287     break;
 288   case scGROW:
 289   case scLAST:
 290         incomplete();
 291   }
 292 }
 293 /*--------------------------------------------------------------------------*/
 294 int TRANSIENT::step_cause()const
 295 {
 296   return ::status.control;
 297 }
 298 /*--------------------------------------------------------------------------*/
 299 void TRANSIENT::first()
 300 {
 301   /* usually, _sim->_time0, time1 == 0, from setup */
 302   assert(_sim->_time0 == _time1);
 303   // assert(_sim->_time0 <= _tstart); // oops?
 304   ::status.review.start();
 305
 306   //_eq.Clear();                                        /* empty the queue */
 307   while (!_sim->_eq.empty()) {untested();
 308     _sim->_eq.pop();
 309   }
 310   _sim->_stepno = 0;
 311
 312   //_time_by_user_request = _sim->_time0 + _tstrobe;    /* set next user step */
 313   //set_step_cause(scUSER);
 314
 315   if (_sim->_time0 < _tstart) {                 // skip until _tstart
 316     set_step_cause(scINITIAL);                          // suppressed
 317     _time_by_user_request = _tstart;                    // set first strobe
 318   }else{                                        // no skip
 319     set_step_cause(scUSER);                             // strobe here
 320     _time_by_user_request = _sim->_time0 + _tstrobe;    // set next strobe
 321   }
 322
 323   ::status.hidden_steps = 0;
 324   ::status.review.stop();
 325 }
 326 /*--------------------------------------------------------------------------*/
 327 #define check_consistency() {                                           \
 328     trace4("", __LINE__, newtime, almost_fixed_time, fixed_time);       \
 329     assert(almost_fixed_time <= fixed_time);                            \
 330     assert(newtime <= fixed_time);                                      \
 331     /*assert(newtime == fixed_time || newtime <= fixed_time -_sim->_dtmin);*/   \
 332     assert(newtime <= almost_fixed_time);                               \
 333     /*assert(newtime == almost_fixed_time || newtime <= almost_fixed_time - _sim->_dtmin);*/ \
 334     assert(newtime > _time1);                                           \
 335     assert(newtime > reftime);                                          \
 336     assert(new_dt > 0.);                                                \
 337     assert(new_dt >= _sim->_dtmin);                                     \
 338     /*assert(new_dt == 0 || new_dt / _sim->_dtmin > 0.999999);*/        \
 339    /* assert(newtime <= _time_by_user_request);                         */ \
 340     /*assert(newtime == _time_by_user_request*/                         \
 341     /*    || newtime < _time_by_user_request - _sim->_dtmin);  */       \
 342   }
 343 #define check_consistency2() {                                          \
 344     assert(newtime > _time1);                                           \
 345     assert(new_dt > 0.);                                                \
 346     assert(new_dt >= _sim->_dtmin);                                             \
 347     assert(newtime <= _time_by_user_request);                           \
 348     /*assert(newtime == _time_by_user_request   */                      \
 349     /*     || newtime < _time_by_user_request - _sim->_dtmin);*/                \
 350   }
 351 /*--------------------------------------------------------------------------*/
 352 /* next: go to next time step
 353  * Set _sim->_time0 to the next time step, store the old one in time1.
 354  * Try several methods.  Take the one that gives the shortest step.
 355  */
 356 bool TRANSIENT::next()
 357 {
 358   ::status.review.start();
 359
 360   double old_dt = _sim->_time0 - _time1;
 361   assert(old_dt >= 0);
 362
 363
 364   double newtime = NEVER;
 365   double new_dt = NEVER;
 366   STEP_CAUSE new_control = scNO_ADVANCE;
 367
 368   if (_sim->_time0 == _time1) {
 369     // initial step -- could be either t==0 or continue
 370     // for the first time, just guess
 371     // make it 100x smaller than expected
 372     new_dt = std::max(_dtmax/100., _sim->_dtmin);
 373     newtime = _sim->_time0 + new_dt;
 374     new_control = scINITIAL;
 375
 376     if (_time_by_ambiguous_event < newtime) {
 377       newtime = _time_by_ambiguous_event;
 378       new_dt = _time_by_ambiguous_event - _time1;
 379       new_control = scAMBEVENT;
 380     }else{
 381     }
 382
 383     // UGLY duplicate
 384     if (up_order(newtime-_sim->_dtmin, _time_by_user_request, newtime+_sim->_dtmin)) {
 385       new_control = scUSER;
 386     }
 387   }else if (!_converged) {
 388     new_dt = old_dt / OPT::trstepshrink;
 389     newtime = _time_by_iteration_count = _time1 + new_dt;
 390     new_control = scITER_R;
 391     //    fprintf(stderr,".");
 392     //    _trace=tDEBUG;
 393   }else{
 394   }
 395   {
 396     double reftime;
 397     if (_accepted) {
 398       reftime = _sim->_time0;
 399     }else{
 400       reftime = _time1;
 401       trace0("rejected");
 402     }
 403     trace2("TRANSIENT::next ", step_cause(), old_dt);
 404     trace3("TRANSIENT::next ", _time1, _sim->_time0, reftime);
 405     trace2("TRANSIENT::next ", _time_by_user_request, _sim->_dtmin);
 406
 407     if (_time_by_user_request < newtime) {
 408       newtime = _time_by_user_request;
 409       new_dt = newtime - reftime;
 410       if (new_dt < _sim->_dtmin) { untested();
 411         // last step handler?
 412         new_dt = _sim->_dtmin;
 413         newtime = reftime + _sim->_dtmin;
 414       }else{untested();
 415       }
 416       new_control = scUSER;
 417     }else{
 418     }
 419     double fixed_time = newtime;
 420     double almost_fixed_time = newtime;
 421     trace2("TRANSIENT::next", _time_by_user_request, newtime);
 422     check_consistency();
 423
 424
 425     // event queue, events that absolutely will happen
 426     // exact time.  NOT ok to move or omit, even by _sim->_dtmin
 427     // some action is associated with it.
 428     if (!_sim->_eq.empty() && _sim->_eq.top() < newtime) {
 429       trace2("TRANSIENT eventq", newtime, _time1);
 430       newtime = _sim->_eq.top();
 431       assert( newtime >= _time1 );
 432       new_dt = newtime - reftime;
 433       if (new_dt < _sim->_dtmin) {untested();
 434         //new_dt = _sim->_dtmin;
 435         //newtime = reftime + new_dt;
 436       }else{
 437       }
 438       new_control = scEVENTQ;
 439       fixed_time = newtime;
 440       almost_fixed_time = newtime;
 441       trace2("checking", reftime, newtime);
 442       check_consistency();
 443     }else if ( !_sim->_eq.empty()  ) {
 444       trace2("TRANSIENT skipping non empty eq", _time1, _sim->_eq.top() );
 445     } else {
 446       trace0("TRANSIENT no events pending");
 447     }
 448     // device events that may not happen
 449     // not sure of exact time.  will be rescheduled if wrong.
 450     // ok to move by _sim->_dtmin.  time is not that accurate anyway.
 451     trace1("next ambevt", _time_by_ambiguous_event);
 452     if (_time_by_ambiguous_event < newtime - _sim->_dtmin) {
 453       if (_time_by_ambiguous_event < _time1 + 2*_sim->_dtmin) {untested();
 454         double mintime = _time1 + 2*_sim->_dtmin;
 455         if (newtime - _sim->_dtmin < mintime) {untested();
 456           newtime = mintime;
 457           new_control = scAMBEVENT;
 458         }else{ untested();
 459         }
 460       }else{
 461         newtime = _time_by_ambiguous_event;
 462         new_control = scAMBEVENT;
 463       }
 464       new_dt = newtime - reftime;
 465       almost_fixed_time = newtime;
 466       check_consistency();
 467     }else{
 468     }
 469
 470     // device error estimates
 471     if (_time_by_error_estimate < newtime - _sim->_dtmin) {
 472       newtime = _time_by_error_estimate;
 473       new_dt = newtime - reftime;
 474       new_control = scTE;
 475       check_consistency();
 476       trace3("TRANSIENT::next err", newtime, new_dt, _time_by_error_estimate);
 477     }else{
 478       trace3("TRANSIENT::next", newtime, new_dt, _time_by_error_estimate);
 479     }
 480
 481     // skip parameter
 482     if (new_dt > _dtmax) {
 483       if (new_dt > _dtmax + _sim->_dtmin) {
 484         new_control = scSKIP;
 485       }else{
 486       }
 487       new_dt = _dtmax;
 488       newtime = reftime + new_dt;
 489       assert(newtime <= _time_by_user_request);
 490       check_consistency();
 491     }else{
 492     }
 493
 494     // converged but with more iterations than we like
 495     if ((new_dt > (old_dt + _sim->_dtmin) * OPT::trstephold)
 496         && _sim->exceeds_iteration_limit(OPT::TRLOW)) { untested();
 497       assert(_accepted);
 498       new_dt = old_dt * OPT::trstephold;
 499       newtime = reftime + new_dt;
 500       new_control = scITER_A;
 501       check_consistency();
 502     }else{
 503     }
 504
 505     // limit growth
 506     if (_sim->analysis_is_tran_dynamic() && new_dt > old_dt * OPT::trstepgrow) { untested();
 507       new_dt = old_dt * OPT::trstepgrow;
 508       newtime = reftime + new_dt;
 509       new_control = scADT;
 510       check_consistency();
 511     }else{
 512     }
 513
 514     // quantize
 515     if (newtime < almost_fixed_time) {
 516       assert(new_dt >= 0);
 517       if (newtime < _sim->_time0) {
 518         assert(reftime == _time1);
 519         assert(reftime < _sim->_time0); // not moving forward
 520         // try to pick a step that will end up repeating the rejected step
 521         // with an integer number of same size steps
 522         double target_dt = _sim->_time0 - reftime;
 523         assert(target_dt > new_dt);
 524         double steps = 1 + floor((target_dt - _sim->_dtmin) / new_dt);
 525         assert(steps > 0);
 526         new_dt = target_dt / steps;
 527         newtime = reftime + new_dt;
 528         check_consistency();
 529       }else if (newtime > reftime + old_dt*.8
 530           && newtime < reftime + old_dt*1.5
 531           && reftime + old_dt <= almost_fixed_time) {
 532         // new_dt is close enough to old_dt.
 533         // use old_dt, to avoid a step change.
 534         old_dt = max(_sim->_dtmin,old_dt);  // eliminate numerical noise to pass assertion
 535         assert(reftime == _sim->_time0); // moving forward
 536         assert(reftime > _time1);
 537         new_dt = old_dt;
 538         newtime = reftime + new_dt;
 539         if (newtime > almost_fixed_time) { untested();
 540           new_control = scAMBEVENT;
 541           newtime = almost_fixed_time;
 542           new_dt = newtime - reftime;
 543         }else{
 544         }
 545         trace4("TRANSIENT::next quantized", new_dt, _sim->_dtmin, _sim->_dtmin-new_dt, new_control );
 546         check_consistency();
 547       }else{
 548         // There will be a step change.
 549         // Try to choose one that we will keep for a while.
 550         // Choose new_dt to be in integer fraction of target_dt.
 551         assert(reftime == _sim->_time0); // moving forward
 552         assert(reftime > _time1);
 553         double target_dt = fixed_time - reftime;
 554         assert(target_dt >= new_dt);
 555         double steps = 1 + floor((target_dt - _sim->_dtmin) / new_dt);
 556         assert(steps > 0);
 557         new_dt = target_dt / steps;
 558         trace3("TRANSIENT::next step change", reftime, new_dt, reftime+new_dt);
 559         newtime = reftime + new_dt;
 560         check_consistency();
 561       }
 562     }else{
 563       assert(newtime == almost_fixed_time);
 564     }
 565
 566     // trap time step too small
 567     if (!_accepted && new_dt < _sim->_dtmin) { untested();
 568       new_dt = _sim->_dtmin;
 569       newtime = reftime + new_dt;
 570       new_control = scSMALL;
 571       check_consistency();
 572     }else{
 573     }
 574
 575     // if all that makes it close to user_requested, make it official
 576     if (up_order(newtime-_sim->_dtmin, _time_by_user_request, newtime+_sim->_dtmin)) {
 577       //newtime = _time_by_user_request;
 578       //new_dt = newtime - reftime;
 579       new_control = scUSER;
 580     }
 581   } // end of else converged
 582
 583   set_step_cause(new_control);
 584
 585   trace4("TRANSIENT::next got it i think", newtime, new_control, newtime-_sim->_time0, _time_by_user_request);
 586
 587   /* check to be sure */
 588   if (newtime < _time1 + _sim->_dtmin) {
 589     /* It's really bad. */
 590     /* Reject the most recent step, back up as much as possible, */
 591     /* and creep along */
 592     trace3("TRANSIENT::next ", newtime, _time1, _sim->_dtmin );
 593     assert(!_accepted);
 594     assert(step_cause() < scREJECT);
 595     assert(step_cause() >= 0);
 596     error(bDANGER,"non-recoverable: " + TR::step_cause[step_cause()] + "\n");
 597     error(bDANGER, "newtime=%e  rejectedtime=%e  oldtime=%e  using=%e\n",
 598           newtime, _sim->_time0, _time1, _time1 + _sim->_dtmin);
 599     newtime = _time1 + _sim->_dtmin;
 600     assert(newtime <= _time_by_user_request);
 601     set_step_cause(scSMALL);
 602     //check_consistency2();
 603     throw Exception("tried everything, still doesn't work, giving up");
 604   }else if (newtime < _sim->_time0) {
 605     /* Reject the most recent step. */
 606     /* We have faith that it will work with a smaller time step. */
 607     assert(!_accepted);
 608     assert(newtime >= _time1 + _sim->_dtmin);
 609     error(bLOG, "backwards time step\n");
 610     error(bLOG, "newtime=%e  rejectedtime=%e  oldtime=%e\n", newtime, _sim->_time0, _time1);
 611     set_step_cause(scREJECT);
 612     _sim->mark_inc_mode_bad();
 613     check_consistency2();
 614 #if 0
 615   }else if (_sim->_time0 >= _time_by_user_request) {
 616     trace3("TRANSIENT::next: already there", _sim->_time0, _time1, newtime );
 617     _time1 = _sim->_time0;
 618     new_dt = 0;
 619     check_consistency2();
 620 #endif
 621   }else if (newtime < _sim->_time0 + _sim->_dtmin) { untested();
 622     /* Another evaluation at the same time. */
 623     /* Keep the most recent step, but creep along. */
 624     assert(newtime > _sim->_time0 - _sim->_dtmin);
 625     error(bDANGER, "zero time step\n");
 626     error(bDANGER, "newtime=%e  rejectedtime=%e delta=%e time1=%e requested=%e dtmin=%e, control=%d\n",
 627         newtime, _sim->_time0, newtime-_sim->_time0, _time1, _time_by_user_request, _sim->_dtmin, step_cause());
 628     if (_accepted) {untested();
 629       _time1 = _sim->_time0;
 630     }else{untested();
 631       assert(_converged);
 632     }
 633     trace3( "TRANSIENT::next:", newtime, _time1, new_dt);
 634     check_consistency2();
 635     newtime = _sim->_time0 + _sim->_dtmin;
 636     if (newtime > _time_by_user_request) { untested();
 637       newtime = _time_by_user_request;
 638       set_step_cause(scUSER);
 639     }else{
 640     }
 641     assert (newtime<=_tstop);
 642     set_step_cause(scZERO);
 643     check_consistency2();
 644   }else{
 645     assert(_accepted);
 646     assert(newtime >= _sim->_time0 + _sim->_dtmin);
 647     /* All is OK.  Moving on. */
 648     /* Keep value of newtime */
 649     _time1 = _sim->_time0;
 650     check_consistency2();
 651   }
 652   _sim->_time0 = newtime;
 653   assert(_sim->_time0 <= _time_by_user_request);
 654
 655   /* advance event queue (maybe) */
 656   /* We already looked at it.  Dump what's on top if we took it. */
 657
 658 #ifndef ALT_CQ
 659
 660   while (!_sim->_eq.empty() && _sim->_eq.top() <= _sim->_time0) {
 661     trace1("eq", _sim->_eq.top());
 662     _sim->_eq.pop();
 663   }
 664   while (!_sim->_eq.empty() && _sim->_eq.top() < _sim->_time0 + _sim->_dtmin) {
 665     trace1("eq-extra", _sim->_eq.top());
 666     _sim->_eq.pop();
 667   }
 668 #endif
 669   // event queue cleanup moved....
 670   //BUG// what if it is later rejected?  It's lost!
 671   // -> why not move to tr_advance? tr_accept? hmmm
 672   //
 673
 674   if(_time1 < _tstop - _sim->_dtmin
 675       && _sim->_time0 > _tstop + _sim->_dtmin ) {
 676     _sim->_time0 = _tstop;
 677     check_consistency2();
 678   } else {
 679     check_consistency2();
 680   }
 681   assert(_sim->_time0 <= _time_by_user_request);
 682
 683   check_consistency2();
 684   ++steps_total_;
 685   ::status.review.stop();
 686   bool ret= _sim->_time0 <= _tstop; // throw away last step if it helps.  + _sim->_dtmin;
 687   ret= _sim->_time0 < _tstop + _sim->_dtmin; // this once worked
 688   if(_accepted && _edge_break){
 689     trace3("accepted edge", _sim->_time0, _time1, _sim->_time0 - _time1);
 690     return false;
 691   }else if(_edge_break){
 692     trace5("unaccepted edge", _sim->_time0, _sim->v0dist_min(), _time1, _time_by_ambiguous_event, STEP_CAUSE(step_cause()));
 693   }
 694
 695   return (ret);
 696 }
 697 /*--------------------------------------------------------------------------*/
 698 bool TRANSIENT::review()
 699 {
 700   ::status.review.start();
 701   _sim->count_iterations(iTOTAL);
 702
 703   TIME_PAIR time_by = CARD_LIST::card_list.tr_review();
 704   _time_by_error_estimate = time_by._error_estimate;
 705   _time_by_ambiguous_event = time_by._event;
 706   _edge_break = false;
 707
 708   double dt0 = _sim->_time0 - _time1;
 709   trace3("review times", _sim->_time0, dt0, _time1);
 710
 711   // fixme: move to a component (?)... maybe later.
 712   if (_sim->_time0 == 0.){
 713   }else if (dt0 == 0. ){
 714   }else if (_sim->_time0 < _tstop * .2){
 715     // incomplete();
 716   }else if (_edge_detect & edEVT){
 717     double dist = 0;
 718     double dt_factor = _sim->v0dist_min(&dist);
 719     double dt_by_edge = dt_factor * dt0;
 720     // if(dt_by_edge < 0. && dtedge1 > 0){ untested();
 721     //   dt_by_edge = (dt_by_edge + dtedge1)*.5;
 722     // }
 723
 724     double time_by_edge = dt_by_edge + _time1;
 725     trace4("edge", time_by_edge, _sim->_time0, dt_by_edge, dist);
 726
 727     if(dist>10){
 728       time_by_edge = NEVER;
 729     }else if(fabs(dt_by_edge) > 10 * fabs(_dt_by_edge1)){
 730       time_by_edge = NEVER;
 731     }else if(_dt_by_edge1 < 0.){
 732       time_by_edge = NEVER;
 733     }else if(_dt_by_edge1 < dt0){
 734       time_by_edge = NEVER;
 735       // ?!
 736     }else if(dt_by_edge < 0. && _dt_by_edge1 > 0. ){
 737       trace3("negative", dt_by_edge, _sim->_time0, _dt_by_edge1);
 738       time_by_edge = _time1 + dt0 * .3;
 739       _edge_break = _edge_detect & edBREAK;
 740     }else if(dt_by_edge < 0. ){ untested();
 741       time_by_edge = _sim->_dt0 * .5 + _time1;
 742     }else if(dt_by_edge < dt0 && _dt_by_edge1 > dt0){
 743       trace4("got edge", dt_by_edge, _sim->_time0, _dt_by_edge1, dt_factor);
 744       _edge_break = _edge_detect & edBREAK; // edge between time0 and time1.
 745                                       // if it gets accepted, we are done
 746 //      time_by_edge = dt_by_edge * .5 + _time1;
 747     }else if(dt_by_edge > dt0){
 748       // everything fine. edge is far ahead.
 749     }else if(dt_by_edge + _sim->_dtmin > dt0 && _dt_by_edge1 > dt0 ){ untested();
 750       trace4("close to edge", dt0, dt_by_edge, _sim->_time0, _dt_by_edge1);
 751       _edge_break = _edge_detect & edBREAK; // edge close enogh
 752     }else if(dt_by_edge < dt0){ untested();
 753       // past edge
 754       time_by_edge = NEVER;
 755     }
 756
 757     if(!(_edge_detect & edEVT)){ untested();
 758     }else if(time_by_edge < _time_by_ambiguous_event){
 759       trace5("setting evt", time_by_edge, _sim->_time0, dt_by_edge, _time_by_ambiguous_event, dt_by_edge / _sim->_dtmin);
 760       _time_by_ambiguous_event = time_by_edge;
 761     }else{
 762     }
 763
 764     _dt_by_edge0 = dt_by_edge;
 765   }
 766
 767   // limit minimum time step
 768   // 2*_sim->_dtmin because _time[1] + _sim->_dtmin might be == _time[0].
 769   if (_time_by_ambiguous_event < _time1 + 2*_sim->_dtmin) {
 770     _time_by_ambiguous_event = _time1 + 2*_sim->_dtmin;
 771   }else{
 772   }
 773   // force advance when time too close to previous
 774   if (std::abs(_time_by_ambiguous_event - _sim->_time0) < 2*_sim->_dtmin) {
 775     _time_by_ambiguous_event = _sim->_time0 + 2*_sim->_dtmin;
 776   }else{
 777   }
 778
 779   if (time_by._error_estimate < _time1 + 2*_sim->_dtmin) {
 780     _time_by_error_estimate = _time1 + 2*_sim->_dtmin;
 781   }else{
 782     _time_by_error_estimate = time_by._error_estimate;
 783   }
 784   if (std::abs(_time_by_error_estimate - _sim->_time0) < 1.1*_sim->_dtmin) {
 785     _time_by_error_estimate = _sim->_time0 + 1.1*_sim->_dtmin;
 786   }else{
 787   }
 788
 789   ::status.review.stop();
 790
 791   return (_time_by_error_estimate > _sim->_time0  &&  _time_by_ambiguous_event > _sim->_time0);
 792 }
 793 /*--------------------------------------------------------------------------*/
 794 void TRANSIENT::accept()
 795 {
 796   ::status.accept.start();
 797   for (unsigned ii = _sim->_lu.size(); ii >= 1; --ii) {
 798       _sim->_nstat[ii].set_discont(disNONE);
 799   }
 800
 801   _dt_by_edge1 = _dt_by_edge0;
 802
 803 #ifdef ALT_CQ_PRE
 804   while (!_sim->_eq.empty() && _sim->_eq.top() < _sim->_time0 + _sim->_dtmin) {itested();
 805     trace1("TRANSIENT::accept eq-pop-extra", _sim->_eq.top());
 806     _sim->_eq.pop();
 807   }
 808 #endif
 809
 810   trace2("TRANSIENT::accept dt0", _sim->_dt0, _sim->_time0 - _time1);
 811   _sim->_dt0 = _sim->_time0 - _time1;
 812   if(_sim->_dt0 <=0) assert (_sim->_stepno == 0);
 813   _sim->set_limit();
 814   if (OPT::traceload) { // traceload == "use queue"
 815     while (!_sim->_acceptq.empty()) {
 816       trace1("TRANSIENT::accept", _sim->_acceptq.back()->long_label());
 817       _sim->_acceptq.back()->tr_accept();
 818       _sim->_acceptq.pop_back();
 819     }
 820   }else{untested();
 821     _sim->_acceptq.clear();
 822     CARD_LIST::card_list.tr_accept();
 823   }
 824 //  tmp hack don't know yet how to fix (always_q_for_accept?)
 825   ADP_LIST::adp_list.tr_accept();
 826
 827   ++steps_accepted_;
 828   if( _sim->analysis_is_tt() || OPT::trage ) {
 829     trace0( "TRANSIENT::accept: done stressing cardlist");
 830     if ( OPT::trage ) {
 831       incomplete();
 832       CARD_LIST::card_list.do_forall( &CARD::do_tt );
 833     }
 834   }
 835 # ifdef ALT_CQ
 836   while (!_sim->_eq.empty() && _sim->_eq.top() <= _sim->_time0) {
 837     trace1("eq", _sim->_eq.top());
 838     _sim->_eq.pop();
 839   }
 840 # endif
 841   _sim->_nstat[0].set_discont(disNONE);
 842   ::status.accept.stop();
 843 }
 844 /*--------------------------------------------------------------------------*/
 845 void TRANSIENT::reject()
 846 {
 847   ::status.accept.start();
 848   _sim->_acceptq.clear();
 849   ++steps_rejected_;
 850   ::status.accept.stop();
 851 }
 852 /*--------------------------------------------------------------------------*/
 853 /*--------------------------------------------------------------------------*/
 854 // vim:ts=8:sw=2:noet: