ThymeSolve.mesa
Copyright © 1985 by Xerox Corporation. All rights reserved.
Last Edited by:
Christian Le Cocq April 29, 1987 10:25:01 am PDT
Barth, November 12, 1986 3:55:21 pm PST
Pradeep Sindhu March 19, 1986 11:45:02 pm PST
Sweetsun Chen, November 24, 1985 4:02:56 pm PST
DIRECTORY
BasicTime USING [Now],
IO USING [char, int, PutF, PutFR, PutRope, real, rope],
Labels USING [Set],
Real USING [Fix, InlineFixC, InlineFixI, Float, Round, 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 {
since ~allInactive, from above we know that source is not Nil
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]);
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);
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.InlineFixC[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.Round[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.InlineFixI[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.Round[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;
handle.msgStream.PutF["Initial step is %d.\n",
IO.int[handle.vars.countsPerStep]];
};
}; -- 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.