lib/e_storag.cc

   1 /*$Id: e_storag.cc,v 26.132 2009/11/24 04:26:37 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  * Base class for storage elements of a circuit
  23  */
  24 //testing=obsolete,script 2006.06.14
  25 #include "e_storag.h"
  26 /*--------------------------------------------------------------------------*/
  27 /* table for selecting local integraton method
  28  * Determines which one wins in a conflict.
  29  * "only" wins over non-only.  local (_method_u) wins over opt.
  30  */
  31 //                     OPT::method    _method_u
  32 METHOD STORAGE::method_select[meNUM_METHODS][meNUM_METHODS] = {
  33   /*vv OPT vv*/
  34   //local>>>EULER,EULERONLY,TRAP,TRAPONLY,GEAR2,GEAR2ONLY,TRAPGEAR,TRAPEULER
  35   /*meUNKNOWN*/
  36   {mTRAPGEAR,mEULER,mEULER,mTRAP, mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER},
  37   /*meEULER*/
  38   {mEULER,   mEULER,mEULER,mTRAP, mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER},
  39   /*meEULERONLY*/
  40   {mEULER,   mEULER,mEULER,mEULER,mTRAP,mEULER,mGEAR,mEULER,   mEULER},
  41   /*meTRAP*/
  42   {mTRAP,    mEULER,mEULER,mTRAP, mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER},
  43   /*meTRAPONLY*/
  44   {mTRAP,    mTRAP, mEULER,mTRAP, mTRAP,mTRAP, mGEAR,mTRAP,    mTRAP},
  45   /*meGEAR*/
  46   {mGEAR,    mEULER,mEULER,mTRAP, mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER},
  47   /*meGEAR2ONLY*/
  48   {mGEAR,    mGEAR, mEULER,mGEAR, mTRAP,mGEAR, mGEAR,mGEAR,    mGEAR},
  49   /*meTRAPGEAR*/
  50   {mTRAPGEAR,mEULER,mEULER,mTRAP, mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER},
  51   /*meTRAPEULER*/
  52   {mTRAPEULER,mEULER,mEULER,mTRAP,mTRAP,mGEAR, mGEAR,mTRAPGEAR,mTRAPEULER}
  53 };
  54 /*--------------------------------------------------------------------------*/
  55 void STORAGE::precalc_last()
  56 {
  57   ELEMENT::precalc_last();
  58
  59   set_converged();
  60   assert(!is_constant()); /* because of integration */
  61
  62   _method_a = method_select[OPT::method][_method_u];
  63   //assert(_loss0 == 0.);
  64   //assert(_loss1 == 0.);
  65   /* m0 and acg are frequency/time dependent and cannot be set here.
  66    * If this is a coupled inductor, there is a subckt, which is expanded
  67    * by the mutual pseudo-element.
  68    * Assigning the values here becomes unnecessary, but harmless.
  69    */
  70 }
  71 /*--------------------------------------------------------------------------*/
  72 void STORAGE::tr_begin()
  73 {
  74   ELEMENT::tr_begin();
  75   _method_a = method_select[OPT::method][_method_u];
  76   for (int i = 0;  i < OPT::_keep_time_steps;  ++i) {
  77     _i[i] = FPOLY1(0., 0., 0.);
  78   }
  79   _m1 = _m0 = CPOLY1(0., 0., 0.);
  80 }
  81 /*--------------------------------------------------------------------------*/
  82 void STORAGE::tr_restore()
  83 {
  84   ELEMENT::tr_restore();
  85   _method_a = method_select[OPT::method][_method_u];
  86 }
  87 /*--------------------------------------------------------------------------*/
  88 void STORAGE::dc_advance()
  89 {
  90   ELEMENT::dc_advance();
  91
  92   for (int i = 1;  i < OPT::_keep_time_steps;  ++i) {
  93     _i[i] = _i[0];
  94   }
  95 }
  96 /*--------------------------------------------------------------------------*/
  97 void STORAGE::tr_advance()
  98 {
  99   ELEMENT::tr_advance();
 100
 101   for (int i=OPT::_keep_time_steps-1; i>0; --i) {
 102     _i[i] = _i[i-1];
 103   }
 104 }
 105 /*--------------------------------------------------------------------------*/
 106 /* tr_needs_eval: check to see if this device needs to be evaluated
 107  * this works, and saves significant time
 108  * but possible errors.
 109  * Errors do not seem significant, but I can't tell without more data.
 110  * To disable:  option nolcbypass
 111  */
 112 bool STORAGE::tr_needs_eval()const
 113 {
 114   //assert(!is_q_for_eval());
 115   return (!OPT::lcbypass
 116           || !converged()
 117           || _sim->is_advance_or_first_iteration()
 118           || !conchk(_y[0].x, tr_input(), OPT::abstol)
 119           || _sim->uic_now());
 120 }
 121 /*--------------------------------------------------------------------------*/
 122 /* differentiate: this is what Spice calls "integrate".
 123  * returns an approximation of the time derivative of _q,
 124  * where _q is an array of states .. charge for capacitors, flux for inductors.
 125  * return value is current for capacitors, volts for inductors.
 126  */
 127 FPOLY1 differentiate(const FPOLY1* q, const FPOLY1* i, double* time, METHOD method)
 128 {
 129   if (CKT_BASE::_sim->analysis_is_static()) {
 130     assert(time[0] == 0.);
 131     return FPOLY1(q[0].x, 0., 0.);
 132   }else if (CKT_BASE::_sim->analysis_is_restore()) {
 133     /* leave it alone to restart from a previous solution */
 134     /* it goes this way to continue a transient analysis */
 135     assert(time[0] > 0);
 136     return i[0];
 137   }else{
 138     assert(CKT_BASE::_sim->analysis_is_tran_dynamic());
 139     if (time[1] == 0) {
 140       method = mEULER;  // Bogus current in previous step.  Force Euler.
 141     }else{
 142     }
 143     double dt = time[0] - time[1];
 144     assert(dt > 0.);
 145     switch (method) {
 146     case mTRAPGEAR:
 147       incomplete();
 148     case mGEAR:
 149       assert(OPT::_keep_time_steps >= 3);
 150       return FPOLY1(q[0].x,
 151                     (3./2.) * (q[0].f0 - q[1].f0) / dt
 152                     - (1./2.) * (q[1].f0 - q[2].f0) / (time[1] - time[2]),
 153                     q[0].f1 * (3./2.) / dt);
 154     case mTRAPEULER:
 155       incomplete();
 156     case mEULER:
 157       assert(OPT::_keep_time_steps >= 2);
 158       return FPOLY1(q[0].x,
 159                     (q[0].f0 - q[1].f0) / dt,
 160                     q[0].f1 / dt);
 161     case mTRAP:
 162       assert(OPT::_keep_time_steps >= 2);
 163       return FPOLY1(q[0].x,
 164                     2 * (q[0].f0 - q[1].f0) / dt - i[1].f0,
 165                     q[0].f1 * 2 / dt);
 166     }
 167     unreachable();
 168     return FPOLY1();
 169   }
 170   unreachable();
 171 }
 172 /*--------------------------------------------------------------------------*/
 173 double STORAGE::tr_c_to_g(double c, double g)const
 174 {
 175   if (_sim->analysis_is_static()) {
 176     assert(_time[0] == 0.);
 177     return 0.;
 178   }else if (_sim->analysis_is_restore()) {itested();
 179     assert(_time[0] > 0);
 180     return g;
 181     // no change, fake
 182   }else{
 183     assert(_sim->analysis_is_tran_dynamic());
 184     METHOD method;
 185     if (_time[1] == 0) {
 186       method = mEULER; // Bogus current in previous step.  Force Euler.
 187     }else{
 188       method = _method_a;
 189     }
 190     g = c / _dt;
 191     switch (method) {
 192     case mTRAPGEAR: incomplete();
 193     case mGEAR:  g *= 3./2.;    break;
 194     case mTRAPEULER: incomplete();
 195     case mEULER: /* g *= 1 */   break;
 196     case mTRAP:  g *= 2;        break;
 197     }
 198     return g;
 199   }
 200 }
 201 /*--------------------------------------------------------------------------*/
 202 TIME_PAIR STORAGE::tr_review()
 203 {
 204   COMPONENT::tr_review(); // skip ELEMENT
 205   if (_method_a == mEULER) {
 206     // Backward Euler, no step control, take it as it comes
 207   }else{
 208     double timestep = tr_review_trunc_error(_y);
 209     _time_by.min_error_estimate(tr_review_check_and_convert(timestep));
 210   }
 211   return _time_by;
 212 }
 213 /*--------------------------------------------------------------------------*/
 214 double STORAGE::tr_probe_num(const std::string& x)const
 215 {
 216   if (Umatch(x, "method ")) {untested();
 217     return static_cast<double>(_method_a);
 218   }else{
 219     return ELEMENT::tr_probe_num(x);
 220   }
 221 }
 222 /*--------------------------------------------------------------------------*/
 223 /*--------------------------------------------------------------------------*/
 224 // vim:ts=8:sw=2:noet: