DIRECTORY BasicTime USING [Now], IO USING [char, int, PutF, PutFR, real, rope], Labels USING [Set], Real USING [Fix, FixC, FixI, Float, RoundLI, SmallestNormalizedNumber], RealFns USING [SqRt], RefText USING [Equal], Rope USING [Equal, Length, ROPE], spGlobals USING [Aborted, argList, argNodes, branchLinkPtr, branchPtr, Canned, CheckPoint, dumpAll, error, findNodeOrBranch, getSignedNumber, Handle, initPlot, initPrint, makeStringNB, maxCountsPerStep, maxLog, modBrPtr, modFuncPtr, next, NextStage, nodePtr, plotFromList, printFromList, PutMsgLine, RefFuncBody, RefModBody, RefC, RefI, RefL, RefR, RefV, ResetCheckPoint, ShowDetails]; 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]= { 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_ (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 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]; 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]); 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 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= 30; 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", ""]; 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=> [initStep, , ]_ 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]; 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]; IF out THEN plotFromList[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]; 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]; 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 ~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 $File: [Cherry]Cedar5.2>System>spSolve.mesa Last Edited by: SChen, May 28, 1985 6:54:27 pm PDT since ~allInactive, from above we know that source is not Nil n.nHist.f0_ IF firstTime THEN 0.0 ELSE 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); handle.msgStream.PutF["Initial step is %d.\n", IO.int[handle.vars.countsPerStep]]; 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. Ê:˜J™1J™2J˜šÏk ˜ Jšœ œ˜Jšœœ&˜.Jšœœ˜Jšœœ=˜GJšœœ˜Jšœœ ˜Jšœœœ˜!Jšœ œò˜J˜—šœ œ˜Jšœ œ+˜@Jšœ ˜J˜—Jšœœ ˜J˜Jš œ œœ œœ˜/Jš œ œœ œœ˜,Jšœ œ œœ˜'Jšœ œ œœ˜&J˜šœœ"˜7J˜ J˜J˜Jšœ œ˜J˜J˜šœœ˜J˜J˜Jšœ œ˜šœ œ˜šœœœ˜!J˜Jšœœ˜0Jšœ˜——šœœ˜˜J˜GJ˜—˜šœœ˜Jšœ=™=J˜#šœœœ˜!J˜!—Jšœ˜J˜J˜0šœœœ˜!J˜—Jšœ˜J˜šœœ˜J˜"J˜—Jšœ˜J˜—Jšœ:˜>—Jšœ˜—J˜Jšœ˜—JšœÏc˜J˜—š Ïn œœœœœ˜PJšœœ˜Jšœœ3œ˜HJ˜ J˜ J˜J˜J˜Jšœœœ˜šœ œœœœœœœ˜GJšœ+˜/—J˜J˜4Jšœœ˜#šœ˜Jšœœ˜ šœœ7˜DJšœœ ˜'—šœœ˜&J˜2—Jšœ œ5˜Fšœ%œœ˜6Jšœœœ˜#Jšœ œ œœ˜'šœœ#˜0Jšœ˜—J˜šœ œ˜Jšœ œ˜šœœ˜J˜J˜šœœ˜J˜@˜ J˜J˜(J˜—Jšœœœ œ˜CJšœœœœ ˜Jš œœœœœ˜šœœ(œ˜;Jšœœ œ ˜0—šœœA˜NJšœœœ œ ˜*—J˜—šœ˜Jšœ ˜ šœ œœ˜J˜Jšœ œ ˜šœœœ#˜?Jšœ$˜(—Jšœ˜J˜ šœœ˜J˜J˜ J˜šœœ˜J˜1J˜IJšœœœ œ˜JJšœœœœ˜IJšœœœœ˜BJšœ˜—J˜Jšœ˜—Jš œœœœœœ˜4Jš œœœœœ˜šœœ(œ˜;Jšœœ œ ˜0—šœ œ œ™&™ J™1šœœœ™-Jšœ™———Jšœ˜šœœ œ˜˜+J˜1Jšœœœœ˜J—Jšœœ(œ˜HJ˜—J˜šœœ3˜@Jšœœœ ˜—J˜—šœ˜J˜šœœ˜J˜J˜šœœ˜˜ J˜J˜.J˜—Jšœœœ œ˜Jšœœœ˜0Jšœ˜—šœœœ˜-Jšœ˜—Jšœ˜—J˜Jšœ˜—J˜Jš œœœœœ˜šœœ(œ˜;Jšœœ œ ˜0—J˜—J˜šœœ6˜CJšœœœ ˜—J˜—Jšœœœ˜J˜šœ ˜šœ˜šœ œ&˜:Jšœ(œ˜2——šœ#œœ˜8Jšœœ˜!Jšœ˜——Jšœž˜ —J˜Jš œœœœœ ˜Išœ˜!šœ;˜;Jšœ˜——Jšœœœœ˜Jšœž˜—Jšœž˜J˜—šœœ˜%J˜ J˜J˜Jšœœ˜ J˜šœ+œœ˜<šœ0˜0Jšœ˜Jšœ˜Jšœ˜—šœ#œœ˜8J˜J˜šœ œ˜šœ<˜šœ ˜˜Jšœœœœ˜J˜#Jšœœ˜1Jšœœ œ˜0J˜—J˜)J˜*J˜*J˜%J˜)J˜&J˜)J˜,J˜&J˜)J˜*J˜'J˜)J˜)J˜'J˜'J˜*J˜*J˜)J˜(J˜(J˜+J˜(J˜+J˜(J˜+J˜1Jšœ-˜-Jšœ˜—Jšœ œ˜ Jšœ˜—Jšœž ˜J˜—šœ œœ˜)J˜(J˜*Jšœœ˜ J˜J˜ J˜šœœ˜šœœ˜Jšœ ˜J˜/J˜J˜—J˜Jšœ˜—J˜šœœ˜Jšœœ ˜"šœœ˜J˜4—J˜Jšœ˜—Jšœž ˜ J˜—šœœœ˜0J˜(J˜*Jšœœ˜J˜J˜ J˜šœœ˜šœœ˜Jšœ ˜J˜-J˜ J˜—J˜Jšœ˜—J˜šœœ˜Jšœœ ˜"šœœ˜J˜/J˜,J˜ —J˜Jšœ˜—Jšœž˜J˜—š œœœœ œ˜NJšœœœ˜J˜(J˜*Jšœœ˜"J˜J˜ šœœ˜šœœ˜Jšœ ˜J˜J˜Jšœœœœ˜šœœœ˜8Jšœœ˜)šœœ˜"Jšœœ˜ Jšœ œ#˜0J˜J˜—J˜—J˜J˜ J˜—J˜Jšœ˜—J˜šœœ˜Jšœœ ˜"šœœ˜J˜/J˜J˜Jšœœœœ˜šœœœ˜8Jšœœ˜)šœœ˜"Jšœœ˜ Jšœ œ#˜0J˜—J˜—J˜J˜J˜—J˜Jšœ˜—Jšœž˜J˜—šŸ œœœ˜,Jšœ œœœ˜.Jšœœ˜ J˜J˜J˜3J˜J˜9šœœ˜3J˜#J˜8Jšœ˜—J˜/Jšœž˜J˜—šŸœœœœ˜6Jšœ œ œœ˜/J˜(J˜*Jšœœ˜J˜šœœ˜šœœ˜J˜6šœœ˜J˜J˜Jšœœœœ˜'šœœ˜#Jšœœ˜6šœœ˜"Jšœœ˜ šœ˜#šœ œ)˜=Jšœ,œ˜6——Jšœ œ"˜/J˜—J˜—J˜J˜J˜—šœ˜J˜J˜'J˜—J˜—J˜Jšœ˜—J˜šœœ˜Jšœœ ˜"J˜/J˜/šœœ˜J˜>J˜J˜—Jšœ-˜1J˜Jšœ˜—Jšœž ˜J˜—šŸ œœœ˜*Jšœ œœœ˜.Jšœœ˜ J˜J˜J˜&Jšœ#œ˜*J˜7Jšœ#œ˜*J˜=J˜Jšœœ˜&J˜9Jšœ4œ˜9Jšœž ˜J˜—š œœœœœ˜>Jšœœ˜ J˜šœ œ.œ˜EJ˜J˜˜-šœœœ˜$Jšœ œx˜‡—Jšœœœ˜!J˜—Jšœ-œ˜8J˜Jšœœœ˜HJ˜˜šœ@˜FJ˜—Jšœœ3˜;J˜—Jš œœœœœ˜šœœ˜˜Jšœœ˜,—˜Jšœœ!œ ˜F—J˜—J˜—Jšœž˜J˜—šœ œ˜"J˜*J˜*J˜J˜J˜šœ6œœ˜KJšœ ˜J˜šŸ œœœ˜/Jš œ œœ,œœ˜tJ˜J˜—J˜J˜šœœ˜Jšœ3˜3Jšœ3˜3JšœD˜GJšœœ˜—Jšœ˜—J˜šœœ˜Jšœ œ ˜šœœ˜J˜J˜J˜0J˜—J˜Jšœ˜—J˜šœœ˜Jšœ œ ˜J˜*J˜Jšœ˜—Jšœž ˜J˜—šœœ˜%J˜%J˜*J˜*J˜J˜J˜šœœ˜J˜ J˜Jšœœ.œ˜MJ˜Jšœ˜—J˜šœœ˜Jšœ œ ˜J˜&J˜Jšœ˜—J˜šœœ˜Jšœ œ ˜J˜*J˜Jšœ˜—Jšœž˜J˜—šœ œœœ˜.Jšœœ˜ Jšœ œ˜Jšœœœ˜Jšœ#œ˜'J˜J˜J˜J˜*J˜,š˜J˜Ošœ-œ˜5J˜9J˜FJšœœ˜ J˜—Jšœ˜ J˜˜KJ˜<—Jšœ+˜-J˜šœœ ˜šœ2œ˜:šœœ˜:˜+Jšœœ˜(——Jšœ$˜$Jšœœœ˜-Jšœ˜—Jšœœ˜ J˜Jšœ ˜ J˜J˜šœœ˜&Jšœ%œ˜)šœœ ˜Jšœœœ˜5—J˜-šœœœ˜J˜7—J˜—šœ˜J˜-Jšœ'œ˜+šœœ˜ Jšœœœ˜@šœ˜Jšœ+œ˜1J˜)Jšœœ˜%J˜J˜—J˜—Jšœ˜ J˜——š˜Jšœ œ˜—J˜Jšœœœ˜J˜J˜J˜,J˜*šœœ˜8J˜ —J˜—Jšœ˜J˜J˜JšœœB˜MJšœœ%˜0šœœ%˜IJšœœ˜—Jšœž˜J˜—šœ œ˜!J˜!J˜*J˜(J˜J˜J˜Jšœœ˜!J˜šœœ˜J˜ Jšœ œ˜J˜Jšœ˜—J˜šœœ˜Jšœ œ ˜J˜J˜Jšœ˜—J˜šœœ˜Jšœ œ ˜J˜+šœœ˜Jšœ!œ˜@š˜Jšœ!œœ˜A——J˜Jšœ˜—Jšœž ˜J˜—šœ œœ œ˜:Jšœœ˜J˜*J˜%J˜J˜J˜J˜KJšœœ˜J˜šœœ˜Jšœœœ˜$šœœœ˜'Jšœ œ˜J˜J˜—J˜Jšœ˜—J˜šœœœ ˜*J˜5Jšœ˜—J˜5J˜>˜Jšœœ˜?—˜Jšœ˜"—šœœ˜%Jšœœ˜J˜1J˜—šœ˜Jšœœ˜ šœ.™.Jšœ#™#—J˜—Jšœž˜J˜—šœœœ˜&J˜ J˜ Jšœœ˜Jšœ œ˜J˜J˜Jšœœœ˜Fšœœ˜!J˜$J˜ Jšœœœ˜DJ˜—š˜J˜"Jšœ œœ"˜DJšœœœ˜HJ˜šœ ˜šœ%œœ˜6J˜ Jš˜——š˜Jšœœœ ˜š˜šœœ˜šœœ˜J˜Jšœ˜——Jšœ˜——Jšœœœœ˜6Jšœ˜—Jšœœœœ˜MJšœœ˜Jšœž ˜ J˜—š œœœœœ˜3Jšœœ˜J˜Jšœœœ˜FJ˜Jšœœœ˜GJšœ&œ˜@šœœ˜#Jšœ&œ˜7Jšœœ˜-J˜%J˜'J˜HJ˜ šœ˜"J˜;—J˜$Jšœ2œœ˜>J˜?J˜$Jšœœœ˜6J˜šœ8˜?J˜.Jšœœœ˜šœœ˜J˜J˜J˜—Jšœ˜—Jšœœ˜Jšœ˜š˜Jšœœ˜—J˜—Jšœœ˜Jšœž˜ J˜—Jšœ˜J˜šœ˜ J™ Jšœ™J™BJ™J™BJ™.J™&J™'™!J™!JšœN™NJ™·Jšœ·™·———…—^šø