Qucs-core
0.0.19
|
00001 /* 00002 * nasolver.cpp - nodal analysis solver class implementation 00003 * 00004 * Copyright (C) 2004, 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 // the types required for qucs library files are defined 00026 // in qucs_typedefs.h, created by configure from 00027 // qucs_typedefs.h.in 00028 #include "qucs_typedefs.h" 00029 00030 #include <stdio.h> 00031 #include <stdlib.h> 00032 #include <string.h> 00033 #include <cmath> 00034 #include <float.h> 00035 #include <assert.h> 00036 #include <limits> 00037 00038 #include "logging.h" 00039 #include "complex.h" 00040 #include "object.h" 00041 #include "node.h" 00042 #include "circuit.h" 00043 #include "vector.h" 00044 #include "dataset.h" 00045 #include "net.h" 00046 #include "analysis.h" 00047 #include "nodelist.h" 00048 #include "nodeset.h" 00049 #include "strlist.h" 00050 #include "tvector.h" 00051 #include "tmatrix.h" 00052 #include "eqnsys.h" 00053 #include "precision.h" 00054 #include "operatingpoint.h" 00055 #include "exception.h" 00056 #include "exceptionstack.h" 00057 #include "nasolver.h" 00058 #include "constants.h" 00059 00060 namespace qucs { 00061 00062 // Constructor creates an unnamed instance of the nasolver class. 00063 template <class nr_type_t> 00064 nasolver<nr_type_t>::nasolver () : analysis () 00065 { 00066 nlist = NULL; 00067 A = C = NULL; 00068 z = x = xprev = zprev = NULL; 00069 reltol = abstol = vntol = 0; 00070 calculate_func = NULL; 00071 convHelper = fixpoint = 0; 00072 eqnAlgo = ALGO_LU_DECOMPOSITION; 00073 updateMatrix = 1; 00074 gMin = srcFactor = 0; 00075 eqns = new eqnsys<nr_type_t> (); 00076 } 00077 00078 // Constructor creates a named instance of the nasolver class. 00079 template <class nr_type_t> 00080 nasolver<nr_type_t>::nasolver (const std::string &n) : analysis (n) 00081 { 00082 nlist = NULL; 00083 A = C = NULL; 00084 z = x = xprev = zprev = NULL; 00085 reltol = abstol = vntol = 0; 00086 calculate_func = NULL; 00087 convHelper = fixpoint = 0; 00088 eqnAlgo = ALGO_LU_DECOMPOSITION; 00089 updateMatrix = 1; 00090 gMin = srcFactor = 0; 00091 eqns = new eqnsys<nr_type_t> (); 00092 } 00093 00094 // Destructor deletes the nasolver class object. 00095 template <class nr_type_t> 00096 nasolver<nr_type_t>::~nasolver () 00097 { 00098 delete nlist; 00099 delete C; 00100 delete A; 00101 delete z; 00102 delete x; 00103 delete xprev; 00104 delete zprev; 00105 delete eqns; 00106 } 00107 00108 /* The copy constructor creates a new instance of the nasolver class 00109 based on the given nasolver object. */ 00110 template <class nr_type_t> 00111 nasolver<nr_type_t>::nasolver (nasolver & o) : analysis (o) 00112 { 00113 nlist = o.nlist ? new nodelist (*(o.nlist)) : NULL; 00114 A = o.A ? new tmatrix<nr_type_t> (*(o.A)) : NULL; 00115 C = o.C ? new tmatrix<nr_type_t> (*(o.C)) : NULL; 00116 z = o.z ? new tvector<nr_type_t> (*(o.z)) : NULL; 00117 x = o.x ? new tvector<nr_type_t> (*(o.x)) : NULL; 00118 xprev = zprev = NULL; 00119 reltol = o.reltol; 00120 abstol = o.abstol; 00121 vntol = o.vntol; 00122 desc = o.desc; 00123 calculate_func = o.calculate_func; 00124 convHelper = o.convHelper; 00125 eqnAlgo = o.eqnAlgo; 00126 updateMatrix = o.updateMatrix; 00127 fixpoint = o.fixpoint; 00128 gMin = o.gMin; 00129 srcFactor = o.srcFactor; 00130 eqns = new eqnsys<nr_type_t> (*(o.eqns)); 00131 solution = nasolution<nr_type_t> (o.solution); 00132 } 00133 00134 /* The function runs the nodal analysis solver once, reports errors if 00135 any and save the results into each circuit. */ 00136 template <class nr_type_t> 00137 int nasolver<nr_type_t>::solve_once (void) 00138 { 00139 qucs::exception * e; 00140 int error = 0, d; 00141 00142 // run the calculation function for each circuit 00143 calculate (); 00144 00145 // generate A matrix and z vector 00146 createMatrix (); 00147 00148 // solve equation system 00149 try_running () 00150 { 00151 runMNA (); 00152 } 00153 // appropriate exception handling 00154 catch_exception () 00155 { 00156 case EXCEPTION_PIVOT: 00157 case EXCEPTION_WRONG_VOLTAGE: 00158 e = new qucs::exception (EXCEPTION_NA_FAILED); 00159 d = top_exception()->getData (); 00160 pop_exception (); 00161 if (d >= countNodes ()) 00162 { 00163 d -= countNodes (); 00164 e->setText ("voltage source `%s' conflicts with some other voltage " 00165 "source", findVoltageSource(d)->getName ()); 00166 } 00167 else 00168 { 00169 e->setText ("circuit admittance matrix in %s solver is singular at " 00170 "node `%s' connected to [%s]", desc.c_str(), nlist->get (d).c_str(), 00171 nlist->getNodeString (d).c_str()); 00172 } 00173 throw_exception (e); 00174 error++; 00175 break; 00176 case EXCEPTION_SINGULAR: 00177 do 00178 { 00179 d = top_exception()->getData (); 00180 pop_exception (); 00181 if (d < countNodes ()) 00182 { 00183 logprint (LOG_ERROR, "WARNING: %s: inserted virtual resistance at " 00184 "node `%s' connected to [%s]\n", getName (), nlist->get (d).c_str(), 00185 nlist->getNodeString (d).c_str()); 00186 } 00187 } 00188 while (top_exception() != NULL && 00189 top_exception()->getCode () == EXCEPTION_SINGULAR); 00190 break; 00191 default: 00192 estack.print (); 00193 break; 00194 } 00195 00196 // save results into circuits 00197 if (!error) saveSolution (); 00198 return error; 00199 } 00200 00201 /* Run this function after the actual solver run and before evaluating 00202 the results. */ 00203 template <class nr_type_t> 00204 void nasolver<nr_type_t>::solve_post (void) 00205 { 00206 delete nlist; 00207 nlist = NULL; 00208 } 00209 00210 /* Run this function before the actual solver. */ 00211 template <class nr_type_t> 00212 void nasolver<nr_type_t>::solve_pre (void) 00213 { 00214 // create node list, enumerate nodes and voltage sources 00215 #if DEBUG 00216 logprint (LOG_STATUS, "NOTIFY: %s: creating node list for %s analysis\n", 00217 getName (), desc.c_str()); 00218 #endif 00219 nlist = new nodelist (subnet); 00220 nlist->assignNodes (); 00221 assignVoltageSources (); 00222 #if DEBUG && 0 00223 nlist->print (); 00224 #endif 00225 00226 // create matrix, solution vector and right hand side vector 00227 int M = countVoltageSources (); 00228 int N = countNodes (); 00229 delete A; 00230 A = new tmatrix<nr_type_t> (M + N); 00231 delete z; 00232 z = new tvector<nr_type_t> (N + M); 00233 delete x; 00234 x = new tvector<nr_type_t> (N + M); 00235 00236 #if DEBUG 00237 logprint (LOG_STATUS, "NOTIFY: %s: solving %s netlist\n", getName (), desc.c_str()); 00238 #endif 00239 } 00240 00241 /* This function goes through the nodeset list of the current netlist 00242 and applies the stored values to the current solution vector. Then 00243 the function saves the solution vector back into the actual 00244 component nodes. */ 00245 template <class nr_type_t> 00246 void nasolver<nr_type_t>::applyNodeset (bool nokeep) 00247 { 00248 if (x == NULL || nlist == NULL) return; 00249 00250 // set each solution to zero 00251 if (nokeep) for (int i = 0; i < x->size (); i++) x->set (i, 0); 00252 00253 // then apply the nodeset itself 00254 for (nodeset * n = subnet->getNodeset (); n; n = n->getNext ()) 00255 { 00256 struct nodelist_t * nl = nlist->getNode (n->getName ()); 00257 if (nl != NULL) 00258 { 00259 x->set (nl->n, n->getValue ()); 00260 } 00261 else 00262 { 00263 logprint (LOG_ERROR, "WARNING: %s: no such node `%s' found, cannot " 00264 "initialize node\n", getName (), n->getName ()); 00265 } 00266 } 00267 if (xprev != NULL) *xprev = *x; 00268 saveSolution (); 00269 } 00270 00271 /* The following function uses the gMin-stepping algorithm in order to 00272 solve the given non-linear netlist by continuous iterations. */ 00273 template <class nr_type_t> 00274 int nasolver<nr_type_t>::solve_nonlinear_continuation_gMin (void) 00275 { 00276 qucs::exception * e; 00277 int convergence, run = 0, MaxIterations, error = 0; 00278 nr_double_t gStep, gPrev; 00279 00280 // fetch simulation properties 00281 MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1; 00282 updateMatrix = 1; 00283 fixpoint = 0; 00284 00285 // initialize the stepper 00286 gPrev = gMin = 0.01; 00287 gStep = gMin / 100; 00288 gMin -= gStep; 00289 00290 do 00291 { 00292 // run solving loop until convergence is reached 00293 run = 0; 00294 do 00295 { 00296 error = solve_once (); 00297 if (!error) 00298 { 00299 // convergence check 00300 convergence = (run > 0) ? checkConvergence () : 0; 00301 savePreviousIteration (); 00302 run++; 00303 } 00304 else break; 00305 } 00306 while (!convergence && run < MaxIterations); 00307 iterations += run; 00308 00309 // not yet converged, so decreased the gMin-step 00310 if (run >= MaxIterations || error) 00311 { 00312 gStep /= 2; 00313 // here the absolute minimum step checker 00314 if (gStep < std::numeric_limits<nr_double_t>::epsilon()) 00315 { 00316 error = 1; 00317 e = new qucs::exception (EXCEPTION_NO_CONVERGENCE); 00318 e->setText ("no convergence in %s analysis after %d gMinStepping " 00319 "iterations", desc.c_str(), iterations); 00320 throw_exception (e); 00321 break; 00322 } 00323 gMin = MAX (gPrev - gStep, 0); 00324 } 00325 // converged, increased the gMin-step 00326 else 00327 { 00328 gPrev = gMin; 00329 gMin = MAX (gMin - gStep, 0); 00330 gStep *= 2; 00331 } 00332 } 00333 // continue until no additional resistances is necessary 00334 while (gPrev > 0); 00335 00336 return error; 00337 } 00338 00339 /* The following function uses the source-stepping algorithm in order 00340 to solve the given non-linear netlist by continuous iterations. */ 00341 template <class nr_type_t> 00342 int nasolver<nr_type_t>::solve_nonlinear_continuation_Source (void) 00343 { 00344 qucs::exception * e; 00345 int convergence, run = 0, MaxIterations, error = 0; 00346 nr_double_t sStep, sPrev; 00347 00348 // fetch simulation properties 00349 MaxIterations = getPropertyInteger ("MaxIter") / 4 + 1; 00350 updateMatrix = 1; 00351 fixpoint = 0; 00352 00353 // initialize the stepper 00354 sPrev = srcFactor = 0; 00355 sStep = 0.01; 00356 srcFactor += sStep; 00357 00358 do 00359 { 00360 // run solving loop until convergence is reached 00361 run = 0; 00362 do 00363 { 00364 subnet->setSrcFactor (srcFactor); 00365 error = solve_once (); 00366 if (!error) 00367 { 00368 // convergence check 00369 convergence = (run > 0) ? checkConvergence () : 0; 00370 savePreviousIteration (); 00371 run++; 00372 } 00373 else break; 00374 } 00375 while (!convergence && run < MaxIterations); 00376 iterations += run; 00377 00378 // not yet converged, so decreased the source-step 00379 if (run >= MaxIterations || error) 00380 { 00381 if (error) 00382 sStep *= 0.1; 00383 else 00384 sStep *= 0.5; 00385 restorePreviousIteration (); 00386 saveSolution (); 00387 // here the absolute minimum step checker 00388 if (sStep < std::numeric_limits<nr_double_t>::epsilon()) 00389 { 00390 error = 1; 00391 e = new qucs::exception (EXCEPTION_NO_CONVERGENCE); 00392 e->setText ("no convergence in %s analysis after %d sourceStepping " 00393 "iterations", desc.c_str(), iterations); 00394 throw_exception (e); 00395 break; 00396 } 00397 srcFactor = std::min (sPrev + sStep, 1.0); 00398 } 00399 // converged, increased the source-step 00400 else if (run < MaxIterations / 4) 00401 { 00402 sPrev = srcFactor; 00403 srcFactor = std::min (srcFactor + sStep, 1.0); 00404 sStep *= 1.5; 00405 } 00406 else 00407 { 00408 srcFactor = std::min (srcFactor + sStep, 1.0); 00409 } 00410 } 00411 // continue until no source factor is necessary 00412 while (sPrev < 1); 00413 00414 subnet->setSrcFactor (1); 00415 return error; 00416 } 00417 00418 /* The function returns an appropriate text representation for the 00419 currently used convergence helper algorithm. */ 00420 template <class nr_type_t> 00421 const char * nasolver<nr_type_t>::getHelperDescription (void) 00422 { 00423 if (convHelper == CONV_Attenuation) 00424 { 00425 return "RHS attenuation"; 00426 } 00427 else if (convHelper == CONV_LineSearch) 00428 { 00429 return "line search"; 00430 } 00431 else if (convHelper == CONV_SteepestDescent) 00432 { 00433 return "steepest descent"; 00434 } 00435 else if (convHelper == CONV_GMinStepping) 00436 { 00437 return "gMin stepping"; 00438 } 00439 else if (convHelper == CONV_SourceStepping) 00440 { 00441 return "source stepping"; 00442 } 00443 return "none"; 00444 } 00445 00446 /* This is the non-linear iterative nodal analysis netlist solver. */ 00447 template <class nr_type_t> 00448 int nasolver<nr_type_t>::solve_nonlinear (void) 00449 { 00450 qucs::exception * e; 00451 int convergence, run = 0, MaxIterations, error = 0; 00452 00453 // fetch simulation properties 00454 MaxIterations = getPropertyInteger ("MaxIter"); 00455 reltol = getPropertyDouble ("reltol"); 00456 abstol = getPropertyDouble ("abstol"); 00457 vntol = getPropertyDouble ("vntol"); 00458 updateMatrix = 1; 00459 00460 if (convHelper == CONV_GMinStepping) 00461 { 00462 // use the alternative non-linear solver solve_nonlinear_continuation_gMin 00463 // instead of the basic solver provided by this function 00464 iterations = 0; 00465 error = solve_nonlinear_continuation_gMin (); 00466 return error; 00467 } 00468 else if (convHelper == CONV_SourceStepping) 00469 { 00470 // use the alternative non-linear solver solve_nonlinear_continuation_Source 00471 // instead of the basic solver provided by this function 00472 iterations = 0; 00473 error = solve_nonlinear_continuation_Source (); 00474 return error; 00475 } 00476 00477 // run solving loop until convergence is reached 00478 do 00479 { 00480 error = solve_once (); 00481 if (!error) 00482 { 00483 // convergence check 00484 convergence = (run > 0) ? checkConvergence () : 0; 00485 savePreviousIteration (); 00486 run++; 00487 // control fixpoint iterations 00488 if (fixpoint) 00489 { 00490 if (convergence && !updateMatrix) 00491 { 00492 updateMatrix = 1; 00493 convergence = 0; 00494 } 00495 else 00496 { 00497 updateMatrix = 0; 00498 } 00499 } 00500 } 00501 else 00502 { 00503 break; 00504 } 00505 } 00506 while (!convergence && 00507 run < MaxIterations * (1 + convHelper ? 1 : 0)); 00508 00509 if (run >= MaxIterations || error) 00510 { 00511 e = new qucs::exception (EXCEPTION_NO_CONVERGENCE); 00512 e->setText ("no convergence in %s analysis after %d iterations", 00513 desc.c_str(), run); 00514 throw_exception (e); 00515 error++; 00516 } 00517 00518 iterations = run; 00519 return error; 00520 } 00521 00522 /* This is the linear nodal analysis netlist solver. */ 00523 template <class nr_type_t> 00524 int nasolver<nr_type_t>::solve_linear (void) 00525 { 00526 updateMatrix = 1; 00527 return solve_once (); 00528 } 00529 00530 /* Applying the MNA (Modified Nodal Analysis) to a circuit with 00531 passive elements and independent current and voltage sources 00532 results in a matrix equation of the form Ax = z. This function 00533 generates the A and z matrix. */ 00534 template <class nr_type_t> 00535 void nasolver<nr_type_t>::createMatrix (void) 00536 { 00537 00538 /* Generate the A matrix. The A matrix consists of four (4) minor 00539 matrices in the form +- -+ 00540 A = | G B | 00541 | C D | 00542 +- -+. 00543 Each of these minor matrices is going to be generated here. */ 00544 if (updateMatrix) 00545 { 00546 createGMatrix (); 00547 createBMatrix (); 00548 createCMatrix (); 00549 createDMatrix (); 00550 } 00551 00552 /* Adjust G matrix if requested. */ 00553 if (convHelper == CONV_GMinStepping) 00554 { 00555 int N = countNodes (); 00556 int M = countVoltageSources (); 00557 for (int n = 0; n < N + M; n++) 00558 { 00559 A->set (n, n, A->get (n, n) + gMin); 00560 } 00561 } 00562 00563 /* Generate the z Matrix. The z Matrix consists of two (2) minor 00564 matrices in the form +- -+ 00565 z = | i | 00566 | e | 00567 +- -+. 00568 Each of these minor matrices is going to be generated here. */ 00569 createZVector (); 00570 } 00571 00572 /* This MatVal() functionality is just helper to get the correct 00573 values from the circuit's matrices. The additional (unused) 00574 argument is used to differentiate between the two possible 00575 types. */ 00576 #define MatVal(x) MatValX (x, (nr_type_t *) 0) 00577 00578 template <class nr_type_t> 00579 nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_complex_t *) 00580 { 00581 return z; 00582 } 00583 00584 template <class nr_type_t> 00585 nr_type_t nasolver<nr_type_t>::MatValX (nr_complex_t z, nr_double_t *) 00586 { 00587 return real (z); 00588 } 00589 00590 /* The B matrix is an MxN matrix with only 0, 1 and -1 elements. Each 00591 location in the matrix corresponds to a particular voltage source 00592 (first dimension) or a node (second dimension). If the positive 00593 terminal of the ith voltage source is connected to node k, then the 00594 element (i,k) in the B matrix is a 1. If the negative terminal of 00595 the ith voltage source is connected to node k, then the element 00596 (i,k) in the B matrix is a -1. Otherwise, elements of the B matrix 00597 are zero. */ 00598 template <class nr_type_t> 00599 void nasolver<nr_type_t>::createBMatrix (void) 00600 { 00601 int N = countNodes (); 00602 int M = countVoltageSources (); 00603 circuit * vs; 00604 struct nodelist_t * n; 00605 nr_type_t val; 00606 00607 // go through each voltage sources (first dimension) 00608 for (int c = 0; c < M; c++) 00609 { 00610 vs = findVoltageSource (c); 00611 // go through each node (second dimension) 00612 for (int r = 0; r < N; r++) 00613 { 00614 val = 0.0; 00615 n = nlist->getNode (r); 00616 for (auto ¤t : *n) 00617 { 00618 // is voltage source connected to node ? 00619 if (current->getCircuit () == vs) 00620 { 00621 val += MatVal (vs->getB (current->getPort (), c)); 00622 } 00623 } 00624 // put value into B matrix 00625 A->set (r, c + N, val); 00626 } 00627 } 00628 } 00629 00630 /* The C matrix is an NxM matrix with only 0, 1 and -1 elements. Each 00631 location in the matrix corresponds to a particular node (first 00632 dimension) or a voltage source (first dimension). If the positive 00633 terminal of the ith voltage source is connected to node k, then the 00634 element (k,i) in the C matrix is a 1. If the negative terminal of 00635 the ith voltage source is connected to node k, then the element 00636 (k,i) in the C matrix is a -1. Otherwise, elements of the C matrix 00637 are zero. */ 00638 template <class nr_type_t> 00639 void nasolver<nr_type_t>::createCMatrix (void) 00640 { 00641 int N = countNodes (); 00642 int M = countVoltageSources (); 00643 circuit * vs; 00644 struct nodelist_t * n; 00645 nr_type_t val; 00646 00647 // go through each voltage sources (second dimension) 00648 for (int r = 0; r < M; r++) 00649 { 00650 vs = findVoltageSource (r); 00651 // go through each node (first dimension) 00652 for (int c = 0; c < N; c++) 00653 { 00654 val = 0.0; 00655 n = nlist->getNode (c); 00656 for (auto ¤t: *n) 00657 { 00658 // is voltage source connected to node ? 00659 if (current->getCircuit () == vs) 00660 { 00661 val += MatVal (vs->getC (r, current->getPort ())); 00662 } 00663 } 00664 // put value into C matrix 00665 A->set (r + N, c, val); 00666 } 00667 } 00668 } 00669 00670 /* The D matrix is an MxM matrix that is composed entirely of zeros. 00671 It can be non-zero if dependent sources are considered. */ 00672 template <class nr_type_t> 00673 void nasolver<nr_type_t>::createDMatrix (void) 00674 { 00675 int M = countVoltageSources (); 00676 int N = countNodes (); 00677 circuit * vsr, * vsc; 00678 nr_type_t val; 00679 for (int r = 0; r < M; r++) 00680 { 00681 vsr = findVoltageSource (r); 00682 for (int c = 0; c < M; c++) 00683 { 00684 vsc = findVoltageSource (c); 00685 val = 0.0; 00686 if (vsr == vsc) 00687 { 00688 val = MatVal (vsr->getD (r, c)); 00689 } 00690 A->set (r + N, c + N, val); 00691 } 00692 } 00693 } 00694 00695 /* The G matrix is an NxN matrix formed in two steps. 00696 1. Each element in the diagonal matrix is equal to the sum of the 00697 conductance of each element connected to the corresponding node. 00698 2. The off diagonal elements are the negative conductance of the 00699 element connected to the pair of corresponding nodes. Therefore a 00700 resistor between nodes 1 and 2 goes into the G matrix at location 00701 (1,2) and location (2,1). If an element is grounded, it will only 00702 have contribute to one entry in the G matrix -- at the appropriate 00703 location on the diagonal. */ 00704 template <class nr_type_t> 00705 void nasolver<nr_type_t>::createGMatrix (void) 00706 { 00707 int pr, pc, N = countNodes (); 00708 nr_type_t g; 00709 struct nodelist_t * nr, * nc; 00710 circuit * ct; 00711 00712 // go through each column of the G matrix 00713 for (int c = 0; c < N; c++) 00714 { 00715 nc = nlist->getNode (c); 00716 // go through each row of the G matrix 00717 for (int r = 0; r < N; r++) 00718 { 00719 nr = nlist->getNode (r); 00720 g = 0.0; 00721 // sum up the conductance of each connected circuit 00722 for (auto & currentnc : *nc) 00723 for (auto & currentnr: *nr) 00724 if (currentnc->getCircuit () == currentnr->getCircuit ()) 00725 { 00726 ct = currentnc->getCircuit (); 00727 pc = currentnc->getPort (); 00728 pr = currentnr->getPort (); 00729 g += MatVal (ct->getY (pr, pc)); 00730 } 00731 // put value into G matrix 00732 A->set (r, c, g); 00733 } 00734 } 00735 } 00736 00737 /* The following function creates the (N+M)x(N+M) noise current 00738 correlation matrix used during the AC noise computations. */ 00739 template <class nr_type_t> 00740 void nasolver<nr_type_t>::createNoiseMatrix (void) 00741 { 00742 int pr, pc, N = countNodes (); 00743 int M = countVoltageSources (); 00744 struct nodelist_t * n; 00745 nr_type_t val; 00746 int r, c, ri, ci; 00747 struct nodelist_t * nr, * nc; 00748 circuit * ct; 00749 00750 // create new Cy matrix if necessary 00751 delete C; 00752 C = new tmatrix<nr_type_t> (N + M); 00753 00754 // go through each column of the Cy matrix 00755 for (c = 0; c < N; c++) 00756 { 00757 nc = nlist->getNode (c); 00758 // go through each row of the Cy matrix 00759 for (r = 0; r < N; r++) 00760 { 00761 nr = nlist->getNode (r); 00762 val = 0.0; 00763 // sum up the noise-correlation of each connected circuit 00764 for (auto & currentnc: *nc) 00765 /* a = 0; a < nc->size(); a++ */ 00766 for (auto ¤tnr : *nr) 00767 /* b = 0; b < nr->size(); b++) */ 00768 if (currentnc->getCircuit () == currentnr->getCircuit ()) 00769 { 00770 ct = currentnc->getCircuit (); 00771 pc = currentnc->getPort (); 00772 pr = currentnr->getPort (); 00773 val += MatVal (ct->getN (pr, pc)); 00774 } 00775 // put value into Cy matrix 00776 C->set (r, c, val); 00777 } 00778 } 00779 00780 // go through each additional voltage source and put coefficients into 00781 // the noise current correlation matrix 00782 circuit * vsr, * vsc; 00783 for (r = 0; r < M; r++) 00784 { 00785 vsr = findVoltageSource (r); 00786 for (c = 0; c < M; c++) 00787 { 00788 vsc = findVoltageSource (c); 00789 val = 0.0; 00790 if (vsr == vsc) 00791 { 00792 ri = vsr->getSize () + r - vsr->getVoltageSource (); 00793 ci = vsc->getSize () + c - vsc->getVoltageSource (); 00794 val = MatVal (vsr->getN (ri, ci)); 00795 } 00796 C->set (r + N, c + N, val); 00797 } 00798 } 00799 00800 // go through each additional voltage source 00801 for (r = 0; r < M; r++) 00802 { 00803 vsr = findVoltageSource (r); 00804 // go through each node 00805 for (c = 0; c < N; c++) 00806 { 00807 val = 0.0; 00808 n = nlist->getNode (c); 00809 for (auto ¤tn: *n) 00810 /*i = 0; i < n->size(); i++ )*/ 00811 { 00812 // is voltage source connected to node ? 00813 if (currentn->getCircuit () == vsr) 00814 { 00815 ri = vsr->getSize () + r - vsr->getVoltageSource (); 00816 ci = currentn->getPort (); 00817 val += MatVal (vsr->getN (ri, ci)); 00818 } 00819 } 00820 // put value into Cy matrix 00821 C->set (r + N, c, val); 00822 } 00823 } 00824 00825 // go through each voltage source 00826 for (c = 0; c < M; c++) 00827 { 00828 vsc = findVoltageSource (c); 00829 // go through each node 00830 for (r = 0; r < N; r++) 00831 { 00832 val = 0.0; 00833 n = nlist->getNode (r); 00834 for (auto & currentn: *n)/*i = 0; i < n->size(); i++)*/ 00835 { 00836 // is voltage source connected to node ? 00837 if (currentn->getCircuit () == vsc) 00838 { 00839 ci = vsc->getSize () + c - vsc->getVoltageSource (); 00840 ri = currentn->getPort (); 00841 val += MatVal (vsc->getN (ri, ci)); 00842 } 00843 } 00844 // put value into Cy matrix 00845 C->set (r, c + N, val); 00846 } 00847 } 00848 00849 } 00850 00851 /* The i matrix is an 1xN matrix with each element of the matrix 00852 corresponding to a particular node. The value of each element of i 00853 is determined by the sum of current sources into the corresponding 00854 node. If there are no current sources connected to the node, the 00855 value is zero. */ 00856 template <class nr_type_t> 00857 void nasolver<nr_type_t>::createIVector (void) 00858 { 00859 int N = countNodes (); 00860 nr_type_t val; 00861 struct nodelist_t * n; 00862 circuit * is; 00863 00864 // go through each node 00865 for (int r = 0; r < N; r++) 00866 { 00867 val = 0.0; 00868 n = nlist->getNode (r); 00869 // go through each circuit connected to the node 00870 for (auto ¤tn: *n)/* int i = 0; i < n->size(); i++)*/ 00871 { 00872 is = currentn->getCircuit (); 00873 // is this a current source ? 00874 if (is->isISource () || is->isNonLinear ()) 00875 { 00876 val += MatVal (is->getI (currentn->getPort ())); 00877 } 00878 } 00879 // put value into i vector 00880 z->set (r, val); 00881 } 00882 } 00883 00884 /* The e matrix is a 1xM matrix with each element of the matrix equal 00885 in value to the corresponding independent voltage source. */ 00886 template <class nr_type_t> 00887 void nasolver<nr_type_t>::createEVector (void) 00888 { 00889 int N = countNodes (); 00890 int M = countVoltageSources (); 00891 nr_type_t val; 00892 circuit * vs; 00893 00894 // go through each voltage source 00895 for (int r = 0; r < M; r++) 00896 { 00897 vs = findVoltageSource (r); 00898 val = MatVal (vs->getE (r)); 00899 // put value into e vector 00900 z->set (r + N, val); 00901 } 00902 } 00903 00904 // The function loads the right hand side vector. 00905 template <class nr_type_t> 00906 void nasolver<nr_type_t>::createZVector (void) 00907 { 00908 createIVector (); 00909 createEVector (); 00910 } 00911 00912 // Returns the number of nodes in the nodelist, excluding the ground node. 00913 template <class nr_type_t> 00914 int nasolver<nr_type_t>::countNodes (void) 00915 { 00916 return nlist->length () - 1; 00917 } 00918 00919 // Returns the node number of the give node name. 00920 template <class nr_type_t> 00921 int nasolver<nr_type_t>::getNodeNr (const std::string &str) 00922 { 00923 return nlist->getNodeNr (str); 00924 } 00925 00926 /* The function returns the assigned node number for the port of the 00927 given circuits. It returns -1 if there is no such node. */ 00928 template <class nr_type_t> 00929 int nasolver<nr_type_t>::findAssignedNode (circuit * c, int port) 00930 { 00931 int N = countNodes (); 00932 struct nodelist_t * n; 00933 for (int r = 0; r < N; r++) 00934 { 00935 n = nlist->getNode (r); 00936 for (auto ¤tn : *n) /*int i = 0; i < n->size(); i++)*/ 00937 if (c == currentn->getCircuit ()) 00938 if (port == currentn->getPort ()) 00939 return r; 00940 } 00941 return -1; 00942 } 00943 00944 // Returns the number of voltage sources in the nodelist. 00945 template <class nr_type_t> 00946 int nasolver<nr_type_t>::countVoltageSources (void) 00947 { 00948 return subnet->getVoltageSources (); 00949 } 00950 00951 /* The function returns the voltage source circuit object 00952 corresponding to the given number. If there is no such voltage 00953 source it returns NULL. */ 00954 template <class nr_type_t> 00955 circuit * nasolver<nr_type_t>::findVoltageSource (int n) 00956 { 00957 circuit * root = subnet->getRoot (); 00958 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) 00959 { 00960 if (n >= c->getVoltageSource () && 00961 n <= c->getVoltageSource () + c->getVoltageSources () - 1) 00962 return c; 00963 } 00964 return NULL; 00965 } 00966 00967 /* The function applies unique voltage source identifiers to each 00968 voltage source (explicit and built in internal ones) in the list of 00969 registered circuits. */ 00970 template <class nr_type_t> 00971 void nasolver<nr_type_t>::assignVoltageSources (void) 00972 { 00973 circuit * root = subnet->getRoot (); 00974 int nSources = 0; 00975 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) 00976 { 00977 if (c->getVoltageSources () > 0) 00978 { 00979 c->setVoltageSource (nSources); 00980 nSources += c->getVoltageSources (); 00981 } 00982 } 00983 subnet->setVoltageSources (nSources); 00984 } 00985 00986 /* The matrix equation Ax = z is solved by x = A^-1*z. The function 00987 applies the operation to the previously generated matrices. */ 00988 template <class nr_type_t> 00989 void nasolver<nr_type_t>::runMNA (void) 00990 { 00991 00992 // just solve the equation system here 00993 eqns->setAlgo (eqnAlgo); 00994 eqns->passEquationSys (updateMatrix ? A : NULL, x, z); 00995 eqns->solve (); 00996 00997 // if damped Newton-Raphson is requested 00998 if (xprev != NULL && top_exception () == NULL) 00999 { 01000 if (convHelper == CONV_Attenuation) 01001 { 01002 applyAttenuation (); 01003 } 01004 else if (convHelper == CONV_LineSearch) 01005 { 01006 lineSearch (); 01007 } 01008 else if (convHelper == CONV_SteepestDescent) 01009 { 01010 steepestDescent (); 01011 } 01012 } 01013 } 01014 01015 /* This function applies a damped Newton-Raphson (limiting scheme) to 01016 the current solution vector in the form x1 = x0 + a * (x1 - x0). This 01017 convergence helper is heuristic and does not ensure global convergence. */ 01018 template <class nr_type_t> 01019 void nasolver<nr_type_t>::applyAttenuation (void) 01020 { 01021 nr_double_t alpha = 1.0, nMax; 01022 01023 // create solution difference vector and find maximum deviation 01024 tvector<nr_type_t> dx = *x - *xprev; 01025 nMax = maxnorm (dx); 01026 01027 // compute appropriate damping factor 01028 if (nMax > 0.0) 01029 { 01030 nr_double_t g = 1.0; 01031 alpha = std::min (0.9, g / nMax); 01032 if (alpha < 0.1) alpha = 0.1; 01033 } 01034 01035 // apply damped solution vector 01036 *x = *xprev + alpha * dx; 01037 } 01038 01039 /* This is damped Newton-Raphson using nested iterations in order to 01040 find a better damping factor. It identifies a damping factor in 01041 the interval [0,1] which minimizes the right hand side vector. The 01042 algorithm actually ensures global convergence but pushes the 01043 solution to local minimums, i.e. where the Jacobian matrix A may be 01044 singular. */ 01045 template <class nr_type_t> 01046 void nasolver<nr_type_t>::lineSearch (void) 01047 { 01048 nr_double_t alpha = 0.5, n, nMin, aprev = 1.0, astep = 0.5, adiff; 01049 int dir = -1; 01050 01051 // compute solution deviation vector 01052 tvector<nr_type_t> dx = *x - *xprev; 01053 nMin = std::numeric_limits<nr_double_t>::max(); 01054 01055 do 01056 { 01057 // apply current damping factor and see what happens 01058 *x = *xprev + alpha * dx; 01059 01060 // recalculate Jacobian and right hand side 01061 saveSolution (); 01062 calculate (); 01063 createZVector (); 01064 01065 // calculate norm of right hand side vector 01066 n = norm (*z); 01067 01068 // TODO: this is not perfect, but usable 01069 astep /= 2; 01070 adiff = fabs (alpha - aprev); 01071 if (adiff > 0.005) 01072 { 01073 aprev = alpha; 01074 if (n < nMin) 01075 { 01076 nMin = n; 01077 if (alpha == 1) dir = -dir; 01078 alpha += astep * dir; 01079 } 01080 else 01081 { 01082 dir = -dir; 01083 alpha += 1.5 * astep * dir; 01084 } 01085 } 01086 } 01087 while (adiff > 0.005); 01088 01089 // apply final damping factor 01090 assert (alpha > 0 && alpha <= 1); 01091 *x = *xprev + alpha * dx; 01092 } 01093 01094 /* The function looks for the optimal gradient for the right hand side 01095 vector using the so-called 'steepest descent' method. Though 01096 better than the one-dimensional linesearch (it doesn't push 01097 iterations into local minimums) it converges painfully slow. */ 01098 template <class nr_type_t> 01099 void nasolver<nr_type_t>::steepestDescent (void) 01100 { 01101 nr_double_t alpha = 1.0, sl, n; 01102 01103 // compute solution deviation vector 01104 tvector<nr_type_t> dx = *x - *xprev; 01105 tvector<nr_type_t> dz = *z - *zprev; 01106 n = norm (*zprev); 01107 01108 do 01109 { 01110 // apply current damping factor and see what happens 01111 *x = *xprev + alpha * dx; 01112 01113 // recalculate Jacobian and right hand side 01114 saveSolution (); 01115 calculate (); 01116 createZVector (); 01117 01118 // check gradient criteria, ThinkME: Is this correct? 01119 dz = *z - *zprev; 01120 sl = real (sum (dz * -dz)); 01121 if (norm (*z) < n + alpha * sl) break; 01122 alpha *= 0.7; 01123 } 01124 while (alpha > 0.001); 01125 01126 // apply final damping factor 01127 *x = *xprev + alpha * dx; 01128 } 01129 01130 /* The function checks whether the iterative algorithm for linearizing 01131 the non-linear components in the network shows convergence. It 01132 returns non-zero if it converges and zero otherwise. */ 01133 template <class nr_type_t> 01134 int nasolver<nr_type_t>::checkConvergence (void) 01135 { 01136 01137 int N = countNodes (); 01138 int M = countVoltageSources (); 01139 nr_double_t v_abs, v_rel, i_abs, i_rel; 01140 int r; 01141 01142 // check the nodal voltage changes against the allowed absolute 01143 // and relative tolerance values 01144 for (r = 0; r < N; r++) 01145 { 01146 v_abs = abs (x->get (r) - xprev->get (r)); 01147 v_rel = abs (x->get (r)); 01148 if (v_abs >= vntol + reltol * v_rel) return 0; 01149 if (!convHelper) 01150 { 01151 i_abs = abs (z->get (r) - zprev->get (r)); 01152 i_rel = abs (z->get (r)); 01153 if (i_abs >= abstol + reltol * i_rel) return 0; 01154 } 01155 } 01156 01157 for (r = 0; r < M; r++) 01158 { 01159 i_abs = abs (x->get (r + N) - xprev->get (r + N)); 01160 i_rel = abs (x->get (r + N)); 01161 if (i_abs >= abstol + reltol * i_rel) return 0; 01162 if (!convHelper) 01163 { 01164 v_abs = abs (z->get (r + N) - zprev->get (r + N)); 01165 v_rel = abs (z->get (r + N)); 01166 if (v_abs >= vntol + reltol * v_rel) return 0; 01167 } 01168 } 01169 return 1; 01170 } 01171 01172 /* The function saves the solution and right hand vector of the previous 01173 iteration. */ 01174 template <class nr_type_t> 01175 void nasolver<nr_type_t>::savePreviousIteration (void) 01176 { 01177 if (xprev != NULL) 01178 *xprev = *x; 01179 else 01180 xprev = new tvector<nr_type_t> (*x); 01181 if (zprev != NULL) 01182 *zprev = *z; 01183 else 01184 zprev = new tvector<nr_type_t> (*z); 01185 } 01186 01187 /* The function restores the solution and right hand vector of the 01188 previous (successful) iteration. */ 01189 template <class nr_type_t> 01190 void nasolver<nr_type_t>::restorePreviousIteration (void) 01191 { 01192 if (xprev != NULL) *x = *xprev; 01193 if (zprev != NULL) *z = *zprev; 01194 } 01195 01196 /* The function restarts the NR iteration for each non-linear 01197 circuit. */ 01198 template <class nr_type_t> 01199 void nasolver<nr_type_t>::restartNR (void) 01200 { 01201 circuit * root = subnet->getRoot (); 01202 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) 01203 { 01204 if (c->isNonLinear ()) c->restartDC (); 01205 } 01206 } 01207 01208 /* This function goes through solution (the x vector) and saves the 01209 node voltages of the last iteration into each non-linear 01210 circuit. */ 01211 template <class nr_type_t> 01212 void nasolver<nr_type_t>::saveNodeVoltages (void) 01213 { 01214 int N = countNodes (); 01215 struct nodelist_t * n; 01216 // save all nodes except reference node 01217 for (int r = 0; r < N; r++) 01218 { 01219 n = nlist->getNode (r); 01220 /* for (int i = 0; i < n->size(); i++)*/ 01221 for(auto ¤tn: *n) 01222 { 01223 currentn->getCircuit()->setV (currentn->getPort (), x->get (r)); 01224 } 01225 } 01226 // save reference node 01227 n = nlist->getNode (-1); 01228 for(auto ¤tn: *n) 01229 currentn->getCircuit()->setV (currentn->getPort (), 0.0); 01230 } 01231 01232 /* This function goes through solution (the x vector) and saves the 01233 branch currents through the voltage sources of the last iteration 01234 into each circuit. */ 01235 template <class nr_type_t> 01236 void nasolver<nr_type_t>::saveBranchCurrents (void) 01237 { 01238 int N = countNodes (); 01239 int M = countVoltageSources (); 01240 circuit * vs; 01241 // save all branch currents of voltage sources 01242 for (int r = 0; r < M; r++) 01243 { 01244 vs = findVoltageSource (r); 01245 vs->setJ (r, x->get (r + N)); 01246 } 01247 } 01248 01249 // The function saves the solution vector into each circuit. 01250 template <class nr_type_t> 01251 void nasolver<nr_type_t>::saveSolution (void) 01252 { 01253 saveNodeVoltages (); 01254 saveBranchCurrents (); 01255 } 01256 01257 // This function stores the solution (node voltages and branch currents). 01258 template <class nr_type_t> 01259 void nasolver<nr_type_t>::storeSolution (void) 01260 { 01261 // cleanup solution previously 01262 solution.clear (); 01263 int r; 01264 int N = countNodes (); 01265 int M = countVoltageSources (); 01266 // store all nodes except reference node 01267 for (r = 0; r < N; r++) 01268 { 01269 struct nodelist_t * n = nlist->getNode (r); 01270 nr_type_t gr = x->get (r); 01271 qucs::naentry<nr_type_t> entry(gr, 0); 01272 solution.insert({{n->name, entry }}); 01273 } 01274 // store all branch currents of voltage sources 01275 for (r = 0; r < M; r++) 01276 { 01277 circuit * vs = findVoltageSource (r); 01278 int vn = r - vs->getVoltageSource () + 1; 01279 nr_type_t xg = x->get (r + N); 01280 qucs::naentry<nr_type_t> entry(xg, vn); 01281 solution.insert({{vs->getName (), entry}}); 01282 } 01283 } 01284 01285 // This function recalls the solution (node voltages and branch currents). 01286 template <class nr_type_t> 01287 void nasolver<nr_type_t>::recallSolution (void) 01288 { 01289 int r; 01290 int N = countNodes (); 01291 int M = countVoltageSources (); 01292 // store all nodes except reference node 01293 for (r = 0; r < N; r++) 01294 { 01295 struct nodelist_t * n = nlist->getNode (r); 01296 auto na = solution.find(n->name); 01297 if (na != solution.end()) 01298 if ((*na).second.current == 0) 01299 x->set (r, (*na).second.value); 01300 } 01301 // store all branch currents of voltage sources 01302 for (r = 0; r < M; r++) 01303 { 01304 circuit * vs = findVoltageSource (r); 01305 int vn = r - vs->getVoltageSource () + 1; 01306 auto na = solution.find(vs->getName ()); 01307 if (na != solution.end()) 01308 if ((*na).second.current == vn) 01309 x->set (r + N, (*na).second.value); 01310 } 01311 } 01312 01313 /* This function saves the results of a single solve() functionality 01314 into the output dataset. */ 01315 template <class nr_type_t> 01316 void nasolver<nr_type_t>::saveResults (const std::string &volts, const std::string &s, 01317 int saveOPs, qucs::vector * f) 01318 { 01319 int N = countNodes (); 01320 int M = countVoltageSources (); 01321 01322 // add node voltage variables 01323 if (!volts.empty()) 01324 { 01325 for (int r = 0; r < N; r++) 01326 { 01327 std::string n = createV (r, volts, saveOPs); 01328 if(!n.empty()) 01329 { 01330 saveVariable (n, x->get (r), f); 01331 } 01332 } 01333 } 01334 01335 // add branch current variables 01336 if (!amps.empty()) 01337 { 01338 for (int r = 0; r < M; r++) 01339 { 01340 std::string n = createI (r, amps, saveOPs); 01341 if (!n.empty()) 01342 { 01343 saveVariable (n, x->get (r + N), f); 01344 } 01345 } 01346 } 01347 01348 // add voltage probe data 01349 if (!volts.empty()) 01350 { 01351 circuit * root = subnet->getRoot (); 01352 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) 01353 { 01354 if (!c->isProbe ()) continue; 01355 if (!c->getSubcircuit().empty() && !(saveOPs & SAVE_ALL)) continue; 01356 if (volts != "vn") 01357 c->saveOperatingPoints (); 01358 std::string n = createOP (c->getName (), volts); 01359 saveVariable (n, nr_complex_t (c->getOperatingPoint ("Vr"), 01360 c->getOperatingPoint ("Vi")), f); 01361 } 01362 } 01363 01364 // save operating points of non-linear circuits if requested 01365 if (saveOPs & SAVE_OPS) 01366 { 01367 circuit * root = subnet->getRoot (); 01368 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) 01369 { 01370 if (!c->isNonLinear ()) continue; 01371 if (!c->getSubcircuit ().empty() && !(saveOPs & SAVE_ALL)) continue; 01372 c->calcOperatingPoints (); 01373 for (auto ops: c->getOperatingPoints ()) 01374 { 01375 operatingpoint &p = ops.second; 01376 std::string n = createOP (c->getName (), p.getName ()); 01377 saveVariable (n, p.getValue (), f); 01378 } 01379 } 01380 } 01381 } 01382 01383 /* Create an appropriate variable name for operating points. The 01384 caller is responsible to free() the returned string. */ 01385 template <class nr_type_t> 01386 std::string nasolver<nr_type_t>::createOP (const std::string &c, const std::string &n) 01387 { 01388 return c+"."+n; 01389 } 01390 01391 /* Creates an appropriate variable name for voltages. The caller is 01392 responsible to free() the returned string. */ 01393 template <class nr_type_t> 01394 std::string nasolver<nr_type_t>::createV (int n, const std::string &volts, int saveOPs) 01395 { 01396 if (nlist->isInternal (n)) 01397 return std::string(); 01398 std::string node = nlist->get (n); 01399 if(node.find('.')!=std::string::npos && !(saveOPs & SAVE_ALL)) 01400 return std::string(); 01401 std::string ret = node+"."+volts; 01402 return ret; 01403 } 01404 01405 /* Create an appropriate variable name for currents. The caller is 01406 responsible to free() the returned string. */ 01407 template <class nr_type_t> 01408 std::string nasolver<nr_type_t>::createI (int n, const std::string &s, int saveOPs) 01409 { 01410 circuit * vs = findVoltageSource (n); 01411 01412 // don't output internal (helper) voltage sources 01413 if (vs->isInternalVoltageSource ()) 01414 return std::string(); 01415 01416 /* save only current through real voltage sources and explicit 01417 current probes */ 01418 if (!vs->isVSource () && !(saveOPs & SAVE_OPS)) 01419 return std::string(); 01420 01421 // don't output subcircuit components if not requested 01422 if (!vs->getSubcircuit ().empty() && !(saveOPs & SAVE_ALL)) 01423 return std::string(); 01424 01425 // create appropriate current name for single/multiple voltage sources 01426 std::string name = vs->getName (); 01427 if (vs->getVoltageSources () > 1) 01428 return name+"."+amps+std::to_string(n - vs->getVoltageSource () + 1); 01429 else 01430 return name+"."+amps; 01431 } 01432 01433 01434 /* Alternaive to countNodes () */ 01435 template <class nr_type_t> 01436 int nasolver<nr_type_t>::getN() 01437 { 01438 return countNodes (); 01439 } 01440 01441 /* Alternative to countVoltageSources () */ 01442 template <class nr_type_t> 01443 int nasolver<nr_type_t>::getM() 01444 { 01445 return countVoltageSources (); 01446 } 01447 01448 } // namespace qucs