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