Qucs-core
0.0.19
|
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