Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <vector>
- #include <string>
- #include <cstdio>
- #include <cmath>
- // set to 1 to make LU factorization show pivots
- #define VERBOSE_LU 0
- // gMin for diodes etc..
- static const double gMin = 1e-12;
- // voltage tolerance
- static const double vTolerance = 5e-5;
- // thermal voltage for diodes/transistors
- static const double vThermal = 0.026;
- static const unsigned maxIter = 200;
- //
- // General overview
- // ----------------
- //
- // Circuits are built from nodes and Components, where nodes are
- // simply positive integers (with 0 designating ground).
- //
- // Every Component has one or more pins connecting to the circuit
- // nodes as well as zero or more internal nets.
- //
- // While we map pins directly to nets here, the separation would
- // be useful if the solver implemented stuff like net-reordering.
- //
- // MNACell represents a single entry in the solution matrix,
- // where we store constants and time-step dependent constants
- // separately, plus collect pointers to dynamic variables.
- //
- // We track enough information here that we only need to stamp once.
- //
- struct MNACell
- {
- double g; // simple values (eg. resistor conductance)
- double gtimed; // time-scaled values (eg. capacitor conductance)
- // pointers to dynamic variables, added in once per solve
- std::vector<double*> gdyn;
- double lu, prelu; // lu-solver values and matrix pre-LU cache
- std::string txt; // text version of MNA for pretty-printing
- void clear()
- {
- g = 0;
- gtimed = 0;
- txt = "";
- }
- void initLU(double stepScale)
- {
- prelu = g + gtimed * stepScale;
- }
- // restore matrix state and update dynamic values
- void updatePre()
- {
- lu = prelu;
- for(int i = 0; i < gdyn.size(); ++i)
- {
- lu += 0[gdyn[i]];
- }
- }
- };
- // this is for keeping track of node information
- // for the purposes of more intelligent plotting
- struct MNANodeInfo
- {
- enum Type
- {
- tVoltage,
- tCurrent,
- tCount
- };
- Type type; // one auto-range per unit-type
- double scale; // scale factor (eg. charge to voltage)
- std::string name; // node name for display
- };
- // Stores A and b for A*x - b = 0, where x is the solution.
- //
- // A is stored as a vector of rows, for easy in-place pivots
- //
- struct MNASystem
- {
- typedef std::vector<MNACell> MNAVector;
- typedef std::vector<MNAVector> MNAMatrix;
- // node names - for output
- std::vector<MNANodeInfo> nodes;
- MNAMatrix A;
- MNAVector b;
- double time;
- void setSize(int n)
- {
- A.resize(n);
- b.resize(n);
- nodes.resize(n);
- for(unsigned i = 0; i < n; ++i)
- {
- b[i].clear();
- A[i].resize(n);
- char buf[16];
- sprintf(buf, "v%d", i);
- nodes[i].name = buf;
- nodes[i].type = MNANodeInfo::tVoltage;
- nodes[i].scale = 1;
- for(unsigned j = 0; j < n; ++j)
- {
- A[i][j].clear();
- }
- }
- time = 0;
- }
- void stampTimed(double g, int r, int c, const std::string & txt)
- {
- A[r][c].gtimed += g;
- A[r][c].txt += txt;
- }
- void stampStatic(double g, int r, int c, const std::string & txt)
- {
- A[r][c].g += g;
- A[r][c].txt += txt;
- }
- };
- struct IComponent
- {
- virtual ~IComponent() {}
- // return the number of pins for this component
- virtual int pinCount() = 0;
- // return a pointer to array of pin locations
- // NOTE: these will eventually be GUI locations to be unified
- virtual int * getPinLocs() = 0;
- // setup pins and calculate the size of the full netlist
- // the Component<> will handle this automatically
- //
- // - netSize is the current size of the netlist
- // - pins is an array of circuits nodes
- //
- virtual void setupNets(int & netSize, int & states, int * pins) = 0;
- // stamp constants into the matrix
- virtual void stamp(MNASystem & m) = 0;
- // this is for allocating state variables
- virtual void setupStates(int & states) {}
- // update state variables, only tagged nodes
- // this is intended for fixed-time compatible
- // testing to make sure we can code-gen stuff
- virtual void update(MNASystem & m) {}
- // return true if we're done - will keep iterating
- // until all the components are happy
- virtual bool newton(MNASystem & m) { return true; }
- // time-step change, for caps to fix their state-variables
- virtual void scaleTime(double told_per_new) {}
- };
- template <int nPins = 0, int nInternalNets = 0>
- struct Component : IComponent
- {
- static const int nNets = nPins + nInternalNets;
- int pinLoc[nPins];
- int nets[nNets];
- int pinCount() { return nPins; }
- int * getPinLocs() { return pinLoc; }
- void setupNets(int & netSize, int & states, int * pins)
- {
- for(int i = 0; i < nPins; ++i)
- {
- nets[i] = pins[i];
- }
- for(int i = 0; i < nInternalNets; ++i)
- {
- nets[nPins + i] = netSize++;
- }
- setupStates(states);
- }
- };
- static const int unitValueOffset = 4;
- static const int unitValueMax = 8;
- static const char * unitValueSuffixes[unitValueMax+1] = {
- "p", "n", "u", "m", "", "k", "M", "G"
- };
- static void formatUnitValue(char * buf, double v, const char * unit)
- {
- int suff = unitValueOffset + int(log(v) / log(10.)) / 3;
- if(v < 1) suff -= 1;
- if(suff < 0) suff = 0;
- if(suff > unitValueMax) suff = unitValueMax;
- double vr = v / pow(10., 3*double(suff - unitValueOffset));
- sprintf(buf, "%.0f%s%s", vr, unitValueSuffixes[suff], unit);
- }
- struct Resistor : Component<2>
- {
- double r;
- Resistor(double r, int l0, int l1) : r(r)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- }
- void stamp(MNASystem & m)
- {
- char txt[16];
- txt[0] = 'R';
- formatUnitValue(txt+1, r, "");
- double g = 1. / r;
- m.stampStatic(+g, nets[0], nets[0], std::string("+") + txt);
- m.stampStatic(-g, nets[0], nets[1], std::string("-") + txt);
- m.stampStatic(-g, nets[1], nets[0], std::string("-") + txt);
- m.stampStatic(+g, nets[1], nets[1], std::string("+") + txt);
- }
- };
- struct Capacitor : Component<2, 1>
- {
- double c;
- double stateVar;
- double voltage;
- Capacitor(double c, int l0, int l1) : c(c)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- stateVar = 0;
- voltage = 0;
- }
- void stamp(MNASystem & m)
- {
- char buf[16];
- formatUnitValue(buf, c, "F");
- // we can use a trick here, to get the capacitor to
- // work on it's own line with direct trapezoidal:
- //
- // | -g*t +g*t +t | v+
- // | +g*t -g*t -t | v-
- // | +2*g -2*g -1 | state
- //
- // the logic with this is that for constant timestep:
- //
- // i1 = g*v1 - s0 , s0 = g*v0 + i0
- // s1 = 2*g*v1 - s0 <-> s0 = 2*g*v1 - s1
- //
- // then if we substitute back:
- // i1 = g*v1 - (2*g*v1 - s1)
- // = s1 - g*v1
- //
- // this way we just need to copy the new state to the
- // next timestep and there's no actual integration needed
- //
- // the "half time-step" error here means that our state
- // is 2*c*v - i/t but we fix this for display in update
- // and correct the current-part on time-step changes
- // trapezoidal needs another factor of two for the g
- // since c*(v1 - v0) = (i1 + i0)/(2*t), where t = 1/T
- double g = 2*c;
- m.stampTimed(+1, nets[0], nets[2], "+t");
- m.stampTimed(-1, nets[1], nets[2], "-t");
- m.stampTimed(-g, nets[0], nets[0], std::string("-t*") + buf);
- m.stampTimed(+g, nets[0], nets[1], std::string("+t*") + buf);
- m.stampTimed(+g, nets[1], nets[0], std::string("+t*") + buf);
- m.stampTimed(-g, nets[1], nets[1], std::string("-t*") + buf);
- m.stampStatic(+2*g, nets[2], nets[0], std::string("+2*") + buf);
- m.stampStatic(-2*g, nets[2], nets[1], std::string("-2*") + buf);
- m.stampStatic(-1, nets[2], nets[2], "-1");
- // see the comment about v:C[%d] below
- sprintf(buf, "q:C:%d,%d", pinLoc[0], pinLoc[1]);
- m.b[nets[2]].gdyn.push_back(&stateVar);
- m.b[nets[2]].txt = buf;
- // this isn't quite right as state stores 2*c*v - i/t
- // however, we'll fix this in updateFull() for display
- sprintf(buf, "v:C:%d,%d", pinLoc[0], pinLoc[1]);
- m.nodes[nets[2]].name = buf;
- m.nodes[nets[2]].scale = 1 / c;
- }
- void update(MNASystem & m)
- {
- stateVar = m.b[nets[2]].lu;
- // solve legit voltage from the pins
- voltage = m.b[nets[0]].lu - m.b[nets[1]].lu;
- // then we can store this for display here
- // since this value won't be used at this point
- m.b[nets[2]].lu = c*voltage;
- }
- void scaleTime(double told_per_new)
- {
- // the state is 2*c*voltage - i/t0
- // so we subtract out the voltage, scale current
- // and then add the voltage back to get new state
- //
- // note that this also works if the old rate is infinite
- // (ie. t0=0) when going from DC analysis to transient
- //
- double qq = 2*c*voltage;
- stateVar = qq + (stateVar - qq)*told_per_new;
- }
- };
- struct Voltage : Component<2, 1>
- {
- double v;
- Voltage(double v, int l0, int l1) : v(v)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- }
- void stamp(MNASystem & m)
- {
- m.stampStatic(-1, nets[0], nets[2], "-1");
- m.stampStatic(+1, nets[1], nets[2], "+1");
- m.stampStatic(+1, nets[2], nets[0], "+1");
- m.stampStatic(-1, nets[2], nets[1], "-1");
- char buf[16];
- sprintf(buf, "%+.2gV", v);
- m.b[nets[2]].g = v;
- m.b[nets[2]].txt = buf;
- sprintf(buf, "i:V(%+.2g):%d,%d", v, pinLoc[0], pinLoc[1]);
- m.nodes[nets[2]].name = buf;
- m.nodes[nets[2]].type = MNANodeInfo::tCurrent;
- }
- };
- // probe a differential voltage
- // also forces this voltage to actually get solved :)
- struct Probe : Component<2, 1>
- {
- Probe(int l0, int l1)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- }
- void stamp(MNASystem & m)
- {
- // vp + vn - vd = 0
- m.stampStatic(+1, nets[2], nets[0], "+1");
- m.stampStatic(-1, nets[2], nets[1], "-1");
- m.stampStatic(-1, nets[2], nets[2], "-1");
- m.nodes[nets[2]].name = "v:probe";
- }
- void update(MNASystem & m)
- {
- // we could do output here :)
- }
- };
- // function voltage generator
- struct Function : Component<2,1>
- {
- typedef double (*FuncPtr)(double t);
- FuncPtr fn;
- double v;
- Function(FuncPtr fn, int l0, int l1) : fn(fn)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- v = fn(0);
- }
- void stamp(MNASystem & m)
- {
- // this is identical to voltage source
- // except voltage is dynanic
- m.stampStatic(-1, nets[0], nets[2], "-1");
- m.stampStatic(+1, nets[1], nets[2], "+1");
- m.stampStatic(+1, nets[2], nets[0], "+1");
- m.stampStatic(-1, nets[2], nets[1], "-1");
- char buf[16];
- m.b[nets[2]].gdyn.push_back(&v);
- sprintf(buf, "Vfn:%d,%d", pinLoc[0], pinLoc[1]);
- m.b[nets[2]].txt = buf;
- sprintf(buf, "i:Vfn:%d,%d", pinLoc[0], pinLoc[1]);
- m.nodes[nets[2]].name = buf;
- m.nodes[nets[2]].type = MNANodeInfo::tCurrent;
- }
- void update(MNASystem & m)
- {
- v = fn(m.time);
- }
- };
- // POD-struct for PN-junction data, for diodes and BJTs
- //
- struct JunctionPN
- {
- // variables
- double geq, ieq, veq;
- // parameters
- double is, nvt, rnvt, vcrit;
- };
- void initJunctionPN(JunctionPN & pn, double is, double n)
- {
- pn.is = is;
- pn.nvt = n * vThermal;
- pn.rnvt = 1 / pn.nvt;
- pn.vcrit = pn.nvt * log(pn.nvt / (pn.is * sqrt(2.)));
- }
- // linearize junction at the specified voltage
- //
- // ideally we could handle series resistance here as well
- // to avoid putting it on a separate node, but not sure how
- // to make that work as it looks like we'd need Lambert-W then
- void linearizeJunctionPN(JunctionPN & pn, double v)
- {
- double e = pn.is * exp(v * pn.rnvt);
- double i = e - pn.is + gMin * v;
- double g = e * pn.rnvt + gMin;
- pn.geq = g;
- pn.ieq = v*g - i;
- pn.veq = v;
- }
- // returns true if junction is good enough
- bool newtonJunctionPN(JunctionPN & pn, double v)
- {
- double dv = v - pn.veq;
- if(fabs(dv) < vTolerance) return true;
- // check critical voltage and adjust voltage if over
- if(v > pn.vcrit)
- {
- // this formula comes from Qucs documentation
- v = pn.veq + pn.nvt*log((std::max)(pn.is, 1+dv*pn.rnvt));
- }
- linearizeJunctionPN(pn, v);
- return false;
- }
- struct Diode : Component<2, 2>
- {
- JunctionPN pn;
- // should make these parameters
- double rs;
- // l0 -->|-- l1 -- parameters default to approx 1N4148
- Diode(int l0, int l1,
- double rs = 10., double is = 35e-12, double n = 1.24)
- : rs(rs)
- {
- pinLoc[0] = l0;
- pinLoc[1] = l1;
- initJunctionPN(pn, is, n);
- // FIXME: move init to some restart routine?
- // initial condition v = 0
- linearizeJunctionPN(pn, 0);
- }
- bool newton(MNASystem & m)
- {
- return newtonJunctionPN(pn, m.b[nets[2]].lu);
- }
- void stamp(MNASystem & m)
- {
- // Diode could be built with 3 extra nodes:
- //
- // | . . . . +1 | V+
- // | . . . . -1 | V-
- // | . . grs -grs -1 | v:D
- // | . . -grs grs+geq . | v:pn = ieq
- // | -1 +1 +1 . . | i:pn
- //
- // Here grs is the 1/rs series conductance.
- //
- // This gives us the junction voltage (v:pn) and
- // current (i:pn) and the composite voltage (v:D).
- //
- // The i:pn row is an ideal transformer connecting
- // the floating diode to the ground-referenced v:D
- // where we connect the series resistance to v:pn
- // that solves the diode equation with Newton.
- //
- // We can then add the 3rd row to the bottom 2 with
- // multipliers 1 and -rs = -1/grs and drop it:
- //
- // | . . . +1 | V+
- // | . . . -1 | V-
- // | . . geq -1 | v:pn = ieq
- // | -1 +1 +1 rs | i:pn
- //
- // Note that only the v:pn row here is non-linear.
- //
- // We could even do away without the separate row for
- // the current, which would lead to the following:
- //
- // | +grs -grs -grs |
- // | -grs +grs +grs |
- // | -grs +grs +grs+geq | = ieq
- //
- // In practice we keep the current row since it's
- // nice to have it as an output anyway.
- //
- m.stampStatic(-1, nets[3], nets[0], "-1");
- m.stampStatic(+1, nets[3], nets[1], "+1");
- m.stampStatic(+1, nets[3], nets[2], "+1");
- m.stampStatic(+1, nets[0], nets[3], "+1");
- m.stampStatic(-1, nets[1], nets[3], "-1");
- m.stampStatic(-1, nets[2], nets[3], "-1");
- m.stampStatic(rs, nets[3], nets[3], "rs:pn");
- m.A[nets[2]][nets[2]].gdyn.push_back(&pn.geq);
- m.A[nets[2]][nets[2]].txt = "gm:D";
- m.b[nets[2]].gdyn.push_back(&pn.ieq);
- char buf[16];
- sprintf(buf, "i0:D:%d,%d", pinLoc[0], pinLoc[1]);
- m.b[nets[2]].txt = buf;
- sprintf(buf, "v:D:%d,%d", pinLoc[0], pinLoc[1]);
- m.nodes[nets[2]].name = buf;
- sprintf(buf, "i:D:%d,%d", pinLoc[0], pinLoc[1]);
- m.nodes[nets[3]].name = buf;
- m.nodes[nets[3]].type = MNANodeInfo::tCurrent;
- }
- };
- struct BJT : Component<3, 4>
- {
- // emitter and collector junctions
- JunctionPN pnC, pnE;
- // forward and reverse alpha
- double af, ar, rsbc, rsbe;
- bool pnp;
- BJT(int b, int c, int e, bool pnp = false) : pnp(pnp)
- {
- pinLoc[0] = b;
- pinLoc[1] = c;
- pinLoc[2] = e;
- // this attempts a 2n3904-style small-signal
- // transistor, although the values are a bit
- // arbitrarily set to "something reasonable"
- // forward and reverse beta
- double bf = 200;
- double br = 20;
- // forward and reverse alpha
- af = bf / (1 + bf);
- ar = br / (1 + br);
- // these are just rb+re and rb+rc
- // this is not necessarily the best way to
- // do anything, but having junction series
- // resistances helps handle degenerate cases
- rsbc = 5.8376+0.0001;
- rsbe = 5.8376+2.65711;
- //
- // the basic rule is that:
- // af * ise = ar * isc = is
- //
- // FIXME: with non-equal ideality factors
- // we can get non-sensical results, why?
- //
- double is = 6.734e-15;
- double n = 1.24;
- initJunctionPN(pnE, is / af, n);
- initJunctionPN(pnC, is / ar, n);
- linearizeJunctionPN(pnE, 0);
- linearizeJunctionPN(pnC, 0);
- }
- bool newton(MNASystem & m)
- {
- return newtonJunctionPN(pnC, m.b[nets[3]].lu)
- & newtonJunctionPN(pnE, m.b[nets[4]].lu);
- }
- void stamp(MNASystem & m)
- {
- // The basic idea here is the same as with diodes
- // except we do it once for each junction.
- //
- // With the transfer currents sourced from the
- // diode currents, NPN then looks like this:
- //
- // 0 | . . . . . 1-ar 1-af | vB
- // 1 | . . . . . -1 +af | vC
- // 2 | . . . . . +ar -1 | vE
- // 3 | . . . gc . -1 . | v:Qbc = ic
- // 4 | . . . . ge . -1 | v:Qbe = ie
- // 5 | -1 +1 . +1 . rsbc . | i:Qbc
- // 6 | -1 . +1 . +1 . rsbe | i:Qbe
- // ------------------------
- // 0 1 2 3 4 5 6
- //
- // For PNP version, we simply flip the junctions
- // by changing signs of (3,5),(5,3) and (4,6),(6,4).
- //
- // Also just like diodes, we have junction series
- // resistances, rather than terminal resistances.
- //
- // This works just as well, but should be kept
- // in mind when fitting particular transistors.
- //
- // diode currents to external base
- m.stampStatic(1-ar, nets[0], nets[5], "1-ar");
- m.stampStatic(1-af, nets[0], nets[6], "1-af");
- // diode currents to external collector and emitter
- m.stampStatic(-1, nets[1], nets[5], "-1");
- m.stampStatic(-1, nets[2], nets[6], "-1");
- // series resistances
- m.stampStatic(rsbc, nets[5], nets[5], "rsbc");
- m.stampStatic(rsbe, nets[6], nets[6], "rsbe");
- // current - junction connections
- // for the PNP case we flip the signs of these
- // to flip the diode junctions wrt. the above
- if(pnp)
- {
- m.stampStatic(-1, nets[5], nets[3], "-1");
- m.stampStatic(+1, nets[3], nets[5], "+1");
- m.stampStatic(-1, nets[6], nets[4], "-1");
- m.stampStatic(+1, nets[4], nets[6], "+1");
- }
- else
- {
- m.stampStatic(+1, nets[5], nets[3], "+1");
- m.stampStatic(-1, nets[3], nets[5], "-1");
- m.stampStatic(+1, nets[6], nets[4], "+1");
- m.stampStatic(-1, nets[4], nets[6], "-1");
- }
- // external voltages to collector current
- m.stampStatic(-1, nets[5], nets[0], "-1");
- m.stampStatic(+1, nets[5], nets[1], "+1");
- // external voltages to emitter current
- m.stampStatic(-1, nets[6], nets[0], "-1");
- m.stampStatic(+1, nets[6], nets[2], "+1");
- // source transfer currents to external pins
- m.stampStatic(+ar, nets[2], nets[5], "+ar");
- m.stampStatic(+af, nets[1], nets[6], "+af");
- char buf[16];
- // dynamic variables
- m.A[nets[3]][nets[3]].gdyn.push_back(&pnC.geq);
- m.A[nets[3]][nets[3]].txt = "gm:Qbc";
- m.b[nets[3]].gdyn.push_back(&pnC.ieq);
- sprintf(buf, "i0:Q:%d,%d,%d:cb", pinLoc[0], pinLoc[1], pinLoc[2]);
- m.b[nets[3]].txt = buf;
- m.A[nets[4]][nets[4]].gdyn.push_back(&pnE.geq);
- m.A[nets[4]][nets[4]].txt = "gm:Qbe";
- m.b[nets[4]].gdyn.push_back(&pnE.ieq);
- sprintf(buf, "i0:Q:%d,%d,%d:eb", pinLoc[0], pinLoc[1], pinLoc[2]);
- m.b[nets[4]].txt = buf;
- sprintf(buf, "v:Q:%d,%d,%d:%s",
- pinLoc[0], pinLoc[1], pinLoc[2], pnp ? "cb" : "bc");
- m.nodes[nets[3]].name = buf;
- sprintf(buf, "v:Q:%d,%d,%d:%s",
- pinLoc[0], pinLoc[1], pinLoc[2], pnp ? "eb" : "be");
- m.nodes[nets[4]].name = buf;
- sprintf(buf, "i:Q:%d,%d,%d:bc", pinLoc[0], pinLoc[1], pinLoc[2]);
- m.nodes[nets[5]].name = buf;
- m.nodes[nets[5]].type = MNANodeInfo::tCurrent;
- m.nodes[nets[5]].scale = 1 - ar;
- sprintf(buf, "i:Q:%d,%d,%d:be", pinLoc[0], pinLoc[1], pinLoc[2]);
- m.nodes[nets[6]].name = buf;
- m.nodes[nets[6]].type = MNANodeInfo::tCurrent;
- m.nodes[nets[6]].scale = 1 - af;
- }
- };
- struct NetList
- {
- typedef std::vector<IComponent*> ComponentList;
- NetList(int nodes) : nets(nodes), states(0)
- {
- }
- void addComponent(IComponent * c)
- {
- // this is a bit "temporary" for now
- c->setupNets(nets, states, c->getPinLocs());
- components.push_back(c);
- }
- void buildSystem()
- {
- system.setSize(nets);
- for(int i = 0; i < components.size(); ++i)
- {
- components[i]->stamp(system);
- }
- printf("Prepare for DC analysis..\n");
- setStepScale(0);
- tStep = 0;
- }
- void dumpMatrix()
- {
- std::vector<int> maxWidth(nets);
- for(int i = 0; i < nets; ++i) maxWidth[i] = 1;
- int nnMax = 1;
- for(int i = 0; i < nets; ++i)
- {
- nnMax = std::max(nnMax, (int)system.nodes[i].name.size());
- for(int j = 0; j < nets; ++j)
- {
- maxWidth[j] = std::max(maxWidth[j],
- (int)system.A[i][j].txt.size());
- }
- }
- char buf[1024];
- for(unsigned i = 0; i < nets; ++i)
- {
- int off = sprintf(buf, "%2d: | ", i);
- for(int j = 0; j < nets; ++j)
- {
- off += sprintf(buf+off,
- " %*s ", maxWidth[j],
- system.A[i][j].txt.size()
- ? system.A[i][j].txt.c_str()
- : ((system.A[i][j].lu==0) ? "." : "#"));
- }
- sprintf(buf+off, " | %-*s = %s",
- nnMax, system.nodes[i].name.c_str(),
- system.b[i].txt.size()
- ? system.b[i].txt.c_str() : (!i ? "ground" : "0"));
- puts(buf);
- }
- }
- void setTimeStep(double tStepSize)
- {
- for(int i = 0; i < components.size(); ++i)
- {
- components[i]->scaleTime(tStep / tStepSize);
- }
- tStep = tStepSize;
- double stepScale = 1. / tStep;
- printf("timeStep changed to %.2g (%.2g Hz)\n", tStep, stepScale);
- setStepScale(stepScale);
- }
- void simulateTick()
- {
- int iter;
- for(iter = 0; iter < maxIter; ++iter)
- {
- // restore matrix state and add dynamic values
- updatePre();
- luFactor();
- luForward();
- luSolve();
- if(newton()) break;
- }
- system.time += tStep;
- update();
- printf(" %02.4f |", system.time);
- int fillPost = 0;
- for(int i = 1; i < nets; ++i)
- {
- printf("\t%+.4e", system.b[i].lu * system.nodes[i].scale);
- for(int j = 1; j < nets; ++j)
- {
- if(system.A[i][j].lu != 0) ++fillPost;
- }
- }
- printf("\t %d iters, LU density: %.1f%%\n",
- iter, 100 * fillPost / ((nets-1.f)*(nets-1.f)));
- }
- void printHeaders()
- {
- printf("\n time: | ");
- for(int i = 1; i < nets; ++i)
- {
- printf("%16s", system.nodes[i].name.c_str());
- }
- printf("\n\n");
- }
- // plotting and such would want to use this
- const MNASystem & getMNA() { return system; }
- protected:
- double tStep;
- int nets, states;
- ComponentList components;
- MNASystem system;
- void update()
- {
- for(int i = 0; i < components.size(); ++i)
- {
- components[i]->update(system);
- }
- }
- // return true if we're done
- bool newton()
- {
- bool done = 1;
- for(int i = 0; i < components.size(); ++i)
- {
- done &= components[i]->newton(system);
- }
- return done;
- }
- void initLU(double stepScale)
- {
- for(int i = 0; i < nets; ++i)
- {
- system.b[i].initLU(stepScale);
- for(int j = 0; j < nets; ++j)
- {
- system.A[i][j].initLU(stepScale);
- }
- }
- }
- void setStepScale(double stepScale)
- {
- // initialize matrix for LU and save it to cache
- initLU(stepScale);
- int fill = 0;
- for(int i = 1; i < nets; ++i)
- {
- for(int j = 1; j < nets; ++j)
- {
- if(system.A[i][j].prelu != 0
- || system.A[i][j].gdyn.size()) ++fill;
- }
- }
- printf("MNA density %.1f%%\n", 100 * fill / ((nets-1.)*(nets-1.)));
- }
- void updatePre()
- {
- for(int i = 0; i < nets; ++i)
- {
- system.b[i].updatePre();
- for(int j = 0; j < nets; ++j)
- {
- system.A[i][j].updatePre();
- }
- }
- }
- void luFactor()
- {
- int p;
- for(p = 1; p < nets; ++p)
- {
- // FIND PIVOT
- {
- int pr = p;
- for(int r = p; r < nets; ++r)
- {
- if(fabs(system.A[r][p].lu)
- > fabs(system.A[pr][p].lu))
- {
- pr = r;
- }
- }
- // swap if necessary
- if(pr != p)
- {
- std::swap(system.A[p], system.A[pr]);
- std::swap(system.b[p], system.b[pr]);
- }
- #if VERBOSE_LU
- printf("pivot %d (from %d): %+.2e\n",
- p, pr, system.A[p][p].lu);
- #endif
- }
- if(0 == system.A[p][p].lu)
- {
- printf("Failed to find a pivot!!");
- return;
- }
- // take reciprocal for D entry
- system.A[p][p].lu = 1 / system.A[p][p].lu;
- // perform reduction on rows below
- for(int r = p+1; r < nets; ++r)
- {
- if(system.A[r][p].lu == 0) continue;
- system.A[r][p].lu *= system.A[p][p].lu;
- for(int c = p+1; c < nets; ++c)
- {
- if(system.A[p][c].lu == 0) continue;
- system.A[r][c].lu -=
- system.A[p][c].lu * system.A[r][p].lu;
- }
- }
- }
- }
- // this does forward substitution for the solution vector
- int luForward()
- {
- int p;
- for(p = 1; p < nets; ++p)
- {
- // perform reduction on rows below
- for(int r = p+1; r < nets; ++r)
- {
- if(system.A[r][p].lu == 0) continue;
- if(system.b[p].lu == 0) continue;
- system.b[r].lu -= system.b[p].lu * system.A[r][p].lu;
- }
- }
- return p;
- }
- // solves nodes backwards from limit-1
- // if solveAll is true, solves all the nodes
- // otherwise if solveNoIter is true, solves until !wantUpdate
- // if both flags are false, solves until !wantIter
- int luSolve()
- {
- for(int r = nets; --r;)
- {
- //printf("solve node %d\n", r);
- for(int s = r+1; s < nets; ++s)
- {
- system.b[r].lu -= system.b[s].lu * system.A[r][s].lu;
- }
- system.b[r].lu *= system.A[r][r].lu;
- }
- return 1;
- }
- };
- double fnGen(double t)
- {
- if(t < 0.0001) return 0;
- return (fmod(2000*t, 1) > (.5+.4*sin(2*acos(-1)*100*t))) ? .25 : -.25;
- }
- int main()
- {
- /*
- BJT test circuit
- this is kinda good at catching problems :)
- */
- NetList * net = new NetList(8);
- // node 1 is +12V
- net->addComponent(new Voltage(+12, 1, 0));
- // bias for base
- net->addComponent(new Resistor(100e3, 2, 1));
- net->addComponent(new Resistor(11e3, 2, 0));
- // input cap and function
- net->addComponent(new Capacitor(.1e-6, 2, 5));
- net->addComponent(new Resistor(1e3, 5, 6));
- net->addComponent(new Function(fnGen, 6, 0));
- // collector resistance
- net->addComponent(new Resistor(10e3, 1, 3));
- // emitter resistance
- net->addComponent(new Resistor(680, 4, 0));
- // bjt on b:2,c:3,e:4
- net->addComponent(new BJT(2,3,4));
- // output emitter follower (direct rail collector?)
- net->addComponent(new BJT(3, 1, 7));
- net->addComponent(new Resistor(100e3, 7, 0));
- net->addComponent(new Probe(7, 0));
- net->buildSystem();
- // get DC solution (optional)
- net->simulateTick();
- net->setTimeStep(1e-4);
- net->printHeaders();
- for(int i = 0; i < 25; ++i)
- {
- net->simulateTick();
- }
- net->printHeaders();
- net->dumpMatrix();
- };
Add Comment
Please, Sign In to add comment