Qucs-core  0.0.19
trsolver.cpp
Go to the documentation of this file.
00001 /*
00002  * trsolver.cpp - transient solver 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 
00025 #if HAVE_CONFIG_H
00026 # include <config.h>
00027 #endif
00028 
00029 #include <stdio.h>
00030 #include <string.h>
00031 #include <float.h>
00032 #include <algorithm>
00033 
00034 #include "compat.h"
00035 #include "object.h"
00036 #include "logging.h"
00037 #include "complex.h"
00038 #include "circuit.h"
00039 #include "sweep.h"
00040 #include "net.h"
00041 #include "netdefs.h"
00042 #include "analysis.h"
00043 #include "nasolver.h"
00044 #include "history.h"
00045 #include "trsolver.h"
00046 #include "transient.h"
00047 #include "exception.h"
00048 #include "exceptionstack.h"
00049 
00050 #define STEPDEBUG   0 // set to zero for release
00051 #define BREAKPOINTS 0 // exact breakpoint calculation
00052 
00053 #define dState 0 // delta T state
00054 #define sState 1 // solution state
00055 
00056 // Macro for the n-th state of the solution vector history.
00057 #define SOL(state) (solution[(int) getState (sState, (state))])
00058 
00059 namespace qucs {
00060 
00061 using namespace transient;
00062 
00063 // Constructor creates an unnamed instance of the trsolver class.
00064 trsolver::trsolver ()
00065     : nasolver<nr_double_t> (), states<nr_double_t> ()
00066 {
00067     swp = NULL;
00068     type = ANALYSIS_TRANSIENT;
00069     setDescription ("transient");
00070     for (int i = 0; i < 8; i++) solution[i] = NULL;
00071     tHistory = NULL;
00072     relaxTSR = false;
00073     initialDC = true;
00074 }
00075 
00076 // Constructor creates a named instance of the trsolver class.
00077 trsolver::trsolver (const std::string &n)
00078     : nasolver<nr_double_t> (n), states<nr_double_t> ()
00079 {
00080     swp = NULL;
00081     type = ANALYSIS_TRANSIENT;
00082     setDescription ("transient");
00083     for (int i = 0; i < 8; i++) solution[i] = NULL;
00084     tHistory = NULL;
00085     relaxTSR = false;
00086     initialDC = true;
00087 }
00088 
00089 // Destructor deletes the trsolver class object.
00090 trsolver::~trsolver ()
00091 {
00092     delete swp;
00093     for (int i = 0; i < 8; i++)
00094     {
00095         if (solution[i] != NULL)
00096         {
00097             delete solution[i];
00098         }
00099     }
00100     delete tHistory;
00101 }
00102 
00103 /* The copy constructor creates a new instance of the trsolver class
00104    based on the given trsolver object. */
00105 trsolver::trsolver (trsolver & o)
00106     : nasolver<nr_double_t> (o), states<nr_double_t> (o)
00107 {
00108     swp = o.swp ? new sweep (*o.swp) : NULL;
00109     for (int i = 0; i < 8; i++) solution[i] = NULL;
00110     tHistory = o.tHistory ? new history (*o.tHistory) : NULL;
00111     relaxTSR = o.relaxTSR;
00112     initialDC = o.initialDC;
00113 }
00114 
00115 // This function creates the time sweep if necessary.
00116 void trsolver::initSteps (void)
00117 {
00118     delete swp;
00119     swp = createSweep ("time");
00120 }
00121 
00122 // Performs the initial DC analysis.
00123 int trsolver::dcAnalysis (void)
00124 {
00125     int error = 0;
00126 
00127     // First calculate a initial state using the non-linear DC analysis.
00128     setDescription ("initial DC");
00129     initDC ();
00130     setCalculation ((calculate_func_t) &calcDC);
00131     solve_pre ();
00132     applyNodeset ();
00133 
00134     // Run the DC solver once.
00135     try_running ()
00136     {
00137         error = solve_nonlinear ();
00138     }
00139     // Appropriate exception handling.
00140     catch_exception ()
00141     {
00142     case EXCEPTION_NO_CONVERGENCE:
00143         pop_exception ();
00144         convHelper = CONV_LineSearch;
00145         logprint (LOG_ERROR, "WARNING: %s: %s analysis failed, using line search "
00146                   "fallback\n", getName (), getDescription ().c_str());
00147         applyNodeset ();
00148         restart ();
00149         error = solve_nonlinear ();
00150         break;
00151     default:
00152         // Otherwise return.
00153         estack.print ();
00154         error++;
00155         break;
00156     }
00157 
00158     // Save the DC solution.
00159     storeSolution ();
00160 
00161     // Cleanup nodal analysis solver.
00162     solve_post ();
00163 
00164     // Really failed to find initial DC solution?
00165     if (error)
00166     {
00167         logprint (LOG_ERROR, "ERROR: %s: %s analysis failed\n",
00168                   getName (), getDescription ().c_str());
00169     }
00170     return error;
00171 }
00172 
00173 /* This is the transient netlist solver.  It prepares the circuit list
00174    for each requested time and solves it then. */
00175 int trsolver::solve (void)
00176 {
00177     nr_double_t time, saveCurrent;
00178     int error = 0, convError = 0;
00179     const char * const solver = getPropertyString ("Solver");
00180     relaxTSR = !strcmp (getPropertyString ("relaxTSR"), "yes") ? true : false;
00181     initialDC = !strcmp (getPropertyString ("initialDC"), "yes") ? true : false;
00182 
00183     runs++;
00184     saveCurrent = current = 0;
00185     stepDelta = -1;
00186     converged = 0;
00187     fixpoint = 0;
00188     statRejected = statSteps = statIterations = statConvergence = 0;
00189 
00190     // Choose a solver.
00191     if (!strcmp (solver, "CroutLU"))
00192         eqnAlgo = ALGO_LU_DECOMPOSITION;
00193     else if (!strcmp (solver, "DoolittleLU"))
00194         eqnAlgo = ALGO_LU_DECOMPOSITION_DOOLITTLE;
00195     else if (!strcmp (solver, "HouseholderQR"))
00196         eqnAlgo = ALGO_QR_DECOMPOSITION;
00197     else if (!strcmp (solver, "HouseholderLQ"))
00198         eqnAlgo = ALGO_QR_DECOMPOSITION_LS;
00199     else if (!strcmp (solver, "GolubSVD"))
00200         eqnAlgo = ALGO_SV_DECOMPOSITION;
00201 
00202     // Perform initial DC analysis.
00203     if (initialDC)
00204     {
00205         error = dcAnalysis ();
00206         if (error)
00207             return -1;
00208     }
00209 
00210     // Initialize transient analysis.
00211     setDescription ("transient");
00212     initTR ();
00213     setCalculation ((calculate_func_t) &calcTR);
00214     solve_pre ();
00215 
00216     // Create time sweep if necessary.
00217     initSteps ();
00218     swp->reset ();
00219 
00220     // Recall the DC solution.
00221     recallSolution ();
00222 
00223     // Apply the nodesets and adjust previous solutions.
00224     applyNodeset (false);
00225     fillSolution (x);
00226 
00227     // Tell integrators to be initialized.
00228     setMode (MODE_INIT);
00229 
00230     int running = 0;
00231     rejected = 0;
00232     delta /= 10;
00233     fillState (dState, delta);
00234     adjustOrder (1);
00235 
00236     // Start to sweep through time.
00237     for (int i = 0; i < swp->getSize (); i++)
00238     {
00239         time = swp->next ();
00240         if (progress) logprogressbar (i, swp->getSize (), 40);
00241 
00242 #if DEBUG && 0
00243         logprint (LOG_STATUS, "NOTIFY: %s: solving netlist for t = %e\n",
00244                   getName (), (double) time);
00245 #endif
00246 
00247         do // while (saveCurrent < time), i.e. until a requested breakpoint is hit
00248         {
00249 #if STEPDEBUG
00250             if (delta == deltaMin)
00251             {
00252                 // the integrator step size has become smaller than the
00253                 // specified allowed minimum, Qucs is unable to solve the circuit
00254                 // while meeting the tolerance conditions
00255                 logprint (LOG_ERROR,
00256                           "WARNING: %s: minimum delta h = %.3e at t = %.3e\n",
00257                           getName (), (double) delta, (double) current);
00258             }
00259 #endif
00260             // updates the integrator coefficients, and updates the array of prev
00261             // 8 deltas with the new delta for this step
00262             updateCoefficients (delta);
00263 
00264             // Run predictor to get a start value for the solution vector for
00265             // the successive iterative corrector process
00266             error += predictor ();
00267 
00268             // restart Newton iteration
00269             if (rejected)
00270             {
00271                 restart ();      // restart non-linear devices
00272                 rejected = 0;
00273             }
00274 
00275             // Run corrector process with appropriate exception handling.
00276             // The corrector iterates through the solutions of the integration
00277             // process until a certain error tolerance has been reached.
00278             try_running () // #defined as:    do {
00279             {
00280                 error += corrector ();
00281             }
00282             catch_exception () // #defined as:   } while (0); if (estack.top ()) switch (estack.top()->getCode ())
00283             {
00284             case EXCEPTION_NO_CONVERGENCE:
00285                 pop_exception ();
00286 
00287                 // step back from the current time value to the previous time
00288                 if (current > 0) current -= delta;
00289                 // Reduce step-size (by half) if failed to converge.
00290                 delta /= 2;
00291                 if (delta <= deltaMin)
00292                 {
00293                     // but do not reduce the step size below a specified minimum
00294                     delta = deltaMin;
00295                     // instead reduce the order of the integration
00296                     adjustOrder (1);
00297                 }
00298                 // step forward to the new current time value
00299                 if (current > 0) current += delta;
00300 
00301                 // Update statistics.
00302                 statRejected++;
00303                 statConvergence++;
00304                 rejected++; // mark the previous step size choice as rejected
00305                 converged = 0;
00306                 error = 0;
00307 
00308                 // Start using damped Newton-Raphson.
00309                 convHelper = CONV_SteepestDescent;
00310                 convError = 2;
00311 #if DEBUG
00312                 logprint (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
00313                           "(no convergence)\n", (double) saveCurrent, (double) delta);
00314 #endif
00315                 break;
00316             default:
00317                 // Otherwise return.
00318                 estack.print ();
00319                 error++;
00320                 break;
00321             }
00322             // return if any errors occured other than convergence failure
00323             if (error) return -1;
00324 
00325             // if the step was rejected, the solution loop is restarted here
00326             if (rejected) continue;
00327 
00328             // check whether Jacobian matrix is still non-singular
00329             if (!A->isFinite ())
00330             {
00331                 logprint (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
00332                           "aborting %s analysis\n", getName (), (double) current,
00333                           getDescription ().c_str());
00334                 return -1;
00335             }
00336 
00337             // Update statistics and no more damped Newton-Raphson.
00338             statIterations += iterations;
00339             if (--convError < 0) convHelper = 0;
00340 
00341             // Now advance in time or not...
00342             if (running > 1)
00343             {
00344                 adjustDelta (time);
00345                 adjustOrder ();
00346             }
00347             else
00348             {
00349                 fillStates ();
00350                 nextStates ();
00351                 rejected = 0;
00352             }
00353 
00354             saveCurrent = current;
00355             current += delta;
00356             running++;
00357             converged++;
00358 
00359             // Tell integrators to be running.
00360             setMode (MODE_NONE);
00361 
00362             // Initialize or update history.
00363             if (running > 1)
00364             {
00365                 updateHistory (saveCurrent);
00366             }
00367             else
00368             {
00369                 initHistory (saveCurrent);
00370             }
00371         }
00372         while (saveCurrent < time); // Hit a requested time point?
00373 
00374         // Save results.
00375 #if STEPDEBUG
00376         logprint (LOG_STATUS, "DEBUG: save point at t = %.3e, h = %.3e\n",
00377                   (double) saveCurrent, (double) delta);
00378 #endif
00379 
00380 #if BREAKPOINTS
00381         saveAllResults (saveCurrent);
00382 #else
00383         saveAllResults (time);
00384 #endif
00385     } // for (int i = 0; i < swp->getSize (); i++)
00386 
00387     solve_post ();
00388     if (progress) logprogressclear (40);
00389     logprint (LOG_STATUS, "NOTIFY: %s: average time-step %g, %d rejections\n",
00390               getName (), (double) (saveCurrent / statSteps), statRejected);
00391     logprint (LOG_STATUS, "NOTIFY: %s: average NR-iterations %g, "
00392               "%d non-convergences\n", getName (),
00393               (double) statIterations / statSteps, statConvergence);
00394 
00395     // cleanup
00396     deinitTR ();
00397     return 0;
00398 }
00399 
00400 // The function initializes the history.
00401 void trsolver::initHistory (nr_double_t t)
00402 {
00403     // initialize time vector
00404     tHistory = new history ();
00405     tHistory->push_back(t);
00406     tHistory->self ();
00407     // initialize circuit histories
00408     nr_double_t age = 0.0;
00409     circuit * root = subnet->getRoot ();
00410     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00411     {
00412         if (c->hasHistory ())
00413         {
00414             c->applyHistory (tHistory);
00415             saveHistory (c);
00416             if (c->getHistoryAge () > age)
00417             {
00418                 age = c->getHistoryAge ();
00419             }
00420         }
00421     }
00422     // set maximum required age for all circuits
00423     tHistory->setAge (age);
00424 }
00425 
00426 /* The following function updates the histories for the circuits which
00427    requested them. */
00428 void trsolver::updateHistory (nr_double_t t)
00429 {
00430     if (t > tHistory->last ())
00431     {
00432         // update time vector
00433         tHistory->push_back (t);
00434         // update circuit histories
00435         circuit * root = subnet->getRoot ();
00436         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00437         {
00438             if (c->hasHistory ()) saveHistory (c);
00439         }
00440         tHistory->drop ();
00441     }
00442 }
00443 
00444 // Stores node voltages and branch currents in the given circuits history.
00445 void trsolver::saveHistory (circuit * c)
00446 {
00447 
00448     int N = countNodes ();
00449     int r, i, s = c->getSize ();
00450 
00451     for (i = 0; i < s; i++)
00452     {
00453         // save node voltages
00454         r = findAssignedNode (c, i);
00455         if (r < 0)
00456             // the node was not found, append a zero to the history
00457             // matching this index
00458             c->appendHistory (i, 0.0);
00459         else
00460             // the node was found, append the voltage value to
00461             // that node's history
00462             c->appendHistory (i, x->get (r));
00463     }
00464 
00465     for (i = 0; i < c->getVoltageSources (); i++)
00466     {
00467         // save branch currents
00468         r = c->getVoltageSource () + i;
00469         c->appendHistory (i + s, x->get (r + N));
00470     }
00471 
00472 }
00473 
00474 /* This function predicts a start value for the solution vector for
00475    the successive iterative corrector process. */
00476 int trsolver::predictor (void)
00477 {
00478     int error = 0;
00479     switch (predType)
00480     {
00481     case INTEGRATOR_GEAR: // explicit GEAR
00482         predictGear ();
00483         break;
00484     case INTEGRATOR_ADAMSBASHFORD: // ADAMS-BASHFORD
00485         predictBashford ();
00486         break;
00487     case INTEGRATOR_EULER: // FORWARD EULER
00488         predictEuler ();
00489         break;
00490     default:
00491         *x = *SOL (1);  // This is too a simple predictor...
00492         break;
00493     }
00494     saveSolution ();
00495     *SOL (0) = *x;
00496     return error;
00497 }
00498 
00499 // Stores the given vector into all the solution vectors.
00500 void trsolver::fillSolution (tvector<nr_double_t> * s)
00501 {
00502     for (int i = 0; i < 8; i++) *SOL (i) = *s;
00503 }
00504 
00505 /* The function predicts the successive solution vector using the
00506    explicit Adams-Bashford integration formula. */
00507 void trsolver::predictBashford (void)
00508 {
00509     int N = countNodes ();
00510     int M = countVoltageSources ();
00511     nr_double_t xn, dd, hn;
00512 
00513     // go through each solution
00514     for (int r = 0; r < N + M; r++)
00515     {
00516         xn = predCoeff[0] * SOL(1)->get (r); // a0 coefficient
00517         for (int o = 1; o <= predOrder; o++)
00518         {
00519             hn = getState (dState, o);         // previous time-step
00520             // divided differences
00521             dd = (SOL(o)->get (r) - SOL(o + 1)->get (r)) / hn;
00522             xn += predCoeff[o] * dd;           // b0, b1, ... coefficients
00523         }
00524         x->set (r, xn);                      // save prediction
00525     }
00526 }
00527 
00528 /* The function predicts the successive solution vector using the
00529    explicit forward Euler integration formula.  Actually this is
00530    Adams-Bashford order 1. */
00531 void trsolver::predictEuler (void)
00532 {
00533     int N = countNodes ();
00534     int M = countVoltageSources ();
00535     nr_double_t xn, dd, hn;
00536 
00537     for (int r = 0; r < N + M; r++)
00538     {
00539         xn = predCoeff[0] * SOL(1)->get (r);
00540         hn = getState (dState, 1);
00541         dd = (SOL(1)->get (r) - SOL(2)->get (r)) / hn;
00542         xn += predCoeff[1] * dd;
00543         x->set (r, xn);
00544     }
00545 }
00546 
00547 /* The function predicts the successive solution vector using the
00548    explicit Gear integration formula. */
00549 void trsolver::predictGear (void)
00550 {
00551     int N = countNodes ();
00552     int M = countVoltageSources ();
00553     nr_double_t xn;
00554 
00555     // go through each solution
00556     for (int r = 0; r < N + M; r++)
00557     {
00558         xn = 0;
00559         for (int o = 0; o <= predOrder; o++)
00560         {
00561             // a0, a1, ... coefficients
00562             xn += predCoeff[o] * SOL(o + 1)->get (r);
00563         }
00564         x->set (r, xn); // save prediction
00565     }
00566 }
00567 
00568 /* The function iterates through the solutions of the integration
00569    process until a certain error tolerance has been reached. */
00570 int trsolver::corrector (void)
00571 {
00572     int error = 0;
00573     error += solve_nonlinear ();
00574     return error;
00575 }
00576 
00577 // The function advances one more time-step.
00578 void trsolver::nextStates (void)
00579 {
00580     circuit * root = subnet->getRoot ();
00581     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00582     {
00583         // for each circuit get the next state
00584         c->nextState ();
00585     }
00586 
00587     *SOL (0) = *x; // save current solution
00588     nextState ();
00589     statSteps++;
00590 }
00591 
00592 /* This function stores the current state of each circuit into all
00593    other states as well.  It is useful for higher order integration
00594    methods in order to initialize the states after the initial
00595    transient solution. */
00596 void trsolver::fillStates (void)
00597 {
00598     circuit * root = subnet->getRoot ();
00599     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00600     {
00601         for (int s = 0; s < c->getStates (); s++)
00602             c->fillState (s, c->getState (s));
00603     }
00604 }
00605 
00606 // The function modifies the circuit lists integrator mode.
00607 void trsolver::setMode (int state)
00608 {
00609     circuit * root = subnet->getRoot ();
00610     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00611         c->setMode (state);
00612 }
00613 
00614 // The function passes the time delta array to the circuit list.
00615 void trsolver::setDelta (void)
00616 {
00617     circuit * root = subnet->getRoot ();
00618     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00619         c->setDelta (deltas);
00620 }
00621 
00622 /* This function tries to adapt the current time-step according to the
00623    global truncation error. */
00624 void trsolver::adjustDelta (nr_double_t t)
00625 {
00626     deltaOld = delta;
00627     delta = checkDelta ();
00628     if (delta > deltaMax) delta = deltaMax;
00629     if (delta < deltaMin) delta = deltaMin;
00630 
00631     // delta correction in order to hit exact breakpoint
00632     int good = 0;
00633     if (!relaxTSR)   // relaxed step raster?
00634     {
00635         if (!statConvergence || converged > 64)   /* Is this a good guess? */
00636         {
00637             // check next breakpoint
00638             if (stepDelta > 0.0)
00639             {
00640                 // restore last valid delta
00641                 delta = stepDelta;
00642                 stepDelta = -1.0;
00643             }
00644             else
00645             {
00646                 if (delta > (t - current) && t > current)
00647                 {
00648                     // save last valid delta and set exact step
00649                     stepDelta = deltaOld;
00650                     delta = t - current;
00651                     good = 1;
00652                 }
00653                 else
00654                 {
00655                     stepDelta = -1.0;
00656                 }
00657             }
00658             if (delta > deltaMax) delta = deltaMax;
00659             if (delta < deltaMin) delta = deltaMin;
00660         }
00661     }
00662 
00663     // usual delta correction
00664     if (delta > 0.9 * deltaOld || good)   // accept current delta
00665     {
00666         nextStates ();
00667         rejected = 0;
00668 #if STEPDEBUG
00669         logprint (LOG_STATUS,
00670                   "DEBUG: delta accepted at t = %.3e, h = %.3e\n",
00671                   (double) current, (double) delta);
00672 #endif
00673     }
00674     else if (deltaOld > delta)   // reject current delta
00675     {
00676         rejected++;
00677         statRejected++;
00678 #if STEPDEBUG
00679         logprint (LOG_STATUS,
00680                   "DEBUG: delta rejected at t = %.3e, h = %.3e\n",
00681                   (double) current, (double) delta);
00682 #endif
00683         if (current > 0) current -= deltaOld;
00684     }
00685     else
00686     {
00687         nextStates ();
00688         rejected = 0;
00689     }
00690 }
00691 
00692 /* The function can be used to increase the current order of the
00693    integration method or to reduce it. */
00694 void trsolver::adjustOrder (int reduce)
00695 {
00696     if ((corrOrder < corrMaxOrder && !rejected) || reduce)
00697     {
00698         if (reduce)
00699         {
00700             corrOrder = 1;
00701         }
00702         else if (!rejected)
00703         {
00704             corrOrder++;
00705         }
00706 
00707         // adjust type and order of corrector and predictor
00708         corrType = correctorType (CMethod, corrOrder);
00709         predType = predictorType (corrType, corrOrder, predOrder);
00710 
00711         // apply new corrector method and order to each circuit
00712         circuit * root = subnet->getRoot ();
00713         for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00714         {
00715             c->setOrder (corrOrder);
00716             setIntegrationMethod (c, corrType);
00717         }
00718     }
00719 }
00720 
00721 /* Goes through the list of circuit objects and runs its calcDC()
00722    function. */
00723 void trsolver::calcDC (trsolver * self)
00724 {
00725     circuit * root = self->getNet()->getRoot ();
00726     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00727     {
00728         c->calcDC ();
00729     }
00730 }
00731 
00732 /* Goes through the list of circuit objects and runs its calcTR()
00733    function. */
00734 void trsolver::calcTR (trsolver * self)
00735 {
00736     circuit * root = self->getNet()->getRoot ();
00737     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00738     {
00739         c->calcTR (self->current);
00740     }
00741 }
00742 
00743 /* Goes through the list of non-linear circuit objects and runs its
00744    restartDC() function. */
00745 void trsolver::restart (void)
00746 {
00747     circuit * root = subnet->getRoot ();
00748     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00749     {
00750         if (c->isNonLinear ()) c->restartDC ();
00751     }
00752 }
00753 
00754 /* Goes through the list of circuit objects and runs its initDC()
00755    function. */
00756 void trsolver::initDC (void)
00757 {
00758     circuit * root = subnet->getRoot ();
00759     for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
00760     {
00761         c->initDC ();
00762     }
00763 }
00764 
00765 /* Goes through the list of circuit objects and runs its initTR()
00766    function. */
00767 void trsolver::initTR (void)
00768 {
00769     const char * const IMethod = getPropertyString ("IntegrationMethod");
00770     nr_double_t start = getPropertyDouble ("Start");
00771     nr_double_t stop = getPropertyDouble ("Stop");
00772     nr_double_t points = getPropertyDouble ("Points");
00773 
00774     // fetch corrector integration method and determine predicor method
00775     corrMaxOrder = getPropertyInteger ("Order");
00776     corrType = CMethod = correctorType (IMethod, corrMaxOrder);
00777     predType = PMethod = predictorType (CMethod, corrMaxOrder, predMaxOrder);
00778     corrOrder = corrMaxOrder;
00779     predOrder = predMaxOrder;
00780 
00781     // initialize step values
00782     delta = getPropertyDouble ("InitialStep");
00783     deltaMin = getPropertyDouble ("MinStep");
00784     deltaMax = getPropertyDouble ("MaxStep");
00785     if (deltaMax == 0.0)
00786         deltaMax = std::min ((stop - start) / (points - 1), stop / 200);
00787     if (deltaMin == 0.0)
00788         deltaMin = NR_TINY * 10 * deltaMax;
00789     if (delta == 0.0)
00790         delta = std::min (stop / 200, deltaMax) / 10;
00791     if (delta < deltaMin) delta = deltaMin;
00792     if (delta > deltaMax) delta = deltaMax;
00793 
00794     // initialize step history
00795     setStates (2);
00796     initStates ();
00797     // initialise the history of states, setting them all to 'delta'
00798     fillState (dState, delta);
00799 
00800     // copy the initialised states to the 'deltas' array
00801     saveState (dState, deltas);
00802     // copy the deltas to all the circuits
00803     setDelta ();
00804     // set the initial corrector and predictor coefficients
00805     calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
00806     calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
00807 
00808     // initialize history of solution vectors (solutions)
00809     for (int i = 0; i < 8; i++)
00810     {
00811         // solution contains the last sets of node voltages and branch
00812         // currents at each of the last 8 'deltas'.
00813         // Note for convenience the definition:
00814         //   #define SOL(state) (solution[(int) getState (sState, (state))])
00815         // is provided and used elsewhere to update the solutions
00816         solution[i] = new tvector<nr_double_t>;
00817         setState (sState, (nr_double_t) i, i);
00818     }
00819 
00820     // tell circuits about the transient analysis
00821     circuit *c, * root = subnet->getRoot ();
00822     for (c = root; c != NULL; c = (circuit *) c->getNext ())
00823         initCircuitTR (c);
00824     // also initialize created circuits
00825     for (c = root; c != NULL; c = (circuit *) c->getPrev ())
00826         initCircuitTR (c);
00827 }
00828 
00829 // This function cleans up some memory used by the transient analysis.
00830 void trsolver::deinitTR (void)
00831 {
00832     // cleanup solutions
00833     for (int i = 0; i < 8; i++)
00834     {
00835         delete solution[i];
00836         solution[i] = NULL;
00837     }
00838     // cleanup history
00839     if (tHistory)
00840     {
00841         delete tHistory;
00842         tHistory = NULL;
00843     }
00844 }
00845 
00846 // The function initialize a single circuit.
00847 void trsolver::initCircuitTR (circuit * c)
00848 {
00849     c->initTR ();
00850     c->initStates ();
00851     c->setCoefficients (corrCoeff);
00852     c->setOrder (corrOrder);
00853     setIntegrationMethod (c, corrType);
00854 }
00855 
00856 /* This function saves the results of a single solve() functionality
00857    (for the given timestamp) into the output dataset. */
00858 void trsolver::saveAllResults (nr_double_t time)
00859 {
00860     qucs::vector * t;
00861     // add current frequency to the dependency of the output dataset
00862     if ((t = data->findDependency ("time")) == NULL)
00863     {
00864       t = new qucs::vector ("time");
00865         data->addDependency (t);
00866     }
00867     if (runs == 1) t->add (time);
00868     saveResults ("Vt", "It", 0, t);
00869 }
00870 
00871 /* This function is meant to adapt the current time-step the transient
00872    analysis advanced.  For the computation of the new time-step the
00873    truncation error depending on the integration method is used. */
00874 nr_double_t trsolver::checkDelta (void)
00875 {
00876     nr_double_t LTEreltol = getPropertyDouble ("LTEreltol");
00877     nr_double_t LTEabstol = getPropertyDouble ("LTEabstol");
00878     nr_double_t LTEfactor = getPropertyDouble ("LTEfactor");
00879     nr_double_t dif, rel, tol, lte, q, n =  std::numeric_limits<nr_double_t>::max();
00880     int N = countNodes ();
00881     int M = countVoltageSources ();
00882 
00883     // cec = corrector error constant
00884     nr_double_t cec = getCorrectorError (corrType, corrOrder);
00885     // pec = predictor error constant
00886     nr_double_t pec = getPredictorError (predType, predOrder);
00887 
00888     // go through each solution
00889     for (int r = 0; r < N + M; r++)
00890     {
00891 
00892         // skip real voltage sources
00893         if (r >= N)
00894         {
00895             if (findVoltageSource(r - N)->isVSource ())
00896                 continue;
00897         }
00898 
00899         dif = x->get (r) - SOL(0)->get (r);
00900         if (std::isfinite (dif) && dif != 0)
00901         {
00902             // use Milne' estimate for the local truncation error
00903             rel = MAX (fabs (x->get (r)), fabs (SOL(0)->get (r)));
00904             tol = LTEreltol * rel + LTEabstol;
00905             lte = LTEfactor * (cec / (pec - cec)) * dif;
00906             q =  delta * exp (log (fabs (tol / lte)) / (corrOrder + 1));
00907             n = std::min (n, q);
00908         }
00909     }
00910 #if STEPDEBUG
00911     logprint (LOG_STATUS, "DEBUG: delta according to local truncation "
00912               "error h = %.3e\n", (double) n);
00913 #endif
00914     delta = std::min ((n > 1.9 * delta) ? 2 * delta : delta, n);
00915     return delta;
00916 }
00917 
00918 // The function updates the integration coefficients.
00919 void trsolver::updateCoefficients (nr_double_t delta)
00920 {
00921     setState (dState, delta);
00922     saveState (dState, deltas);
00923     calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
00924     calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
00925 }
00926 
00927 // properties
00928 PROP_REQ [] =
00929 {
00930     {
00931         "Type", PROP_STR, { PROP_NO_VAL, "lin" },
00932         PROP_RNG_STR2 ("lin", "log")
00933     },
00934     { "Start", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00935     { "Stop", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE },
00936     { "Points", PROP_INT, { 10, PROP_NO_STR }, PROP_MIN_VAL (2) },
00937     PROP_NO_PROP
00938 };
00939 PROP_OPT [] =
00940 {
00941     {
00942         "IntegrationMethod", PROP_STR, { PROP_NO_VAL, "Trapezoidal" },
00943         PROP_RNG_STR4 ("Euler", "Trapezoidal", "Gear", "AdamsMoulton")
00944     },
00945     { "Order", PROP_INT, { 2, PROP_NO_STR }, PROP_RNGII (1, 6) },
00946     { "InitialStep", PROP_REAL, { 1e-9, PROP_NO_STR }, PROP_POS_RANGE },
00947     { "MinStep", PROP_REAL, { 1e-16, PROP_NO_STR }, PROP_POS_RANGE },
00948     { "MaxStep", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
00949     { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
00950     { "abstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
00951     { "vntol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
00952     { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
00953     { "LTEabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
00954     { "LTEreltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
00955     { "LTEfactor", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 16) },
00956     { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00957     { "Solver", PROP_STR, { PROP_NO_VAL, "CroutLU" }, PROP_RNG_SOL },
00958     { "relaxTSR", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
00959     { "initialDC", PROP_STR, { PROP_NO_VAL, "yes" }, PROP_RNG_YESNO },
00960     PROP_NO_PROP
00961 };
00962 struct define_t trsolver::anadef =
00963     { "TR", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
00964 
00965 } // namespace qucs