lib/s__solve.cc

   1 /*                              -*- 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  * solve one step of a transient or dc analysis
  23  */
  24 #include "e_cardlist.h"
  25 #include "u_status.h"
  26 #include "e_node.h"
  27 #include "s__.h"
  28 #include "io_matrix.h"
  29 /*--------------------------------------------------------------------------*/
  30 //      bool    SIM::solve(int,int);
  31 //      void    SIM::finish_building_evalq();
  32 //      void    SIM::advance_time();
  33 //      void    SIM::set_flags();
  34 //      void    SIM::clear_arrays();
  35 //      void    SIM::evaluate_models();
  36 //      void    SIM::set_damp();
  37 //      void    SIM::load_matrix();
  38 //      void    SIM::solve_equations();
  39 /*--------------------------------------------------------------------------*/
  40 static bool converged = false;
  41 /*--------------------------------------------------------------------------*/
  42 bool SIM::solve(OPT::ITL itl, TRACE trace)
  43 {
  44   trace1("SIM::solve", trace);
  45   converged = false;
  46   int convergedcount = 0;
  47
  48   _sim->reset_iteration_counter(iSTEP);
  49   advance_time();
  50
  51   _sim->_damp = OPT::dampmax; // default 1.0
  52
  53   do{
  54     if (( trace & (tMATRIX-1) ) >= tITERATION) {
  55       trace1("SIM::solve results", _sim->sim_mode());
  56       SIM::print_results(static_cast<double>(-_sim->iteration_number()));
  57       // Hack: Added SIM:: / print_results has been overloade but puts out the wrong data
  58     }
  59     set_flags();
  60     clear_arrays();
  61     finish_building_evalq();
  62
  63     _sim->count_iterations(iPRINTSTEP);
  64     _sim->count_iterations(iSTEP);
  65     _sim->count_iterations(_sim->_mode);
  66     _sim->count_iterations(iTOTAL);
  67
  68     evaluate_models();
  69
  70     if (converged) {
  71       if (_sim->_limiting) {
  72         error(bDEBUG, "converged beyond limit, resetting limit\n");
  73         _sim->set_limit();
  74         convergedcount = 0;
  75       }else{
  76         ++convergedcount;
  77       }
  78     }else{
  79       convergedcount = 0;
  80     }
  81     if (convergedcount <= OPT::itermin) {
  82       converged = false;
  83     }
  84
  85     if (!converged || !OPT::fbbypass || _sim->_damp < .99) {
  86       set_damp();
  87       load_matrix();
  88       assert(_sim->_loadq.empty());
  89       trace0("solve_equations");
  90       try{
  91         solve_equations(trace);
  92       }catch(Exception e){
  93        error(bWARNING, "%s\n", e.message().c_str());
  94        return false;
  95       }
  96
  97       trace0("solve_equations done");
  98     }
  99   }while (!converged && !_sim->exceeds_iteration_limit(itl));
 100
 101   //for(unsigned int y=0; y<_sim->_loadq.size(); y++) {
 102     //untested0( ("after solve loadq " + _sim->_loadq[y]->long_label()).c_str() );
 103   //}
 104   return converged;
 105 }
 106 /*--------------------------------------------------------------------------*/
 107 bool SIM::solve_with_homotopy(OPT::ITL itl, TRACE trace)
 108 {
 109   solve(itl, trace);
 110   trace3("plain", _sim->_iter[iSTEP], OPT::gmin, converged);
 111   if (!converged && OPT::itl[OPT::SSTEP] > 0) {
 112     int save_itermin = OPT::itermin;
 113     OPT::itermin = 0;
 114     double save_gmin = OPT::gmin;
 115     OPT::gmin = .9;
 116     while (_sim->_iter[iPRINTSTEP] < OPT::itl[OPT::SSTEP] && OPT::gmin > save_gmin) {
 117       //CARD_LIST::card_list.precalc();
 118       _sim->set_inc_mode_no();
 119       trace2("again... ", _sim->_iter[iSTEP], OPT::gmin);
 120       solve(itl, trace);
 121       if (!converged) {
 122         trace2("fail", _sim->_iter[iSTEP], OPT::gmin);
 123         if (OPT::gmin*OPT::shortckt > 1){
 124           error(bDANGER, "ramping failed at gmin=%E\n", OPT::gmin);
 125           break;
 126         }
 127         OPT::gmin *= 3.5;
 128       }else{
 129         trace2("success", _sim->_iter[iSTEP], OPT::gmin);
 130         OPT::gmin /= 4;
 131       }
 132     }
 133     OPT::itermin = save_itermin;
 134     OPT::gmin = save_gmin;
 135     //CARD_LIST::card_list.precalc();
 136     solve(itl, trace);
 137     if (!converged) {
 138       trace2("final fail", _sim->_iter[iSTEP], OPT::gmin);
 139     }else{
 140       trace2("final success", _sim->_iter[iSTEP], OPT::gmin);
 141     }
 142   }else{
 143   }
 144   return converged;
 145 }
 146 /*--------------------------------------------------------------------------*/
 147 /* finish_building_evalq
 148  * This function scans the circuit to queue for evaluation anything
 149  * that is relevant that the devices missed themselves.
 150  * For now, it scans the whole circuit, but this will change.
 151  * This is slow, but the result is still faster than not having a queue.
 152  * The plan is .. every node knows whether it needs eval or not, and
 153  * only those nodes needing eval will be scanned.
 154  * Its purpose is to catch nodes that wake up after being dormant
 155  */
 156 void SIM::finish_building_evalq(void)
 157 {
 158   ::status.queue.start();
 159   CARD_LIST::card_list.tr_queue_eval();
 160   ::status.queue.stop();
 161 }
 162 /*--------------------------------------------------------------------------*/
 163 void SIM::advance_time(void)
 164 {
 165   static double last_iter_time;
 166   trace2("SIM::advance_time()", _sim->_time0, last_iter_time);
 167   ::status.advance.start();
 168   if (_sim->_time0 > 0) {
 169     _sim->_nstat[0].set_discont(disNONE);
 170     if (_sim->_time0 > last_iter_time) {        /* moving forward */
 171       notstd::copy_n(_sim->_v0, _sim->_total_nodes+1, _sim->_vt1);
 172       CARD_LIST::card_list.tr_advance();
 173     }else{
 174       /* moving backward */
 175       /* don't save voltages.  They're wrong! */
 176       /* instead, restore a clean start for iteration */
 177       notstd::copy_n(_sim->_vt1, _sim->_total_nodes+1, _sim->_v0);
 178       CARD_LIST::card_list.tr_regress();
 179     }
 180   }else{
 181     trace1("SIM::advance_time dc_adv.", _sim->_time0);
 182     CARD_LIST::card_list.dc_advance();
 183   }
 184   last_iter_time = _sim->_time0;
 185   ::status.advance.stop();
 186 }
 187 /* last_iter_time is initially 0 by C definition.
 188  * On subsequent runs it will start with an arbitrary positive value.
 189  * _sim->_time0 starts at either 0 or the ending time of the last run.
 190  * In either case, (time0 > last_iter_time) is false on the first step.
 191  * This correctly results in "don't save voltages..."
 192  */
 193 /*--------------------------------------------------------------------------*/
 194 void SIM::set_flags()
 195 {
 196   _sim->_limiting = false;
 197   _sim->_fulldamp = false;
 198
 199   if (OPT::incmode == false) {
 200     _sim->set_inc_mode_no();
 201   }else if (_sim->inc_mode_is_bad()) {
 202     _sim->set_inc_mode_no();
 203   }else if (_sim->is_iteration_number(OPT::itl[OPT::TRLOW])) {
 204     _sim->set_inc_mode_no();
 205   }else if (_sim->is_iteration_number(0)) {
 206     // leave it as is
 207   }else{
 208     _sim->set_inc_mode_yes();
 209   }
 210
 211   _sim->_bypass_ok =
 212     (is_step_rejected()  ||  _sim->_damp < OPT::dampmax*OPT::dampmax)
 213     ? false : bool(OPT::bypass);
 214 }
 215 /*--------------------------------------------------------------------------*/
 216 void SIM::clear_arrays(void)
 217 {
 218   trace0("SIM::clear_arrays");
 219   if (!_sim->is_inc_mode()) {                   /* Clear working array */
 220     _sim->_aa.zero();
 221     _sim->_aa.dezero(OPT::gmin);                /* gmin fudge */
 222     trace2("SIM::clear_arrays ", _sim->_aa.size(), hp(_sim->_i) );
 223     assert(_sim->_i);
 224     assert(_sim->_aa.size()<=_sim->_total_nodes);
 225     std::fill_n(_sim->_i, _sim->_aa.size()+1, 0);
 226   }
 227   trace0("loadq clear");
 228   _sim->_loadq.clear();
 229 }
 230 /*--------------------------------------------------------------------------*/
 231 void SIM::evaluate_models()
 232 {
 233   ::status.evaluate.start();
 234   _sim->_nstat[0].set_discont(disNONE);
 235   if (OPT::bypass) {
 236     converged = true;
 237     std::swap(_sim->_evalq, _sim->_evalq_uc);
 238     while (!_sim->_evalq->empty()) {
 239       converged &= _sim->_evalq->front()->do_tr();
 240       _sim->_evalq->pop_front();
 241     }
 242   }else{
 243     _sim->_evalq_uc->clear();
 244     converged = CARD_LIST::card_list.do_tr();
 245   }
 246   while (!_sim->_late_evalq.empty()) { //BUG// encapsulation violation
 247     converged &= _sim->_late_evalq.front()->do_tr_last();
 248     _sim->_late_evalq.pop_front();
 249   }
 250   ::status.evaluate.stop();
 251 }
 252 /*--------------------------------------------------------------------------*/
 253 void SIM::set_damp()
 254 {
 255   if (_sim->is_second_iteration() && !converged && OPT::dampstrategy&dsINIT) {
 256     _sim->_damp = OPT::dampmin;
 257   }else if (_sim->is_first_iteration()  ||  converged) {
 258     _sim->_damp = OPT::dampmax;
 259   }else if (_sim->_fulldamp) {
 260     _sim->_damp = OPT::dampmin;
 261   }else{
 262     _sim->_damp = OPT::dampmax;
 263   }
 264   trace1("SIM::set_damp", _sim->_damp);
 265 }
 266 /*--------------------------------------------------------------------------*/
 267 void SIM::load_matrix()
 268 {
 269   ::status.load.start();
 270   if (OPT::traceload && _sim->is_inc_mode()) {
 271     trace0("SIM::load_matrix loading some");
 272     while (!_sim->_loadq.empty()) {
 273       _sim->_loadq.back()->tr_load();
 274       _sim->_loadq.pop_back();
 275     }
 276   }else{
 277     trace0("SIM::load_matrix loading all :( ");
 278     _sim->_loadq.clear();
 279     CARD_LIST::card_list.tr_load();
 280   }
 281   ::status.load.stop();
 282 }
 283 /*--------------------------------------------------------------------------*/
 284 #ifndef NDEBUG
 285 static double* v0;
 286 unsigned init_vectors;
 287 #endif
 288 void SIM::solve_equations(TRACE trace)
 289 {
 290 #ifndef NDEBUG
 291   if (init_vectors!=_sim->_total_nodes+1) {
 292     delete[] v0;
 293     v0 = new double[_sim->_total_nodes+1];
 294     init_vectors = _sim->_total_nodes+1;
 295   }
 296 #endif
 297
 298   trace1("SIM::solve_equations", trace);
 299   if(trace & tMATRIX /* && printhere */) {
 300     _out << "\n--- aa -------\n" <<  _sim->_aa << "\n--- i ----\n";
 301     for (unsigned i=0; i<_sim->_total_nodes; ++i){
 302       _out << _sim->_i[i+1] << ", ";
 303     }
 304     _out << "\n-----\n";
 305   }
 306
 307   ::status.lud.start();
 308   _sim->_lu.lu_decomp(_sim->_aa, bool(OPT::lubypass && _sim->is_inc_mode()));
 309   ::status.lud.stop();
 310
 311 #ifndef NDEBUG
 312   notstd::copy_n(_sim->_v0, _sim->_total_nodes+1, v0);
 313   std::string exnumbers("");
 314 #endif
 315   ::status.back.start();
 316   _sim->_lu.fbsub(_sim->_v0, _sim->_i, _sim->_v0);
 317   ::status.back.stop();
 318
 319   _sim->_dxm = 0;
 320   for(unsigned i=0;i<=_sim->_total_nodes; ++i){
 321     if (!is_number(_sim->_v0[i])) {
 322       assert(i); // cannot happen in gnd
 323 #ifdef NDEBUG
 324       unreachable();
 325 #else
 326       _sim->_v0[i] = 0;
 327       exnumbers += " " + to_string(i);
 328     }
 329     _sim->_dxm = std::max(_sim->_dxm, fabs(_sim->_v0[i] - v0[i]));
 330   }
 331   {
 332     if(exnumbers!=""){
 333       throw(Exception("cannot solve linear equation, nan in positions" + exnumbers));
 334 #endif
 335     }
 336   }
 337
 338   if (_sim->_nstat) {
 339     // mixed mode
 340     for (unsigned ii = _sim->_lu.size(); ii >= 1; --ii) {
 341 //      _sim->_nstat[ii].set_discont(disNONE);
 342       _sim->_nstat[ii].set_a_iter();
 343 //      assert(!_sim->_nstat[ii].discont());
 344     }
 345   }else{ untested();
 346     // pure analog
 347   }
 348 }
 349 /*--------------------------------------------------------------------------*/
 350 /*--------------------------------------------------------------------------*/
 351 // vim:ts=8:sw=2:noet: