Branch data Line data Source code
1 : : /*
2 : : * e_trsolver.cpp - external transient solver interface 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 : : /** \file ecvs.h
26 : : * \brief The externally controlled transient solver implementation file.
27 : : *
28 : : */
29 : :
30 : : /**
31 : : * \ingroup QucsInterface
32 : : */
33 : :
34 : : #if HAVE_CONFIG_H
35 : : # include <config.h>
36 : : #endif
37 : :
38 : : #include <stdio.h>
39 : : #include <string.h>
40 : : #include <cmath>
41 : : #include <float.h>
42 : :
43 : : #include "compat.h"
44 : : #include "object.h"
45 : : #include "logging.h"
46 : : #include "complex.h"
47 : : #include "circuit.h"
48 : : #include "sweep.h"
49 : : #include "net.h"
50 : : #include "netdefs.h"
51 : : #include "analysis.h"
52 : : #include "nasolver.h"
53 : : #include "history.h"
54 : : #include "transient.h"
55 : : #include "exception.h"
56 : : #include "exceptionstack.h"
57 : : #include "environment.h"
58 : : #include "e_trsolver.h"
59 : : #include "component_id.h"
60 : : #include "ecvs.h"
61 : :
62 : : #define STEPDEBUG 0 // set to zero for release
63 : : #define BREAKPOINTS 0 // exact breakpoint calculation
64 : :
65 : : #ifndef dState
66 : : #define dState 0 // delta T state
67 : : #endif
68 : :
69 : : #ifndef sState
70 : : #define sState 1 // solution state
71 : : #endif
72 : :
73 : : // Macro for the n-th state of the solution vector history.
74 : : #ifndef SOL
75 : : #define SOL(state) (solution[(int) getState (sState, (state))])
76 : : #endif
77 : :
78 : : namespace qucs {
79 : :
80 : : using namespace transient;
81 : :
82 : : // Constructor creates an unnamed instance of the trsolver class.
83 : 0 : e_trsolver::e_trsolver ()
84 [ # # ][ # # ]: 0 : : trsolver ()
85 : : {
86 : : //swp = NULL;
87 : 0 : type = ANALYSIS_E_TRANSIENT;
88 : : //setDescription ("m-transient");
89 : : //for (int i = 0; i < 8; i++) solution[i] = NULL;
90 : : //tHistory = NULL;
91 : : //relaxTSR = false;
92 : : //initialDC = true;
93 : :
94 : : // initialise the message function pointer to
95 : : // point to the PrintWarningMsg function
96 : 0 : messagefcn = &logprint;
97 : : #if DEBUG
98 : : //loginit ();
99 : : // produce an actual log file
100 : : // file_error = file_status = fopen("e_trsolver.log", "w+");
101 : : // precinit ();
102 : : // ::srand (::time (NULL));
103 : : #endif // DEBUG
104 : 0 : }
105 : :
106 : : // Constructor creates a named instance of the e_trsolver class.
107 : 0 : e_trsolver::e_trsolver (char * n)
108 [ # # ][ # # ]: 0 : : trsolver (n)
109 : : {
110 : : //swp = NULL;
111 : 0 : type = ANALYSIS_E_TRANSIENT;
112 : : //setDescription ("m-transient");
113 : : //for (int i = 0; i < 8; i++) solution[i] = NULL;
114 : : //tHistory = NULL;
115 : : //relaxTSR = false;
116 : : //initialDC = true;
117 : :
118 : : // initialise the message function pointer to
119 : : // point to the PrintWarningMsg function
120 : 0 : messagefcn = &logprint;
121 : 0 : }
122 : :
123 : : // Destructor deletes the e_trsolver class object.
124 : 0 : e_trsolver::~e_trsolver ()
125 : : {
126 : :
127 [ # # ][ # # ]: 0 : solve_post ();
128 [ # # ][ # # ]: 0 : if (progress) logprogressclear (40);
[ # # ][ # # ]
129 : :
130 [ # # ][ # # ]: 0 : deinitTR ();
131 : :
132 [ # # ][ # # ]: 0 : if (swp) delete swp;
[ # # ][ # # ]
[ # # ][ # # ]
133 : :
134 [ # # ][ # # ]: 0 : for (int i = 0; i < 8; i++)
135 : : {
136 [ # # ][ # # ]: 0 : if (solution[i] != NULL)
137 : : {
138 [ # # ][ # # ]: 0 : delete solution[i];
139 : : }
140 [ # # ][ # # ]: 0 : if (lastsolution[i] != NULL)
141 : : {
142 [ # # ][ # # ]: 0 : delete lastsolution[i];
143 : : }
144 : : }
145 [ # # ][ # # ]: 0 : if (tHistory) delete tHistory;
[ # # ][ # # ]
146 : :
147 [ # # ][ # # ]: 0 : }
148 : :
149 : : /* The copy constructor creates a new instance of the e_trsolver class
150 : : based on the given e_trsolver object. */
151 : 0 : e_trsolver::e_trsolver (e_trsolver & o)
152 [ # # ][ # # ]: 0 : : trsolver (o)
153 : : {
154 [ # # ][ # # ]: 0 : swp = o.swp ? new sweep (*o.swp) : NULL;
[ # # ][ # # ]
[ # # ][ # # ]
155 [ # # ][ # # ]: 0 : for (int i = 0; i < 8; i++) solution[i] = NULL;
156 [ # # ][ # # ]: 0 : tHistory = o.tHistory ? new history (*o.tHistory) : NULL;
[ # # ][ # # ]
[ # # ][ # # ]
157 : 0 : relaxTSR = o.relaxTSR;
158 : 0 : initialDC = o.initialDC;
159 : 0 : }
160 : :
161 : 0 : void e_trsolver::debug()
162 : : {
163 : 0 : circuit * root = subnet->getRoot ();
164 : :
165 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
166 : : {
167 : 0 : messagefcn (0, c->getName() );
168 : :
169 [ # # ]: 0 : if (c->getSubcircuit ()) {
170 : 0 : printf ("subcircuit Name %s\n", c->getSubcircuit ());
171 : : }
172 : : }
173 : :
174 : : //netlist_list();
175 : 0 : }
176 : :
177 : :
178 : : /* This is the initialisation function for the slaved transient
179 : : netlist solver. It prepares the circuit for simulation. */
180 : 0 : int e_trsolver::init (nr_double_t start, nr_double_t firstdelta, int mode)
181 : : {
182 : : // run the equation solver for the netlist
183 : 0 : this->getEnv()->runSolver();
184 : :
185 : 0 : int error = 0;
186 : 0 : const char * const solver = getPropertyString ("Solver");
187 : 0 : relaxTSR = !strcmp (getPropertyString ("relaxTSR"), "yes") ? true : false;
188 : 0 : initialDC = !strcmp (getPropertyString ("initialDC"), "yes") ? true : false;
189 : : // fetch simulation properties
190 : 0 : MaxIterations = getPropertyInteger ("MaxIter");
191 : 0 : reltol = getPropertyDouble ("reltol");
192 : 0 : abstol = getPropertyDouble ("abstol");
193 : 0 : vntol = getPropertyDouble ("vntol");
194 : :
195 : 0 : runs++;
196 : 0 : saveCurrent = current = 0;
197 : 0 : stepDelta = -1;
198 : 0 : converged = 0;
199 : 0 : fixpoint = 0;
200 : 0 : lastsynctime = 0.0;
201 : 0 : statRejected = statSteps = statIterations = statConvergence = 0;
202 : :
203 : : // Choose a solver.
204 [ # # ]: 0 : if (!strcmp (solver, "CroutLU"))
205 : 0 : eqnAlgo = ALGO_LU_DECOMPOSITION;
206 [ # # ]: 0 : else if (!strcmp (solver, "DoolittleLU"))
207 : 0 : eqnAlgo = ALGO_LU_DECOMPOSITION_DOOLITTLE;
208 [ # # ]: 0 : else if (!strcmp (solver, "HouseholderQR"))
209 : 0 : eqnAlgo = ALGO_QR_DECOMPOSITION;
210 [ # # ]: 0 : else if (!strcmp (solver, "HouseholderLQ"))
211 : 0 : eqnAlgo = ALGO_QR_DECOMPOSITION_LS;
212 [ # # ]: 0 : else if (!strcmp (solver, "GolubSVD"))
213 : 0 : eqnAlgo = ALGO_SV_DECOMPOSITION;
214 : :
215 : : // Perform initial DC analysis.
216 [ # # ]: 0 : if (initialDC)
217 : : {
218 : 0 : error = dcAnalysis ();
219 [ # # ]: 0 : if (error)
220 : 0 : return -1;
221 : : }
222 : :
223 : : // Initialize transient analysis.
224 : 0 : setDescription ("transient");
225 : 0 : initETR (start, firstdelta, mode);
226 : 0 : setCalculation ((calculate_func_t) &calcTR);
227 : 0 : solve_pre ();
228 : :
229 : : // Recall the DC solution.
230 : 0 : recallSolution ();
231 : :
232 : : // Apply the nodesets and adjust previous solutions.
233 : 0 : applyNodeset (false);
234 : 0 : fillSolution (x);
235 : 0 : fillLastSolution (x);
236 : :
237 : : // Tell integrators to be initialized.
238 : 0 : setMode (MODE_INIT);
239 : :
240 : : // reset the circuit status to not running so the histories
241 : : // etc will be reinitialised on the first time step
242 : 0 : running = 0;
243 : 0 : rejected = 0;
244 [ # # ]: 0 : if (mode == ETR_MODE_ASYNC)
245 : : {
246 : 0 : delta /= 10;
247 : :
248 : : }
249 [ # # ]: 0 : else if (mode == ETR_MODE_SYNC)
250 : : {
251 : : //lastsynctime = start - delta;
252 : : }
253 : : else
254 : : {
255 [ # # ]: 0 : qucs::exception * e = new qucs::exception (EXCEPTION_UNKNOWN_ETR_MODE);
256 : 0 : e->setText ("Unknown ETR mode.");
257 : 0 : throw_exception (e);
258 : 0 : return -2;
259 : : }
260 : 0 : fillState (dState, delta);
261 : 0 : adjustOrder (1);
262 : :
263 : 0 : storeHistoryAges ();
264 : :
265 : 0 : return 0;
266 : :
267 : : }
268 : :
269 : : /* Stores all the initial history lengths requested by the circuit
270 : : elements for later use (to make sure we don't set the histories
271 : : to be less than these initial requested values) */
272 : 0 : void e_trsolver::storeHistoryAges (void)
273 : : {
274 : 0 : circuit * root = subnet->getRoot ();
275 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
276 : : {
277 : : // get the a
278 [ # # ]: 0 : if (c->hasHistory ())
279 : : {
280 [ # # ]: 0 : initialhistages.push_back (c->getHistoryAge ());
281 : : }
282 : : }
283 : 0 : }
284 : :
285 : 0 : void e_trsolver::fillLastSolution (tvector<nr_double_t> * s)
286 : : {
287 [ # # ]: 0 : for (int i = 0; i < 8; i++) * lastsolution[(int) getState (sState, (i))] = *s;
288 : 0 : }
289 : :
290 : : /* Goes through the list of circuit objects and runs its initTR()
291 : : function. */
292 : 0 : void e_trsolver::initETR (nr_double_t start, nr_double_t firstdelta, int mode)
293 : : {
294 : 0 : const char * const IMethod = getPropertyString ("IntegrationMethod");
295 : : //nr_double_t start = getPropertyDouble ("Start");
296 : : //nr_double_t stop = start + firstdelta;
297 : : //nr_double_t points = 1.0;
298 : :
299 : : // fetch corrector integration method and determine predicor method
300 : 0 : corrMaxOrder = getPropertyInteger ("Order");
301 : 0 : corrType = CMethod = correctorType (IMethod, corrMaxOrder);
302 : 0 : predType = PMethod = predictorType (CMethod, corrMaxOrder, predMaxOrder);
303 : 0 : corrOrder = corrMaxOrder;
304 : 0 : predOrder = predMaxOrder;
305 : :
306 : : // initialize step values
307 [ # # ]: 0 : if (mode == ETR_MODE_ASYNC){
308 : 0 : delta = getPropertyDouble ("InitialStep");
309 : 0 : deltaMin = getPropertyDouble ("MinStep");
310 : 0 : deltaMax = getPropertyDouble ("MaxStep");
311 [ # # ]: 0 : if (deltaMax == 0.0)
312 : 0 : deltaMax = firstdelta; // MIN ((stop - start) / (points - 1), stop / 200);
313 [ # # ]: 0 : if (deltaMin == 0.0)
314 : 0 : deltaMin = NR_TINY * 10 * deltaMax;
315 [ # # ]: 0 : if (delta == 0.0)
316 : 0 : delta = firstdelta; // MIN (stop / 200, deltaMax) / 10;
317 [ # # ]: 0 : if (delta < deltaMin) delta = deltaMin;
318 [ # # ]: 0 : if (delta > deltaMax) delta = deltaMax;
319 : : }
320 [ # # ]: 0 : else if (mode == ETR_MODE_SYNC)
321 : : {
322 : 0 : delta = firstdelta;
323 : 0 : deltaMin = NR_TINY * 10;
324 : 0 : deltaMax = std::numeric_limits<nr_double_t>::max() / 10;
325 : : }
326 : :
327 : : // initialize step history
328 : 0 : setStates (2);
329 : 0 : initStates ();
330 : : // initialise the history of states, setting them all to 'delta'
331 : 0 : fillState (dState, delta);
332 : :
333 : : // copy the initialised states to the 'deltas' array
334 : 0 : saveState (dState, deltas);
335 : : // copy the deltas to all the circuits
336 : 0 : setDelta ();
337 : : // set the initial corrector and predictor coefficients
338 : 0 : calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
339 : 0 : calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
340 : :
341 : : // initialize history of solution vectors (solutions)
342 [ # # ]: 0 : for (int i = 0; i < 8; i++)
343 : : {
344 : : // solution contains the last sets of node voltages and branch
345 : : // currents at each of the last 8 'deltas'.
346 : : // Note for convenience the definition:
347 : : // #define SOL(state) (solution[(int) getState (sState, (state))])
348 : : // is provided and used elsewhere to update the solutions
349 : 0 : solution[i] = new tvector<nr_double_t>;
350 : 0 : setState (sState, (nr_double_t) i, i);
351 : 0 : lastsolution[i] = new tvector<nr_double_t>;
352 : : }
353 : :
354 : : // Initialise history tracking for asynchronous solvers
355 : : // See acceptstep_async and rejectstep_async for more
356 : : // information
357 : 0 : lastasynctime = start;
358 : 0 : saveState (dState, lastdeltas);
359 : 0 : lastdelta = delta;
360 : :
361 : : // tell circuit elements about the transient analysis
362 : 0 : circuit *c, * root = subnet->getRoot ();
363 [ # # ]: 0 : for (c = root; c != NULL; c = (circuit *) c->getNext ())
364 : 0 : initCircuitTR (c);
365 : : // also initialize the created circuit elements
366 [ # # ]: 0 : for (c = root; c != NULL; c = (circuit *) c->getPrev ())
367 : 0 : initCircuitTR (c);
368 : 0 : }
369 : :
370 : 0 : void e_trsolver::printx()
371 : : {
372 : : char buf [1024];
373 : :
374 [ # # ]: 0 : for (int r = 0; r < x->getSize(); r++) {
375 : 0 : buf[0] = '\0';
376 : : //sprintf (buf, "%+.2e%+.2ei", (double) real (x->get (r)), (double) imag (x->get (r)));
377 : :
378 [ # # ]: 0 : if (r == 2)
379 : : {
380 : : // print just the currents
381 : : // SOL (0)->get->(r)
382 : : // solution[0]->get(r);
383 : : sprintf (buf, "%f\t%18.17f\t%6.5f %6.5f %6.5f %6.5f %6.5f %6.5f %6.5f %6.5f",
384 : : current,
385 : : (double) real (x->get (r)),
386 : 0 : solution[0]->get(r) ,
387 : 0 : solution[1]->get(r) ,
388 : 0 : solution[2]->get(r) ,
389 : 0 : solution[3]->get(r) ,
390 : 0 : solution[4]->get(r) ,
391 : 0 : solution[5]->get(r) ,
392 : 0 : solution[6]->get(r) ,
393 [ # # ][ # # ]: 0 : solution[7]->get(r) );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
394 : :
395 [ # # ]: 0 : messagefcn(0, buf);
396 : : }
397 : : }
398 : 0 : }
399 : :
400 : : /* synchronous step solver for external ode routine
401 : : *
402 : : * This function solves the circuit for a single time delta provided
403 : : * by an external source. Convergence issues etc. are expected to
404 : : * be handled by the external solver, as it is in full control of the
405 : : * time stepping.
406 : : */
407 : 0 : int e_trsolver::stepsolve_sync(nr_double_t synctime)
408 : : {
409 : :
410 : 0 : int error = 0;
411 : 0 : convError = 0;
412 : :
413 : 0 : time = synctime;
414 : : // update the interpolation time of any externally controlled
415 : : // components which require it.
416 : 0 : updateExternalInterpTime(time);
417 : :
418 : : // copy the externally chosen time step to delta
419 : 0 : delta = time - lastsynctime;
420 : :
421 : : // get the current solution time
422 : : //current += delta;
423 : :
424 : : // updates the integrator coefficients, and updates the array of prev
425 : : // 8 deltas with the new delta for this step
426 : 0 : updateCoefficients (delta);
427 : :
428 : : // Run predictor to get a start value for the solution vector for
429 : : // the successive iterative corrector process
430 : 0 : error += predictor ();
431 : :
432 : : // restart Newton iteration
433 : 0 : restart (); // restart non-linear devices
434 : :
435 : : // Attempt to solve the circuit with the given delta
436 : : try_running () // #defined as: do {
437 : : {
438 : : //error += solve_nonlinear_step ();
439 : 0 : error += corrector ();
440 : : }
441 [ # # ][ # # ]: 0 : catch_exception () // #defined as: } while (0); if (estack.top ()) switch (estack.top()->getCode ())
442 : : {
443 : : case EXCEPTION_NO_CONVERGENCE:
444 : 0 : pop_exception ();
445 : :
446 : : // Retry using damped Newton-Raphson.
447 : 0 : this->convHelper = CONV_SteepestDescent;
448 : 0 : convError = 2;
449 : : #if DEBUG
450 : : messagefcn (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
451 : 0 : "(no convergence)\n", (double) saveCurrent, (double) delta);
452 : : #endif
453 : :
454 : : try_running () // #defined as: do {
455 : : {
456 : : // error += solve_nonlinear_step ();
457 : 0 : error += solve_nonlinear ();
458 : : }
459 [ # # ][ # # ]: 0 : catch_exception () // #defined as: } while (0); if (estack.top ()) switch (estack.top()->getCode ())
460 : : {
461 : : case EXCEPTION_NO_CONVERGENCE:
462 : 0 : pop_exception ();
463 : :
464 : : // Update statistics.
465 : 0 : statRejected++;
466 : 0 : statConvergence++;
467 : 0 : rejected++;
468 : 0 : converged = 0;
469 : 0 : error = 0;
470 : :
471 : 0 : break;
472 : : default:
473 : : // Otherwise return.
474 : 0 : estack.print ();
475 : 0 : error++;
476 : 0 : break;
477 : : }
478 : :
479 : : // Update statistics and no more damped Newton-Raphson.
480 : : // statIterations += iterations;
481 : : // if (--convError < 0) this->convHelper = 0;
482 : :
483 : :
484 : 0 : break;
485 : : default:
486 : : // Otherwise return.
487 : 0 : estack.print ();
488 : 0 : error++;
489 : 0 : break;
490 : : }
491 : : // if there was an error other than non-convergence, return -1
492 [ # # ]: 0 : if (error) return -1;
493 : :
494 : : // check whether Jacobian matrix is still non-singular
495 [ # # ]: 0 : if (!A->isFinite ())
496 : : {
497 : : // messagefcn (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
498 : : // "aborting %s analysis\n", getName (), (double) current,
499 : : // getDescription ());
500 : 0 : return -1;
501 : : }
502 : :
503 : 0 : return 0;
504 : : }
505 : :
506 : : // accept a time step into the solution history
507 : 0 : void e_trsolver::acceptstep_sync()
508 : : {
509 : 0 : statIterations += iterations;
510 [ # # ]: 0 : if (--convError < 0) convHelper = 0;
511 : :
512 : : // Now advance in time or not...
513 [ # # ]: 0 : if (running > 1)
514 : : {
515 : 0 : adjustDelta_sync (current);
516 : : // deltaOld = delta;
517 : : // stepDelta = deltaOld;
518 : : // nextStates ();
519 : : // rejected = 0;
520 : 0 : adjustOrder ();
521 : : }
522 : : else
523 : : {
524 : 0 : fillStates ();
525 : 0 : nextStates ();
526 : 0 : rejected = 0;
527 : : }
528 : :
529 : 0 : saveCurrent = current;
530 : 0 : current += delta;
531 : 0 : running++;
532 : 0 : converged++;
533 : :
534 : : // Tell integrators to be running.
535 : 0 : setMode (MODE_NONE);
536 : :
537 : : // Initialize or update history.
538 [ # # ]: 0 : if (running > 1)
539 : : {
540 : : // update the solution history with the new results
541 : 0 : updateHistory (current);
542 : : }
543 : : else
544 : : {
545 : : // we have just solved the first transient state
546 : 0 : initHistory (current);
547 : : }
548 : :
549 : : // store the current time
550 : 0 : lastsynctime = current;
551 : 0 : }
552 : :
553 : : /* This function tries to adapt the current time-step according to the
554 : : global truncation error. */
555 : 0 : void e_trsolver::adjustDelta_sync (nr_double_t t)
556 : : {
557 : 0 : deltaOld = delta;
558 : :
559 : : // makes a new delta based on truncation error
560 : : // delta = checkDelta ();
561 : :
562 [ # # ]: 0 : if (delta > deltaMax)
563 : : {
564 : 0 : delta = deltaMax;
565 : : }
566 : :
567 [ # # ]: 0 : if (delta < deltaMin)
568 : : {
569 : 0 : delta = deltaMin;
570 : : }
571 : :
572 : : // delta correction in order to hit exact breakpoint
573 : 0 : int good = 0;
574 : : // if (!relaxTSR) // relaxed step raster?
575 : : // {
576 : : // if (!statConvergence || converged > 64) /* Is this a good guess? */
577 : : // {
578 : : // // check next breakpoint
579 : : // if (stepDelta > 0.0)
580 : : // {
581 : : // // restore last valid delta
582 : : // delta = stepDelta;
583 : : // stepDelta = -1.0;
584 : : // }
585 : : // else
586 : : // {
587 : : // if (delta > (t - current) && t > current)
588 : : // {
589 : : // // save last valid delta and set exact step
590 : : // stepDelta = deltaOld;
591 : : // delta = t - current;
592 : : // good = 1;
593 : : // }
594 : : // else
595 : : // {
596 : : // stepDelta = -1.0;
597 : : // }
598 : : // }
599 : : // if (delta > deltaMax) delta = deltaMax;
600 : : // if (delta < deltaMin) delta = deltaMin;
601 : : // }
602 : : // }
603 : :
604 : 0 : stepDelta = -1;
605 : 0 : good = 1;
606 : :
607 : : // usual delta correction
608 [ # # ][ # # ]: 0 : if (delta > 0.9 * deltaOld || good) // accept current delta
609 : : {
610 : 0 : nextStates ();
611 : 0 : rejected = 0;
612 : : #if STEPDEBUG
613 : : logprint (LOG_STATUS,
614 : : "DEBUG: delta accepted at t = %.3e, h = %.3e\n",
615 : : (double) current, (double) delta);
616 : : #endif
617 : : }
618 [ # # ]: 0 : else if (deltaOld > delta) // reject current delta
619 : : {
620 : 0 : rejected++;
621 : 0 : statRejected++;
622 : : #if STEPDEBUG
623 : : logprint (LOG_STATUS,
624 : : "DEBUG: delta rejected at t = %.3e, h = %.3e\n",
625 : : (double) current, (double) delta);
626 : : #endif
627 [ # # ]: 0 : if (current > 0) current -= deltaOld;
628 : : }
629 : : else
630 : : {
631 : 0 : nextStates ();
632 : 0 : rejected = 0;
633 : : }
634 : 0 : }
635 : :
636 : : // asynchronous step solver
637 : 0 : int e_trsolver::stepsolve_async(nr_double_t steptime)
638 : : {
639 : : // Start to sweep through time.
640 : 0 : int error = 0;
641 : 0 : convError = 0;
642 : :
643 : 0 : time = steptime;
644 : : // update the interpolation time of any externally controlled
645 : : // components which require it.
646 : 0 : updateExternalInterpTime(time);
647 : : // make the stored histories for all ircuits that have
648 : : // requested them at least as long as the next major time
649 : : // step so we can reject the step later if needed and
650 : : // restore all the histories to their previous state
651 : 0 : updateHistoryAges (time - lastasynctime);
652 : :
653 : : //delta = (steptime - time) / 10;
654 : : //if (progress) logprogressbar (i, swp->getSize (), 40);
655 : : #if DEBUG && 0
656 : : messagefcn (LOG_STATUS, "NOTIFY: %s: solving netlist for t = %e\n",
657 : : getName (), (double) time);
658 : : #endif
659 : :
660 [ # # ]: 0 : do
661 : : {
662 : : #if STEPDEBUG
663 : : if (delta == deltaMin)
664 : : {
665 : : messagefcn (LOG_ERROR,
666 : : "WARNING: %s: minimum delta h = %.3e at t = %.3e\n",
667 : : getName (), (double) delta, (double) current);
668 : : }
669 : : #endif
670 : : // update the integration coefficients
671 : 0 : updateCoefficients (delta);
672 : :
673 : : // Run predictor to get a start value for the solution vector for
674 : : // the successive iterative corrector process
675 : 0 : error += predictor ();
676 : :
677 : : // restart Newton iteration
678 [ # # ]: 0 : if (rejected)
679 : : {
680 : 0 : restart (); // restart non-linear devices
681 : 0 : rejected = 0;
682 : : }
683 : :
684 : : // Run corrector process with appropriate exception handling.
685 : : // The corrector iterates through the solutions of the integration
686 : : // process until a certain error tolerance has been reached.
687 : : try_running () // #defined as: do {
688 : : {
689 : 0 : error += corrector ();
690 : : }
691 [ # # ][ # # ]: 0 : catch_exception () // #defined as: } while (0); if (estack.top ()) switch (estack.top()->getCode ())
692 : : {
693 : : case EXCEPTION_NO_CONVERGENCE:
694 : 0 : pop_exception ();
695 : :
696 : : // Reduce step-size (by half) if failed to converge.
697 [ # # ]: 0 : if (current > 0) current -= delta;
698 : 0 : delta /= 2;
699 [ # # ]: 0 : if (delta <= deltaMin)
700 : : {
701 : 0 : delta = deltaMin;
702 : 0 : adjustOrder (1);
703 : : }
704 [ # # ]: 0 : if (current > 0) current += delta;
705 : :
706 : : // Update statistics.
707 : 0 : statRejected++;
708 : 0 : statConvergence++;
709 : 0 : rejected++;
710 : 0 : converged = 0;
711 : 0 : error = 0;
712 : :
713 : : // Start using damped Newton-Raphson.
714 : 0 : convHelper = CONV_SteepestDescent;
715 : 0 : convError = 2;
716 : : #if DEBUG
717 : : messagefcn (LOG_ERROR, "WARNING: delta rejected at t = %.3e, h = %.3e "
718 : 0 : "(no convergence)\n", (double) saveCurrent, (double) delta);
719 : : #endif
720 : 0 : break;
721 : : default:
722 : : // Otherwise return.
723 : 0 : estack.print ();
724 : 0 : error++;
725 : 0 : break;
726 : : }
727 [ # # ]: 0 : if (error) return -1;
728 [ # # ]: 0 : if (rejected) continue;
729 : :
730 : : // check whether Jacobian matrix is still non-singular
731 [ # # ]: 0 : if (!A->isFinite ())
732 : : {
733 : : messagefcn (LOG_ERROR, "ERROR: %s: Jacobian singular at t = %.3e, "
734 : : "aborting %s analysis\n", getName (), (double) current,
735 : 0 : getDescription ());
736 : 0 : return -1;
737 : : }
738 : :
739 : : // Update statistics and no more damped Newton-Raphson.
740 : 0 : statIterations += iterations;
741 [ # # ]: 0 : if (--convError < 0) convHelper = 0;
742 : :
743 : : // Now advance in time or not...
744 [ # # ]: 0 : if (running > 1)
745 : : {
746 : 0 : adjustDelta (time);
747 : 0 : adjustOrder ();
748 : : }
749 : : else
750 : : {
751 : 0 : fillStates ();
752 : 0 : nextStates ();
753 : 0 : rejected = 0;
754 : : }
755 : :
756 : 0 : saveCurrent = current;
757 : 0 : current += delta;
758 : 0 : running++;
759 : 0 : converged++;
760 : :
761 : : // Tell integrators to be running.
762 : 0 : setMode (MODE_NONE);
763 : :
764 : : // Initialize or update history.
765 [ # # ]: 0 : if (running > 1)
766 : : {
767 : 0 : updateHistory (saveCurrent);
768 : : }
769 : : else
770 : : {
771 : 0 : initHistory (saveCurrent);
772 : : }
773 : : }
774 : : while (saveCurrent < time); // Hit a requested time point?
775 : :
776 : 0 : return 0;
777 : : }
778 : :
779 : : // accept an asynchronous step into the solution history
780 : 0 : void e_trsolver::acceptstep_async(void)
781 : : {
782 : : // copy the solution in case we wish to step back to this
783 : : // point later
784 : 0 : copySolution (solution, lastsolution);
785 : :
786 : : // Store the time
787 : 0 : lastasynctime = time;
788 : :
789 : : // Store the deltas history
790 : 0 : saveState (dState, lastdeltas);
791 : :
792 : : // Store the current delta
793 : 0 : lastdelta = delta;
794 : 0 : }
795 : :
796 : : // reject the last asynchronous step and restore the history
797 : : // of solutions to the last major step
798 : 0 : void e_trsolver::rejectstep_async(void)
799 : : {
800 : : // restore the solution (node voltages and branch currents) from
801 : : // the previously stored solution
802 : 0 : copySolution (lastsolution, solution);
803 : :
804 : : // Restore the circuit histories to their previous states
805 : 0 : truncateHistory (lastasynctime);
806 : :
807 : : // Restore the time deltas
808 : 0 : inputState (dState, lastdeltas);
809 : :
810 [ # # ]: 0 : for (int i = 0; i < 8; i++)
811 : : {
812 : 0 : deltas[i] = lastdeltas[i];
813 : : }
814 : :
815 : 0 : delta = lastdelta;
816 : :
817 : : // copy the deltas to all the circuit elements
818 : 0 : setDelta ();
819 : :
820 : : // reset the corrector and predictor coefficients
821 : 0 : calcCorrectorCoeff (corrType, corrOrder, corrCoeff, deltas);
822 : 0 : calcPredictorCoeff (predType, predOrder, predCoeff, deltas);
823 : 0 : }
824 : :
825 : : /* Makes a copy of a set of solution vectors */
826 : 0 : void e_trsolver::copySolution (tvector<nr_double_t> * src[8], tvector<nr_double_t> * dest[8])
827 : : {
828 [ # # ]: 0 : for (int i = 0; i < 8; i++)
829 : : {
830 : : // check sizes are the same
831 [ # # ]: 0 : assert (src[i]->getSize () == dest[i]->getSize ());
832 : : // copy over the data values
833 [ # # ]: 0 : for (int j = 0; j < src[i]->getSize (); j++)
834 : : {
835 : 0 : dest[i]->set (j, src[i]->get (j));
836 : : }
837 : : }
838 : 0 : }
839 : :
840 : 0 : void e_trsolver::updateHistoryAges (nr_double_t newage)
841 : : {
842 : 0 : int i = 0;
843 : 0 : circuit * root = subnet->getRoot ();
844 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
845 : : {
846 : : // set the history length to retain to be at least
847 : : // the length of the supplied age
848 [ # # ]: 0 : if (c->hasHistory ())
849 : : {
850 : 0 : c->setHistoryAge (std::max (initialhistages[i], newage));
851 : 0 : i++;
852 : : }
853 : : }
854 : 0 : }
855 : :
856 : : //int e_trsolver::finish()
857 : : //{
858 : : // solve_post ();
859 : : //
860 : : // if (progress) logprogressclear (40);
861 : : //
862 : : // messagefcn (LOG_STATUS, "NOTIFY: %s: average time-step %g, %d rejections\n",
863 : : // getName (), (double) (saveCurrent / statSteps), statRejected);
864 : : // messagefcn (LOG_STATUS, "NOTIFY: %s: average NR-iterations %g, "
865 : : // "%d non-convergences\n", getName (),
866 : : // (double) statIterations / statSteps, statConvergence);
867 : : //
868 : : // // cleanup
869 : : // return 0;
870 : : //}
871 : :
872 : :
873 : 0 : void e_trsolver::getsolution (double * lastsol)
874 : : {
875 : 0 : int N = countNodes ();
876 : 0 : int M = countVoltageSources ();
877 : :
878 : : // copy solution
879 [ # # ]: 0 : for (int r = 0; r < N + M; r++)
880 : : {
881 : 0 : lastsol[r] = real(x->get(r));
882 : : }
883 : 0 : }
884 : :
885 : : /* getNodeV attempts to get the voltage of the node with a
886 : : given name. returns -1 if the node name was not found */
887 : 0 : int e_trsolver::getNodeV (char * label, nr_double_t& nodeV)
888 : : {
889 : 0 : int r = nlist->getNodeNr (label);
890 : :
891 [ # # ]: 0 : if (r == -1)
892 : : {
893 : 0 : return r;
894 : : }
895 : : else
896 : : {
897 : 0 : nodeV = x->get(r);
898 : 0 : return 0;
899 : : }
900 : : }
901 : :
902 : : /* Get the voltage reported by a voltage probe */
903 : 0 : int e_trsolver::getVProbeV (char * probename, nr_double_t& probeV)
904 : : {
905 : : // string to hold the full name of the circuit
906 [ # # ]: 0 : std::string fullname;
907 : :
908 : : // check for NULL name
909 [ # # ]: 0 : if (probename)
910 : : {
911 : 0 : circuit * root = subnet->getRoot ();
912 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
913 : : {
914 [ # # ]: 0 : if (c->getType () == CIR_VPROBE) {
915 : :
916 : 0 : fullname.clear ();
917 : :
918 : : // Subcircuit elements top level name is the
919 : : // subcircuit type (the base name of the subcircuit
920 : : // file)
921 [ # # ]: 0 : if (c->getSubcircuit ())
922 : : {
923 [ # # ]: 0 : fullname.append (c->getSubcircuit ());
924 [ # # ]: 0 : fullname.append (".");
925 : : }
926 : :
927 : : // append the user supplied name to search for
928 [ # # ]: 0 : fullname.append (probename);
929 : :
930 : : // Check if it is the desired voltage source
931 [ # # ]: 0 : if (strcmp (fullname.c_str(), c->getName ()) == 0)
932 : : {
933 : : // Saves the real and imaginary voltages in the probe to the
934 : : // named variables Vr and Vi
935 [ # # ]: 0 : c->saveOperatingPoints ();
936 : : // We are only interested in the real part for transient
937 : : // analysis
938 [ # # ]: 0 : probeV = c->getOperatingPoint ("Vr");
939 : :
940 : 0 : return 0;
941 : : }
942 : : }
943 : : }
944 : : }
945 : 0 : return -1;
946 : : }
947 : :
948 : : /* Get the current reported by a current probe */
949 : 0 : int e_trsolver::getIProbeI (char * probename, nr_double_t& probeI)
950 : : {
951 : : // string to hold the full name of the circuit
952 [ # # ]: 0 : std::string fullname;
953 : :
954 : : // check for NULL name
955 [ # # ]: 0 : if (probename)
956 : : {
957 : 0 : circuit * root = subnet->getRoot ();
958 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
959 : : {
960 [ # # ]: 0 : if (c->getType () == CIR_IPROBE) {
961 : :
962 : 0 : fullname.clear ();
963 : :
964 : : // Subcircuit elements top level name is the
965 : : // subcircuit type (the base name of the subcircuit
966 : : // file)
967 [ # # ]: 0 : if (c->getSubcircuit ())
968 : : {
969 [ # # ]: 0 : fullname.append (c->getSubcircuit ());
970 [ # # ]: 0 : fullname.append (".");
971 : : }
972 : :
973 : : // append the user supplied name to search for
974 [ # # ]: 0 : fullname.append (probename);
975 : :
976 : : // Check if it is the desired voltage source
977 [ # # ]: 0 : if (strcmp (fullname.c_str(), c->getName ()) == 0)
978 : : {
979 : : // Get the current reported by the probe
980 [ # # ][ # # ]: 0 : probeI = real (x->get (c->getVoltageSource () + getN ()));
[ # # ]
981 : :
982 : 0 : return 0;
983 : : }
984 : : }
985 : : }
986 : : }
987 : 0 : return -1;
988 : : }
989 : :
990 : 0 : int e_trsolver::setECVSVoltage(char * ecvsname, nr_double_t V)
991 : : {
992 : : // string to hold the full name of the circuit
993 [ # # ]: 0 : std::string fullname;
994 : :
995 : : // check for NULL name
996 [ # # ]: 0 : if (ecvsname)
997 : : {
998 : 0 : circuit * root = subnet->getRoot ();
999 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1000 : : {
1001 : : // examine only ECVS elements
1002 [ # # ]: 0 : if (c->getType () == CIR_ECVS) {
1003 : :
1004 : 0 : fullname.clear ();
1005 : :
1006 : : // Subcircuit elements top level name is the
1007 : : // subcircuit type (the base name of the subcircuit
1008 : : // file)
1009 [ # # ]: 0 : if (c->getSubcircuit ())
1010 : : {
1011 [ # # ]: 0 : fullname.append (c->getSubcircuit ());
1012 [ # # ]: 0 : fullname.append (".");
1013 : : }
1014 : :
1015 : : // append the user supplied name to search for
1016 [ # # ]: 0 : fullname.append (ecvsname);
1017 : :
1018 : : // Check if it is the desired voltage source
1019 [ # # ]: 0 : if (strcmp (fullname.c_str(), c->getName ()) == 0)
1020 : : {
1021 : : // Set the voltage to the desired value
1022 [ # # ]: 0 : c->setProperty("U", V);
1023 : 0 : return 0;
1024 : : }
1025 : : }
1026 : : }
1027 : : }
1028 : 0 : return -1;
1029 : : }
1030 : :
1031 : 0 : void e_trsolver::updateExternalInterpTime(nr_double_t t)
1032 : : {
1033 : 0 : circuit * root = subnet->getRoot ();
1034 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1035 : : {
1036 : : // examine only external elements that have interpolation,
1037 : : // such as ECVS elements
1038 [ # # ]: 0 : if (c->getType () == CIR_ECVS) {
1039 : : // Set the voltage to the desired value
1040 : 0 : c->setProperty ("Tnext", t);
1041 [ # # ][ # # ]: 0 : if (tHistory != NULL && tHistory->getSize () > 0)
[ # # ]
1042 : : {
1043 : 0 : c->setHistoryAge ( t - tHistory->last () + 0.1 * (t - tHistory->last ()) );
1044 : : }
1045 : : }
1046 : : }
1047 : 0 : }
1048 : :
1049 : : /* The following function removed stored times newer than a specified time
1050 : : from all the circuit element histories */
1051 : 0 : void e_trsolver::truncateHistory (nr_double_t t)
1052 : : {
1053 : : // truncate all the circuit element histories
1054 : 0 : circuit * root = subnet->getRoot ();
1055 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ())
1056 : : {
1057 [ # # ]: 0 : if (c->hasHistory ()) c->truncateHistory (t);
1058 : : }
1059 : 0 : }
1060 : :
1061 : 0 : int e_trsolver::getJacRows()
1062 : : {
1063 : 0 : return A->getRows();
1064 : : }
1065 : :
1066 : 0 : int e_trsolver::getJacCols()
1067 : : {
1068 : 0 : return A->getCols();
1069 : : }
1070 : :
1071 : 0 : void e_trsolver::getJacData(int r, int c, nr_double_t& data)
1072 : : {
1073 : 0 : data = A->get(r,c);
1074 : 0 : }
1075 : :
1076 : : // properties
1077 : : PROP_REQ [] =
1078 : : {
1079 : : PROP_NO_PROP
1080 : : };
1081 : : PROP_OPT [] =
1082 : : {
1083 : : {
1084 : : "IntegrationMethod", PROP_STR, { PROP_NO_VAL, "Trapezoidal" },
1085 : : PROP_RNG_STR4 ("Euler", "Trapezoidal", "Gear", "AdamsMoulton")
1086 : : },
1087 : : { "Order", PROP_INT, { 2, PROP_NO_STR }, PROP_RNGII (1, 6) },
1088 : : { "InitialStep", PROP_REAL, { 1e-9, PROP_NO_STR }, PROP_POS_RANGE },
1089 : : { "MinStep", PROP_REAL, { 1e-16, PROP_NO_STR }, PROP_POS_RANGE },
1090 : : { "MaxStep", PROP_REAL, { 0, PROP_NO_STR }, PROP_POS_RANGE },
1091 : : { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
1092 : : { "abstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
1093 : : { "vntol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
1094 : : { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
1095 : : { "LTEabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
1096 : : { "LTEreltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
1097 : : { "LTEfactor", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 16) },
1098 : : { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
1099 : : { "Solver", PROP_STR, { PROP_NO_VAL, "CroutLU" }, PROP_RNG_SOL },
1100 : : { "relaxTSR", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
1101 : : { "initialDC", PROP_STR, { PROP_NO_VAL, "yes" }, PROP_RNG_YESNO },
1102 : : PROP_NO_PROP
1103 : : };
1104 : : struct define_t e_trsolver::anadef =
1105 : : { "ETR", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
1106 : :
1107 : :
1108 : : } // namespace qucs
|
