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{itested();
 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         if(reftime <= _time1){ unreachable();
 553           // this is a weird bug.
 554           // FIX LATER
 555         }else{
 556         }
 557         double target_dt = fixed_time - reftime;
 558         assert(target_dt >= new_dt);
 559         double steps = 1 + floor((target_dt - _sim->_dtmin) / new_dt);
 560         assert(steps > 0);
 561         new_dt = target_dt / steps;
 562         trace3("TRANSIENT::next step change", reftime, new_dt, reftime+new_dt);
 563         newtime = reftime + new_dt;
 564         check_consistency();
 565       }
 566     }else{
 567       assert(newtime == almost_fixed_time);
 568     }
 569
 570     // trap time step too small
 571     if (!_accepted && new_dt < _sim->_dtmin) { untested();
 572       new_dt = _sim->_dtmin;
 573       newtime = reftime + new_dt;
 574       new_control = scSMALL;
 575       check_consistency();
 576     }else{
 577     }
 578
 579     // if all that makes it close to user_requested, make it official
 580     if (up_order(newtime-_sim->_dtmin, _time_by_user_request, newtime+_sim->_dtmin)) {
 581       //newtime = _time_by_user_request;
 582       //new_dt = newtime - reftime;
 583       new_control = scUSER;
 584     }
 585   } // end of else converged
 586
 587   set_step_cause(new_control);
 588
 589   trace4("TRANSIENT::next got it i think", newtime, new_control, newtime-_sim->_time0, _time_by_user_request);
 590
 591   /* check to be sure */
 592   if (newtime < _time1 + _sim->_dtmin) {
 593     /* It's really bad. */
 594     /* Reject the most recent step, back up as much as possible, */
 595     /* and creep along */
 596     trace3("TRANSIENT::next ", newtime, _time1, _sim->_dtmin );
 597     assert(!_accepted);
 598     assert(step_cause() < scREJECT);
 599     assert(step_cause() >= 0);
 600     error(bDANGER,"non-recoverable: " + TR::step_cause[step_cause()] + "\n");
 601     error(bDANGER, "newtime=%e  rejectedtime=%e  oldtime=%e  using=%e\n",
 602           newtime, _sim->_time0, _time1, _time1 + _sim->_dtmin);
 603     newtime = _time1 + _sim->_dtmin;
 604     assert(newtime <= _time_by_user_request);
 605     set_step_cause(scSMALL);
 606     //check_consistency2();
 607     throw Exception("tried everything, still doesn't work, giving up");
 608   }else if (newtime < _sim->_time0) {
 609     /* Reject the most recent step. */
 610     /* We have faith that it will work with a smaller time step. */
 611     assert(!_accepted);
 612     assert(newtime >= _time1 + _sim->_dtmin);
 613     error(bLOG, "backwards time step\n");
 614     error(bLOG, "newtime=%e  rejectedtime=%e  oldtime=%e\n", newtime, _sim->_time0, _time1);
 615     set_step_cause(scREJECT);
 616     _sim->mark_inc_mode_bad();
 617     check_consistency2();
 618 #if 0
 619   }else if (_sim->_time0 >= _time_by_user_request) {
 620     trace3("TRANSIENT::next: already there", _sim->_time0, _time1, newtime );
 621     _time1 = _sim->_time0;
 622     new_dt = 0;
 623     check_consistency2();
 624 #endif
 625   }else if (newtime < _sim->_time0 + _sim->_dtmin) { untested();
 626     /* Another evaluation at the same time. */
 627     /* Keep the most recent step, but creep along. */
 628     assert(newtime > _sim->_time0 - _sim->_dtmin);
 629     error(bDANGER, "zero time step\n");
 630     error(bDANGER, "newtime=%e  rejectedtime=%e delta=%e time1=%e requested=%e dtmin=%e, control=%d\n",
 631         newtime, _sim->_time0, newtime-_sim->_time0, _time1, _time_by_user_request, _sim->_dtmin, step_cause());
 632     if (_accepted) {untested();
 633       _time1 = _sim->_time0;
 634     }else{untested();
 635       assert(_converged);
 636     }
 637     trace3( "TRANSIENT::next:", newtime, _time1, new_dt);
 638     check_consistency2();
 639     newtime = _sim->_time0 + _sim->_dtmin;
 640     if (newtime > _time_by_user_request) { untested();
 641       newtime = _time_by_user_request;
 642       set_step_cause(scUSER);
 643     }else{
 644     }
 645     assert (newtime<=_tstop);
 646     set_step_cause(scZERO);
 647     check_consistency2();
 648   }else{
 649     assert(_accepted);
 650     assert(newtime >= _sim->_time0 + _sim->_dtmin);
 651     /* All is OK.  Moving on. */
 652     /* Keep value of newtime */
 653     _time1 = _sim->_time0;
 654     check_consistency2();
 655   }
 656   _sim->_time0 = newtime;
 657   assert(_sim->_time0 <= _time_by_user_request);
 658
 659   /* advance event queue (maybe) */
 660   /* We already looked at it.  Dump what's on top if we took it. */
 661
 662 #ifndef ALT_CQ
 663
 664   while (!_sim->_eq.empty() && _sim->_eq.top() <= _sim->_time0) {
 665     trace1("eq", _sim->_eq.top());
 666     _sim->_eq.pop();
 667   }
 668   while (!_sim->_eq.empty() && _sim->_eq.top() < _sim->_time0 + _sim->_dtmin) {
 669     trace1("eq-extra", _sim->_eq.top());
 670     _sim->_eq.pop();
 671   }
 672 #endif
 673   // event queue cleanup moved....
 674   //BUG// what if it is later rejected?  It's lost!
 675   // -> why not move to tr_advance? tr_accept? hmmm
 676   //
 677
 678   if(_time1 < _tstop - _sim->_dtmin
 679       && _sim->_time0 > _tstop + _sim->_dtmin ) {
 680     _sim->_time0 = _tstop;
 681     check_consistency2();
 682   } else {
 683     check_consistency2();
 684   }
 685   assert(_sim->_time0 <= _time_by_user_request);
 686
 687   check_consistency2();
 688   ++steps_total_;
 689   ::status.review.stop();
 690   bool ret= _sim->_time0 <= _tstop; // throw away last step if it helps.  + _sim->_dtmin;
 691   ret= _sim->_time0 < _tstop + _sim->_dtmin; // this once worked
 692   if(_accepted && _edge_break){
 693     trace3("accepted edge", _sim->_time0, _time1, _sim->_time0 - _time1);
 694     return false;
 695   }else if(_edge_break){
 696     trace5("unaccepted edge", _sim->_time0, _sim->v0dist_min(), _time1, _time_by_ambiguous_event, STEP_CAUSE(step_cause()));
 697   }
 698
 699   return (ret);
 700 }
 701 /*--------------------------------------------------------------------------*/
 702 bool TRANSIENT::review()
 703 {
 704   ::status.review.start();
 705   _sim->count_iterations(iTOTAL);
 706
 707   TIME_PAIR time_by = CARD_LIST::card_list.tr_review();
 708   _time_by_error_estimate = time_by._error_estimate;
 709   _time_by_ambiguous_event = time_by._event;
 710   _edge_break = false;
 711
 712   double dt0 = _sim->_time0 - _time1;
 713   trace3("review times", _sim->_time0, dt0, _time1);
 714
 715   // fixme: move to a component (?)... maybe later.
 716   if (_sim->_time0 == 0.){
 717   }else if (dt0 == 0. ){
 718   }else if (_sim->_time0 < _tstop * .2){
 719     // incomplete();
 720   }else if (_edge_detect & edEVT){
 721     double dist = 0;
 722     double dt_factor = _sim->v0dist_min(&dist);
 723     double dt_by_edge = dt_factor * dt0;
 724     // if(dt_by_edge < 0. && dtedge1 > 0){ untested();
 725     //   dt_by_edge = (dt_by_edge + dtedge1)*.5;
 726     // }
 727
 728     double time_by_edge = dt_by_edge + _time1;
 729     trace4("edge", time_by_edge, _sim->_time0, dt_by_edge, dist);
 730
 731     if(dist>10){
 732       time_by_edge = NEVER;
 733     }else if(fabs(dt_by_edge) > 10 * fabs(_dt_by_edge1)){
 734       time_by_edge = NEVER;
 735     }else if(_dt_by_edge1 < 0.){
 736       time_by_edge = NEVER;
 737     }else if(_dt_by_edge1 < dt0){
 738       time_by_edge = NEVER;
 739       // ?!
 740     }else if(dt_by_edge < 0. && _dt_by_edge1 > 0. ){
 741       trace3("negative", dt_by_edge, _sim->_time0, _dt_by_edge1);
 742       time_by_edge = _time1 + dt0 * .3;
 743       _edge_break = _edge_detect & edBREAK;
 744     }else if(dt_by_edge < 0. ){ untested();
 745       time_by_edge = _sim->_dt0 * .5 + _time1;
 746     }else if(dt_by_edge < dt0 && _dt_by_edge1 > dt0){
 747       trace4("got edge", dt_by_edge, _sim->_time0, _dt_by_edge1, dt_factor);
 748       _edge_break = _edge_detect & edBREAK; // edge between time0 and time1.
 749                                       // if it gets accepted, we are done
 750 //      time_by_edge = dt_by_edge * .5 + _time1;
 751     }else if(dt_by_edge > dt0){
 752       // everything fine. edge is far ahead.
 753     }else if(dt_by_edge + _sim->_dtmin > dt0 && _dt_by_edge1 > dt0 ){ untested();
 754       trace4("close to edge", dt0, dt_by_edge, _sim->_time0, _dt_by_edge1);
 755       _edge_break = _edge_detect & edBREAK; // edge close enogh
 756     }else if(dt_by_edge < dt0){ untested();
 757       // past edge
 758       time_by_edge = NEVER;
 759     }
 760
 761     if(!(_edge_detect & edEVT)){ untested();
 762     }else if(time_by_edge < _time_by_ambiguous_event){
 763       trace5("setting evt", time_by_edge, _sim->_time0, dt_by_edge, _time_by_ambiguous_event, dt_by_edge / _sim->_dtmin);
 764       _time_by_ambiguous_event = time_by_edge;
 765     }else{
 766     }
 767
 768     _dt_by_edge0 = dt_by_edge;
 769   }
 770
 771   // limit minimum time step
 772   // 2*_sim->_dtmin because _time[1] + _sim->_dtmin might be == _time[0].
 773   if (_time_by_ambiguous_event < _time1 + 2*_sim->_dtmin) {
 774     _time_by_ambiguous_event = _time1 + 2*_sim->_dtmin;
 775   }else{
 776   }
 777   // force advance when time too close to previous
 778   if (std::abs(_time_by_ambiguous_event - _sim->_time0) < 2*_sim->_dtmin) {
 779     _time_by_ambiguous_event = _sim->_time0 + 2*_sim->_dtmin;
 780   }else{
 781   }
 782
 783   if (time_by._error_estimate < _time1 + 2*_sim->_dtmin) {
 784     _time_by_error_estimate = _time1 + 2*_sim->_dtmin;
 785   }else{
 786     _time_by_error_estimate = time_by._error_estimate;
 787   }
 788   if (std::abs(_time_by_error_estimate - _sim->_time0) < 1.1*_sim->_dtmin) {
 789     _time_by_error_estimate = _sim->_time0 + 1.1*_sim->_dtmin;
 790   }else{
 791   }
 792
 793   ::status.review.stop();
 794
 795   return (_time_by_error_estimate > _sim->_time0  &&  _time_by_ambiguous_event > _sim->_time0);
 796 }
 797 /*--------------------------------------------------------------------------*/
 798 void TRANSIENT::accept()
 799 {
 800   ::status.accept.start();
 801   for (unsigned ii = _sim->_lu.size(); ii >= 1; --ii) {
 802       _sim->_nstat[ii].set_discont(disNONE);
 803   }
 804
 805   _dt_by_edge1 = _dt_by_edge0;
 806
 807 #ifdef ALT_CQ_PRE
 808   while (!_sim->_eq.empty() && _sim->_eq.top() < _sim->_time0 + _sim->_dtmin) {itested();
 809     trace1("TRANSIENT::accept eq-pop-extra", _sim->_eq.top());
 810     _sim->_eq.pop();
 811   }
 812 #endif
 813
 814   trace2("TRANSIENT::accept dt0", _sim->_dt0, _sim->_time0 - _time1);
 815   _sim->_dt0 = _sim->_time0 - _time1;
 816   if(_sim->_dt0 <=0) assert (_sim->_stepno == 0);
 817   _sim->set_limit();
 818   if (OPT::traceload) { // traceload == "use queue"
 819     while (!_sim->_acceptq.empty()) {
 820       trace1("TRANSIENT::accept", _sim->_acceptq.back()->long_label());
 821       _sim->_acceptq.back()->tr_accept();
 822       _sim->_acceptq.pop_back();
 823     }
 824   }else{untested();
 825     _sim->_acceptq.clear();
 826     CARD_LIST::card_list.tr_accept();
 827   }
 828 //  tmp hack don't know yet how to fix (always_q_for_accept?)
 829   ADP_LIST::adp_list.tr_accept();
 830
 831   ++steps_accepted_;
 832   if( _sim->analysis_is_tt() || OPT::trage ) {
 833     trace0( "TRANSIENT::accept: done stressing cardlist");
 834     if ( OPT::trage ) {
 835       incomplete();
 836       CARD_LIST::card_list.do_forall( &CARD::do_tt );
 837     }
 838   }
 839 # ifdef ALT_CQ
 840   while (!_sim->_eq.empty() && _sim->_eq.top() <= _sim->_time0) {
 841     trace1("eq", _sim->_eq.top());
 842     _sim->_eq.pop();
 843   }
 844 # endif
 845   _sim->_nstat[0].set_discont(disNONE);
 846   ::status.accept.stop();
 847 }
 848 /*--------------------------------------------------------------------------*/
 849 void TRANSIENT::reject()
 850 {
 851   ::status.accept.start();
 852   _sim->_acceptq.clear();
 853   ++steps_rejected_;
 854   ::status.accept.stop();
 855 }
 856 /*--------------------------------------------------------------------------*/
 857 /*--------------------------------------------------------------------------*/
 858 // vim:ts=8:sw=2:noet: