Qucs-core  0.0.19
e_trsolver.cpp
Go to the documentation of this file.
00001 /*
00002  * e_trsolver.cpp - external transient solver interface class implementation
00003  *
00004  * Copyright (C) 2004, 2005, 2006, 2007, 2009 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00034 #if HAVE_CONFIG_H
00035 # include <config.h>
00036 #endif
00037 
00038 #include <stdio.h>
00039 #include <string.h>
00040 #include <cmath>
00041 #include <float.h>
00042 
00043 #include "compat.h"
00044 #include "object.h"
00045 #include "logging.h"
00046 #include "complex.h"
00047 #include "circuit.h"
00048 #include "sweep.h"
00049 #include "net.h"
00050 #include "netdefs.h"
00051 #include "analysis.h"
00052 #include "nasolver.h"
00053 #include "history.h"
00054 #include "transient.h"
00055 #include "exception.h"
00056 #include "exceptionstack.h"
00057 #include "environment.h"
00058 #include "e_trsolver.h"
00059 #include "component_id.h"
00060 #include "ecvs.h"
00061 
00062 #define STEPDEBUG   0 // set to zero for release
00063 #define BREAKPOINTS 0 // exact breakpoint calculation
00064 
00065 #ifndef dState
00066 #define dState 0 // delta T state
00067 #endif
00068 
00069 #ifndef sState
00070 #define sState 1 // solution state
00071 #endif
00072 
00073 // Macro for the n-th state of the solution vector history.
00074 #ifndef SOL
00075 #define SOL(state) (solution[(int) getState (sState, (state))])
00076 #endif
00077 
00078 namespace qucs {
00079 
00080 using namespace transient;
00081 
00082 // Constructor creates an unnamed instance of the trsolver class.
00083 e_trsolver::e_trsolver ()
00084     : trsolver ()
00085 {
00086     //swp = NULL;
00087     type = ANALYSIS_E_TRANSIENT;
00088     //setDescription ("m-transient");
00089     //for (int i = 0; i < 8; i++) solution[i] = NULL;
00090     //tHistory = NULL;
00091     //relaxTSR = false;
00092     //initialDC = true;
00093 
00094     // initialise the message function pointer to
00095     // point to the PrintWarningMsg function
00096     messagefcn = &logprint;
00097 #if DEBUG
00098     //loginit ();
00099     // produce an actual log file
00100 //    file_error = file_status = fopen("e_trsolver.log", "w+");
00101 //    precinit ();
00102 //    ::srand (::time (NULL));
00103 #endif // DEBUG
00104 }
00105 
00106 // Constructor creates a named instance of the e_trsolver class.
00107 e_trsolver::e_trsolver (char * n)
00108     : trsolver (n)
00109 {
00110     //swp = NULL;
00111     type = ANALYSIS_E_TRANSIENT;
00112     //setDescription ("m-transient");
00113     //for (int i = 0; i < 8; i++) solution[i] = NULL;
00114     //tHistory = NULL;
00115     //relaxTSR = false;
00116     //initialDC = true;
00117 
00118     // initialise the message function pointer to
00119     // point to the PrintWarningMsg function
00120     messagefcn = &logprint;
00121 }
00122 
00123 // Destructor deletes the e_trsolver class object.
00124 e_trsolver::~e_trsolver ()
00125 {
00126 
00127     solve_post ();
00128     if (progress) logprogressclear (40);
00129 
00130     deinitTR ();
00131 
00132     delete swp;
00133 
00134     for (int i = 0; i < 8; i++)
00135     {
00136         if (solution[i] != NULL)
00137         {
00138             delete solution[i];
00139         }
00140         if (lastsolution[i] != NULL)
00141         {
00142             delete lastsolution[i];
00143         }
00144     }
00145     delete tHistory;
00146 
00147 }
00148 
00149 /* The copy constructor creates a new instance of the e_trsolver class
00150    based on the given e_trsolver object. */
00151 e_trsolver::e_trsolver (e_trsolver & o)
00152     : trsolver (o)
00153 {
00154     swp = o.swp ? new sweep (*o.swp) : NULL;
00155     for (int i = 0; i < 8; i++) solution[i] = NULL;
00156     tHistory = o.tHistory ? new history (*o.tHistory) : NULL;
00157     relaxTSR = o.relaxTSR;
00158     initialDC = o.initialDC;
00159 }
00160 
00161 void e_trsolver::debug()
00162 {
00163     circuit * root = subnet->getRoot ();
00164 
00165     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00166     {
00167         messagefcn (0, c->getName() );
00168 
00169         if (!c->getSubcircuit ().empty()) {
00170           printf ("subcircuit Name %s\n", c->getSubcircuit ().c_str());
00171         }
00172     }
00173 
00174     //netlist_list();
00175 }
00176 
00177 
00178 /* This is the initialisation function for the slaved transient
00179    netlist solver. It prepares the circuit for simulation. */
00180 int e_trsolver::init (nr_double_t start, nr_double_t firstdelta, int mode)
00181 {
00182     // run the equation solver for the netlist
00183     this->getEnv()->runSolver();
00184 
00185     int error = 0;
00186     const char * const solver = getPropertyString ("Solver");
00187     relaxTSR = !strcmp (getPropertyString ("relaxTSR"), "yes") ? true : false;
00188     initialDC = !strcmp (getPropertyString ("initialDC"), "yes") ? true : false;
00189     // fetch simulation properties
00190     MaxIterations = getPropertyInteger ("MaxIter");
00191     reltol = getPropertyDouble ("reltol");
00192     abstol = getPropertyDouble ("abstol");
00193     vntol = getPropertyDouble ("vntol");
00194 
00195     runs++;
00196     saveCurrent = current = 0;
00197     stepDelta = -1;
00198     converged = 0;
00199     fixpoint = 0;
00200     lastsynctime = 0.0;
00201     statRejected = statSteps = statIterations = statConvergence = 0;
00202 
00203     // Choose a solver.
00204     if (!strcmp (solver, "CroutLU"))
00205         eqnAlgo = ALGO_LU_DECOMPOSITION;
00206     else if (!strcmp (solver, "DoolittleLU"))
00207         eqnAlgo = ALGO_LU_DECOMPOSITION_DOOLITTLE;
00208     else if (!strcmp (solver, "HouseholderQR"))
00209         eqnAlgo = ALGO_QR_DECOMPOSITION;
00210     else if (!strcmp (solver, "HouseholderLQ"))
00211         eqnAlgo = ALGO_QR_DECOMPOSITION_LS;
00212     else if (!strcmp (solver, "GolubSVD"))
00213         eqnAlgo = ALGO_SV_DECOMPOSITION;
00214 
00215     // Perform initial DC analysis.
00216     if (initialDC)
00217     {
00218         error = dcAnalysis ();
00219         if (error)
00220             return -1;
00221     }
00222 
00223     // Initialize transient analysis.
00224     setDescription ("transient");
00225     initETR (start, firstdelta, mode);
00226     setCalculation ((calculate_func_t) &calcTR);
00227     solve_pre ();
00228 
00229     // Recall the DC solution.
00230     recallSolution ();
00231 
00232     // Apply the nodesets and adjust previous solutions.
00233     applyNodeset (false);
00234     fillSolution (x);
00235     fillLastSolution (x);
00236 
00237     // Tell integrators to be initialized.
00238     setMode (MODE_INIT);
00239 
00240     // reset the circuit status to not running so the histories
00241     // etc will be reinitialised on the first time step
00242     running = 0;
00243     rejected = 0;
00244     if (mode == ETR_MODE_ASYNC)
00245     {
00246         delta /= 10;
00247 
00248     }
00249     else if (mode == ETR_MODE_SYNC)
00250     {
00251         //lastsynctime = start - delta;
00252     }
00253     else
00254     {
00255         qucs::exception * e = new qucs::exception (EXCEPTION_UNKNOWN_ETR_MODE);
00256         e->setText ("Unknown ETR mode.");
00257         throw_exception (e);
00258         return -2;
00259     }
00260     fillState (dState, delta);
00261     adjustOrder (1);
00262 
00263     storeHistoryAges ();
00264 
00265     return 0;
00266 
00267 }
00268 
00269 /* Stores all the initial history lengths requested by the circuit
00270    elements for later use (to make sure we don't set the histories
00271    to be less than these initial requested values) */
00272 void e_trsolver::storeHistoryAges (void)
00273 {
00274     circuit * root = subnet->getRoot ();
00275     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00276     {
00277         // get the a
00278         if (c->hasHistory ())
00279         {
00280             initialhistages.push_back (c->getHistoryAge ());
00281         }
00282     }
00283 }
00284 
00285 void e_trsolver::fillLastSolution (tvector<nr_double_t> * s)
00286 {
00287     for (int i = 0; i < 8; i++) * lastsolution[(int) getState (sState, (i))] = *s;
00288 }
00289 
00290 /* Goes through the list of circuit objects and runs its initTR()
00291    function. */
00292 void e_trsolver::initETR (nr_double_t start, nr_double_t firstdelta, int mode)
00293 {
00294     const char * const IMethod = getPropertyString ("IntegrationMethod");
00295     //nr_double_t start = getPropertyDouble ("Start");
00296     //nr_double_t stop = start + firstdelta;
00297     //nr_double_t points = 1.0;
00298 
00299     // fetch corrector integration method and determine predicor method
00300     corrMaxOrder = getPropertyInteger ("Order");
00301     corrType = CMethod = correctorType (IMethod, corrMaxOrder);
00302     predType = PMethod = predictorType (CMethod, corrMaxOrder, predMaxOrder);
00303     corrOrder = corrMaxOrder;
00304     predOrder = predMaxOrder;
00305 
00306     // initialize step values
00307     if (mode == ETR_MODE_ASYNC){
00308         delta = getPropertyDouble ("InitialStep");
00309         deltaMin = getPropertyDouble ("MinStep");
00310         deltaMax = getPropertyDouble ("MaxStep");
00311         if (deltaMax == 0.0)
00312             deltaMax = firstdelta; // MIN ((stop - start) / (points - 1), stop / 200);
00313         if (deltaMin == 0.0)
00314             deltaMin = NR_TINY * 10 * deltaMax;
00315         if (delta == 0.0)
00316             delta = firstdelta; // MIN (stop / 200, deltaMax) / 10;
00317         if (delta < deltaMin) delta = deltaMin;
00318         if (delta > deltaMax) delta = deltaMax;
00319     }
00320     else if (mode == ETR_MODE_SYNC)
00321     {
00322         delta = firstdelta;
00323         deltaMin = NR_TINY * 10;
00324         deltaMax =  std::numeric_limits<nr_double_t>::max() / 10;
00325     }
00326 
00327     // initialize step history
00328     setStates (2);
00329     initStates ();
00330     // initialise the history of states, setting them all to 'delta'
00331     fillState (dState, delta);
00332 
00333     // copy the initialised states to the 'deltas' array
00334     saveState (dState, deltas);
00335     // copy the deltas to all the circuits
00336     setDelta ();
00337     // set the initial corrector and predictor coefficients
00338     calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
00339     calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
00340 
00341     // initialize history of solution vectors (solutions)
00342     for (int i = 0; i < 8; i++)
00343     {
00344         // solution contains the last sets of node voltages and branch
00345         // currents at each of the last 8 'deltas'.
00346         // Note for convenience the definition:
00347         //   #define SOL(state) (solution[(int) getState (sState, (state))])
00348         // is provided and used elsewhere to update the solutions
00349         solution[i] = new tvector<nr_double_t>;
00350         setState (sState, (nr_double_t) i, i);
00351         lastsolution[i] = new tvector<nr_double_t>;
00352     }
00353 
00354     // Initialise history tracking for asynchronous solvers
00355     // See acceptstep_async and rejectstep_async for more
00356     // information
00357     lastasynctime = start;
00358     saveState (dState, lastdeltas);
00359     lastdelta = delta;
00360 
00361     // tell circuit elements about the transient analysis
00362     circuit *c, * root = subnet->getRoot ();
00363     for (c = root; c != NULL; c = (circuit *) c->getNext ())
00364         initCircuitTR (c);
00365     // also initialize the created circuit elements
00366     for (c = root; c != NULL; c = (circuit *) c->getPrev ())
00367         initCircuitTR (c);
00368 }
00369 
00370 void e_trsolver::printx()
00371 {
00372     char buf [1024];
00373 
00374     for (int r = 0; r < x->size(); r++) {
00375         buf[0] = '\0';
00376         //sprintf (buf, "%+.2e%+.2ei", (double) real (x->get (r)), (double) imag (x->get (r)));
00377 
00378         if (r == 2)
00379         {
00380             // print just the currents
00381 //            SOL (0)->get->(r)
00382 //            solution[0]->get(r);
00383             sprintf (buf, "%f\t%18.17f\t%6.5f %6.5f %6.5f %6.5f %6.5f %6.5f %6.5f %6.5f",
00384                      current,
00385                      (double) real (x->get (r)),
00386                      solution[0]->get(r) ,
00387                      solution[1]->get(r) ,
00388                      solution[2]->get(r) ,
00389                      solution[3]->get(r) ,
00390                      solution[4]->get(r) ,
00391                      solution[5]->get(r) ,
00392                      solution[6]->get(r) ,
00393                      solution[7]->get(r) );
00394 
00395             messagefcn(0, buf);
00396         }
00397     }
00398 }
00399 
00400 /* synchronous step solver for external ode routine
00401  *
00402  * This function solves the circuit for a single time delta provided
00403  * by an external source. Convergence issues etc. are expected to
00404  * be handled by the external solver, as it is in full control of the
00405  * time stepping.
00406  */
00407 int e_trsolver::stepsolve_sync(nr_double_t synctime)
00408 {
00409 
00410     int error = 0;
00411     convError = 0;
00412 
00413     time = synctime;
00414     // update the interpolation time of any externally controlled
00415     // components which require it.
00416     updateExternalInterpTime(time);
00417 
00418     // copy the externally chosen time step to delta
00419     delta = time - lastsynctime;
00420 
00421     // get the current solution time
00422     //current += delta;
00423 
00424     // updates the integrator coefficients, and updates the array of prev
00425     // 8 deltas with the new delta for this step
00426     updateCoefficients (delta);
00427 
00428     // Run predictor to get a start value for the solution vector for
00429     // the successive iterative corrector process
00430     error += predictor ();
00431 
00432     // restart Newton iteration
00433     restart (); // restart non-linear devices
00434 
00435     // Attempt to solve the circuit with the given delta
00436     try_running () // #defined as:    do {
00437     {
00438         //error += solve_nonlinear_step ();
00439         error += corrector ();
00440     }
00441     catch_exception () // #defined as:   } while (0); if (estack.top ()) switch (estack.top()->getCode ())
00442     {
00443     case EXCEPTION_NO_CONVERGENCE:
00444         pop_exception ();
00445 
00446         // Retry using damped Newton-Raphson.
00447         this->convHelper = CONV_SteepestDescent;
00448         convError = 2;
00449 #if DEBUG
00450         messagefcn (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
00451                   "(no convergence)\n", (double) saveCurrent, (double) delta);
00452 #endif
00453 
00454         try_running () // #defined as:    do {
00455         {
00456 //            error += solve_nonlinear_step ();
00457             error += solve_nonlinear ();
00458         }
00459         catch_exception () // #defined as:   } while (0); if (estack.top ()) switch (estack.top()->getCode ())
00460         {
00461         case EXCEPTION_NO_CONVERGENCE:
00462             pop_exception ();
00463 
00464             // Update statistics.
00465             statRejected++;
00466             statConvergence++;
00467             rejected++;
00468             converged = 0;
00469             error = 0;
00470 
00471             break;
00472         default:
00473             // Otherwise return.
00474             estack.print ();
00475             error++;
00476             break;
00477         }
00478 
00479         // Update statistics and no more damped Newton-Raphson.
00480 //        statIterations += iterations;
00481 //        if (--convError < 0) this->convHelper = 0;
00482 
00483 
00484         break;
00485     default:
00486         // Otherwise return.
00487         estack.print ();
00488         error++;
00489         break;
00490     }
00491     // if there was an error other than non-convergence, return -1
00492     if (error) return -1;
00493 
00494     // check whether Jacobian matrix is still non-singular
00495     if (!A->isFinite ())
00496     {
00497 //        messagefcn (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
00498 //                  "aborting %s analysis\n", getName (), (double) current,
00499 //                  getDescription ());
00500         return -1;
00501     }
00502 
00503     return 0;
00504 }
00505 
00506 // accept a time step into the solution history
00507 void e_trsolver::acceptstep_sync()
00508 {
00509     statIterations += iterations;
00510     if (--convError < 0) convHelper = 0;
00511 
00512     // Now advance in time or not...
00513     if (running > 1)
00514     {
00515         adjustDelta_sync (current);
00516 //        deltaOld = delta;
00517 //        stepDelta = deltaOld;
00518 //        nextStates ();
00519 //        rejected = 0;
00520         adjustOrder ();
00521     }
00522     else
00523     {
00524         fillStates ();
00525         nextStates ();
00526         rejected = 0;
00527     }
00528 
00529     saveCurrent = current;
00530     current += delta;
00531     running++;
00532     converged++;
00533 
00534     // Tell integrators to be running.
00535     setMode (MODE_NONE);
00536 
00537     // Initialize or update history.
00538     if (running > 1)
00539     {
00540         // update the solution history with the new results
00541         updateHistory (current);
00542     }
00543     else
00544     {
00545         // we have just solved the first transient state
00546         initHistory (current);
00547     }
00548 
00549     // store the current time
00550     lastsynctime = current;
00551 }
00552 
00553 /* This function tries to adapt the current time-step according to the
00554    global truncation error. */
00555 void e_trsolver::adjustDelta_sync (nr_double_t t)
00556 {
00557     deltaOld = delta;
00558 
00559     // makes a new delta based on truncation error
00560 //    delta = checkDelta ();
00561 
00562     if (delta > deltaMax)
00563     {
00564         delta = deltaMax;
00565     }
00566 
00567     if (delta < deltaMin)
00568     {
00569         delta = deltaMin;
00570     }
00571 
00572     // delta correction in order to hit exact breakpoint
00573     int good = 0;
00574 //    if (!relaxTSR)   // relaxed step raster?
00575 //    {
00576 //        if (!statConvergence || converged > 64)   /* Is this a good guess? */
00577 //        {
00578 //            // check next breakpoint
00579 //            if (stepDelta > 0.0)
00580 //            {
00581 //                // restore last valid delta
00582 //                delta = stepDelta;
00583 //                stepDelta = -1.0;
00584 //            }
00585 //            else
00586 //            {
00587 //                if (delta > (t - current) && t > current)
00588 //                {
00589 //                    // save last valid delta and set exact step
00590 //                    stepDelta = deltaOld;
00591 //                    delta = t - current;
00592 //                    good = 1;
00593 //                }
00594 //                else
00595 //                {
00596 //                    stepDelta = -1.0;
00597 //                }
00598 //            }
00599 //            if (delta > deltaMax) delta = deltaMax;
00600 //            if (delta < deltaMin) delta = deltaMin;
00601 //        }
00602 //    }
00603 
00604     stepDelta = -1;
00605     good = 1;
00606 
00607     // usual delta correction
00608     if (delta > 0.9 * deltaOld || good)   // accept current delta
00609     {
00610         nextStates ();
00611         rejected = 0;
00612 #if STEPDEBUG
00613         logprint (LOG_STATUS,
00614                   "DEBUG: delta accepted at t = %.3e, h = %.3e\n",
00615                   (double) current, (double) delta);
00616 #endif
00617     }
00618     else if (deltaOld > delta)   // reject current delta
00619     {
00620         rejected++;
00621         statRejected++;
00622 #if STEPDEBUG
00623         logprint (LOG_STATUS,
00624                   "DEBUG: delta rejected at t = %.3e, h = %.3e\n",
00625                   (double) current, (double) delta);
00626 #endif
00627         if (current > 0) current -= deltaOld;
00628     }
00629     else
00630     {
00631         nextStates ();
00632         rejected = 0;
00633     }
00634 }
00635 
00636 // asynchronous step solver
00637 int e_trsolver::stepsolve_async(nr_double_t steptime)
00638 {
00639     // Start to sweep through time.
00640     int error = 0;
00641     convError = 0;
00642 
00643     time = steptime;
00644     // update the interpolation time of any externally controlled
00645     // components which require it.
00646     updateExternalInterpTime(time);
00647     // make the stored histories for all ircuits that have
00648     // requested them at least as long as the next major time
00649     // step so we can reject the step later if needed and
00650     // restore all the histories to their previous state
00651     updateHistoryAges (time - lastasynctime);
00652 
00653     //delta = (steptime - time) / 10;
00654     //if (progress) logprogressbar (i, swp->getSize (), 40);
00655 #if DEBUG && 0
00656     messagefcn (LOG_STATUS, "NOTIFY: %s: solving netlist for t = %e\n",
00657               getName (), (double) time);
00658 #endif
00659 
00660     do
00661     {
00662 #if STEPDEBUG
00663         if (delta == deltaMin)
00664         {
00665             messagefcn (LOG_ERROR,
00666                       "WARNING: %s: minimum delta h = %.3e at t = %.3e\n",
00667                       getName (), (double) delta, (double) current);
00668         }
00669 #endif
00670         // update the integration coefficients
00671         updateCoefficients (delta);
00672 
00673         // Run predictor to get a start value for the solution vector for
00674         // the successive iterative corrector process
00675         error += predictor ();
00676 
00677         // restart Newton iteration
00678         if (rejected)
00679         {
00680             restart ();      // restart non-linear devices
00681             rejected = 0;
00682         }
00683 
00684         // Run corrector process with appropriate exception handling.
00685         // The corrector iterates through the solutions of the integration
00686         // process until a certain error tolerance has been reached.
00687         try_running () // #defined as:    do {
00688         {
00689             error += corrector ();
00690         }
00691         catch_exception () // #defined as:   } while (0); if (estack.top ()) switch (estack.top()->getCode ())
00692         {
00693         case EXCEPTION_NO_CONVERGENCE:
00694             pop_exception ();
00695 
00696             // Reduce step-size (by half) if failed to converge.
00697             if (current > 0) current -= delta;
00698             delta /= 2;
00699             if (delta <= deltaMin)
00700             {
00701                 delta = deltaMin;
00702                 adjustOrder (1);
00703             }
00704             if (current > 0) current += delta;
00705 
00706             // Update statistics.
00707             statRejected++;
00708             statConvergence++;
00709             rejected++;
00710             converged = 0;
00711             error = 0;
00712 
00713             // Start using damped Newton-Raphson.
00714             convHelper = CONV_SteepestDescent;
00715             convError = 2;
00716 #if DEBUG
00717             messagefcn (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
00718                       "(no convergence)\n", (double) saveCurrent, (double) delta);
00719 #endif
00720             break;
00721         default:
00722             // Otherwise return.
00723             estack.print ();
00724             error++;
00725             break;
00726         }
00727         if (error) return -1;
00728         if (rejected) continue;
00729 
00730         // check whether Jacobian matrix is still non-singular
00731         if (!A->isFinite ())
00732         {
00733             messagefcn (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
00734                       "aborting %s analysis\n", getName (), (double) current,
00735                         getDescription ().c_str());
00736             return -1;
00737         }
00738 
00739         // Update statistics and no more damped Newton-Raphson.
00740         statIterations += iterations;
00741         if (--convError < 0) convHelper = 0;
00742 
00743         // Now advance in time or not...
00744         if (running > 1)
00745         {
00746             adjustDelta (time);
00747             adjustOrder ();
00748         }
00749         else
00750         {
00751             fillStates ();
00752             nextStates ();
00753             rejected = 0;
00754         }
00755 
00756         saveCurrent = current;
00757         current += delta;
00758         running++;
00759         converged++;
00760 
00761         // Tell integrators to be running.
00762         setMode (MODE_NONE);
00763 
00764         // Initialize or update history.
00765         if (running > 1)
00766         {
00767             updateHistory (saveCurrent);
00768         }
00769         else
00770         {
00771             initHistory (saveCurrent);
00772         }
00773     }
00774     while (saveCurrent < time); // Hit a requested time point?
00775 
00776     return 0;
00777 }
00778 
00779 // accept an asynchronous step into the solution history
00780 void e_trsolver::acceptstep_async(void)
00781 {
00782     // copy the solution in case we wish to step back to this
00783     // point later
00784     copySolution (solution, lastsolution);
00785 
00786     // Store the time
00787     lastasynctime = time;
00788 
00789     // Store the deltas history
00790     saveState (dState, lastdeltas);
00791 
00792     // Store the current delta
00793     lastdelta = delta;
00794 }
00795 
00796 // reject the last asynchronous step and restore the history
00797 // of solutions to the last major step
00798 void e_trsolver::rejectstep_async(void)
00799 {
00800     // restore the solution (node voltages and branch currents) from
00801     // the previously stored solution
00802     copySolution (lastsolution, solution);
00803 
00804     // Restore the circuit histories to their previous states
00805     truncateHistory (lastasynctime);
00806 
00807     // Restore the time deltas
00808     inputState (dState, lastdeltas);
00809 
00810     for (int i = 0; i < 8; i++)
00811     {
00812         deltas[i] = lastdeltas[i];
00813     }
00814 
00815     delta = lastdelta;
00816 
00817     // copy the deltas to all the circuit elements
00818     setDelta ();
00819 
00820     // reset the corrector and predictor coefficients
00821     calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
00822     calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
00823 }
00824 
00825 /* Makes a copy of a set of solution vectors */
00826 void e_trsolver::copySolution (tvector<nr_double_t> * src[8], tvector<nr_double_t> * dest[8])
00827 {
00828     for (int i = 0; i < 8; i++)
00829     {
00830         // check sizes are the same
00831         assert (src[i]->size () == dest[i]->size ());
00832         // copy over the data values
00833         for (int j = 0; j < src[i]->size (); j++)
00834         {
00835             dest[i]->set (j, src[i]->get (j));
00836         }
00837     }
00838 }
00839 
00840 void e_trsolver::updateHistoryAges (nr_double_t newage)
00841 {
00842     int i = 0;
00843     circuit * root = subnet->getRoot ();
00844     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00845     {
00846         // set the history length to retain to be at least
00847         // the length of the supplied age
00848         if (c->hasHistory ())
00849         {
00850             c->setHistoryAge (std::max (initialhistages[i], newage));
00851             i++;
00852         }
00853     }
00854 }
00855 
00856 //int e_trsolver::finish()
00857 //{
00858 //    solve_post ();
00859 //
00860 //    if (progress) logprogressclear (40);
00861 //
00862 //    messagefcn (LOG_STATUS, "NOTIFY: %s: average time-step %g, %d rejections\n",
00863 //              getName (), (double) (saveCurrent / statSteps), statRejected);
00864 //    messagefcn (LOG_STATUS, "NOTIFY: %s: average NR-iterations %g, "
00865 //              "%d non-convergences\n", getName (),
00866 //              (double) statIterations / statSteps, statConvergence);
00867 //
00868 //    // cleanup
00869 //    return 0;
00870 //}
00871 
00872 
00873 void e_trsolver::getsolution (double * lastsol)
00874 {
00875     int N = countNodes ();
00876     int M = countVoltageSources ();
00877 
00878     // copy solution
00879     for (int r = 0; r < N + M; r++)
00880     {
00881         lastsol[r]  = real(x->get(r));
00882     }
00883 }
00884 
00885 /* getNodeV attempts to get the voltage of the node with a
00886    given name. returns -1 if the node name was not found */
00887 int e_trsolver::getNodeV (char * label, nr_double_t& nodeV)
00888 {
00889     int r = nlist->getNodeNr (label);
00890 
00891     if (r == -1)
00892     {
00893       return r;
00894     }
00895     else
00896     {
00897       nodeV = x->get(r);
00898       return 0;
00899     }
00900 }
00901 
00902 /* Get the voltage reported by a voltage probe */
00903 int e_trsolver::getVProbeV (char * probename, nr_double_t& probeV)
00904 {
00905     // string to hold the full name of the circuit
00906     std::string fullname;
00907 
00908     // check for NULL name
00909     if (probename)
00910     {
00911         circuit * root = subnet->getRoot ();
00912         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00913         {
00914             if (c->getType () == CIR_VPROBE) {
00915 
00916                 fullname.clear ();
00917 
00918                 // Subcircuit elements top level name is the
00919                 // subcircuit type (the base name of the subcircuit
00920                 // file)
00921                 if (!c->getSubcircuit ().empty())
00922                 {
00923                     fullname.append (c->getSubcircuit ());
00924                     fullname.append (".");
00925                 }
00926 
00927                 // append the user supplied name to search for
00928                 fullname.append (probename);
00929 
00930                 // Check if it is the desired voltage source
00931                 if (strcmp (fullname.c_str(), c->getName ()) == 0)
00932                 {
00933                     // Saves the real and imaginary voltages in the probe to the
00934                     // named variables Vr and Vi
00935                     c->saveOperatingPoints ();
00936                     // We are only interested in the real part for transient
00937                     // analysis
00938                     probeV = c->getOperatingPoint ("Vr");
00939 
00940                     return 0;
00941                 }
00942             }
00943         }
00944     }
00945     return -1;
00946 }
00947 
00948 /* Get the current reported by a current probe */
00949 int e_trsolver::getIProbeI (char * probename, nr_double_t& probeI)
00950 {
00951     // string to hold the full name of the circuit
00952     std::string fullname;
00953 
00954     // check for NULL name
00955     if (probename)
00956     {
00957         circuit * root = subnet->getRoot ();
00958         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00959         {
00960             if (c->getType () == CIR_IPROBE) {
00961 
00962                 fullname.clear ();
00963 
00964                 // Subcircuit elements top level name is the
00965                 // subcircuit type (the base name of the subcircuit
00966                 // file)
00967                 if (!c->getSubcircuit ().empty())
00968                 {
00969                     fullname.append (c->getSubcircuit ());
00970                     fullname.append (".");
00971                 }
00972 
00973                 // append the user supplied name to search for
00974                 fullname.append (probename);
00975 
00976                 // Check if it is the desired voltage source
00977                 if (strcmp (fullname.c_str(), c->getName ()) == 0)
00978                 {
00979                     // Get the current reported by the probe
00980                     probeI = real (x->get (c->getVoltageSource () + getN ()));
00981 
00982                     return 0;
00983                 }
00984             }
00985         }
00986     }
00987     return -1;
00988 }
00989 
00990 int e_trsolver::setECVSVoltage(char * ecvsname, nr_double_t V)
00991 {
00992     // string to hold the full name of the circuit
00993     std::string fullname;
00994 
00995     // check for NULL name
00996     if (ecvsname)
00997     {
00998         circuit * root = subnet->getRoot ();
00999         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01000         {
01001             // examine only ECVS elements
01002             if (c->getType () == CIR_ECVS) {
01003 
01004                 fullname.clear ();
01005 
01006                 // Subcircuit elements top level name is the
01007                 // subcircuit type (the base name of the subcircuit
01008                 // file)
01009                 if (!c->getSubcircuit ().empty())
01010                 {
01011                     fullname.append (c->getSubcircuit ());
01012                     fullname.append (".");
01013                 }
01014 
01015                 // append the user supplied name to search for
01016                 fullname.append (ecvsname);
01017 
01018                 // Check if it is the desired voltage source
01019                 if (strcmp (fullname.c_str(), c->getName ()) == 0)
01020                 {
01021                     // Set the voltage to the desired value
01022                     c->setProperty("U", V);
01023                     return 0;
01024                 }
01025             }
01026         }
01027     }
01028     return -1;
01029 }
01030 
01031 void e_trsolver::updateExternalInterpTime(nr_double_t t)
01032 {
01033     circuit * root = subnet->getRoot ();
01034     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01035     {
01036         // examine only external elements that have interpolation,
01037         // such as ECVS elements
01038         if (c->getType () == CIR_ECVS) {
01039             // Set the voltage to the desired value
01040             c->setProperty ("Tnext", t);
01041             if (tHistory != NULL && tHistory->size () > 0)
01042             {
01043               c->setHistoryAge ( t - tHistory->last () + 0.1 * (t - tHistory->last ()) );
01044             }
01045         }
01046     }
01047 }
01048 
01049 /* The following function removed stored times newer than a specified time
01050    from all the circuit element histories */
01051 void e_trsolver::truncateHistory (nr_double_t t)
01052 {
01053     // truncate all the circuit element histories
01054     circuit * root = subnet->getRoot ();
01055     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
01056     {
01057         if (c->hasHistory ()) c->truncateHistory (t);
01058     }
01059 }
01060 
01061 int e_trsolver::getJacRows()
01062 {
01063     return A->getRows();
01064 }
01065 
01066 int e_trsolver::getJacCols()
01067 {
01068     return A->getCols();
01069 }
01070 
01071 void e_trsolver::getJacData(int r, int c, nr_double_t& data)
01072 {
01073     data = A->get(r,c);
01074 }
01075 
01076 // properties
01077 PROP_REQ [] =
01078 {
01079     PROP_NO_PROP
01080 };
01081 PROP_OPT [] =
01082 {
01083     {
01084         "IntegrationMethod", PROP_STR, { PROP_NO_VAL, "Trapezoidal" },
01085         PROP_RNG_STR4 ("Euler", "Trapezoidal", "Gear", "AdamsMoulton")
01086     },
01087     { "Order", PROP_INT, { 2, PROP_NO_STR }, PROP_RNGII (1, 6) },
01088     { "InitialStep", PROP_REAL, { 1e-9, PROP_NO_STR }, PROP_POS_RANGE },
01089     { "MinStep", PROP_REAL, { 1e-16, PROP_NO_STR }, PROP_POS_RANGE },
01090     { "MaxStep", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
01091     { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
01092     { "abstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
01093     { "vntol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
01094     { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
01095     { "LTEabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
01096     { "LTEreltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
01097     { "LTEfactor", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 16) },
01098     { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
01099     { "Solver", PROP_STR, { PROP_NO_VAL, "CroutLU" }, PROP_RNG_SOL },
01100     { "relaxTSR", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
01101     { "initialDC", PROP_STR, { PROP_NO_VAL, "yes" }, PROP_RNG_YESNO },
01102     PROP_NO_PROP
01103 };
01104 struct define_t e_trsolver::anadef =
01105     { "ETR", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
01106 
01107 
01108 } // namespace qucs