Branch data Line data Source code
1 : : /*
2 : : * trsolver.cpp - transient solver class implementation
3 : : *
4 : : * Copyright (C) 2004, 2005, 2006, 2007, 2009 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <string.h>
31 : : #include <math.h>
32 : : #include <float.h>
33 : :
34 : : #include "compat.h"
35 : : #include "object.h"
36 : : #include "logging.h"
37 : : #include "complex.h"
38 : : #include "circuit.h"
39 : : #include "sweep.h"
40 : : #include "net.h"
41 : : #include "netdefs.h"
42 : : #include "analysis.h"
43 : : #include "nasolver.h"
44 : : #include "history.h"
45 : : #include "trsolver.h"
46 : : #include "transient.h"
47 : : #include "exception.h"
48 : : #include "exceptionstack.h"
49 : :
50 : : #define STEPDEBUG 0 // set to zero for release
51 : : #define BREAKPOINTS 0 // exact breakpoint calculation
52 : :
53 : : #define dState 0 // delta T state
54 : : #define sState 1 // solution state
55 : :
56 : : // Macro for the n-th state of the solution vector history.
57 : : #define SOL(state) (solution[(int) getState (sState, (state))])
58 : :
59 : : namespace qucs {
60 : :
61 : : using namespace transient;
62 : :
63 : : // Constructor creates an unnamed instance of the trsolver class.
64 : 38 : trsolver::trsolver ()
65 : 38 : : nasolver<nr_double_t> (), states<nr_double_t> ()
66 : : {
67 : 38 : swp = NULL;
68 : 38 : type = ANALYSIS_TRANSIENT;
69 : 38 : setDescription ("transient");
70 [ + + ][ # # ]: 342 : for (int i = 0; i < 8; i++) solution[i] = NULL;
71 : 38 : tHistory = NULL;
72 : 38 : relaxTSR = false;
73 : 38 : initialDC = true;
74 : 38 : }
75 : :
76 : : // Constructor creates a named instance of the trsolver class.
77 : 0 : trsolver::trsolver (char * n)
78 : 0 : : nasolver<nr_double_t> (n), states<nr_double_t> ()
79 : : {
80 : 0 : swp = NULL;
81 : 0 : type = ANALYSIS_TRANSIENT;
82 : 0 : setDescription ("transient");
83 [ # # ][ # # ]: 0 : for (int i = 0; i < 8; i++) solution[i] = NULL;
84 : 0 : tHistory = NULL;
85 : 0 : relaxTSR = false;
86 : 0 : initialDC = true;
87 : 0 : }
88 : :
89 : : // Destructor deletes the trsolver class object.
90 : 38 : trsolver::~trsolver ()
91 : : {
92 [ + - ][ + - ]: 38 : if (swp) delete swp;
[ + - ][ # # ]
[ # # ][ # # ]
93 [ + + ][ # # ]: 342 : for (int i = 0; i < 8; i++)
94 : : {
95 [ - + ][ # # ]: 304 : if (solution[i] != NULL)
96 : : {
97 [ # # ][ # # ]: 0 : delete solution[i];
98 : : }
99 : : }
100 [ - + ][ # # ]: 38 : if (tHistory) delete tHistory;
[ # # ][ # # ]
101 [ - + ][ # # ]: 76 : }
102 : :
103 : : /* The copy constructor creates a new instance of the trsolver class
104 : : based on the given trsolver object. */
105 : 0 : trsolver::trsolver (trsolver & o)
106 : 0 : : nasolver<nr_double_t> (o), states<nr_double_t> (o)
107 : : {
108 [ # # ][ # # ]: 0 : swp = o.swp ? new sweep (*o.swp) : NULL;
[ # # ][ # # ]
[ # # ][ # # ]
109 [ # # ][ # # ]: 0 : for (int i = 0; i < 8; i++) solution[i] = NULL;
110 [ # # ][ # # ]: 0 : tHistory = o.tHistory ? new history (*o.tHistory) : NULL;
[ # # ][ # # ]
[ # # ][ # # ]
111 : 0 : relaxTSR = o.relaxTSR;
112 : 0 : initialDC = o.initialDC;
113 : 0 : }
114 : :
115 : : // This function creates the time sweep if necessary.
116 : 65 : void trsolver::initSteps (void)
117 : : {
118 [ + + ][ + - ]: 65 : if (swp != NULL) delete swp;
119 : 65 : swp = createSweep ("time");
120 : 65 : }
121 : :
122 : : // Performs the initial DC analysis.
123 : 64 : int trsolver::dcAnalysis (void)
124 : : {
125 : 64 : int error = 0;
126 : :
127 : : // First calculate a initial state using the non-linear DC analysis.
128 : 64 : setDescription ("initial DC");
129 : 64 : initDC ();
130 : 64 : setCalculation ((calculate_func_t) &calcDC);
131 : 64 : solve_pre ();
132 : 64 : applyNodeset ();
133 : :
134 : : // Run the DC solver once.
135 : : try_running ()
136 : : {
137 : 64 : error = solve_nonlinear ();
138 : : }
139 : : // Appropriate exception handling.
140 [ - + ][ # # ]: 64 : catch_exception ()
141 : : {
142 : : case EXCEPTION_NO_CONVERGENCE:
143 : 0 : pop_exception ();
144 : 0 : convHelper = CONV_LineSearch;
145 : : logprint (LOG_ERROR, "WARNING: %s: %s analysis failed, using line search "
146 : 0 : "fallback\n", getName (), getDescription ());
147 : 0 : applyNodeset ();
148 : 0 : restart ();
149 : 0 : error = solve_nonlinear ();
150 : 0 : break;
151 : : default:
152 : : // Otherwise return.
153 : 0 : estack.print ();
154 : 0 : error++;
155 : 0 : break;
156 : : }
157 : :
158 : : // Save the DC solution.
159 : 64 : storeSolution ();
160 : :
161 : : // Cleanup nodal analysis solver.
162 : 64 : solve_post ();
163 : :
164 : : // Really failed to find initial DC solution?
165 [ - + ]: 64 : if (error)
166 : : {
167 : : logprint (LOG_ERROR, "ERROR: %s: %s analysis failed\n",
168 : 0 : getName (), getDescription ());
169 : : }
170 : 64 : return error;
171 : : }
172 : :
173 : : /* This is the transient netlist solver. It prepares the circuit list
174 : : for each requested time and solves it then. */
175 : 65 : int trsolver::solve (void)
176 : : {
177 : : nr_double_t time, saveCurrent;
178 : 65 : int error = 0, convError = 0;
179 : 65 : const char * const solver = getPropertyString ("Solver");
180 : 65 : relaxTSR = !strcmp (getPropertyString ("relaxTSR"), "yes") ? true : false;
181 : 65 : initialDC = !strcmp (getPropertyString ("initialDC"), "yes") ? true : false;
182 : :
183 : 65 : runs++;
184 : 65 : saveCurrent = current = 0;
185 : 65 : stepDelta = -1;
186 : 65 : converged = 0;
187 : 65 : fixpoint = 0;
188 : 65 : statRejected = statSteps = statIterations = statConvergence = 0;
189 : :
190 : : // Choose a solver.
191 [ + - ]: 65 : if (!strcmp (solver, "CroutLU"))
192 : 65 : eqnAlgo = ALGO_LU_DECOMPOSITION;
193 [ # # ]: 0 : else if (!strcmp (solver, "DoolittleLU"))
194 : 0 : eqnAlgo = ALGO_LU_DECOMPOSITION_DOOLITTLE;
195 [ # # ]: 0 : else if (!strcmp (solver, "HouseholderQR"))
196 : 0 : eqnAlgo = ALGO_QR_DECOMPOSITION;
197 [ # # ]: 0 : else if (!strcmp (solver, "HouseholderLQ"))
198 : 0 : eqnAlgo = ALGO_QR_DECOMPOSITION_LS;
199 [ # # ]: 0 : else if (!strcmp (solver, "GolubSVD"))
200 : 0 : eqnAlgo = ALGO_SV_DECOMPOSITION;
201 : :
202 : : // Perform initial DC analysis.
203 [ + + ]: 65 : if (initialDC)
204 : : {
205 : 64 : error = dcAnalysis ();
206 [ - + ]: 64 : if (error)
207 : 0 : return -1;
208 : : }
209 : :
210 : : // Initialize transient analysis.
211 : 65 : setDescription ("transient");
212 : 65 : initTR ();
213 : 65 : setCalculation ((calculate_func_t) &calcTR);
214 : 65 : solve_pre ();
215 : :
216 : : // Create time sweep if necessary.
217 : 65 : initSteps ();
218 : 65 : swp->reset ();
219 : :
220 : : // Recall the DC solution.
221 : 65 : recallSolution ();
222 : :
223 : : // Apply the nodesets and adjust previous solutions.
224 : 65 : applyNodeset (false);
225 : 65 : fillSolution (x);
226 : :
227 : : // Tell integrators to be initialized.
228 : 65 : setMode (MODE_INIT);
229 : :
230 : 65 : int running = 0;
231 : 65 : rejected = 0;
232 : 65 : delta /= 10;
233 : 65 : fillState (dState, delta);
234 : 65 : adjustOrder (1);
235 : :
236 : : // Start to sweep through time.
237 [ + + ]: 43242 : for (int i = 0; i < swp->getSize (); i++)
238 : : {
239 : 43177 : time = swp->next ();
240 [ + + ]: 43177 : if (progress) logprogressbar (i, swp->getSize (), 40);
241 : :
242 : : #if DEBUG && 0
243 : : logprint (LOG_STATUS, "NOTIFY: %s: solving netlist for t = %e\n",
244 : : getName (), (double) time);
245 : : #endif
246 : :
247 [ + + ]: 214485 : do // while (saveCurrent < time), i.e. until a requested breakpoint is hit
248 : : {
249 : : #if STEPDEBUG
250 : : if (delta == deltaMin)
251 : : {
252 : : // the integrator step size has become smaller than the
253 : : // specified allowed minimum, Qucs is unable to solve the circuit
254 : : // while meeting the tolerance conditions
255 : : logprint (LOG_ERROR,
256 : : "WARNING: %s: minimum delta h = %.3e at t = %.3e\n",
257 : : getName (), (double) delta, (double) current);
258 : : }
259 : : #endif
260 : : // updates the integrator coefficients, and updates the array of prev
261 : : // 8 deltas with the new delta for this step
262 : 214485 : updateCoefficients (delta);
263 : :
264 : : // Run predictor to get a start value for the solution vector for
265 : : // the successive iterative corrector process
266 : 214485 : error += predictor ();
267 : :
268 : : // restart Newton iteration
269 [ + + ]: 214485 : if (rejected)
270 : : {
271 : 22030 : restart (); // restart non-linear devices
272 : 22030 : rejected = 0;
273 : : }
274 : :
275 : : // Run corrector process with appropriate exception handling.
276 : : // The corrector iterates through the solutions of the integration
277 : : // process until a certain error tolerance has been reached.
278 : : try_running () // #defined as: do {
279 : : {
280 : 214485 : error += corrector ();
281 : : }
282 [ + + ][ + - ]: 214485 : catch_exception () // #defined as: } while (0); if (estack.top ()) switch (estack.top()->getCode ())
283 : : {
284 : : case EXCEPTION_NO_CONVERGENCE:
285 : 65 : pop_exception ();
286 : :
287 : : // step back from the current time value to the previous time
288 [ + + ]: 65 : if (current > 0) current -= delta;
289 : : // Reduce step-size (by half) if failed to converge.
290 : 65 : delta /= 2;
291 [ + + ]: 65 : if (delta <= deltaMin)
292 : : {
293 : : // but do not reduce the step size below a specified minimum
294 : 14 : delta = deltaMin;
295 : : // instead reduce the order of the integration
296 : 14 : adjustOrder (1);
297 : : }
298 : : // step forward to the new current time value
299 [ + + ]: 65 : if (current > 0) current += delta;
300 : :
301 : : // Update statistics.
302 : 65 : statRejected++;
303 : 65 : statConvergence++;
304 : 65 : rejected++; // mark the previous step size choice as rejected
305 : 65 : converged = 0;
306 : 65 : error = 0;
307 : :
308 : : // Start using damped Newton-Raphson.
309 : 65 : convHelper = CONV_SteepestDescent;
310 : 65 : convError = 2;
311 : : #if DEBUG
312 : : logprint (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
313 : 65 : "(no convergence)\n", (double) saveCurrent, (double) delta);
314 : : #endif
315 : 65 : break;
316 : : default:
317 : : // Otherwise return.
318 : 0 : estack.print ();
319 : 0 : error++;
320 : 65 : break;
321 : : }
322 : : // return if any errors occured other than convergence failure
323 [ - + ]: 214485 : if (error) return -1;
324 : :
325 : : // if the step was rejected, the solution loop is restarted here
326 [ + + ]: 214485 : if (rejected) continue;
327 : :
328 : : // check whether Jacobian matrix is still non-singular
329 [ - + ]: 214420 : if (!A->isFinite ())
330 : : {
331 : : logprint (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
332 : : "aborting %s analysis\n", getName (), (double) current,
333 : 0 : getDescription ());
334 : 0 : return -1;
335 : : }
336 : :
337 : : // Update statistics and no more damped Newton-Raphson.
338 : 214420 : statIterations += iterations;
339 [ + + ]: 214420 : if (--convError < 0) convHelper = 0;
340 : :
341 : : // Now advance in time or not...
342 [ + + ]: 214420 : if (running > 1)
343 : : {
344 : 214290 : adjustDelta (time);
345 : 214290 : adjustOrder ();
346 : : }
347 : : else
348 : : {
349 : 130 : fillStates ();
350 : 130 : nextStates ();
351 : 130 : rejected = 0;
352 : : }
353 : :
354 : 214420 : saveCurrent = current;
355 : 214420 : current += delta;
356 : 214420 : running++;
357 : 214420 : converged++;
358 : :
359 : : // Tell integrators to be running.
360 : 214420 : setMode (MODE_NONE);
361 : :
362 : : // Initialize or update history.
363 [ + + ]: 214420 : if (running > 1)
364 : : {
365 : 214355 : updateHistory (saveCurrent);
366 : : }
367 : : else
368 : : {
369 : 65 : initHistory (saveCurrent);
370 : : }
371 : : }
372 : : while (saveCurrent < time); // Hit a requested time point?
373 : :
374 : : // Save results.
375 : : #if STEPDEBUG
376 : : logprint (LOG_STATUS, "DEBUG: save point at t = %.3e, h = %.3e\n",
377 : : (double) saveCurrent, (double) delta);
378 : : #endif
379 : :
380 : : #if BREAKPOINTS
381 : : saveAllResults (saveCurrent);
382 : : #else
383 : 43177 : saveAllResults (time);
384 : : #endif
385 : : } // for (int i = 0; i < swp->getSize (); i++)
386 : :
387 : 65 : solve_post ();
388 [ + + ]: 65 : if (progress) logprogressclear (40);
389 : : logprint (LOG_STATUS, "NOTIFY: %s: average time-step %g, %d rejections\n",
390 : 65 : getName (), (double) (saveCurrent / statSteps), statRejected);
391 : : logprint (LOG_STATUS, "NOTIFY: %s: average NR-iterations %g, "
392 : : "%d non-convergences\n", getName (),
393 : 65 : (double) statIterations / statSteps, statConvergence);
394 : :
395 : : // cleanup
396 : 65 : deinitTR ();
397 : 65 : return 0;
398 : : }
399 : :
400 : : // The function initializes the history.
401 : 65 : void trsolver::initHistory (nr_double_t t)
402 : : {
403 : : // initialize time vector
404 [ + - ]: 65 : tHistory = new history ();
405 : 65 : tHistory->append (t);
406 : 65 : tHistory->self ();
407 : : // initialize circuit histories
408 : 65 : nr_double_t age = 0.0;
409 : 65 : circuit * root = subnet->getRoot ();
410 [ + + ]: 805 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
411 : : {
412 [ - + ]: 740 : if (c->hasHistory ())
413 : : {
414 : 0 : c->applyHistory (tHistory);
415 : 0 : saveHistory (c);
416 [ # # ]: 0 : if (c->getHistoryAge () > age)
417 : : {
418 : 0 : age = c->getHistoryAge ();
419 : : }
420 : : }
421 : : }
422 : : // set maximum required age for all circuits
423 : 65 : tHistory->setAge (age);
424 : 65 : }
425 : :
426 : : /* The following function updates the histories for the circuits which
427 : : requested them. */
428 : 214355 : void trsolver::updateHistory (nr_double_t t)
429 : : {
430 [ + + ]: 214355 : if (t > tHistory->last ())
431 : : {
432 : : // update time vector
433 : 192404 : tHistory->append (t);
434 : : // update circuit histories
435 : 192404 : circuit * root = subnet->getRoot ();
436 [ + + ]: 1981771 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
437 : : {
438 [ - + ]: 1789367 : if (c->hasHistory ()) saveHistory (c);
439 : : }
440 : 192404 : tHistory->drop ();
441 : : }
442 : 214355 : }
443 : :
444 : : // Stores node voltages and branch currents in the given circuits history.
445 : 0 : void trsolver::saveHistory (circuit * c)
446 : : {
447 : :
448 : 0 : int N = countNodes ();
449 : 0 : int r, i, s = c->getSize ();
450 : :
451 [ # # ]: 0 : for (i = 0; i < s; i++)
452 : : {
453 : : // save node voltages
454 : 0 : r = findAssignedNode (c, i);
455 [ # # ]: 0 : if (r < 0)
456 : : // the node was not found, append a zero to the history
457 : : // matching this index
458 : 0 : c->appendHistory (i, 0.0);
459 : : else
460 : : // the node was found, append the voltage value to
461 : : // that node's history
462 : 0 : c->appendHistory (i, x->get (r));
463 : : }
464 : :
465 [ # # ]: 0 : for (i = 0; i < c->getVoltageSources (); i++)
466 : : {
467 : : // save branch currents
468 : 0 : r = c->getVoltageSource () + i;
469 : 0 : c->appendHistory (i + s, x->get (r + N));
470 : : }
471 : :
472 : 0 : }
473 : :
474 : : /* This function predicts a start value for the solution vector for
475 : : the successive iterative corrector process. */
476 : 214485 : int trsolver::predictor (void)
477 : : {
478 : 214485 : int error = 0;
479 [ + + + - ]: 214485 : switch (predType)
480 : : {
481 : : case INTEGRATOR_GEAR: // explicit GEAR
482 : 1671 : predictGear ();
483 : 1671 : break;
484 : : case INTEGRATOR_ADAMSBASHFORD: // ADAMS-BASHFORD
485 : 212616 : predictBashford ();
486 : 212616 : break;
487 : : case INTEGRATOR_EULER: // FORWARD EULER
488 : 198 : predictEuler ();
489 : 198 : break;
490 : : default:
491 : 0 : *x = *SOL (1); // This is too a simple predictor...
492 : 0 : break;
493 : : }
494 : 214485 : saveSolution ();
495 : 214485 : *SOL (0) = *x;
496 : 214485 : return error;
497 : : }
498 : :
499 : : // Stores the given vector into all the solution vectors.
500 : 65 : void trsolver::fillSolution (tvector<nr_double_t> * s)
501 : : {
502 [ + + ]: 585 : for (int i = 0; i < 8; i++) *SOL (i) = *s;
503 : 65 : }
504 : :
505 : : /* The function predicts the successive solution vector using the
506 : : explicit Adams-Bashford integration formula. */
507 : 212616 : void trsolver::predictBashford (void)
508 : : {
509 : 212616 : int N = countNodes ();
510 : 212616 : int M = countVoltageSources ();
511 : : nr_double_t xn, dd, hn;
512 : :
513 : : // go through each solution
514 [ + + ]: 2097318 : for (int r = 0; r < N + M; r++)
515 : : {
516 : 1884702 : xn = predCoeff[0] * SOL(1)->get (r); // a0 coefficient
517 [ + + ]: 5654106 : for (int o = 1; o <= predOrder; o++)
518 : : {
519 : 3769404 : hn = getState (dState, o); // previous time-step
520 : : // divided differences
521 : 3769404 : dd = (SOL(o)->get (r) - SOL(o + 1)->get (r)) / hn;
522 : 3769404 : xn += predCoeff[o] * dd; // b0, b1, ... coefficients
523 : : }
524 : 1884702 : x->set (r, xn); // save prediction
525 : : }
526 : 212616 : }
527 : :
528 : : /* The function predicts the successive solution vector using the
529 : : explicit forward Euler integration formula. Actually this is
530 : : Adams-Bashford order 1. */
531 : 198 : void trsolver::predictEuler (void)
532 : : {
533 : 198 : int N = countNodes ();
534 : 198 : int M = countVoltageSources ();
535 : : nr_double_t xn, dd, hn;
536 : :
537 [ + + ]: 2352 : for (int r = 0; r < N + M; r++)
538 : : {
539 : 2154 : xn = predCoeff[0] * SOL(1)->get (r);
540 : 2154 : hn = getState (dState, 1);
541 : 2154 : dd = (SOL(1)->get (r) - SOL(2)->get (r)) / hn;
542 : 2154 : xn += predCoeff[1] * dd;
543 : 2154 : x->set (r, xn);
544 : : }
545 : 198 : }
546 : :
547 : : /* The function predicts the successive solution vector using the
548 : : explicit Gear integration formula. */
549 : 1671 : void trsolver::predictGear (void)
550 : : {
551 : 1671 : int N = countNodes ();
552 : 1671 : int M = countVoltageSources ();
553 : : nr_double_t xn;
554 : :
555 : : // go through each solution
556 [ + + ]: 25065 : for (int r = 0; r < N + M; r++)
557 : : {
558 : 23394 : xn = 0;
559 [ + + ]: 70182 : for (int o = 0; o <= predOrder; o++)
560 : : {
561 : : // a0, a1, ... coefficients
562 : 46788 : xn += predCoeff[o] * SOL(o + 1)->get (r);
563 : : }
564 : 23394 : x->set (r, xn); // save prediction
565 : : }
566 : 1671 : }
567 : :
568 : : /* The function iterates through the solutions of the integration
569 : : process until a certain error tolerance has been reached. */
570 : 214485 : int trsolver::corrector (void)
571 : : {
572 : 214485 : int error = 0;
573 : 214485 : error += solve_nonlinear ();
574 : 214485 : return error;
575 : : }
576 : :
577 : : // The function advances one more time-step.
578 : 192455 : void trsolver::nextStates (void)
579 : : {
580 : 192455 : circuit * root = subnet->getRoot ();
581 [ + + ]: 1982427 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
582 : : {
583 : : // for each circuit get the next state
584 : 1789972 : c->nextState ();
585 : : }
586 : :
587 : 192455 : *SOL (0) = *x; // save current solution
588 : 192455 : nextState ();
589 : 192455 : statSteps++;
590 : 192455 : }
591 : :
592 : : /* This function stores the current state of each circuit into all
593 : : other states as well. It is useful for higher order integration
594 : : methods in order to initialize the states after the initial
595 : : transient solution. */
596 : 130 : void trsolver::fillStates (void)
597 : : {
598 : 130 : circuit * root = subnet->getRoot ();
599 [ + + ]: 1610 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
600 : : {
601 [ + + ]: 2708 : for (int s = 0; s < c->getStates (); s++)
602 : 1228 : c->fillState (s, c->getState (s));
603 : : }
604 : 130 : }
605 : :
606 : : // The function modifies the circuit lists integrator mode.
607 : 214485 : void trsolver::setMode (int state)
608 : : {
609 : 214485 : circuit * root = subnet->getRoot ();
610 [ + + ]: 2199981 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
611 : 1985496 : c->setMode (state);
612 : 214485 : }
613 : :
614 : : // The function passes the time delta array to the circuit list.
615 : 65 : void trsolver::setDelta (void)
616 : : {
617 : 65 : circuit * root = subnet->getRoot ();
618 [ + + ]: 805 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
619 : 740 : c->setDelta (deltas);
620 : 65 : }
621 : :
622 : : /* This function tries to adapt the current time-step according to the
623 : : global truncation error. */
624 : 214290 : void trsolver::adjustDelta (nr_double_t t)
625 : : {
626 : 214290 : deltaOld = delta;
627 : 214290 : delta = checkDelta ();
628 [ + + ]: 214290 : if (delta > deltaMax) delta = deltaMax;
629 [ + + ]: 214290 : if (delta < deltaMin) delta = deltaMin;
630 : :
631 : : // delta correction in order to hit exact breakpoint
632 : 214290 : int good = 0;
633 [ + + ]: 214290 : if (!relaxTSR) // relaxed step raster?
634 : : {
635 [ + + ][ + + ]: 212681 : if (!statConvergence || converged > 64) /* Is this a good guess? */
636 : : {
637 : : // check next breakpoint
638 [ + + ]: 212424 : if (stepDelta > 0.0)
639 : : {
640 : : // restore last valid delta
641 : 59546 : delta = stepDelta;
642 : 59546 : stepDelta = -1.0;
643 : : }
644 : : else
645 : : {
646 [ + + ][ + + ]: 152878 : if (delta > (t - current) && t > current)
647 : : {
648 : : // save last valid delta and set exact step
649 : 59546 : stepDelta = deltaOld;
650 : 59546 : delta = t - current;
651 : 59546 : good = 1;
652 : : }
653 : : else
654 : : {
655 : 93332 : stepDelta = -1.0;
656 : : }
657 : : }
658 [ - + ]: 212424 : if (delta > deltaMax) delta = deltaMax;
659 [ + + ]: 212424 : if (delta < deltaMin) delta = deltaMin;
660 : : }
661 : : }
662 : :
663 : : // usual delta correction
664 [ + + ][ + + ]: 214290 : if (delta > 0.9 * deltaOld || good) // accept current delta
665 : : {
666 : 192325 : nextStates ();
667 : 192325 : rejected = 0;
668 : : #if STEPDEBUG
669 : : logprint (LOG_STATUS,
670 : : "DEBUG: delta accepted at t = %.3e, h = %.3e\n",
671 : : (double) current, (double) delta);
672 : : #endif
673 : : }
674 [ + - ]: 21965 : else if (deltaOld > delta) // reject current delta
675 : : {
676 : 21965 : rejected++;
677 : 21965 : statRejected++;
678 : : #if STEPDEBUG
679 : : logprint (LOG_STATUS,
680 : : "DEBUG: delta rejected at t = %.3e, h = %.3e\n",
681 : : (double) current, (double) delta);
682 : : #endif
683 [ + - ]: 21965 : if (current > 0) current -= deltaOld;
684 : : }
685 : : else
686 : : {
687 : 0 : nextStates ();
688 : 0 : rejected = 0;
689 : : }
690 : 214290 : }
691 : :
692 : : /* The function can be used to increase the current order of the
693 : : integration method or to reduce it. */
694 : 214369 : void trsolver::adjustOrder (int reduce)
695 : : {
696 [ + + ][ + + ]: 214369 : if ((corrOrder < corrMaxOrder && !rejected) || reduce)
[ + + ]
697 : : {
698 [ + + ]: 143 : if (reduce)
699 : : {
700 : 79 : corrOrder = 1;
701 : : }
702 [ + - ]: 64 : else if (!rejected)
703 : : {
704 : 64 : corrOrder++;
705 : : }
706 : :
707 : : // adjust type and order of corrector and predictor
708 : 143 : corrType = correctorType (CMethod, corrOrder);
709 : 143 : predType = predictorType (corrType, corrOrder, predOrder);
710 : :
711 : : // apply new corrector method and order to each circuit
712 : 143 : circuit * root = subnet->getRoot ();
713 [ + + ]: 1805 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
714 : : {
715 : 1662 : c->setOrder (corrOrder);
716 : 1662 : setIntegrationMethod (c, corrType);
717 : : }
718 : : }
719 : 214369 : }
720 : :
721 : : /* Goes through the list of circuit objects and runs its calcDC()
722 : : function. */
723 : 267 : void trsolver::calcDC (trsolver * self)
724 : : {
725 : 267 : circuit * root = self->getNet()->getRoot ();
726 [ + + ]: 4525 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
727 : : {
728 : 4258 : c->calcDC ();
729 : : }
730 : 267 : }
731 : :
732 : : /* Goes through the list of circuit objects and runs its calcTR()
733 : : function. */
734 : 483576 : void trsolver::calcTR (trsolver * self)
735 : : {
736 : 483576 : circuit * root = self->getNet()->getRoot ();
737 [ + + ]: 5148245 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
738 : : {
739 : 4664669 : c->calcTR (self->current);
740 : : }
741 : 483576 : }
742 : :
743 : : /* Goes through the list of non-linear circuit objects and runs its
744 : : restartDC() function. */
745 : 22030 : void trsolver::restart (void)
746 : : {
747 : 22030 : circuit * root = subnet->getRoot ();
748 [ + + ]: 217710 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
749 : : {
750 [ + + ]: 195680 : if (c->isNonLinear ()) c->restartDC ();
751 : : }
752 : 22030 : }
753 : :
754 : : /* Goes through the list of circuit objects and runs its initDC()
755 : : function. */
756 : 64 : void trsolver::initDC (void)
757 : : {
758 : 64 : circuit * root = subnet->getRoot ();
759 [ + + ]: 766 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
760 : : {
761 : 702 : c->initDC ();
762 : : }
763 : 64 : }
764 : :
765 : : /* Goes through the list of circuit objects and runs its initTR()
766 : : function. */
767 : 65 : void trsolver::initTR (void)
768 : : {
769 : 65 : const char * const IMethod = getPropertyString ("IntegrationMethod");
770 : 65 : nr_double_t start = getPropertyDouble ("Start");
771 : 65 : nr_double_t stop = getPropertyDouble ("Stop");
772 : 65 : nr_double_t points = getPropertyDouble ("Points");
773 : :
774 : : // fetch corrector integration method and determine predicor method
775 : 65 : corrMaxOrder = getPropertyInteger ("Order");
776 : 65 : corrType = CMethod = correctorType (IMethod, corrMaxOrder);
777 : 65 : predType = PMethod = predictorType (CMethod, corrMaxOrder, predMaxOrder);
778 : 65 : corrOrder = corrMaxOrder;
779 : 65 : predOrder = predMaxOrder;
780 : :
781 : : // initialize step values
782 : 65 : delta = getPropertyDouble ("InitialStep");
783 : 65 : deltaMin = getPropertyDouble ("MinStep");
784 : 65 : deltaMax = getPropertyDouble ("MaxStep");
785 [ + + ]: 65 : if (deltaMax == 0.0)
786 [ + + ]: 62 : deltaMax = MIN ((stop - start) / (points - 1), stop / 200);
787 [ - + ]: 65 : if (deltaMin == 0.0)
788 : 0 : deltaMin = NR_TINY * 10 * deltaMax;
789 [ - + ]: 65 : if (delta == 0.0)
790 [ # # ]: 0 : delta = MIN (stop / 200, deltaMax) / 10;
791 [ + + ]: 65 : if (delta < deltaMin) delta = deltaMin;
792 [ + + ]: 65 : if (delta > deltaMax) delta = deltaMax;
793 : :
794 : : // initialize step history
795 : 65 : setStates (2);
796 : 65 : initStates ();
797 : : // initialise the history of states, setting them all to 'delta'
798 : 65 : fillState (dState, delta);
799 : :
800 : : // copy the initialised states to the 'deltas' array
801 : 65 : saveState (dState, deltas);
802 : : // copy the deltas to all the circuits
803 : 65 : setDelta ();
804 : : // set the initial corrector and predictor coefficients
805 : 65 : calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
806 : 65 : calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
807 : :
808 : : // initialize history of solution vectors (solutions)
809 [ + + ]: 585 : for (int i = 0; i < 8; i++)
810 : : {
811 : : // solution contains the last sets of node voltages and branch
812 : : // currents at each of the last 8 'deltas'.
813 : : // Note for convenience the definition:
814 : : // #define SOL(state) (solution[(int) getState (sState, (state))])
815 : : // is provided and used elsewhere to update the solutions
816 : 520 : solution[i] = new tvector<nr_double_t>;
817 : 520 : setState (sState, (nr_double_t) i, i);
818 : : }
819 : :
820 : : // tell circuits about the transient analysis
821 : 65 : circuit *c, * root = subnet->getRoot ();
822 [ + + ]: 805 : for (c = root; c != NULL; c = (circuit *) c->getNext ())
823 : 740 : initCircuitTR (c);
824 : : // also initialize created circuits
825 [ + + ]: 130 : for (c = root; c != NULL; c = (circuit *) c->getPrev ())
826 : 65 : initCircuitTR (c);
827 : 65 : }
828 : :
829 : : // This function cleans up some memory used by the transient analysis.
830 : 65 : void trsolver::deinitTR (void)
831 : : {
832 : : // cleanup solutions
833 [ + + ]: 585 : for (int i = 0; i < 8; i++)
834 : : {
835 [ + - ]: 520 : delete solution[i];
836 : 520 : solution[i] = NULL;
837 : : }
838 : : // cleanup history
839 [ + - ]: 65 : if (tHistory)
840 : : {
841 [ + - ]: 65 : delete tHistory;
842 : 65 : tHistory = NULL;
843 : : }
844 : 65 : }
845 : :
846 : : // The function initialize a single circuit.
847 : 805 : void trsolver::initCircuitTR (circuit * c)
848 : : {
849 : 805 : c->initTR ();
850 : 805 : c->initStates ();
851 : 805 : c->setCoefficients (corrCoeff);
852 : 805 : c->setOrder (corrOrder);
853 : 805 : setIntegrationMethod (c, corrType);
854 : 805 : }
855 : :
856 : : /* This function saves the results of a single solve() functionality
857 : : (for the given timestamp) into the output dataset. */
858 : 43177 : void trsolver::saveAllResults (nr_double_t time)
859 : : {
860 : : qucs::vector * t;
861 : : // add current frequency to the dependency of the output dataset
862 [ + + ]: 43177 : if ((t = data->findDependency ("time")) == NULL)
863 : : {
864 [ + - ]: 38 : t = new qucs::vector ("time");
865 : 38 : data->addDependency (t);
866 : : }
867 [ + + ][ + - ]: 43177 : if (runs == 1) t->add (time);
868 : 43177 : saveResults ("Vt", "It", 0, t);
869 : 43177 : }
870 : :
871 : : /* This function is meant to adapt the current time-step the transient
872 : : analysis advanced. For the computation of the new time-step the
873 : : truncation error depending on the integration method is used. */
874 : 214290 : nr_double_t trsolver::checkDelta (void)
875 : : {
876 : 214290 : nr_double_t LTEreltol = getPropertyDouble ("LTEreltol");
877 : 214290 : nr_double_t LTEabstol = getPropertyDouble ("LTEabstol");
878 : 214290 : nr_double_t LTEfactor = getPropertyDouble ("LTEfactor");
879 : 214290 : nr_double_t dif, rel, tol, lte, q, n = std::numeric_limits<nr_double_t>::max();
880 : 214290 : int N = countNodes ();
881 : 214290 : int M = countVoltageSources ();
882 : :
883 : : // cec = corrector error constant
884 : 214290 : nr_double_t cec = getCorrectorError (corrType, corrOrder);
885 : : // pec = predictor error constant
886 : 214290 : nr_double_t pec = getPredictorError (predType, predOrder);
887 : :
888 : : // go through each solution
889 [ + + ]: 2122231 : for (int r = 0; r < N + M; r++)
890 : : {
891 : :
892 : : // skip real voltage sources
893 [ + + ]: 1907941 : if (r >= N)
894 : : {
895 [ + + ]: 798433 : if (findVoltageSource(r - N)->isVSource ())
896 : 421843 : continue;
897 : : }
898 : :
899 : 1486098 : dif = x->get (r) - SOL(0)->get (r);
900 [ + - ][ + + ]: 1486098 : if (std::isfinite (dif) && dif != 0)
[ + + ]
901 : : {
902 : : // use Milne' estimate for the local truncation error
903 [ + + ]: 1210811 : rel = MAX (fabs (x->get (r)), fabs (SOL(0)->get (r)));
904 : 1210811 : tol = LTEreltol * rel + LTEabstol;
905 : 1210811 : lte = LTEfactor * (cec / (pec - cec)) * dif;
906 : 1210811 : q = delta * exp (log (fabs (tol / lte)) / (corrOrder + 1));
907 [ + + ]: 1210811 : n = MIN (n, q);
908 : : }
909 : : }
910 : : #if STEPDEBUG
911 : : logprint (LOG_STATUS, "DEBUG: delta according to local truncation "
912 : : "error h = %.3e\n", (double) n);
913 : : #endif
914 [ + + ][ + + ]: 214290 : delta = MIN ((n > 1.9 * delta) ? 2 * delta : delta, n);
[ + + ]
915 : 214290 : return delta;
916 : : }
917 : :
918 : : // The function updates the integration coefficients.
919 : 214485 : void trsolver::updateCoefficients (nr_double_t delta)
920 : : {
921 : 214485 : setState (dState, delta);
922 : 214485 : saveState (dState, deltas);
923 : 214485 : calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
924 : 214485 : calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
925 : 214485 : }
926 : :
927 : : // properties
928 : : PROP_REQ [] =
929 : : {
930 : : {
931 : : "Type", PROP_STR, { PROP_NO_VAL, "lin" },
932 : : PROP_RNG_STR2 ("lin", "log")
933 : : },
934 : : { "Start", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
935 : : { "Stop", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE },
936 : : { "Points", PROP_INT, { 10, PROP_NO_STR }, PROP_MIN_VAL (2) },
937 : : PROP_NO_PROP
938 : : };
939 : : PROP_OPT [] =
940 : : {
941 : : {
942 : : "IntegrationMethod", PROP_STR, { PROP_NO_VAL, "Trapezoidal" },
943 : : PROP_RNG_STR4 ("Euler", "Trapezoidal", "Gear", "AdamsMoulton")
944 : : },
945 : : { "Order", PROP_INT, { 2, PROP_NO_STR }, PROP_RNGII (1, 6) },
946 : : { "InitialStep", PROP_REAL, { 1e-9, PROP_NO_STR }, PROP_POS_RANGE },
947 : : { "MinStep", PROP_REAL, { 1e-16, PROP_NO_STR }, PROP_POS_RANGE },
948 : : { "MaxStep", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
949 : : { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
950 : : { "abstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
951 : : { "vntol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
952 : : { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
953 : : { "LTEabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
954 : : { "LTEreltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
955 : : { "LTEfactor", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 16) },
956 : : { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
957 : : { "Solver", PROP_STR, { PROP_NO_VAL, "CroutLU" }, PROP_RNG_SOL },
958 : : { "relaxTSR", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
959 : : { "initialDC", PROP_STR, { PROP_NO_VAL, "yes" }, PROP_RNG_YESNO },
960 : : PROP_NO_PROP
961 : : };
962 : : struct define_t trsolver::anadef =
963 : : { "TR", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
964 : :
965 : : } // namespace qucs
|
