<> <> <> <> <> <> <> 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], ThymeGlobals 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]; ThymeSolve: CEDAR PROGRAM IMPORTS BasicTime, IO, Labels, Real, RealFns, RefText, ThymeGlobals EXPORTS ThymeGlobals = BEGIN OPEN ThymeGlobals; 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[ThymeGlobals.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[ThymeGlobals.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[ThymeGlobals.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.tMax]; IF handle.saveAllData THEN InitPP[handle, 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.