Qucs-core
0.0.19
|
00001 /* 00002 * hbsolver.cpp - harmonic balance solver class implementation 00003 * 00004 * Copyright (C) 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 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include<algorithm> 00030 00031 #include <stdio.h> 00032 00033 #include "object.h" 00034 #include "logging.h" 00035 #include "complex.h" 00036 #include "vector.h" 00037 #include "node.h" 00038 #include "circuit.h" 00039 #include "component_id.h" 00040 #include "net.h" 00041 #include "netdefs.h" 00042 #include "strlist.h" 00043 #include "ptrlist.h" 00044 #include "tvector.h" 00045 #include "tmatrix.h" 00046 #include "eqnsys.h" 00047 #include "analysis.h" 00048 #include "dataset.h" 00049 #include "fourier.h" 00050 #include "hbsolver.h" 00051 00052 #define HB_DEBUG 0 00053 00054 namespace qucs { 00055 00056 using namespace fourier; 00057 00058 // Constructor creates an unnamed instance of the hbsolver class. 00059 hbsolver::hbsolver () : analysis () { 00060 type = ANALYSIS_HBALANCE; 00061 frequency = 0; 00062 nlnodes = lnnodes = banodes = nanodes = NULL; 00063 Y = Z = A = NULL; 00064 NA = YV = JQ = JG = JF = NULL; 00065 OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL; 00066 vs = x = NULL; 00067 runs = 0; 00068 ndfreqs = NULL; 00069 } 00070 00071 // Constructor creates a named instance of the hbsolver class. 00072 hbsolver::hbsolver (char * n) : analysis (n) { 00073 type = ANALYSIS_HBALANCE; 00074 frequency = 0; 00075 nlnodes = lnnodes = banodes = nanodes = NULL; 00076 Y = Z = A = NULL; 00077 NA = YV = JQ = JG = JF = NULL; 00078 OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL; 00079 vs = x = NULL; 00080 runs = 0; 00081 ndfreqs = NULL; 00082 } 00083 00084 // Destructor deletes the hbsolver class object. 00085 hbsolver::~hbsolver () { 00086 // delete nodelists 00087 delete nlnodes; 00088 delete lnnodes; 00089 delete banodes; 00090 delete nanodes; 00091 00092 // delete temporary matrices 00093 delete A; 00094 delete Z; 00095 delete Y; 00096 00097 // delete matrices 00098 delete NA; 00099 delete YV; 00100 delete JQ; 00101 delete JG; 00102 delete JF; 00103 00104 // delete vectors 00105 delete IC; 00106 delete IS; 00107 delete FV; 00108 delete IL; 00109 delete IN; 00110 delete IG; 00111 delete FQ; 00112 delete VS; 00113 delete VP; 00114 delete vs; 00115 delete OM; 00116 delete IR; 00117 delete QR; 00118 delete RH; 00119 00120 delete x; 00121 delete[] ndfreqs; 00122 } 00123 00124 /* The copy constructor creates a new instance of the hbsolver class 00125 based on the given hbsolver object. */ 00126 hbsolver::hbsolver (hbsolver & o) : analysis (o) { 00127 frequency = o.frequency; 00128 negfreqs = o.negfreqs; 00129 posfreqs = o.posfreqs; 00130 nlnodes = o.nlnodes; 00131 lnnodes = o.lnnodes; 00132 banodes = o.banodes; 00133 nanodes = o.nanodes; 00134 Y = Z = A = NULL; 00135 NA = YV = JQ = JG = JF = NULL; 00136 OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL; 00137 vs = x = NULL; 00138 runs = o.runs; 00139 ndfreqs = NULL; 00140 } 00141 00142 #define VS_(r) (*VS) (r) 00143 #define OM_(r) (*OM) (r) 00144 00145 /* This is the HB netlist solver. It prepares the circuit list and 00146 solves it then. */ 00147 int hbsolver::solve (void) { 00148 00149 int iterations = 0, done = 0; 00150 int MaxIterations = getPropertyInteger ("MaxIter"); 00151 00152 // collect different parts of the circuit 00153 splitCircuits (); 00154 00155 // create frequency array 00156 collectFrequencies (); 00157 00158 // find interconnects between the linear and non-linear subcircuit 00159 getNodeLists (); 00160 00161 // prepares the linear part --> 0 = IC + [YV] * VS 00162 prepareLinear (); 00163 00164 runs++; 00165 logprint (LOG_STATUS, "NOTIFY: %s: solving for %d frequencies\n", 00166 getName (), lnfreqs); 00167 00168 if (nbanodes > 0) { 00169 00170 // start balancing 00171 logprint (LOG_STATUS, "NOTIFY: %s: balancing at %d nodes\n", getName (), 00172 nbanodes); 00173 00174 // prepares the non-linear part 00175 prepareNonLinear (); 00176 00177 #if HB_DEBUG 00178 fprintf (stderr, "YV -- transY in f:\n"); YV->print (); 00179 fprintf (stderr, "IC -- constant current in f:\n"); IC->print (); 00180 #endif 00181 00182 // start iteration 00183 do { 00184 iterations++; 00185 00186 #if HB_DEBUG 00187 fprintf (stderr, "\n -- iteration step: %d\n", iterations); 00188 fprintf (stderr, "vs -- voltage in t:\n"); vs->print (); 00189 #endif 00190 00191 // evaluate component functionality and fill matrices and vectors 00192 loadMatrices (); 00193 00194 #if HB_DEBUG 00195 fprintf (stderr, "FQ -- charge in t:\n"); FQ->print (); 00196 fprintf (stderr, "IG -- current in t:\n"); IG->print (); 00197 #endif 00198 00199 // currents into frequency domain 00200 VectorFFT (IG); 00201 00202 // charges into frequency domain 00203 VectorFFT (FQ); 00204 00205 // right hand side currents and charges into the frequency domain 00206 VectorFFT (IR); 00207 VectorFFT (QR); 00208 00209 #if HB_DEBUG 00210 fprintf (stderr, "VS -- voltage in f:\n"); VS->print (); 00211 fprintf (stderr, "FQ -- charge in f:\n"); FQ->print (); 00212 fprintf (stderr, "IG -- current in f:\n"); IG->print (); 00213 fprintf (stderr, "IR -- corrected Jacobi current in f:\n"); IR->print (); 00214 #endif 00215 00216 // solve HB equation --> FV = IC + [YV] * VS + j[O] * FQ + IG 00217 solveHB (); 00218 00219 #if HB_DEBUG 00220 fprintf (stderr, "FV -- error vector F(V) in f:\n"); FV->print (); 00221 fprintf (stderr, "IL -- linear currents in f:\n"); IL->print (); 00222 fprintf (stderr, "IN -- non-linear currents in f:\n"); IN->print (); 00223 fprintf (stderr, "RH -- right-hand side currents in f:\n"); RH->print (); 00224 #endif 00225 00226 // termination criteria met 00227 if (iterations > 1 && checkBalance ()) { 00228 done = 1; 00229 break; 00230 } 00231 00232 #if HB_DEBUG 00233 fprintf (stderr, "JG -- G-Jacobian in t:\n"); JG->print (); 00234 fprintf (stderr, "JQ -- C-Jacobian in t:\n"); JQ->print (); 00235 #endif 00236 00237 // G-Jacobian into frequency domain 00238 MatrixFFT (JG); 00239 00240 // C-Jacobian into frequency domain 00241 MatrixFFT (JQ); 00242 00243 #if HB_DEBUG 00244 fprintf (stderr, "JQ -- dQ/dV C-Jacobian in f:\n"); JQ->print (); 00245 fprintf (stderr, "JG -- dI/dV G-Jacobian in f:\n"); JG->print (); 00246 #endif 00247 00248 // calculate Jacobian --> JF = [YV] + j[O] * JQ + JG 00249 calcJacobian (); 00250 00251 #if HB_DEBUG 00252 fprintf (stderr, "JF -- full Jacobian in f:\n"); JF->print (); 00253 #endif 00254 00255 // solve equation system --> JF * VS(n+1) = JF * VS(n) - FV 00256 solveVoltages (); 00257 00258 #if HB_DEBUG 00259 fprintf (stderr, "VS -- next voltage in f:\n"); VS->print (); 00260 #endif 00261 00262 // inverse FFT of frequency domain voltage vector VS(n+1) 00263 VectorIFFT (vs); 00264 } 00265 // check termination criteria (balanced frequency domain currents) 00266 while (!done && iterations < MaxIterations); 00267 00268 if (iterations >= MaxIterations) { 00269 qucs::exception * e = new qucs::exception (EXCEPTION_NO_CONVERGENCE); 00270 e->setText ("no convergence in %s analysis after %d iterations", 00271 getName (), iterations); 00272 throw_exception (e); 00273 logprint (LOG_ERROR, "%s: no convergence after %d iterations\n", 00274 getName (), iterations); 00275 } 00276 else { 00277 logprint (LOG_STATUS, "%s: convergence reached after %d iterations\n", 00278 getName (), iterations); 00279 } 00280 } 00281 else { 00282 // no balancing necessary 00283 logprint (LOG_STATUS, "NOTIFY: %s: no balancing necessary\n", getName ()); 00284 } 00285 00286 // print exception stack 00287 estack.print (); 00288 00289 // apply AC analysis to the complete network in order to obtain the 00290 // final results 00291 finalSolution (); 00292 00293 // save results into dataset 00294 saveResults (); 00295 return 0; 00296 } 00297 00298 /* Goes through the list of circuit objects and runs its calcHB() 00299 function. */ 00300 void hbsolver::calc (hbsolver * self) { 00301 circuit * root = self->getNet()->getRoot (); 00302 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) { 00303 c->calcHB (self->frequency); 00304 } 00305 } 00306 00307 /* Goes through the list of circuit objects and runs its initHB() 00308 function. */ 00309 void hbsolver::initHB (void) { 00310 circuit * root = subnet->getRoot (); 00311 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) { 00312 c->initHB (); 00313 } 00314 } 00315 00316 /* Goes through the list of circuit objects and runs its initDC() 00317 function. */ 00318 void hbsolver::initDC (void) { 00319 circuit * root = subnet->getRoot (); 00320 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) { 00321 c->initDC (); 00322 } 00323 } 00324 00325 // Returns true if circuit is a HB source. 00326 bool hbsolver::isExcitation (circuit * c) { 00327 int type = c->getType (); 00328 if (type == CIR_PAC || type == CIR_VAC || type == CIR_VDC) 00329 return true; 00330 return false; 00331 } 00332 00333 // Expands the frequency array using the given frequency and the order. 00334 void hbsolver::expandFrequencies (nr_double_t f, int n) { 00335 auto nfreqs = negfreqs; 00336 auto pfreqs = posfreqs; 00337 int i, k, len = nfreqs.size (); 00338 negfreqs.clear (); 00339 posfreqs.clear (); 00340 if (len > 0) { 00341 // frequency expansion for full frequency sets 00342 for (i = 0; i <= n + 1; i++) { 00343 for (k = 0; k < len; k++) { 00344 negfreqs.push_back (i * f + nfreqs[k]); 00345 } 00346 } 00347 for (i = -n; i < 0; i++) { 00348 for (k = 0; k < len; k++) { 00349 negfreqs.push_back (i * f + nfreqs[k]); 00350 } 00351 } 00352 for (i = 0; i <= 2 * n + 1; i++) { 00353 for (k = 0; k < len; k++) { 00354 posfreqs.push_back (i * f + pfreqs[k]); 00355 } 00356 } 00357 } 00358 else { 00359 // first frequency 00360 for (i = 0; i <= n + 1; i++) negfreqs.push_back (i * f); 00361 for (i = -n; i < 0; i++) negfreqs.push_back (i * f); 00362 for (i = 0; i <= 2 * n + 1; i++) posfreqs.push_back (i * f); 00363 } 00364 } 00365 00366 // Calculates an order fulfilling the "power of two" requirement. 00367 int hbsolver::calcOrder (int n) { 00368 int o, order = n * 2; // current order + DC + negative freqencies 00369 for (o = 1; o < order; o <<= 1) ; // a power of 2 00370 return o / 2 - 1; 00371 } 00372 00373 /* The function computes the harmonic frequencies excited in the 00374 circuit list depending on the maximum number of harmonics per 00375 exitation and saves its results into the 'negfreqs' vector. */ 00376 void hbsolver::collectFrequencies (void) { 00377 00378 // initialization 00379 negfreqs.clear (); 00380 posfreqs.clear (); 00381 rfreqs.clear (); 00382 dfreqs.clear (); 00383 delete[] ndfreqs; 00384 00385 // obtain order 00386 int i, n = calcOrder (getPropertyInteger ("n")); 00387 00388 // expand frequencies for each exitation 00389 nr_double_t f; 00390 for (auto * c : excitations) { 00391 if (c->getType () != CIR_VDC) { // no extra DC sources 00392 if ((f = c->getPropertyDouble ("f")) != 0.0) { 00393 const auto epsilon = std::numeric_limits<nr_double_t>::epsilon(); 00394 auto found = 00395 std::find_if(dfreqs.cbegin(),dfreqs.cend(), 00396 [f,epsilon](decltype(*dfreqs.cbegin()) x) { 00397 return (std::abs(x-f) < epsilon); 00398 }) 00399 ; 00400 if (found == dfreqs.cend()) { // no double frequencies 00401 dfreqs.push_back (f); 00402 expandFrequencies (f, n); 00403 } 00404 } 00405 } 00406 } 00407 00408 // no excitations 00409 if (negfreqs.size () == 0) { 00410 // use specified frequency 00411 f = getPropertyDouble ("f"); 00412 dfreqs.push_back (f); 00413 expandFrequencies (f, n); 00414 } 00415 00416 // build frequency dimension lengths 00417 ndfreqs = new int[dfreqs.size ()]; 00418 for (i = 0; i < dfreqs.size (); i++) { 00419 ndfreqs[i] = (n + 1) * 2; 00420 } 00421 00422 #if HB_DEBUG 00423 fprintf (stderr, "%d frequencies: [ ", negfreqs.getSize ()); 00424 for (i = 0; i < negfreqs.getSize (); i++) { 00425 fprintf (stderr, "%g ", (double) real (negfreqs.get (i))); 00426 } 00427 fprintf (stderr, "]\n"); 00428 #endif /* HB_DEBUG */ 00429 00430 // build list of positive frequencies including DC 00431 for (n = 0; n < negfreqs.size (); n++) { 00432 if ((f = negfreqs[n]) < 0.0) continue; 00433 rfreqs.push_back (f); 00434 } 00435 lnfreqs = rfreqs.size (); 00436 nlfreqs = negfreqs.size (); 00437 00438 // pre-calculate the j[O] vector 00439 OM = new tvector<nr_complex_t> (nlfreqs); 00440 for (n = i = 0; n < nlfreqs; n++, i++) 00441 OM_(n) = nr_complex_t (0, 2 * pi * negfreqs[i]); 00442 } 00443 00444 // Split netlist into excitation, linear and non-linear part. 00445 void hbsolver::splitCircuits (void) { 00446 circuit * root = subnet->getRoot (); 00447 for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) { 00448 if (c->isNonLinear ()) { 00449 // non-linear part 00450 nolcircuits.push_front(c); 00451 } 00452 else if (isExcitation (c)) { 00453 // get sinusoidal sources 00454 excitations.push_front(c); 00455 } 00456 else if (c->getType () != CIR_GROUND) { 00457 // linear part 00458 lincircuits.push_front (c); 00459 } 00460 } 00461 } 00462 00463 // Creates a nodelist for the given circuit list. 00464 strlist * hbsolver::circuitNodes (ptrlist<circuit> circuits) { 00465 strlist * nodes = new strlist (); 00466 for (auto * c : circuits) { 00467 for (int i = 0; i < c->getSize (); i++) { 00468 const char * n = c->getNode(i)->getName (); 00469 if (strcmp (n, "gnd")) { 00470 if (!nodes->contains (n)) nodes->add (n); 00471 } 00472 } 00473 } 00474 return nodes; 00475 } 00476 00477 // Obtain node lists for linear and non-linear part. 00478 void hbsolver::getNodeLists (void) { 00479 // non-linear nodes 00480 nlnodes = circuitNodes (nolcircuits); 00481 // linear nodes 00482 lnnodes = circuitNodes (lincircuits); 00483 // excitation nodes 00484 exnodes = circuitNodes (excitations); 00485 00486 #if DEBUG && 0 00487 fprintf (stderr, "nonlinear nodes: [ %s ]\n", nlnodes->toString ()); 00488 fprintf (stderr, " linear nodes: [ %s ]\n", lnnodes->toString ()); 00489 #endif 00490 00491 // organization of the nodes for the MNA: 00492 // -------------------------------------- 00493 // 1.) balanced nodes: all connected to at least one non-linear device 00494 // 2.a) the excitation nodes 00495 // 2.b) all linear nodes not contained in 1. and 2.a 00496 // 2.c) additional gyrators of linear nodes (built-in voltage sources) 00497 // please note: excitation nodes also in 2.b; 1. and 2.a are 'ports' 00498 00499 nanodes = new strlist (*nlnodes); // list 1. 00500 strlistiterator it; 00501 // add excitation nodes; list 2.a 00502 for (it = strlistiterator (exnodes); *it; ++it) 00503 nanodes->append (*it); 00504 // add linear nodes; list 2.b 00505 for (it = strlistiterator (lnnodes); *it; ++it) { 00506 if (!nanodes->contains (*it)) 00507 nanodes->append (*it); 00508 } 00509 00510 banodes = new strlist (*nlnodes); 00511 #if DEBUG && 0 00512 fprintf (stderr, " balanced nodes: [ %s ]\n", banodes->toString ()); 00513 fprintf (stderr, " exnodes nodes: [ %s ]\n", exnodes->toString ()); 00514 fprintf (stderr, " nanodes nodes: [ %s ]\n", nanodes->toString ()); 00515 #endif 00516 } 00517 00518 /* The function applies unique voltage source identifiers to each 00519 built in internal voltage source in the list of circuits. */ 00520 int hbsolver::assignVoltageSources (ptrlist<circuit> circuits) { 00521 int sources = 0; 00522 for (auto *c: circuits) { 00523 if (c->getVoltageSources () > 0) { 00524 c->setVoltageSource (sources); 00525 sources += c->getVoltageSources (); 00526 } 00527 } 00528 return sources; 00529 } 00530 00531 /* The function assigns unique node identifiers to each circuit node. */ 00532 int hbsolver::assignNodes (ptrlist<circuit> circuits, strlist * nodes, 00533 int offset) { 00534 // through all collected nodes 00535 for (int nr = 0; nr < nodes->length (); nr++) { 00536 char * nn = nodes->get (nr); 00537 // through all circuits 00538 for (auto *c : circuits) { 00539 // assign current identifier if any of the circuit nodes equals 00540 // the current one 00541 for (int i = 0; i < c->getSize (); i++) { 00542 node * n = c->getNode (i); 00543 if (!strcmp (n->getName (), nn)) n->setNode (offset + nr + 1); 00544 } 00545 } 00546 } 00547 return nodes->length (); 00548 } 00549 00550 // Prepares the linear operations. 00551 void hbsolver::prepareLinear (void) { 00552 for (auto *lc : lincircuits) 00553 lc->initHB (); 00554 nlnvsrcs = assignVoltageSources (lincircuits); 00555 nnlvsrcs = excitations.size (); 00556 nnanodes = nanodes->length (); 00557 nexnodes = exnodes->length (); 00558 nbanodes = banodes->length (); 00559 assignNodes (lincircuits, nanodes); 00560 assignNodes (excitations, nanodes); 00561 createMatrixLinearA (); 00562 createMatrixLinearY (); 00563 calcConstantCurrent (); 00564 } 00565 00566 /* The function creates the complex linear network MNA matrix. It 00567 contains the MNA entries for all linear components for each 00568 requested frequency. */ 00569 void hbsolver::createMatrixLinearA (void) { 00570 int M = nlnvsrcs; 00571 int N = nnanodes; 00572 int f = 0; 00573 nr_double_t freq; 00574 00575 // create new MNA matrix 00576 A = new tmatrix<nr_complex_t> ((N + M) * lnfreqs); 00577 00578 // through each frequency 00579 for (int i = 0; i < rfreqs.size (); i++) { 00580 freq = rfreqs[i]; 00581 // calculate components' MNA matrix for the given frequency 00582 for (auto *lc : lincircuits) 00583 lc->calcHB (freq); 00584 // fill in all matrix entries for the given frequency 00585 fillMatrixLinearA (A, f++); 00586 } 00587 00588 // save a copy of the original MNA matrix 00589 NA = new tmatrix<nr_complex_t> (*A); 00590 } 00591 00592 // some definitions for the linear matrix filler 00593 #undef A_ 00594 #undef B_ 00595 #define A_(r,c) (*A) ((r)*lnfreqs+f,(c)*lnfreqs+f) 00596 #define G_(r,c) A_(r,c) 00597 #define B_(r,c) A_(r,c+N) 00598 #define C_(r,c) A_(r+N,c) 00599 #define D_(r,c) A_(r+N,c+N) 00600 00601 /* This function fills in the MNA matrix entries into the A matrix for 00602 a given frequency index. */ 00603 void hbsolver::fillMatrixLinearA (tmatrix<nr_complex_t> * A, int f) { 00604 int N = nnanodes; 00605 00606 // through each linear circuit 00607 for (auto *cir : lincircuits) { 00608 int s = cir->getSize (); 00609 int nr, nc, r, c, v; 00610 00611 // apply G-matrix entries 00612 for (r = 0; r < s; r++) { 00613 if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue; 00614 for (c = 0; c < s; c++) { 00615 if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue; 00616 G_(nr, nc) += cir->getY (r, c); 00617 } 00618 } 00619 00620 // augmented part -- built in voltage sources 00621 if ((v = cir->getVoltageSources ()) > 0) { 00622 00623 // apply B-matrix entries 00624 for (r = 0; r < s; r++) { 00625 if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue; 00626 for (c = 0; c < v; c++) { 00627 nc = cir->getVoltageSource () + c; 00628 B_(nr, nc) += cir->getB (r, nc); 00629 } 00630 } 00631 00632 // apply C-matrix entries 00633 for (r = 0; r < v; r++) { 00634 nr = cir->getVoltageSource () + r; 00635 for (c = 0; c < s; c++) { 00636 if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue; 00637 C_(nr, nc) += cir->getC (nr, c); 00638 } 00639 } 00640 00641 // apply D-matrix entries 00642 for (r = 0; r < v; r++) { 00643 nr = cir->getVoltageSource () + r; 00644 for (c = 0; c < v; c++) { 00645 nc = cir->getVoltageSource () + c; 00646 D_(nr, nc) += cir->getD (nr, nc); 00647 } 00648 } 00649 } 00650 } 00651 } 00652 00653 // The function inverts the given matrix A into the matrix H. 00654 void hbsolver::invertMatrix (tmatrix<nr_complex_t> * A, 00655 tmatrix<nr_complex_t> * H) { 00656 eqnsys<nr_complex_t> eqns; 00657 int N = A->getCols (); 00658 tvector<nr_complex_t> * x = new tvector<nr_complex_t> (N); 00659 tvector<nr_complex_t> * z = new tvector<nr_complex_t> (N); 00660 00661 try_running () { 00662 // create LU decomposition of the A matrix 00663 eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT); 00664 eqns.passEquationSys (A, x, z); 00665 eqns.solve (); 00666 } 00667 // appropriate exception handling 00668 catch_exception () { 00669 case EXCEPTION_PIVOT: 00670 default: 00671 logprint (LOG_ERROR, "WARNING: %s: during TI inversion\n", getName ()); 00672 estack.print (); 00673 } 00674 00675 // use the LU decomposition to obtain the inverse H 00676 eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT); 00677 for (int c = 0; c < N; c++) { 00678 z->set (0.0); 00679 z->set (c, 1.0); 00680 eqns.passEquationSys (A, x, z); 00681 eqns.solve (); 00682 for (int r = 0; r < N; r++) H->set (r, c, x->get (r)); 00683 } 00684 delete x; 00685 delete z; 00686 } 00687 00688 // Some defines for matrix element access. 00689 #define V_(r) (*V) (r) 00690 #define I_(r) (*I) (r) 00691 #undef A_ 00692 #define A_(r,c) (*A) (r,c) 00693 00694 #define Z_(r,c) (*Z) (r,c) 00695 #define Y_(r,c) (*Y) (r,c) 00696 00697 #define ZVU_(r,c) Z_(r,c) 00698 #define ZVL_(r,c) Z_((r)*lnfreqs+f+sn,c) 00699 #define ZCU_(r,c) Z_(r,(c)*lnfreqs+f+sn) 00700 #define ZCL_(r,c) Z_((r)*lnfreqs+f+sn,(c)*lnfreqs+f+sn) 00701 00702 #define YV_(r,c) (*YV) (r,c) 00703 #define NA_(r,c) (*NA) (r,c) 00704 #define JF_(r,c) (*JF) (r,c) 00705 00706 /* The following function performs the following steps: 00707 1. form the MNA matrix A including all nodes (linear, non-linear and 00708 excitations) 00709 2. compute the variable transimpedance matrix entries for the nodes 00710 to be balanced 00711 3. compute the constant transimpedance matrix entries for the constant 00712 current vector caused by the excitations 00713 4. invert this overall transimpedance matrix 00714 5. extract the variable transadmittance matrix entries 00715 */ 00716 void hbsolver::createMatrixLinearY (void) { 00717 int M = nlnvsrcs; 00718 int N = nnanodes; 00719 int c, r, f; 00720 00721 // size of overall MNA matrix 00722 int sa = (N + M) * lnfreqs; 00723 int sv = nbanodes; 00724 int se = nnlvsrcs; 00725 int sy = sv + se; 00726 00727 // allocate new transimpedance matrix 00728 Z = new tmatrix<nr_complex_t> (sy * lnfreqs); 00729 00730 // prepare equation system 00731 eqnsys<nr_complex_t> eqns; 00732 tvector<nr_complex_t> * V; 00733 tvector<nr_complex_t> * I; 00734 00735 // 1. create variable transimpedance matrix entries relating 00736 // voltages at the balanced nodes to the currents through these 00737 // nodes into the non-linear part 00738 int sn = sv * lnfreqs; 00739 V = new tvector<nr_complex_t> (sa); 00740 I = new tvector<nr_complex_t> (sa); 00741 00742 // connect a 100 Ohm resistor (to ground) to balanced node in the MNA matrix 00743 for (c = 0; c < sv * lnfreqs; c++) A_(c, c) += 0.01; 00744 00745 // connect a 100 Ohm resistor (in parallel) to each excitation 00746 for (auto *vs : excitations) { 00747 // get positive and negative node 00748 int pnode = vs->getNode(NODE_1)->getNode (); 00749 int nnode = vs->getNode(NODE_2)->getNode (); 00750 for (f = 0; f < lnfreqs; f++) { // for each frequency 00751 int pn = (pnode - 1) * lnfreqs + f; 00752 int nn = (nnode - 1) * lnfreqs + f; 00753 if (pnode) A_(pn, pn) += 0.01; 00754 if (nnode) A_(nn, nn) += 0.01; 00755 if (pnode && nnode) { 00756 A_(pn, nn) -= 0.01; 00757 A_(nn, pn) -= 0.01; 00758 } 00759 } 00760 } 00761 00762 // LU decompose the MNA matrix 00763 try_running () { 00764 eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT); 00765 eqns.passEquationSys (A, V, I); 00766 eqns.solve (); 00767 } 00768 // appropriate exception handling 00769 catch_exception () { 00770 case EXCEPTION_PIVOT: 00771 default: 00772 logprint (LOG_ERROR, "WARNING: %s: during A factorization\n", getName ()); 00773 estack.print (); 00774 } 00775 00776 // aquire variable transimpedance matrix entries 00777 eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT); 00778 for (c = 0; c < sn; c++) { 00779 I->set (0.0); 00780 I_(c) = 1.0; 00781 eqns.passEquationSys (A, V, I); 00782 eqns.solve (); 00783 // ZV | .. 00784 // ---+--- 00785 // .. | .. 00786 for (r = 0; r < sn; r++) ZVU_(r, c) = V_(r); 00787 // .. | .. 00788 // ---+--- 00789 // ZV | .. 00790 r = 0; 00791 for (auto ite = excitations.begin(); ite != excitations.end(); ++ite, r++) { 00792 // lower part entries 00793 for (f = 0; f < lnfreqs; f++) { 00794 ZVL_(r, c) = excitationZ (V, *ite, f); 00795 } 00796 } 00797 } 00798 00799 // create constant transimpedance matrix entries relating the 00800 // source voltages to the interconnection currents 00801 int vsrc = 0; 00802 // aquire constant transadmittance matrix entries 00803 for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) { 00804 circuit * vs = *it; 00805 // get positive and negative node 00806 int pnode = vs->getNode(NODE_1)->getNode (); 00807 int nnode = vs->getNode(NODE_2)->getNode (); 00808 for (f = 0; f < lnfreqs; f++) { // for each frequency 00809 int pn = (pnode - 1) * lnfreqs + f; 00810 int nn = (nnode - 1) * lnfreqs + f; 00811 I->set (0.0); 00812 if (pnode) I_(pn) = +1.0; 00813 if (nnode) I_(nn) = -1.0; 00814 eqns.passEquationSys (A, V, I); 00815 eqns.solve (); 00816 // .. | ZC 00817 // ---+--- 00818 // .. | .. 00819 for (r = 0; r < sn; r++) { 00820 // upper part of the entries 00821 ZCU_(r, vsrc) = V_(r); 00822 } 00823 // .. | .. 00824 // ---+--- 00825 // .. | ZC 00826 r = 0; 00827 for (auto ite = excitations.begin(); ite != excitations.end(); ++ite, r++) { 00828 // lower part entries 00829 ZCL_(r, vsrc) = excitationZ (V, *ite, f); 00830 } 00831 } 00832 } 00833 delete I; 00834 delete V; 00835 00836 // allocate new transadmittance matrix 00837 Y = new tmatrix<nr_complex_t> (sy * lnfreqs); 00838 00839 // invert the Z matrix to a Y matrix 00840 invertMatrix (Z, Y); 00841 00842 // substract the 100 Ohm resistor 00843 for (c = 0; c < sy * lnfreqs; c++) Y_(c, c) -= 0.01; 00844 00845 // extract the variable transadmittance matrix 00846 YV = new tmatrix<nr_complex_t> (sv * nlfreqs); 00847 00848 // variable transadmittance matrix must be continued conjugately 00849 *YV = expandMatrix (*Y, sv); 00850 00851 // delete overall temporary MNA matrix 00852 delete A; A = NULL; 00853 // delete transimpedance matrix 00854 delete Z; Z = NULL; 00855 } 00856 00857 /* Little helper function obtaining a transimpedance value for the 00858 given voltage source (excitation) and for a given frequency 00859 index. */ 00860 nr_complex_t hbsolver::excitationZ (tvector<nr_complex_t> * V, circuit * vs, 00861 int f) { 00862 // get positive and negative node 00863 int pnode = vs->getNode(NODE_1)->getNode (); 00864 int nnode = vs->getNode(NODE_2)->getNode (); 00865 nr_complex_t z = 0.0; 00866 if (pnode) z += V_((pnode - 1) * lnfreqs + f); 00867 if (nnode) z -= V_((nnode - 1) * lnfreqs + f); 00868 return z; 00869 } 00870 00871 /* This function computes the constant current vectors using the 00872 voltage of the excitations and the transadmittance matrix 00873 entries. */ 00874 void hbsolver::calcConstantCurrent (void) { 00875 int se = nnlvsrcs * lnfreqs; 00876 int sn = nbanodes * lnfreqs; 00877 int r, c, vsrc = 0; 00878 00879 // collect excitation voltages 00880 tvector<nr_complex_t> VC (se); 00881 for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) { 00882 circuit * vs = *it; 00883 vs->initHB (); 00884 vs->setVoltageSource (0); 00885 for (int f = 0; f < rfreqs.size (); f++) { // for each frequency 00886 nr_double_t freq = rfreqs[f]; 00887 vs->calcHB (freq); 00888 VC (vsrc * lnfreqs + f) = vs->getE (VSRC_1); 00889 } 00890 } 00891 00892 // compute constant current vector for balanced nodes 00893 IC = new tvector<nr_complex_t> (sn); 00894 // .. | YC * VC 00895 // ---+--- 00896 // .. | .. 00897 for (r = 0; r < sn; r++) { 00898 nr_complex_t i = 0.0; 00899 for (c = 0; c < se; c++) { 00900 i += Y_(r, c + sn) * VC (c); 00901 } 00902 int f = r % lnfreqs; 00903 if (f != 0 && f != lnfreqs - 1) i /= 2; 00904 IC->set (r, i); 00905 } 00906 // expand the constant current conjugate 00907 *IC = expandVector (*IC, nbanodes); 00908 00909 // compute constant current vector for sources itself 00910 IS = new tvector<nr_complex_t> (se); 00911 // .. | .. 00912 // ---+--- 00913 // .. | YC * VC 00914 for (r = 0; r < se; r++) { 00915 nr_complex_t i = 0.0; 00916 for (c = 0; c < se; c++) { 00917 i += Y_(r + sn, c + sn) * VC (c); 00918 } 00919 IS->set (r, i); 00920 } 00921 00922 // delete overall transadmittance matrix 00923 delete Y; Y = NULL; 00924 } 00925 00926 /* Checks whether currents through the interconnects of the linear and 00927 non-linear subcircuit (in the frequency domain) are equal. */ 00928 int hbsolver::checkBalance (void) { 00929 nr_double_t iabstol = getPropertyDouble ("iabstol"); 00930 nr_double_t vabstol = getPropertyDouble ("vabstol"); 00931 nr_double_t reltol = getPropertyDouble ("reltol"); 00932 int n, len = FV->size (); 00933 for (n = 0; n < len; n++) { 00934 // check iteration voltages 00935 nr_double_t v_abs = abs (VS->get (n) - VP->get (n)); 00936 nr_double_t v_rel = abs (VS->get (n)); 00937 if (v_abs >= vabstol + reltol * v_rel) return 0; 00938 // check balanced currents 00939 nr_complex_t il = IL->get (n); 00940 nr_complex_t in = IN->get (n); 00941 if (il != in) { 00942 nr_double_t i_abs = abs (il + in); 00943 nr_double_t i_rel = abs ((il + in) / (il - in)); 00944 if (i_abs >= iabstol && 2.0 * i_rel >= reltol) return 0; 00945 } 00946 } 00947 return 1; 00948 } 00949 00950 // some definitions for the non-linear matrix filler 00951 #undef G_ 00952 #undef C_ 00953 #define G_(r,c) (*jg) ((r)*nlfreqs+f,(c)*nlfreqs+f) 00954 #define C_(r,c) (*jq) ((r)*nlfreqs+f,(c)*nlfreqs+f) 00955 #undef FI_ 00956 #undef FQ_ 00957 #define FI_(r) (*ig) ((r)*nlfreqs+f) 00958 #define FQ_(r) (*fq) ((r)*nlfreqs+f) 00959 #define IR_(r) (*ir) ((r)*nlfreqs+f) 00960 #define QR_(r) (*qr) ((r)*nlfreqs+f) 00961 00962 /* This function fills in the matrix and vector entries for the 00963 non-linear HB equations for a given frequency index. */ 00964 void hbsolver::fillMatrixNonLinear (tmatrix<nr_complex_t> * jg, 00965 tmatrix<nr_complex_t> * jq, 00966 tvector<nr_complex_t> * ig, 00967 tvector<nr_complex_t> * fq, 00968 tvector<nr_complex_t> * ir, 00969 tvector<nr_complex_t> * qr, 00970 int f) { 00971 // through each linear circuit 00972 for (auto *cir: nolcircuits) { 00973 int s = cir->getSize (); 00974 int nr, nc, r, c; 00975 00976 for (r = 0; r < s; r++) { 00977 if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue; 00978 // apply G- and C-matrix entries 00979 for (c = 0; c < s; c++) { 00980 if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue; 00981 G_(nr, nc) += cir->getY (r, c); 00982 C_(nr, nc) += cir->getQV (r, c); 00983 } 00984 // apply I- and Q-vector entries 00985 FI_(nr) -= cir->getI (r); 00986 FQ_(nr) -= cir->getQ (r); 00987 // ThinkME: positive or negative? 00988 IR_(nr) += cir->getGV (r) + cir->getI (r); 00989 QR_(nr) += cir->getCV (r) + cir->getQ (r); 00990 } 00991 } 00992 } 00993 00994 /* The function initializes the non-linear part of the HB. */ 00995 void hbsolver::prepareNonLinear (void) { 00996 int N = nbanodes; 00997 00998 // allocate matrices and vectors 00999 if (FQ == NULL) { 01000 FQ = new tvector<nr_complex_t> (N * nlfreqs); 01001 } 01002 if (IG == NULL) { 01003 IG = new tvector<nr_complex_t> (N * nlfreqs); 01004 } 01005 if (IR == NULL) { 01006 IR = new tvector<nr_complex_t> (N * nlfreqs); 01007 } 01008 if (QR == NULL) { 01009 QR = new tvector<nr_complex_t> (N * nlfreqs); 01010 } 01011 if (JG == NULL) { 01012 JG = new tmatrix<nr_complex_t> (N * nlfreqs); 01013 } 01014 if (JQ == NULL) { 01015 JQ = new tmatrix<nr_complex_t> (N * nlfreqs); 01016 } 01017 if (JF == NULL) { 01018 JF = new tmatrix<nr_complex_t> (N * nlfreqs); 01019 } 01020 01021 // voltage vector in frequency and time domain 01022 if (VS == NULL) { 01023 VS = new tvector<nr_complex_t> (N * nlfreqs); 01024 } 01025 if (vs == NULL) { 01026 vs = new tvector<nr_complex_t> (N * nlfreqs); 01027 } 01028 if (VP == NULL) { 01029 VP = new tvector<nr_complex_t> (N * nlfreqs); 01030 } 01031 01032 // error vector 01033 if (FV == NULL) { 01034 FV = new tvector<nr_complex_t> (N * nlfreqs); 01035 } 01036 // right hand side vector 01037 if (RH == NULL) { 01038 RH = new tvector<nr_complex_t> (N * nlfreqs); 01039 } 01040 01041 // linear and non-linear current vector 01042 if (IL == NULL) { 01043 IL = new tvector<nr_complex_t> (N * nlfreqs); 01044 } 01045 if (IN == NULL) { 01046 IN = new tvector<nr_complex_t> (N * nlfreqs); 01047 } 01048 01049 // assign nodes 01050 assignNodes (nolcircuits, nanodes); 01051 01052 // initialize circuits 01053 for (auto *cir : nolcircuits) { 01054 cir->initHB (nlfreqs); 01055 } 01056 } 01057 01058 /* Saves the node voltages of the given circuit and for the given 01059 frequency entry into the circuit voltage vector. */ 01060 void hbsolver::saveNodeVoltages (circuit * cir, int f) { 01061 int r, nr, s = cir->getSize (); 01062 for (r = 0; r < s; r++) { 01063 if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue; 01064 // apply V-vector entries 01065 cir->setV (r, real (vs->get (nr * nlfreqs + f))); 01066 } 01067 } 01068 01069 /* The function saves voltages into non-linear circuits, runs each 01070 non-linear components' HB calculator for each frequency and applies 01071 the matrix and vector entries appropriately. */ 01072 void hbsolver::loadMatrices (void) { 01073 // clear matrices and vectors before 01074 IG->set (0.0); 01075 FQ->set (0.0); 01076 IR->set (0.0); 01077 QR->set (0.0); 01078 JG->set (0.0); 01079 JQ->set (0.0); 01080 // through each frequency 01081 for (int f = 0; f < nlfreqs; f++) { 01082 // calculate components' HB matrices and vector for the given frequency 01083 for (auto *cir : nolcircuits) { 01084 saveNodeVoltages (cir, f); // node voltages 01085 cir->calcHB (f); // HB calculator 01086 } 01087 // fill in all matrix entries for the given frequency 01088 fillMatrixNonLinear (JG, JQ, IG, FQ, IR, QR, f); 01089 } 01090 } 01091 01092 /* The following function transforms a vector using a Fast Fourier 01093 Transformation from the time domain to the frequency domain. 01094 \todo rewrite ugly sould die 01095 */ 01096 void hbsolver::VectorFFT (tvector<nr_complex_t> * V, int isign) { 01097 int i, k, r; 01098 int n = nlfreqs; 01099 int nd = dfreqs.size (); 01100 int nodes = V->size () / n; 01101 nr_double_t * d = (double *)V->getData (); 01102 01103 if (nd == 1) { 01104 // for each node a single 1d-FFT 01105 for (k = i = 0; i < nodes; i++, k += 2 * n) { 01106 nr_double_t * dst = &d[k]; 01107 _fft_1d (dst, n, isign); 01108 if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= n; 01109 } 01110 } 01111 else { 01112 // for each node a single nd-FFT 01113 for (k = i = 0; i < nodes; i++, k += 2 * n) { 01114 nr_double_t * dst = &d[k]; 01115 _fft_nd (dst, ndfreqs, nd, isign); 01116 if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= ndfreqs[0]; 01117 } 01118 } 01119 } 01120 01121 /* The following function transforms a vector using an Inverse Fast 01122 Fourier Transformation from the frequency domain to the domain 01123 time. */ 01124 void hbsolver::VectorIFFT (tvector<nr_complex_t> * V, int isign) { 01125 VectorFFT (V, -isign); 01126 } 01127 01128 /* The following function transforms a matrix using a Fast Fourier 01129 Transformation from the time domain to the frequency domain. */ 01130 void hbsolver::MatrixFFT (tmatrix<nr_complex_t> * M) { 01131 01132 #if THE_SLO_ALGO 01133 // each column FFT 01134 for (int c = 0; c < M->getCols (); c++) { 01135 tvector<nr_complex_t> V = M->getCol (c); 01136 VectorFFT (&V); 01137 M->setCol (c, V); 01138 } 01139 01140 // each row IFFT 01141 for (int r = 0; r < M->getRows (); r++) { 01142 tvector<nr_complex_t> V = M->getRow (r); 01143 VectorIFFT (&V); 01144 M->setRow (r, V); 01145 } 01146 #else 01147 int c, r, nc, nr; 01148 // for each non-linear node block 01149 for (nc = c = 0; c < nbanodes; c++, nc += nlfreqs) { 01150 for (nr = r = 0; r < nbanodes; r++, nr += nlfreqs) { 01151 tvector<nr_complex_t> V (nlfreqs); 01152 int fr, fc, fi; 01153 // transform the sub-diagonal only 01154 for (fc = 0; fc < nlfreqs; fc++) V (fc) = M->get (nr + fc, nc + fc); 01155 VectorFFT (&V); 01156 // fill in resulting sub-matrix for the node 01157 for (fc = 0; fc < nlfreqs; fc++) { 01158 for (fi = nlfreqs - 1 - fc, fr = 0; fr < nlfreqs; fr++) { 01159 if (++fi >= nlfreqs) fi = 0; 01160 M->set (nr + fr, nc + fc, V (fi)); 01161 } 01162 } 01163 } 01164 } 01165 #endif 01166 } 01167 01168 /* This function solves the actual HB equation in the frequency domain. 01169 F(V) = IC + [YV] * VS + j[O] * FQ + IG -> 0 01170 Care must be taken with indexing here: In the frequency domain only 01171 real positive frequencies are used and computed, but in the time 01172 domain we have more values because of the periodic continuation in 01173 the frequency domain. 01174 RHS = j[O] * CV + GV - (IC + j[O] * FQ + IG) 01175 Also the right hand side of the equation system for the new voltage 01176 vector is computed here. */ 01177 void hbsolver::solveHB (void) { 01178 // for each non-linear node 01179 for (int r = 0; r < nbanodes * nlfreqs; ) { 01180 // for each frequency 01181 for (int f = 0; f < nlfreqs; f++, r++) { 01182 nr_complex_t il = 0.0, in = 0.0, ir = 0.0; 01183 // constant current vector due to sources 01184 il += IC->get (r); 01185 // part 1 of right hand side vector 01186 ir -= il; 01187 // transadmittance matrix multiplied by voltage vector 01188 for (int c = 0; c < nbanodes * nlfreqs; c++) { 01189 il += YV_(r, c) * VS_(c); 01190 } 01191 // charge vector 01192 in += OM_(f) * FQ->get (r); 01193 // current vector 01194 in += IG->get (r); 01195 // part 2, 3 and 4 of right hand side vector 01196 ir += IR->get (r); 01197 ir += OM_(f) * QR->get (r); 01198 // put values into result vectors 01199 RH->set (r, ir); 01200 FV->set (r, il + in); 01201 IL->set (r, il); 01202 IN->set (r, in); 01203 } 01204 } 01205 } 01206 01207 /* The function calculates the full Jacobian JF = [YV] + j[O] * JQ + JG */ 01208 void hbsolver::calcJacobian (void) { 01209 int c, r, fc, fr, rt, ct; 01210 /* add admittances of capacitance matrix JQ and non-linear 01211 admittances matrix JG into complete Jacobian JF */ 01212 for (c = 0; c < nbanodes; c++) { 01213 ct = c * nlfreqs; 01214 for (fc = 0; fc < nlfreqs; fc++, ct++) { 01215 for (r = 0; r < nbanodes; r++) { 01216 rt = r * nlfreqs; 01217 for (fr = 0; fr < nlfreqs; fr++, rt++) { 01218 JF_(rt, ct) = JG->get (rt, ct) + JQ->get (rt, ct) * OM_(fr); 01219 } 01220 } 01221 } 01222 } 01223 *JF += *YV; // add linear admittance matrix 01224 } 01225 01226 /* The function expands the given vector in the frequency domain to 01227 make it a real valued signal in the time domain. */ 01228 tvector<nr_complex_t> hbsolver::expandVector (tvector<nr_complex_t> V, 01229 int nodes) { 01230 tvector<nr_complex_t> res (nodes * nlfreqs); 01231 int r, ff, rf, rt; 01232 for (r = 0; r < nodes; r++) { 01233 rt = r * nlfreqs; 01234 rf = r * lnfreqs; 01235 // copy first part of vector 01236 for (ff = 0; ff < lnfreqs; ff++, rf++, rt++) { 01237 res (rt) = V (rf); 01238 } 01239 // continue vector conjugated 01240 for (rf -= 2; ff < nlfreqs; ff++, rf--, rt++) { 01241 res (rt) = conj (V (rf)); 01242 } 01243 } 01244 return res; 01245 } 01246 01247 /* The function expands the given matrix in the frequency domain to 01248 make it a real valued signal in the time domain. */ 01249 tmatrix<nr_complex_t> hbsolver::expandMatrix (tmatrix<nr_complex_t> M, 01250 int nodes) { 01251 tmatrix<nr_complex_t> res (nodes * nlfreqs); 01252 int r, c, rf, rt, cf, ct, ff; 01253 for (r = 0; r < nodes; r++) { 01254 for (c = 0; c < nodes; c++) { 01255 rf = r * lnfreqs; 01256 rt = r * nlfreqs; 01257 cf = c * lnfreqs; 01258 ct = c * nlfreqs; 01259 // copy first part of diagonal 01260 for (ff = 0; ff < lnfreqs; ff++, cf++, ct++, rf++, rt++) { 01261 res (rt, ct) = M (rf, cf); 01262 } 01263 // continue diagonal conjugated 01264 for (cf -= 2, rf -= 2; ff < nlfreqs; ff++, cf--, ct++, rf--, rt++) { 01265 res (rt, ct) = conj (M (rf, cf)); 01266 } 01267 } 01268 } 01269 return res; 01270 } 01271 01272 /* This function solves the equation system 01273 JF * VS(n+1) = JF * VS(n) - FV 01274 in order to obtains a new voltage vector in the frequency domain. */ 01275 void hbsolver::solveVoltages (void) { 01276 // save previous iteration voltage 01277 *VP = *VS; 01278 01279 // setup equation system 01280 eqnsys<nr_complex_t> eqns; 01281 try_running () { 01282 // use LU decomposition for solving 01283 eqns.setAlgo (ALGO_LU_DECOMPOSITION); 01284 eqns.passEquationSys (JF, VS, RH); 01285 eqns.solve (); 01286 } 01287 // appropriate exception handling 01288 catch_exception () { 01289 case EXCEPTION_PIVOT: 01290 default: 01291 logprint (LOG_ERROR, "WARNING: %s: during NR iteration\n", getName ()); 01292 estack.print (); 01293 } 01294 01295 // save new voltages in time domain vector 01296 *vs = *VS; 01297 } 01298 01299 /* The following function extends the existing linear MNA matrix to 01300 contain the additional rows and columns for the excitation voltage 01301 sources. */ 01302 tmatrix<nr_complex_t> hbsolver::extendMatrixLinear (tmatrix<nr_complex_t> M, 01303 int nodes) { 01304 int no = M.getCols (); 01305 tmatrix<nr_complex_t> res (no + nodes * lnfreqs); 01306 // copy the existing part 01307 for (int r = 0; r < no; r++) { 01308 for (int c = 0; c < no; c++) { 01309 res (r, c) = M (r, c); 01310 } 01311 } 01312 return res; 01313 } 01314 01315 /* The function fills in the missing MNA entries for the excitation 01316 voltage sources into the extended rows and columns as well as the 01317 actual voltage values into the right hand side vector. */ 01318 void hbsolver::fillMatrixLinearExtended (tmatrix<nr_complex_t> * Y, 01319 tvector<nr_complex_t> * I) { 01320 // through each excitation source 01321 int of = lnfreqs * (nlnvsrcs + nnanodes); 01322 int sc = of; 01323 01324 for (auto *vs : excitations) { 01325 // get positive and negative node 01326 int pnode = vs->getNode(NODE_1)->getNode (); 01327 int nnode = vs->getNode(NODE_2)->getNode (); 01328 for (int f = 0; f < lnfreqs; f++, sc++) { // for each frequency 01329 // fill right hand side vector 01330 nr_double_t freq = rfreqs[f]; 01331 vs->calcHB (freq); 01332 I_(sc) = vs->getE (VSRC_1); 01333 // fill MNA entries 01334 int pn = (pnode - 1) * lnfreqs + f; 01335 int nn = (nnode - 1) * lnfreqs + f; 01336 if (pnode) { 01337 Y_(pn, sc) = +1.0; 01338 Y_(sc, pn) = +1.0; 01339 } 01340 if (nnode) { 01341 Y_(nn, sc) = -1.0; 01342 Y_(sc, nn) = -1.0; 01343 } 01344 } 01345 } 01346 } 01347 01348 /* The function calculates and saves the final solution. */ 01349 void hbsolver::finalSolution (void) { 01350 01351 // extend the linear MNA matrix 01352 *NA = extendMatrixLinear (*NA, nnlvsrcs); 01353 01354 int S = NA->getCols (); 01355 int N = nnanodes * lnfreqs; 01356 01357 // right hand side vector 01358 tvector<nr_complex_t> * I = new tvector<nr_complex_t> (S); 01359 // temporary solution 01360 tvector<nr_complex_t> * V = new tvector<nr_complex_t> (S); 01361 // final solution 01362 x = new tvector<nr_complex_t> (N); 01363 01364 // fill in missing MNA entries 01365 fillMatrixLinearExtended (NA, I); 01366 01367 // put currents through balanced nodes into right hand side 01368 for (int n = 0; n < nbanodes; n++) { 01369 for (int f = 0; f < lnfreqs; f++) { 01370 nr_complex_t i = IL->get (n * nlfreqs + f); 01371 if (f != 0 && f != lnfreqs - 1) i *= 2; 01372 I_(n * lnfreqs + f) = i; 01373 } 01374 } 01375 01376 // use QR decomposition for the final solution 01377 try_running () { 01378 eqnsys<nr_complex_t> eqns; 01379 eqns.setAlgo (ALGO_LU_DECOMPOSITION); 01380 eqns.passEquationSys (NA, V, I); 01381 eqns.solve (); 01382 } 01383 // appropriate exception handling 01384 catch_exception () { 01385 case EXCEPTION_PIVOT: 01386 default: 01387 logprint (LOG_ERROR, "WARNING: %s: during final AC analysis\n", 01388 getName ()); 01389 estack.print (); 01390 } 01391 for (int i = 0; i < N; i++) x->set (i, V_(i)); 01392 } 01393 01394 // Saves simulation results. 01395 void hbsolver::saveResults (void) { 01396 vector * f; 01397 // add current frequency to the dependency of the output dataset 01398 if ((f = data->findDependency ("hbfrequency")) == NULL) { 01399 f = new vector ("hbfrequency"); 01400 data->addDependency (f); 01401 } 01402 // save frequency vector 01403 if (runs == 1) { 01404 for (int i = 0; i < lnfreqs; i++) f->add (rfreqs[i]); 01405 } 01406 // save node voltage vectors 01407 int nanode = 0; 01408 for (strlistiterator it (nanodes); *it; ++it, nanode++) { 01409 int l = strlen (*it); 01410 char * n = (char *) malloc (l + 4); 01411 sprintf (n, "%s.Vb", *it); 01412 for (int i = 0; i < lnfreqs; i++) { 01413 saveVariable (n, x->get (i + nanode * lnfreqs), f); 01414 } 01415 } 01416 } 01417 01418 // properties 01419 PROP_REQ [] = { 01420 { "n", PROP_INT, { 1, PROP_NO_STR }, PROP_MIN_VAL (1) }, 01421 PROP_NO_PROP }; 01422 PROP_OPT [] = { 01423 { "f", PROP_REAL, { 1e9, PROP_NO_STR }, PROP_POS_RANGEX }, 01424 { "iabstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I }, 01425 { "vabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I }, 01426 { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I }, 01427 { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) }, 01428 PROP_NO_PROP }; 01429 struct define_t hbsolver::anadef = 01430 { "HB", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF }; 01431 01432 } // namespace qucs 01433 01434 /* TODO list for HB Solver: 01435 - Take care about nodes with non-linear components only. 01436 - AC Power Sources (extra Z and open loop voltage). 01437 - Current sources. 01438 - Balancing of multi-dimensional non-linear networks. 01439 - Sources directly connected to non-linear components and no other 01440 linear component (insert short). 01441 - Bug: With capacitors at hand there is voltage convergence but no 01442 current balancing. 01443 - Output currents through voltage sources. 01444 */