<> <> <> <> <> <> DIRECTORY BasicTime USING [Now], IO USING [char, int, PutF, PutFR, PutRope, real, rope], Labels USING [Set], Real USING [Fix, FixC, FixI, Float, RoundLI, SmallestNormalizedNumber], RealFns USING [SqRt], RefText USING [Equal], Rope USING [Equal, FromRefText, Length, ROPE], spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, InitPP, initPrint, makeStringNB, maxCountsPerStep, maxLog, modBrPtr, modFuncPtr, next, NextStage, nodePtr, plotFromList, printFromList, PutMsgLine, RefFuncBody, RefModBody, RefC, RefI, RefL, RefR, RefV, ResetCheckPoint, ShowDetails, UpdatePP]; spSolve: CEDAR PROGRAM IMPORTS BasicTime, IO, Labels, Real, RealFns, RefText, spGlobals EXPORTS spGlobals= BEGIN OPEN spGlobals; Retreat: PUBLIC SIGNAL[cause: Rope.ROPE]= CODE; Failure: PUBLIC SIGNAL[errorNum: NAT]= CODE; solveError: SIGNAL[s: Rope.ROPE]= CODE; solveTest: SIGNAL[s: Rope.ROPE]= CODE; updateFunctions: PROC[handle: Handle, f: modFuncPtr, firstTime: BOOL _ FALSE]= { mb: modBrPtr; source: argNodes; args, res: argList; allInactive: BOOL; handle.vars.modelsSkipped_ 0; UNTIL f=NIL DO source_ f.arguments; args_ f.argVector; allInactive_ TRUE; IF source # NIL THEN FOR i: NAT IN [0..source.size) DO args[i]_ source[i].nHist.y; allInactive_ ((NOT firstTime) AND allInactive AND source[i].marked); ENDLOOP; WITH f.body SELECT FROM fcn: RefFuncBody => { fcn.branch.comVal_ fcn.functionProc[handle.vars.t, args, f.parmVector]; }; mod: RefModBody => IF ~allInactive THEN { <> oldArgs: argList_ mod.oldArgVector; FOR i: NAT IN [0..source.size) DO oldArgs[i]_ source[i].nHist.oldy; ENDLOOP; res_ mod.modelResults; mod.modelProc[args, oldArgs, f.parmVector, res]; FOR i: NAT IN [0..source.size) DO -- what the hell is this doing here? source[i].nHist.y_ args[i]; ENDLOOP; mb_ mod.modelBranches; UNTIL mb=NIL DO mb.b.comVal_ res[mb.b.modelIndex]; mb_ mb.nextBranch; ENDLOOP; } ELSE handle.vars.modelsSkipped_ handle.vars.modelsSkipped + 1; ENDCASE; f_ f.nextMFPtr; ENDLOOP; }; -- updateFunctions GSiteration: PROC[handle: Handle, dT: REAL, print, firstTime, mark: BOOL_ FALSE] RETURNS[nIter: INT_ 0]= { E, newE, dvdt, I, newI, cond, csum, otherV, maxI, denom, tempDvDt: REAL; n: nodePtr; b: branchPtr; vCur: branchPtr; vCurBody: RefV; bLink: branchLinkPtr; plus, nodeOK, OK: BOOL; realMaxIter: INT = IF handle.vars.maxIter >= LAST[INT]/2 THEN LAST[INT] ELSE handle.vars.maxIter + handle.vars.maxIter; updateFunctions[handle, handle.vars.intNodeModFunc, firstTime]; IF print THEN printNodeEqn[handle]; UNTIL nIter > realMaxIter DO OK_ TRUE; IF print THEN handle.msgStream.PutF["t= %10.4f, iteration %2d --\n", IO.real[handle.vars.t], IO.int[nIter]]; IF handle.vars.otherModFunc # NIL THEN updateFunctions[handle, handle.vars.otherModFunc]; IF firstTime THEN updateFunctions[handle, handle.vars.intNodeModFunc]; FOR n_ handle.vars.nodeList, n.nextNode UNTIL n=NIL DO IF n=handle.vars.gndNode THEN LOOP; IF nIter > 1 AND n.converged THEN LOOP; IF print THEN handle.msgStream.PutF["Node %g: ", IO.rope[n.nodeName.name]]; bLink_ n.branches; IF n.integrate THEN { maxI_ csum_ I_ 0.0; UNTIL bLink=NIL DO b _ bLink.branch; plus _ bLink.pos; WITH b.body SELECT FROM r: RefR => newI_ (bLink.otherNode.nHist.y - n.nHist.y)/b.comVal; c: RefC => { csum_ csum + b.comVal; newI_ b.comVal*bLink.otherNode.nHist.f0; }; l: RefL => newI_ (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal; v: RefV => newI_ IF plus THEN -v.vsCurrent ELSE v.vsCurrent; i: RefI => newI_ IF plus THEN -b.comVal ELSE b.comVal ENDCASE; maxI_ MAX[maxI, ABS[newI]]; I_ I + newI; bLink_ bLink.nextLink; ENDLOOP; E_ n.nHist.f0; tempDvDt_ I/csum; dvdt_ n.nHist.f0_ IF ABS[tempDvDt] < handle.vars.ignorDvDtBelow THEN 0.0 ELSE tempDvDt; IF mark THEN n.marked_ (ABS[I] <= handle.vars.dvdtTrunc*maxI); denom_ MAX[ABS[E], ABS[dvdt]]; nodeOK_ IF denom <= Real.SmallestNormalizedNumber THEN TRUE ELSE (handle.vars.tol >= ABS[(dvdt - E)/denom]); IF print THEN handle.msgStream.PutF["Old v'= %9.2f New v'= %9.2f C= %9.2f.\n", IO.real[E], IO.real[dvdt], IO.real[csum]]; } ELSE { E_ n.nHist.y; IF n.curPtr # NIL THEN { vCur_ n.curPtr; vCurBody_ NARROW[vCur.body]; newE_ IF vCur.posNode=n THEN vCur.negNode.nHist.y + vCur.comVal ELSE vCur.posNode.nHist.y - vCur.comVal; I_ vCurBody.vsCurrent; newI_ 0.0; UNTIL bLink=NIL DO b _ bLink.branch; otherV_ bLink.otherNode.nHist.y; plus _ bLink.pos; WITH b.body SELECT FROM r: RefR => newI_ newI + (otherV - newE)/b.comVal; c: RefC => newI_ newI + b.comVal*(bLink.otherNode.nHist.f0 - n.nHist.f0); l: RefL => newI_ newI + (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal; v: RefV => newI_ IF plus THEN newI - v.vsCurrent ELSE newI + v.vsCurrent; i: RefI => newI_ IF plus THEN newI - b.comVal ELSE newI + b.comVal ENDCASE; bLink_ bLink.nextLink ENDLOOP; newI_ IF n=vCur.posNode THEN I + newI ELSE I - newI; denom_ MAX[ABS[I], ABS[newI]]; nodeOK_ IF denom <= Real.SmallestNormalizedNumber THEN TRUE ELSE (handle.vars.tol >= ABS[(newI - I)/denom]); <> <> <<(1.0-handle.vars.squish)*(newE-n.nHist.oldy)/dT +>> <<(IF vCur.posNode=n THEN vCur.negNode.nHist.f0>> <> n.nHist.f0_ 0.0; IF NOT firstTime THEN { tempDvDt _ handle.vars.squish*n.nHist.f1 + (1.0-handle.vars.squish)*(newE-n.nHist.oldy)/dT + (IF vCur.posNode=n THEN vCur.negNode.nHist.f0 ELSE vCur.posNode.nHist.f0); IF ABS[tempDvDt] > handle.vars.ignorDvDtBelow THEN n.nHist.f0_ tempDvDt; }; vCurBody.vsCurrent_ newI; IF print THEN handle.msgStream.PutF["Old I= %9.2f New I= %9.2f", IO.real[I], IO.real[newI]]; } ELSE { cond_ newE_ 0.0; UNTIL bLink=NIL DO b _ bLink.branch; plus _ bLink.pos; WITH b.body SELECT FROM r: RefR => { cond_ cond + 1.0/b.comVal; newE_ newE + bLink.otherNode.nHist.y/b.comVal; }; l: RefL => newE_ newE + (IF plus THEN -l.iHist.y ELSE l.iHist.y)/b.comVal; v: RefV => newE_ IF plus THEN newE - v.vsCurrent ELSE newE + v.vsCurrent; i: RefI => newE_ IF plus THEN newE - b.comVal ELSE newE + b.comVal ENDCASE; bLink_ bLink.nextLink ENDLOOP; newE_ newE/cond; denom_ MAX[ABS[E], ABS[newE]]; nodeOK_ IF denom <= Real.SmallestNormalizedNumber THEN TRUE ELSE (handle.vars.tol >= ABS[(newE - E)/denom]); }; n.nHist.y_ newE; IF print THEN handle.msgStream.PutF["Old E= %9.2f New E= %9.2f.\n", IO.real[E], IO.real[newE]]; }; OK_ OK AND nodeOK; n.converged_ nodeOK; IF ~nodeOK THEN IF nIter = realMaxIter THEN SIGNAL solveError[IO.PutFR["Node %g GS did not converge.", IO.rope[spGlobals.makeStringNB[handle, n, NIL]]]]; FOR bLink_ n.branches, bLink.nextLink UNTIL bLink=NIL DO bLink.otherNode.converged_ FALSE; ENDLOOP; ENDLOOP; -- n nIter_ nIter + 1; IF handle.vars.canned AND NOT handle.vars.forcedDump THEN SIGNAL Aborted; IF nIter=handle.vars.maxIter THEN handle.msgStream.PutF["t= %9.3f, iteration count exceeded", IO.real[handle.vars.t]]; IF OK THEN EXIT; ENDLOOP; -- nIter }; -- GSiteration printBranches: PROC [handle: Handle, n: nodePtr] = { FOR bLink: branchLinkPtr _ n.branches, bLink.nextLink UNTIL bLink=NIL DO branch: branchPtr _ bLink.branch; WITH branch.body SELECT FROM r: RefR => handle.msgStream.PutRope["resistor"]; c: RefC => handle.msgStream.PutRope["capacitor"]; l: RefL => handle.msgStream.PutRope["inductor"]; v: RefV => handle.msgStream.PutRope["voltage source"]; i: RefI => handle.msgStream.PutRope["current source"]; ENDCASE; handle.msgStream.PutF[", %g, %g\n", IO.real[branch.comVal], IO.rope[bLink.otherNode.nodeName.name]]; ENDLOOP; }; printNodeEqn: PROC[handle: Handle]= { n: nodePtr; bLink: branchLinkPtr; branch: branchPtr; plus: BOOL; FOR n_ handle.vars.intNodeList, n.nextIntNode UNTIL n=NIL DO handle.msgStream.PutF["\n%g -- %13.6f %13.6f\n", IO.rope[n.nodeName.name], IO.real[n.nHist.y], IO.real[n.nHist.oldy]]; FOR bLink_ n.branches, bLink.nextLink UNTIL bLink=NIL DO branch_ bLink.branch; plus _ bLink.pos; WITH branch.body SELECT FROM r: RefR => handle.msgStream.PutF["(%13.6f - %13.6f)/%13.6f", IO.real[bLink.otherNode.nHist.y], IO.real[n.nHist.y], IO.real[branch.comVal] ]; c: RefC => handle.msgStream.PutF["%13.6fx%13.6f", IO.real[branch.comVal], IO.real[bLink.otherNode.nHist.f0] ]; l: RefL => IF plus THEN handle.msgStream.PutF["-1x%13.6f/%13.6f", IO.real[l.iHist.y], IO.real[branch.comVal]] ELSE handle.msgStream.PutF["%13.6f/%13.6f", IO.real[l.iHist.y], IO.real[branch.comVal]]; v: RefV => IF plus THEN handle.msgStream.PutF["-1x%13.6f", IO.real[v.vsCurrent]] ELSE handle.msgStream.PutF["%13.6f", IO.real[v.vsCurrent]]; i: RefI => IF plus THEN handle.msgStream.PutF["-1x%13.6f", IO.real[branch.comVal]] ELSE handle.msgStream.PutF["%13.6f",IO.real[branch.comVal]]; ENDCASE; handle.msgStream.PutF[" %s(%13.6f, %13.6f)\n", IO.rope[bLink.otherNode.nodeName.name], IO.real[bLink.otherNode.nHist.y], IO.real[bLink.otherNode.nHist.f0]]; ENDLOOP; ENDLOOP; }; -- printNodeEqn runParms: PROC[handle: Handle]= { OPEN handle.vars; parmValue: PROC[booleanOrNot: {boolean, notBoolean}] RETURNS[r: REAL, b: BOOL, n: NAT]= { IF booleanOrNot=boolean THEN IF item=name THEN { IF RefText.Equal[newString, "FALSE", FALSE] THEN b_ FALSE ELSE IF RefText.Equal[newString, "TRUE", FALSE] THEN b_ TRUE ELSE error[handle, 719]; next[handle] } ELSE error[handle, 704] ELSE { r_ getSignedNumber[handle]; IF r >= 0.0 AND r < Real.Float[LAST[NAT]] THEN n_ Real.FixC[r]; -- else? }; }; -- parmValue nParms: NAT= 36; parmNames: ARRAY[0..nParms) OF REF TEXT= ["maxiter", "printiter", "initstep", "squish", "tol", "dt", "tmax", "printstep", "rungekutta", "tmin", "printall", "amiter", "amup", "inttol", "amdown", "vmin", "vmax", "dvdtmax", "dvdtfac", "tbreak", "rkerr", "amerr", "amdelay", "floor", "minratio", "rkfac", "retratio", "ignorDvDtBelow", "showRetreat", "SaveAllData", "tUnit", "iUnit", "vUnit", "ymin", "ymax", ""]; nameIndex: NAT; tempR: REAL; WHILE item=name DO FOR i: NAT IN [0..nParms) DO nameIndex_ i; IF RefText.Equal[newString, parmNames[nameIndex], FALSE] THEN EXIT; ENDLOOP; next[handle]; IF item # leftArrow THEN error[handle, 703] ELSE next[handle]; SELECT nameIndex FROM 0=> { limit: INT_ LAST[INT]; [tempR, , ]_ parmValue[notBoolean]; IF tempR >= Real.Float[limit] THEN maxIter_ limit ELSE IF tempR > 0 THEN maxIter_ Real.Fix[tempR]; }; 1=> [ , printIter, ]_ parmValue[boolean]; 2=> [initStep, , ]_ parmValue[notBoolean]; 3=> [squish, , ]_ parmValue[notBoolean]; 4=> [tol, , ]_ parmValue[notBoolean]; 5=> [printdT, , ]_ parmValue[notBoolean]; 6=> [tMax, , ]_ parmValue[notBoolean]; 7=> [ , printStep, ]_ parmValue[boolean]; 8=> [ , doRungeKutta, ]_ parmValue[boolean]; 9=> [tMin, , ]_ parmValue[notBoolean]; 10=> [ , printAll, ]_ parmValue[boolean]; 11=> [ , , AMiter]_ parmValue[notBoolean]; 12=> [AMup, , ]_ parmValue[notBoolean]; 13=> [intTol, , ]_ parmValue[notBoolean]; 14=> [AMdown, , ]_ parmValue[notBoolean]; 15=> [Vmin, , ]_ parmValue[notBoolean]; 16=> [Vmax, , ]_ parmValue[notBoolean]; 17=> [dvdtmax, , ]_ parmValue[notBoolean]; 18=> [dvdtFac, , ]_ parmValue[notBoolean]; 19=> [tBreak, , ]_ parmValue[notBoolean]; 20=> [RKerr, , ]_ parmValue[notBoolean]; 21=> [AMerr, , ]_ parmValue[notBoolean]; 22=> [ , , AMdelay]_ parmValue[notBoolean]; 23=> [floor, , ]_ parmValue[notBoolean]; 24=> [minRatio, , ]_ parmValue[notBoolean]; 25=> [RKfac, , ]_ parmValue[notBoolean]; 26=> [retRatio, , ]_ parmValue[notBoolean]; 27=> [ignorDvDtBelow, , ]_ parmValue[notBoolean]; 28=> [ , showRetreat, ]_ parmValue[boolean]; 29=> [ , saveAllData, ]_ parmValue[boolean]; 30=> [tUnit, , ]_ parmValue[notBoolean]; 31=> [iUnit, , ]_ parmValue[notBoolean]; 32=> [vUnit, , ]_ parmValue[notBoolean]; 33=> [yMin, , ]_ parmValue[notBoolean]; 34 => [yMax, , ]_ parmValue[notBoolean]; ENDCASE=> error[handle, 790]; IF item=comma THEN next[handle]; ENDLOOP; }; -- runParms predict: PROC[handle: Handle, h: REAL]= { nodes: nodePtr_ handle.vars.intNodeList; inds: branchPtr_ handle.vars.inductorList; incr: REAL; h_ h/24.0; UNTIL nodes=NIL DO IF ~nodes.marked THEN { OPEN nodes.nHist; incr_ h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4); y_ oldy + incr; }; nodes_ nodes.nextIntNode; ENDLOOP; UNTIL inds=NIL DO indsBody: RefL_ NARROW[inds.body]; {OPEN indsBody.iHist; y_ oldy + h*(55.0*f1 - 59.0*f2 + 37.0*f3 - 9.0*f4)}; inds_ indsBody.nextInductor; ENDLOOP; }; -- predict firstCorrector: PROC[handle: Handle, h: REAL]= { nodes: nodePtr_ handle.vars.intNodeList; inds: branchPtr_ handle.vars.inductorList; f: REAL; h_ h/24.0; UNTIL nodes=NIL DO IF ~nodes.marked THEN { OPEN nodes.nHist; y_ oldy + h*(9.0*f0 + 19.0*f1 - 5.0*f2 + f3); oldf_ f0; }; nodes_ nodes.nextIntNode; ENDLOOP; UNTIL inds=NIL DO indsBody: RefL_ NARROW[inds.body]; {OPEN indsBody.iHist; f_ inds.posNode.nHist.y - inds.negNode.nHist.y; y_ oldy + h*(9.0*f + 19.0*f1 - 5.0*f2 + f3); oldf_ f}; inds_ indsBody.nextInductor; ENDLOOP; }; -- firstCorrector iteratingCorrector: PROC[handle: Handle, h: REAL] RETURNS[newRatio: REAL_ 4.0, OK: BOOL_ TRUE]= { nodes: nodePtr_ handle.vars.intNodeList; inds: branchPtr_ handle.vars.inductorList; f, deltay, newy, err, denom: REAL; h_ 3.0*h/8.0; UNTIL nodes=NIL DO IF ~nodes.marked THEN { OPEN nodes.nHist; deltay_ h*(f0 - oldf); newy_ y + deltay; denom_ MAX[ABS[newy], ABS[y]]; IF denom > handle.vars.floor AND deltay > 1.0E-30 THEN { err_ handle.vars.AMerr*ABS[deltay/denom]; IF handle.vars.intTol < err THEN { OK_ FALSE; newRatio_ MIN[newRatio, handle.vars.intTol/err]; handle.vars.worstNode_ nodes; }; }; y_ newy; oldf_ f0; }; nodes_ nodes.nextIntNode; ENDLOOP; UNTIL inds=NIL DO indsBody: RefL_ NARROW[inds.body]; { OPEN indsBody.iHist; f_ inds.posNode.nHist.y - inds.negNode.nHist.y; deltay_ h*(f - oldf); newy_ y + deltay; denom_ MAX[ABS[newy], ABS[y]]; IF denom > handle.vars.floor AND deltay > 1.0E-30 THEN { err_ handle.vars.AMerr*ABS[deltay/denom]; IF handle.vars.intTol < err THEN { OK_ FALSE; newRatio_ MIN[newRatio, handle.vars.intTol/err]; }; }; y_ newy; oldf_ f; }; inds_ indsBody.nextInductor; ENDLOOP; }; -- iteratingCorrector AdamsMoulton: PROC[handle: Handle, dT: REAL] RETURNS[newRatio: REAL, nGS: INT, ok: BOOL]= { nAM: NAT; predict[handle, dT]; nGS_ GSiteration[handle, dT, handle.vars.printAll]; firstCorrector[handle, dT]; nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS; FOR nAM_ 2, nAM + 1 UNTIL nAM=handle.vars.AMiter DO []_ iteratingCorrector[handle, dT]; nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS ENDLOOP; [newRatio, ok]_ iteratingCorrector[handle, dT]; }; -- AdamsMoulton RKupdate: PROC[handle: Handle, c, h: REAL, last: BOOL] RETURNS[newRatio: REAL_ 1.0, ok: BOOL_ TRUE]= { nodes: nodePtr_ handle.vars.intNodeList; inds: branchPtr_ handle.vars.inductorList; k, incr, y, err, denom: REAL; UNTIL nodes=NIL DO IF ~nodes.marked THEN { nodes.nHist.sumk_ c*nodes.nHist.f0 + nodes.nHist.sumk; IF last THEN { incr_ h*nodes.nHist.sumk; y_ nodes.nHist.oldy + incr; denom_ MAX[ABS[y], ABS[nodes.nHist.y]]; IF denom > handle.vars.floor THEN { err_ handle.vars.RKerr*ABS[(y - nodes.nHist.y)/denom]; IF handle.vars.intTol < err THEN { ok_ FALSE; IF handle.vars.countsPerStep=1 THEN SIGNAL solveError[IO.PutFR["\nNode %g convergence problems.", IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]]]]; newRatio_ MIN[newRatio, handle.vars.intTol/err] }; }; nodes.nHist.y_ y; nodes.nHist.sumk_ 0.0; } ELSE { incr_ h*nodes.nHist.f0; nodes.nHist.y_ nodes.nHist.oldy + incr; }; }; nodes_ nodes.nextIntNode; ENDLOOP; UNTIL inds=NIL DO indsBody: RefL_ NARROW[inds.body]; k_ inds.posNode.nHist.y - inds.negNode.nHist.y; indsBody.iHist.sumk_ c*k + indsBody.iHist.sumk; IF last THEN { indsBody.iHist.y_ indsBody.iHist.oldy + h*indsBody.iHist.sumk; indsBody.iHist.sumk_ 0.0 } ELSE indsBody.iHist.y_ indsBody.iHist.oldy + h*k; inds_ indsBody.nextInductor; ENDLOOP; }; -- RKupdate RungeKutta: PROC[handle: Handle, dT: REAL] RETURNS[newRatio: REAL, nGS: NAT, ok: BOOL]= { time: REAL; time_ handle.vars.t; handle.vars.t_ handle.vars.t - 0.5*dT; [] _ RKupdate[handle, 1.0, 0.5*dT, FALSE]; nGS_ GSiteration[handle, 0.5*dT, handle.vars.printAll]; [] _ RKupdate[handle, 2.0, 0.5*dT, FALSE]; nGS_ GSiteration[handle, 0.5*dT, handle.vars.printAll] + nGS; handle.vars.t_ time; [] _ RKupdate[handle, 2.0, dT, FALSE]; nGS_ GSiteration[handle, dT, handle.vars.printAll] + nGS; [newRatio, ok]_ RKupdate[handle, 1.0, 0.1666667*dT, TRUE] }; -- RungeKutta changeStepSize: PROC[handle: Handle, ratio: REAL, RK: BOOL]= { int: CHAR; IF ratio < 1.0 OR handle.vars.countsPerStep < maxCountsPerStep THEN { handle.vars.numGoodSteps_ 1; handle.vars.worstNodeLog[handle.vars.curLog]_ IF handle.vars.worstNode # NIL THEN [visited: FALSE, node: handle.vars.worstNode, t: handle.vars.t, v: handle.vars.worstNode.nHist.y, dvdt: handle.vars.worstNode.nHist.f0] ELSE [FALSE, NIL, 0.0, 0.0, 0.0]; handle.vars.curLog_ (handle.vars.curLog + 1) MOD maxLog; IF handle.vars.countsPerStep=1 THEN SIGNAL solveError["Step too small"]; handle.vars.countsPerStep_ IF ratio > Real.Float[maxCountsPerStep/handle.vars.countsPerStep] THEN maxCountsPerStep ELSE MAX[1, Real.RoundLI[handle.vars.countsPerStep*ratio]]; int_ IF RK THEN 'R ELSE 'A; IF ShowDetails[handle] THEN { Labels.Set[handle.time, IO.PutFR["%12.6e", IO.real[handle.vars.t]]]; Labels.Set[handle.step, IO.PutFR["%10d %g", IO.int[handle.vars.countsPerStep], IO.char[int]]]; }; }; }; -- changeStepSize saveState: PROC[handle: Handle]= { inds: branchPtr_ handle.vars.inductorList; volts: branchPtr_ handle.vars.vSourceList; indsBody: RefL; voltsBody: RefV; FOR nodes: nodePtr_ handle.vars.nodeList, nodes.nextNode UNTIL nodes=NIL DO OPEN nodes.nHist; SolveTrouble: PROC [explanation: Rope.ROPE] = { SIGNAL solveTest[IO.PutFR["Node %g %g", IO.rope[spGlobals.makeStringNB[handle, nodes, NIL]], IO.rope[explanation]]]; }; oldy_ y; f4_ f3; f3_ f2; f2_ f1; f1_ f0; SELECT TRUE FROM y < handle.vars.Vmin => SolveTrouble["below Vmin"]; y > handle.vars.Vmax => SolveTrouble["above Vmax"]; ABS[f0] > handle.vars.dvdtmax => SolveTrouble["change rate too large"]; ENDCASE => NULL; ENDLOOP; UNTIL inds=NIL DO indsBody_ NARROW[inds.body]; { OPEN indsBody.iHist; oldy_ y; f4_ f3; f3_ f2; f2_ f1; f1_ inds.posNode.nHist.y - inds.negNode.nHist.y; }; inds_ indsBody.nextInductor; ENDLOOP; UNTIL volts=NIL DO voltsBody_ NARROW[volts.body]; voltsBody.oldCurrent_ voltsBody.vsCurrent; volts_ voltsBody.nextvSource; ENDLOOP; }; -- saveState restoreState: PROC[handle: Handle]= { nodes: nodePtr_ handle.vars.nodeList; inds: branchPtr_ handle.vars.inductorList; volts: branchPtr_ handle.vars.vSourceList; indsBody: RefL; voltsBody: RefV; UNTIL nodes=NIL DO nodes.nHist.y_ nodes.nHist.oldy; nodes.nHist.f0_ nodes.nHist.f1; IF ABS[nodes.nHist.f0] < handle.vars.ignorDvDtBelow THEN nodes.nHist.f0_ 0.0; nodes_ nodes.nextNode; ENDLOOP; UNTIL inds=NIL DO indsBody_ NARROW[inds.body]; indsBody.iHist.y_ indsBody.iHist.oldy; inds_ indsBody.nextInductor; ENDLOOP; UNTIL volts=NIL DO voltsBody_ NARROW[volts.body]; voltsBody.vsCurrent_ voltsBody.oldCurrent; volts_ voltsBody.nextvSource; ENDLOOP; }; -- restoreState integration: PROC[handle: Handle, RK: BOOL]= { nGS: INT; out, ok: BOOL_ FALSE; ratio, oldTime, DT: REAL; oldLowStepcount, oldHighStepCount: INT; saveState[handle]; oldTime_ handle.vars.t; oldLowStepcount_ handle.vars.lowStepCount; oldHighStepCount_ handle.vars.highStepCount; DO handle.vars.lowStepCount_ handle.vars.lowStepCount + handle.vars.countsPerStep; IF handle.vars.lowStepCount > maxCountsPerStep THEN { handle.vars.highStepCount_ handle.vars.highStepCount + 1; handle.vars.lowStepCount_ handle.vars.lowStepCount - maxCountsPerStep; out_ TRUE; } ELSE out_ handle.vars.printStep; handle.vars.t_ handle.vars.tMin + handle.vars.dT*handle.vars.lowStepCount + (handle.vars.dT*maxCountsPerStep)*handle.vars.highStepCount; DT_ handle.vars.dT*handle.vars.countsPerStep; { ENABLE Retreat=> { IF ~Rope.Equal[cause, handle.vars.lastRetreatCause] THEN { IF Rope.Length[cause] > 0 AND handle.vars.showRetreat THEN handle.msgStream.PutF["%g at t= %12.6e.\n", IO.rope[cause], IO.real[handle.vars.t]]; handle.vars.lastRetreatCause_ cause; IF handle.vars.countsPerStep = 1 THEN REJECT; }; ok_ FALSE; ratio_ handle.vars.retRatio; GOTO FlushIt; }; IF handle.vars.numGoodSteps < 4 THEN { [ratio, nGS, ok]_ RungeKutta[handle, DT]; IF ok THEN nGS_ nGS + GSiteration[handle, DT, handle.vars.printAll,, TRUE]; ratio_ RealFns.SqRt[ratio]*handle.vars.RKfac; IF ok AND ~RK THEN handle.vars.numGoodSteps_ 1 + handle.vars.numGoodSteps; } ELSE { handle.vars.AMcount_ 1 + handle.vars.AMcount; [ratio, nGS, ok]_ AdamsMoulton[handle, DT]; IF ok THEN { nGS_ nGS + GSiteration[handle, DT, handle.vars.printAll,, TRUE]; IF ratio > handle.vars.AMup AND handle.vars.AMcount >= handle.vars.AMdelay THEN { ratio_ RealFns.SqRt[RealFns.SqRt[ratio]]; changeStepSize[handle, ratio, FALSE]; handle.vars.AMcount_ 0; }; } ELSE ratio_ RealFns.SqRt[ratio]; }; EXITS FlushIt=> NULL; }; IF ok THEN EXIT; restoreState[handle]; handle.vars.t_ oldTime; handle.vars.highStepCount_ oldHighStepCount; handle.vars.lowStepCount_ oldLowStepcount; changeStepSize[handle, MIN[ratio, handle.vars.minRatio], (handle.vars.numGoodSteps < 4)]; handle.vars.AMcount_ 0 ENDLOOP; NextStage[handle]; IF out THEN { printFromList[handle, nGS, handle.vars.t, handle.vars.printIter]; plotFromList[handle, handle.vars.t]; IF handle.saveAllData THEN UpdatePP[handle, handle.vars.t]; }; WHILE handle.vars.tBreak > 0.0 AND handle.vars.t >= handle.vars.tBreak DO oldTime_ handle.vars.t ENDLOOP; }; -- integration zeroVars: PROC[handle: Handle]= { n: nodePtr_ handle.vars.nodeList; inds: branchPtr_ handle.vars.inductorList; vSs: branchPtr_ handle.vars.vSourceList; indsBody: RefL; vSsBody: RefV; handle.vars.gndNode.marked_ TRUE; UNTIL n=NIL DO n.nHist_ []; n.marked_ FALSE; n_ n.nextNode; ENDLOOP; UNTIL inds=NIL DO indsBody_ NARROW[inds.body]; indsBody.iHist_ []; inds_ indsBody.nextInductor; ENDLOOP; UNTIL vSs=NIL DO vSsBody_ NARROW[vSs.body]; vSsBody.vsCurrent_ vSsBody.oldCurrent_ 0.0; IF vSs.controller=NIL THEN IF vSs.posNode=handle.vars.gndNode THEN vSs.negNode.marked_ TRUE ELSE IF vSs.negNode=handle.vars.gndNode THEN vSs.posNode.marked_ TRUE; vSs_ vSsBody.nextvSource; ENDLOOP; }; -- zeroVars getStepSize: PROC[handle: Handle] RETURNS[stepOK: BOOL]= { vMax, dvdtMax, frac: REAL_ 0.0; inds: branchPtr_ handle.vars.inductorList; nodes: nodePtr_ handle.vars.nodeList; handle.vars.lowStepCount_ 0; handle.vars.highStepCount_ 0; handle.vars.numberOfSteps_ Real.FixI[handle.vars.tMax/handle.vars.printdT]; handle.vars.worstNode_ NIL; UNTIL nodes=NIL DO vMax_ MAX[ABS[nodes.nHist.y], vMax]; IF ABS[nodes.nHist.f0] > dvdtMax THEN { dvdtMax_ ABS[nodes.nHist.f0]; handle.vars.worstNode_ nodes; }; nodes_ nodes.nextNode; ENDLOOP; frac_ IF dvdtMax > 0.0 AND vMax > 0.0 THEN handle.vars.intTol*vMax/(dvdtMax*handle.vars.printdT) ELSE handle.vars.initStep; handle.vars.dT_ handle.vars.printdT/maxCountsPerStep; handle.vars.dvdtTrunc_ handle.vars.dvdtFac*handle.vars.intTol; handle.vars.countsPerStep_ Real.RoundLI[maxCountsPerStep*MIN[frac, handle.vars.initStep]]; handle.vars.countsPerStep_ MAX[handle.vars.countsPerStep, 1]; IF handle.vars.countsPerStep=1 THEN { stepOK_ FALSE; PutMsgLine[handle, "Initial step is too small!"]; } ELSE { stepOK_ TRUE; <> <> }; }; -- getStepSize setICs: PUBLIC PROC[handle: Handle]= { n: nodePtr; b: branchPtr; v: REAL; allNodes: BOOL; zeroVars[handle]; IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle]; IF handle.vars.item=number THEN { handle.vars.tMin_ handle.vars.value; next[handle]; IF handle.vars.item=comma THEN next[handle] ELSE error[handle, 702]; }; DO allNodes_ (handle.vars.item=star); IF allNodes THEN next[handle] ELSE [n, b]_ findNodeOrBranch[handle]; IF handle.vars.item=leftArrow THEN next[handle] ELSE error[handle, 703]; v_ getSignedNumber[handle]; IF allNodes THEN FOR n_ handle.vars.nodeList, n.nextNode UNTIL n=NIL DO n.nHist.y_ v; ENDLOOP ELSE IF n # NIL THEN n.nHist.y_ v ELSE IF b # NIL THEN WITH b.body SELECT FROM l: RefL => l.iHist.y_ v; ENDCASE=> error[handle, 721] ELSE error[handle, 720]; IF handle.vars.item=comma THEN next[handle] ELSE EXIT; ENDLOOP; IF handle.vars.item # rightB THEN error[handle, 701, TRUE] ELSE next[handle]; handle.vars.icsSet_ TRUE; }; -- setICs runIt: PUBLIC PROC[handle: Handle] RETURNS[REAL]= { n: INT; IF handle.vars.item # leftB THEN error[handle, 700] ELSE next[handle]; handle.vars.title _ "Post-Processing Graph"; IF handle.vars.item = string THEN { handle.vars.title _ Rope.FromRefText[handle.vars.newString]; next[handle]; IF handle.vars.item = comma THEN next[handle] ELSE error[handle, 603] }; runParms[handle]; IF handle.vars.item # rightB THEN error[handle, 701] ELSE next[handle]; IF handle.vars.tMax <= handle.vars.tMin THEN error[handle, 731]; IF handle.vars.errCount <= 0 THEN { ENABLE Failure => {error[handle, errorNum]; GOTO quit}; IF ~handle.vars.icsSet THEN zeroVars[handle]; handle.vars.simTime_ BasicTime.Now[]; initPrint[handle, handle.vars.simTime]; initPlot[handle, handle.vars.tMin, handle.vars.tMin + handle.vars.tMax]; IF handle.saveAllData THEN InitPP[handle, handle.vars.tMin, handle.vars.tMin + handle.vars.tMax]; handle.vars.t_ handle.vars.tMin; IF handle.vars.printdT <= 0.0 THEN handle.vars.printdT_ handle.vars.tMax/handle.vars.dfltInts; handle.vars.dT_ handle.vars.printdT; n_ GSiteration[handle, 0.0, handle.vars.printAll, TRUE, TRUE]; printFromList[handle, n, handle.vars.t, handle.vars.printIter]; plotFromList[handle, handle.vars.t]; IF handle.saveAllData THEN UpdatePP[handle, handle.vars.t]; IF ~getStepSize[handle] THEN RETURN[handle.vars.tMin]; handle.vars.numGoodSteps_ 1; WHILE handle.vars.highStepCount <= handle.vars.numberOfSteps DO integration[handle, handle.vars.doRungeKutta]; IF Canned[handle] THEN EXIT; IF CheckPoint[handle] THEN { ResetCheckPoint[handle]; dumpAll[handle, handle.vars.t]; }; ENDLOOP; handle.vars.icsSet_ FALSE; RETURN[handle.vars.t]; EXITS quit=> RETURN[handle.vars.t]; } ELSE RETURN[handle.vars.tMin] }; -- runIt END. CHANGE LOG Wilhelm, April 27, 1982 4:04 PM Barth, 7-May-82 10:32:52 PDT Chen, April 22, 1983 10:56 AM. Changed dvdtmax from 1E13 to 1E15. Barth, July 11, 1983 1:42 PM Chen, February 12, 1984 8:13 PM modified to support oldArgVector. Chen, June 12, 1984 6:10:21 pm PDT, cedarized. Barth, February 8, 1985 3:17:34 pm PST McCreight, April 4, 1985 1:57:23 pm PST Chen, May 9, 1985 6:40:28 pm PDT: 1. RealOps.SqRt => RealFns.SqRt; 2. moved lastRetreatCause from a global variable to be one on the handle.vars. 3. in order to get around the concurrency problem in floating point software, added a new run-parameter ignorDvDtBelow to ignor dvdt when its magnitude is smaller than ignorDvDtBelow. 4. added another new run-parameter, showRetreat, so that if one really wants to see the reasons of retreats during "integration", he may do so by setting it true in the run statement. Chen, July 22, 1985 8:09:51 pm PDT, => Cedar6.0. Chen, November 24, 1985 3:21:27 pm PST, added post-processing capability. Barth, November 12, 1986 3:53:53 pm PST, changed GSIteration to pass firstTime to updateFunctions and changed updateFunctions to always evaluate the models the first time. Controlled branches were not being initialized if both of the nodes to which they were connected were constant nodes.